text:
function [A,B,C,G] = ospls(X,Y,nhidden,opts)

% OSPLS  Sparse orthogonalized partial least squares where the elastic
% net is implemented using either a native implementation or the glmnet package
%
% The native implementation supports arbitrary alpha (ridge penalty) and
% lambda (lasso penalty) values and supports a coupling matrix for alpha to
% implement smoothing. It requires a fixed value for lambda.
%
% The glmnet implementation is much faster but uses alpha to mix between
% ridge and lasso while lambda determines that amount of regularization. It
% also does not support a coupling matrix. It learns the optimal lambda
% using inner cross-validation.
%
% X: ninput x nsamples input data matrix
% Y: noutput x nsamples output data matrix
% nhidden: number of components [1]
% opts: glmnet parameters
%
% A: noutput x nhidden weight matrix
% B: ninput x nhidden weight matrix
% C: 1 x nhidden bias vector
% G: output parameters for debugging purposes

% Parse inputs

if nargin < 3,
    nhidden = 1;
end

ninput = size(X,1);
noutput = size(Y,1);
nsamples = size(X,2);

if nhidden > 1,

    G = cell(1,nhidden);

    % Run nhidden times sequentially

    R = Y;
    A = zeros(noutput,nhidden);
    B = zeros(ninput,nhidden);
    C = zeros(1,nhidden);
    for i=1:nhidden,
        [A(:,i),B(:,i),C(i),G{i}] = ospls(X,R,1,opts);
        if i < nhidden,    % deflate
            Z = B(:,i)'*X + C(i); % hidden activations
            R = R - A(:,i)*Z; % Y activations after deflation
        end
        if opts.verbose,
            fprintf('done %d out of %d; proportion selected %g\n',i,nhidden,nnz(B(:,i))/ninput);
        end
    end
else

    % Run for a single hidden unit

    B = zeros(ninput,1);
    C = 0;
    iter = 0;
    maxiter = 100;
    tol = nhidden*noutput*(1e-10);

    % Initialize A to first principal component of Y

    optseig.disp = 0;
    if nsamples < noutput,
        [d1,d2] = eigs(Y'*Y,[],1,'LM',optseig);
        A = Y*d1;
        A = A/sqrt(A'*A);
    else
        [A,d2] = eigs(Y*Y',[],1,'LM',optseig);
    end

    glmnetopts = glmnetSet;

    g = ft_mv_glmnet('method',opts.method,'family','gaussian','validator',ft_mv_crossvalidator('nfolds',5,'metric','correlation'));

    FN = fieldnames(glmnetopts);
    for c=1:length(FN)
      g.(FN{c}) = glmnetopts.(FN{c});
    end
    FN = fieldnames(opts);
    for c=1:length(FN)
      if isfield(glmnetopts,FN{c})
        g.(FN{c}) = opts.(FN{c});
      end
    end

    Aold = A;
    while iter < maxiter,

        if opts.verbose > 1,
            fprintf('   now starting iteration %d\n',iter+1);
        end

        Z = A'*Y;    % reconstruct Z given A from output Y

        % Use glmnet elastic net code to fit reconstructed Z given input X
        try

          f = g.train(X',Z');

          B = f.weights(1:(end-1));
          C = f.weights(end);

        catch

          % we end up in the catch block when the number of included
          % variables pmax has been exceeded; this should not happen
          warning(lasterr);

          B = zeros(ninput,1);
          C = 0;

          % this implies that deflation has no effect anymore

        end

        Z = B'*X + C;   % reconstruct Z given B and C

        % Find optimal A under constraint A'*A = 1

        Syz = Y*Z'/nsamples;
        denom = sqrt(Syz'*Syz);
        if denom,
            A = Syz/denom;
        else
            A = Aold;
        end
        if any(isnan(A(:))),    % check...
            error('nans!!!\n');
        end

        if sumsqr(A - Aold) < tol,
            if opts.verbose > 1,
                fprintf('   done!!\n',iter+1);
            end
		    iter = maxiter;
        else
	        iter = iter + 1;
            Aold = A;
        end
    end

    G = f;

end

author: fieldtrip

repo: fieldtrip


text:
function y = naca4_symmetric ( t, c, x )

%*****************************************************************************80
%
%% NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eastman Jacobs, Kenneth Ward, Robert Pinkerton,
%    "The characteristics of 78 related airfoil sections from tests in
%    the variable-density wind tunnel",
%    NACA Report 460, 1933.
%
%  Parameters:
%
%    Input, real T, the maximum relative thickness.
%
%    Input, real C, the chord length.
%
%    Input, real X(*), points along the chord length.
%    0.0 <= X(*) <= C.
%
%    Output, real Y(*), for each value of X, the corresponding value of Y
%    so that (X,Y) is on the upper wing surface, and (X,-Y) is on the
%    lower wing surface.
%
  y = 5.0 * t * c * ( ...
    0.2969 * sqrt ( x / c ) ...
    + (((( ...
      - 0.1015 ) .* ( x / c ) ...
      + 0.2843 ) .* ( x / c ) ...
      - 0.3516 ) .* ( x / c ) ...
      - 0.1260 ) .* ( x / c ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout=bss_decomp_tvgain(varargin)

% decompose an estimated source into target/interference/noise/artefacts components, assuming the admissible distortion is a time-varying gain.
%
% Usage:
%
% [s_target,e_interf[,e_noise],e_artif]=bss_decomp_tvgain(se,index,S[,N],tvshape,tvstep)
%
% Input:
%   - se: row vector of length T containing the estimated source,
%   - index: points which component of S se has to be compared to,
%   - S: n x T matrix containing the original sources,
%   - N: m x T matrix containing the noise on the obseravtions (if any).
%   - tvshape : row vector of length V at most T containing the shape of the elementary
%     allowed time variations of the gain
%   - tvstep  : hop size (in number of samples) between two consecutive
%     variations of the gain
%
% Output:
%   - s_target: row vector of length T containing the target source(s)
%   contribution,
%   - e_interf: row vector of length T containing the interferences
%   contribution,
%   - e_noise: row vector of length T containing the noise contribution (if
%   any),
%   - e_artif: row vector of length T containing the artifacts
%   contribution.
%
% Developers:  - Cedric Fevotte (cf269@cam.ac.uk) - Emmanuel Vincent
% (vincent@ircam.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

switch nargin
    case 5
        [varargout{1},varargout{2},varargout{3}]=bss_decomp_tvfilt(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5},0);
    case 6
        [varargout{1},varargout{2},varargout{3},varargout{4}]=bss_decomp_filt(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5},varargin{6},0);
    otherwise
        disp('Wrong number of arguments.')
end

author: posenhuang

repo: deeplearningsourceseparation


text:
% Downloaded from MatLab FileExchange as part of the findjobj submission by
% Yair M Altman at
%   https://www.mathworks.com/matlabcentral/fileexchange/14317-findjobj-find-java-handles-of-matlab-graphic-objects
% Thus this falls under the same open license as described by Yair in
% license.txt of his download.
% Programmed and Copyright by Yair M. Altman: altmany(at)gmail.com
% $Revision: 1.1 $  $Date: 2020/02/07 07:00:36 $
function jControl = findjobj_fast(hControl, jContainer)
    try jControl = hControl.getTable; return, catch, end  % fast bail-out for old uitables
    try jControl = hControl.JavaFrame.getGUIDEView; return, catch, end  % bail-out for HG2 matlab.ui.container.Panel
    oldWarn = warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame');
    if nargin < 2 || isempty(jContainer)
        % Use a HG2 matlab.ui.container.Panel jContainer if the control's parent is a uipanel
        try
            hParent = get(hControl,'Parent');
        catch
            % Probably indicates an invalid/deleted/empty handle
            jControl = [];
            return
        end
        try jContainer = hParent.JavaFrame.getGUIDEView; catch, jContainer = []; end
    end
    if isempty(jContainer)
        hFig = ancestor(hControl,'figure');
        jf = get(hFig, 'JavaFrame');
        jContainer = jf.getFigurePanelContainer.getComponent(0);
    end
    warning(oldWarn);
    jControl = [];
    counter = 20;  % 2018-09-21 speedup (100 x 0.001 => 20 x 0.005) - Martin Lehmann suggestion on FEX 2016-06-07
    specialTooltipStr = '!@#$%^&*';
    try  % Fix for R2018b suggested by Eddie (FEX comment 2018-09-19)
        tooltipPropName = 'TooltipString';
        oldTooltip = get(hControl,tooltipPropName);
        set(hControl,tooltipPropName,specialTooltipStr);
    catch
        tooltipPropName = 'Tooltip';
        oldTooltip = get(hControl,tooltipPropName);
        set(hControl,tooltipPropName,specialTooltipStr);
    end
    while isempty(jControl) && counter>0
        counter = counter - 1;
        pause(0.005);
        jControl = findTooltipIn(jContainer, specialTooltipStr);
    end
    set(hControl,tooltipPropName,oldTooltip);
    try jControl.setToolTipText(oldTooltip); catch, end
    try jControl = jControl.getParent.getView.getParent.getParent; catch, end  % return JScrollPane if exists
end

function jControl = findTooltipIn(jContainer, specialTooltipStr)
    try
        jControl = [];  % Fix suggested by H. Koch 11/4/2017
        tooltipStr = jContainer.getToolTipText;
        %if strcmp(char(tooltipStr),specialTooltipStr)
        if ~isempty(tooltipStr) && tooltipStr.startsWith(specialTooltipStr)  % a bit faster
            jControl = jContainer;
        else
            for idx = 1 : jContainer.getComponentCount
                jControl = findTooltipIn(jContainer.getComponent(idx-1), specialTooltipStr);
                if ~isempty(jControl), return; end
            end
        end
    catch
        % ignore
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function acc = ContinuousAccuracy(obj, pattern, predrange, unit)
% Calculate forced choice accuracy for unit increases in continuous
% predictions.  Requires units to be ranked ordered from 1:end.  Will work
% for Gianaros or Pain datasets.  Accuracies are not penalized for missing
% cases.  Probably best to run this on single subjects and then aggregate
% accuracies across subjects.
%
% :Usage:
% ::
%
%     acc = ContinuousAccuracy(obj, pattern, unit)
%
% :Inputs:
%
%   **obj:**
%        fmri_data() object with data stacked by
%        increasing levels of prediction.  Make sure
%        obj.Y includes the training labels
%
%   **pattern:**
%        fmri_data() object with weight pattern
%
%   **predrange:**
%        specify the range of predictions (e.g., 1:5)
%
%   **unit:**
%        specify the unit increase in prediction (e.g., 1 or 2)
%
% :Outputs:
%
%   **acc:**
%        accuracy of prediction for specified units
%
% :Examples:
% ::
%
%    acc = ContinuousAccuracy(dat, pine, 1:5, 1)
%
% ..
%    Original version: Copyright Luke Chang 12/2013
% ..

pexp = apply_mask(obj, pattern, 'pattern_expression', 'ignore_missing'); %Calculate pattern expression

%Create pairwise matrix and find positive values for lower triangle
%-need to force it to be 5 and fill in missing with NaNs
for i = predrange
    for j = predrange
        if ~any(obj.Y==i) || ~any(obj.Y==j)
            a(i,j) = nan;
        else
            a(i,j) = pexp(obj.Y == i) - pexp(obj.Y == j);
        end
    end
end
lt = double(tril(a) > 0);
lt(isnan(a)) = nan;

%remove by a factor of unit
i = 1;
while i < max(predrange)
    lt(i:i+unit-1,i) = 0;
    i = i + 1;
end

n = length(predrange)-unit;
miss = tril(isnan(lt)); %calculate lower triangle accuracy
miss(logical(eye(length(predrange))))=0; %remove any nans on diagonal
total = (n*(n+1)/2) - sum(sum(miss));  %subtract out any nans from total
acc = sum(nansum(lt))/total;


text:
function ind = contributingVoxels(pos)
%CONTRIBUTINGVOXELS  find indices of voxels not on image border
%
%   usage :
%   IND = contributingVoxels(POS)
%   where POS is a number between 1 and 27 indicating position of tile with
%   respect to image borders
%   and IND is an array of indices, giving only voxels which are not
%   located on the border of image.
%

% get position wrt edges of image
% if pk==1 -> config is on lower bounds in direction k
% if pk==2 -> config is in the middle of the 2 bounds
% if pk==3 -> config is on greater bound in direction k
pz = floor((pos-1)/9)+1;
pos2 = pos-9*(pz-1);
py = floor((pos2-1)/3)+1;
px = pos2-3*(py-1);

% Flags are set to one for contributions which should be computed
% If configuration is on the edge of image, the contributions of some
% vertices are not computed
flags = ones(1,8);
if px==1; flags([1 3 5 7]) = 0; end
if px==3; flags([2 4 6 8]) = 0; end
if py==1; flags([1 2 5 6]) = 0; end
if py==3; flags([3 4 7 8]) = 0; end
if pz==1; flags([1 2 3 4]) = 0; end
if pz==3; flags([5 6 7 8]) = 0; end

% get indices of pixels whose contribution will be computed
ind = find(flags);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [view ROI] = createROIFromRange(view, range, varargin)
% createROIFromRange
%   Generates an ROI from a range of values given by the user.  Optional
%   ability to generate histogram with feedback re: data within range.
%
% Usage:
%   [view ROI] = createROIFromRange(view, range, varargin)
%
%   View:
%       Volume/gray struct.
%
%   Range:
%       String encoding desired values.  If empty, will prompt for string.
%           FORMAT: OPAREN MINNUMBER , MAXNUMBER CPAREN
%               OPAREN - ( for exlusive, [ for inclusive
%               MINNUMBER - int/double of min, empty for no min
%               MAXNUMBER - int/double of max, empty for no max
%               CPAREN - ) for exclusive, ] for inclusive
%
%           Examples:
%               (3, 4] - from 3 exclusive to 4 inclusive
%               [2,) - from 2 inclusive to max
%               [,-1] - from min to -1 inclusive
%
%               Note:
%                   inclusive/exclusive does not effect empty nums (always
%                   treated as 'inclusive' on the min/max value)
%
%   Varargin Options:
%       'showhist' - Display histogram of values in range (true/false)
%           * DEFAULT false
%       'addflag' - Add ROI to view
%           * DEFAULT true
%       'roiname' - String for new name of ROI
%           * DEFAULT named by input range string
%
% Output:
%   View:
%       Volume/gray struct with ROI added
%
%   ROI:
%       ROI struct
%
% [renobowen@gmail.com 2010]
%
    if (isempty(view.map))
        ROI = [];
        error('No parameter map loaded.');
        return;
    end

    if (isempty(range))
        range = inputdlg('Enter range (see help CreateROIFromRange for instructions):',...
            'Enter Range', 1, {'[,]'});
        range = range{1};
    end

    [isOpenInclusive num1 num2 isCloseInclusive] = ParseRange(range);
    if (isempty(isOpenInclusive))
        ROI = [];
        return;
    end

    curScan = view.curScan;

    % Defaults
    showHist = false;
    addFlag = true;
    roiname = [];
    for i = 1:2:length(varargin)
        switch lower(varargin{i})
            case {'showhist'}
                showHist = varargin{i + 1};
            case {'addflag'}
                addFlag = varargin{i + 1};
            case {'roiname'}
                roiname = varargin{i + 1};
            otherwise
                fprintf(1, 'Unrecognized option: ''%s''', varargin{i});
        end
    end

    if (isempty(num1)), num1 = min(view.map{curScan}); isOpenInclusive = 1; end
    if (isempty(num2)), num2 = max(view.map{curScan}); isCloseInclusive = 1; end

    if (isOpenInclusive)
        openInds = find(view.map{curScan} >= num1);
    else
        openInds = find(view.map{curScan} > num1);
    end

    if (isCloseInclusive)
        closeInds = find(view.map{curScan} <= num2);
    else
        closeInds = find(view.map{curScan} < num2);
    end

    inds = intersect(openInds, closeInds);

    coords = view.coords(:, inds);

    ROI.coords = coords;
    ROI.color = 'w';
    if isempty(roiname)
        ROI.name = sprintf('range%s',range);
    else
        ROI.name = roiname;
    end
    ROI.viewType = view.viewType;

    if (addFlag), view = addROI(view, ROI); view.ui.showROIs = -1; end
    if (showHist), ShowRangeHistogram(view.map{curScan}(inds)); end
end

function [isOpenInclusive num1 num2 isCloseInclusive] = ParseRange(range)
% [isOpenInclusive num1 num2 isCloseInclusive] = parseRange(range)
%   Parse a range statement using regular expressions, returning properly
%   digested tokens to be used in selecting the values.
%
    isCloseInclusive = 0;
    num1 = 0; num2 = 0;

    % Regular expressions to parse range
    number = '((-)?(\d)+(\.(\d)+)?)';
    optNumber = [number '?'];
    openParen = '(\(|\[)';
    closeParen = '(\)|\])';
    comma = '(\,)';
    whitespace = '((\s)*)';
    rangeExp = ['(?<openParen>' openParen ')' whitespace '(?<num1>' optNumber ')' ...
        whitespace comma whitespace '(?<num2>' optNumber ')' whitespace '(?<closeParen>' closeParen ')'];
    result = regexp(range, rangeExp, 'names');

    if (isempty(result))
        isOpenInclusive = [];
        error('Malformed range input, please see the following help text for instructions.');
        return;
    end

    isOpenInclusive = strcmp(result.openParen, '[');
    isCloseInclusive = strcmp(result.closeParen, ']');
    num1 = str2num(result.num1);
    num2 = str2num(result.num2);

end

function ShowRangeHistogram(data)
    h = figure;
    hist(data);
    axes = get(h, 'CurrentAxes');

    yLim = get(axes, 'YLim');
    yRange = yLim(2) - yLim(1);
    yPos = [(yLim(2) - yRange * .05), (yLim(2) - yRange * .1)];

    xLim = get(axes, 'XLim');
    xRange = xLim(2) - xLim(1);
    xPos = ones(1, 2) * (xLim(1) + xRange * .05);

    text(xPos(1), yPos(1), ['min: ' num2str(min(data))]);
    text(xPos(2), yPos(2), ['max: ' num2str(max(data))]);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function msk=get_lines_executed(obj)
% get a mask indicating which lines have been executed
%
% msk=get_lines_executed(obj)
%
% Input:
%   obj                 MOcovMFile instance
%
% Output:
%   msk                 Nx1 logical mask, with msk(k)==true indicating
%                       that the k-th line has been executed at least once

    msk=get_lines_executed_count(obj)>0;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compare different clustering methods of event sequences
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

options.N = 100; % the number of sequences
options.Nmax = 100; % the maximum number of events per sequence
options.Tmax = 50; % the maximum size of time window
options.tstep = 0.1;
options.dt = 0.1;
options.M = 250;
options.GenerationNum = 10;
D = 3; % the dimension of Hawkes processes
K = 2; % the number of clusters
nTest = 5;
nSeg = 5;
nNum = options.N/nSeg;

disp('Approximate simulation of Hawkes processes via branching process')
% the first cluster contains the event sequences generated by a Hawkes process with exponential kernel
disp('Simple exponential kernel')
para1.kernel = 'exp';
para1.landmark = 0;
L = length(para1.landmark);
para1.mu = rand(D,1)/D;
para1.A = zeros(D, D, L);
for l = 1:L
    para1.A(:,:,l) = (0.7^l)*rand(D);
end
para1.A = 0.5 * para1.A./max(abs(eig(para1.A)));
para1.w = 0.5;
Seqs1 = Simulation_Branch_HP(para1, options);

% the second cluster contains the event sequences generated by a Hawkes process with gaussian kernel
disp('Complicated gaussian kernel')
para2.kernel = 'gauss';
para2.landmark = 0:3:12;
L = length(para2.landmark);
para2.mu = rand(D,1)/D;
para2.A = zeros(D, D, L);
for l = 1:L
    para2.A(:,:,l) = (0.9^l)*rand(D);
end
para2.A = 0.25 * para2.A./max(abs(eig(sum(para2.A,3))));
para2.A = reshape(para2.A, [D, L, D]);
para2.w = 1;
Seqs2 = Simulation_Branch_HP(para2, options);

SeqsMix = [Seqs1,Seqs2];

% ground truth: the similarity matrix of event sequences.
% GT(i,j)=1 if the sequence i and the sequence j belong to a same cluster
GT = blkdiag(ones(options.N), ones(options.N));

% Learning a mixture model of Hawkes processes
% initialize
model = Initialization_Cluster_Basis(SeqsMix, 2);
alg.outer = 8;
alg.rho = 0.1;
alg.inner = 5;
alg.thres = 1e-5;
alg.Tmax = [];
model = Learning_Cluster_Basis(SeqsMix, model, alg);
[~, labels1] = max(model.R, [],2);
% estimate similarity matrix
Est1 = zeros(length(labels1));
for i = 1:length(labels1)
    for j = 1:length(labels1)
        if labels1(i)==labels1(j)
            Est1(i,j) = 1;
        end
    end
end

% Learning a distance metric of marked point processes
configure.M = [options.Tmax; D];
configure.id = {1,2,[1,2]};
configure.tau = 1;
configure.W = 5;
configure.epoch = 15;
configure.lr = 1e-4;
[configure, obj] = Estimate_Weight(configure, SeqsMix);
Dis = zeros(length(SeqsMix));
tic
for n = 1:length(SeqsMix)-1
    Xn = [SeqsMix(n).Time; SeqsMix(n).Mark];
    for m = n+1:length(SeqsMix)
        Xm = [SeqsMix(m).Time; SeqsMix(m).Mark];
        [Dis(n,m), ~] = DistanceSum_MPP(Xn, Xm, configure);
    end
    toc
end
Dis = Dis + Dis';
% estimate similarity matrix
SimilarMat = exp(-Dis.^2./(2*var(Dis(:))));

% visualize clustering results
figure
subplot(131)
imshow(GT,[])
title('Ground truth of clusters')
subplot(132)
imshow(Est1,[])
title('Mixture of Hawkes processess')
subplot(133)
imshow(SimilarMat,[])
title('Distance metrics of MPP')

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function [L,Lc,Lli] = toFramePlucker(C,Li)

% TOFRAMEPLUCKER  Transform plucker line to a given frame.
%   L = TOFRAMEPLUCKER(C,Li) expresses in frame C=[t;q] the line Li
%   originally expressed in global frame.
%
%   The formula for transformation is
%
%       L = [ Rt       [it]_x*Rt ]
%           [ zeros(3)        Rt ] * Li
%
%   where
%       Rt    = q2R(iq)
%      [it]_x = hat(it)
%       iq    = q2qc(q)
%       it    = t2it(t,q)
%       t     = C(1:3)
%       q     = C(4:7)
%
%   [L,Lc,Lli] = TOFRAMEPLUCKER(...) returns the Jacobians wrt C
%   and Li.
%
%   See also FROMFRAMEPLUCKER, TXP, HAT, CROSS, Q2R, Q2QC, T2IT.

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

% extract input vectors ai and bi
ai = Li(1:3);
bi = Li(4:6);

t = C.t;
q = C.q;

if nargout == 1

    % get it and Rt of inverse transform
    it = t2it(t,q);
    Rt = C.Rt;

    % compute a and b
    b  = Rt*bi;
    a  = Rt*ai + hat(it)*b;

    % Build output
    L = [a;b];

else

    % get it and iq of inverse transform
    [it,ITt,ITq]     = t2it(t,q);

    % compute a and b
    [b,Bq,Bbi]       = Rtp(q,bi);
    [txb,TXBit,TXBb] = crossJ(it,b);
    [Ra,RAq,RAai]    = Rtp(q,ai);

    a = Ra + txb;

    % Jacobians from the chain rule
    At  = TXBit*ITt;
    Aq  = RAq + TXBit*ITq + TXBb*Bq;
    Aai = RAai;
    Abi = TXBb*Bbi;

    % Build outputs
    L   = [a;b];

    Lc  = [At Aq;zeros(3) Bq];
    Lli = [Aai Abi;zeros(3) Bbi];

end

return

%%

syms a b c d x y z real
syms L1 L2 L3 L4 L5 L6 real

q = [a;b;c;d];
t = [x;y;z];
C = [t;q];
Li = [L1;L2;L3;L4;L5;L6];

[L,Lc,Lli] = toFramePlucker(C,Li);

simplify(Lc  - jacobian(L,C))
simplify(Lli - jacobian(L,Li))

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [F10, F01] = retrieveQ1001(level, o, H, L)
%------------------------------------------------------------------------------
%
% This function extracts gridfunction F10 and F01 from H.
% F10 is of colour '10' and F01 is of colour '01' such that together they
% constitute a quincunx gridfunction.
%
%                H   H   H   H
%              V   V   V   V   V       V ~ F10
%                H   H   H   H         H ~ F01
%              V   V   V   V   V
%                H   H   H   H
%              V   V   V   V   V
%                H   H   H   H
%              V   V   V   V   V
%                H   H   H   H
%       (domain of quincunx gridfunction)
%
% F10 & F01 are supposed to be uniquely determined by the integer L (level)
% and character o describing its type.
% Function retrieveQ1001 is the counterpart of storeQ1001.
% This function is a two-dimensional lifting scheme utility.
%
% F10 = 2D gridfunction of coefficients of colour '10' extracted from H.
%
% F01 = 2D gridfunction of coefficients of colour '01' extracted from H.
%
% level = integer designated as the level of both F10 and F01.
%
% o = character, should be either 'a' or 'd',
%     describing the type of both F10 and F01:
%     'a' relates to approximation (coefficients) and
%     'd' relates to detail (coefficients)
%
% L = 2D integer array of bookkeeping, see function storeR.
%
% H = 1D array that functions as storage (heap). The coefficients of both F10
%     and F01 are extracted from H as a result of calling retrieveQ1001.
%
% See also: retrieveQ1001, retrieveQ, retrieveQ0110.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: October 17, 2002.
% (c) 1999-2002 Stichting CWI, Amsterdam.
%------------------------------------------------------------------------------
F10 = retrieveQ(level, '10', o, H, L);
F01 = retrieveQ(level, '01', o, H, L);
if isempty(F10) || isempty(F01)
  disp([' retrieveQ1001 - ERROR at type ' o ' with level ' num2str(level)]);
  error(' retrieveQ1001 - empty result ')
else
  [n10, m10]=size(F10);
  [n01, m01]=size(F01);
  if m01 > m10
    disp([' retrieveQ1001 - ERROR at type ' o ' with level ' num2str(level)]);
    error(' retrieveQ1001 - dimensions do not match for m01 > m10 ')
  end
  if n10 > n01
    disp([' retrieveQ1001 - ERROR at type ' o ' with level ' num2str(level)]);
    error(' retrieveQ1001 - dimensions do not match for n10 > n01 ')
  end
  if m10 > m01+1
    disp([' retrieveQ1001 - ERROR at type ' o ' with level ' num2str(level)]);
    error(' retrieveQ1001 - dimensions do not match for m10 > m01+1 ')
  end
  if n01 > n10+1
    disp([' retrieveQ1001 - ERROR at type ' o ' with level ' num2str(level)]);
    error(' retrieveQ1001 - dimensions do not match for n01 > n10+1 ')
  end
%
% m01 <= m10 <= m01+1 is satisfied
% n10 <= n01 <= n10+1 is satisfied
%
end
%------------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function [ftAllNew,transMdl] = ftTrans_mida(ftAll,domainFtAll,target,...
	maLabeled,param)
%Maximum Independence Domain Adaptation (MIDA)
%
%	A feaure-level transfer learning (domain adaptation) algorithm which
% augments the features then learns a domain-invariant subspace.
% Application scope:
%	+ two or multiple discrete domains
%	+ continuous distributional change
%	+ labeled or unlabeled or partially labeled source domain
%	+ labeled or unlabeled or partially labeled target domain
%	+ label type: classification or regression
%	There is no need to distinguish which domain a sample is from. The domain
% information is contained in the variable domainFtAll, which contains the
% domain features of all samples.
%	Domain features indicate the background of samples. For example, if
% the samples are from md discrete domains, then domainFtAll can be a
% n-by-md matrix, domainFtAll(i,j)=1 if sample i is from the j'th domain, 0
% else. If the samples have a time order and their distribution changes
% continuously along with time, then domainFtAll can be n-by-1 and
% domainFtAll(i) is the "time" of sample i. Multiple background info can be
% integrated into domainFtAll in the similar way. See the ref for details.
%
% ftAll:	All samples in all domains. n-by-m matrix, n is the number of
%	samples, m is the dimension of features.
% domainFtAll:	Domain features of all samples. n-by-md matrix, md is the
%	dimension of domain features.
% target:	When some samples in any domain are labeled, their labels can
%	be provided in this variable to enhance the discriminative power of the
%	the learned features. ntr-by-1 matrix, ntr is the number of labeled
%	samples. Classification problems should use class indices as labels,
%	i.e. target(i)=j if the i'th labeled sample is from the j'th class.
% maLabeled:	Mask for the labeled samples. n-by-1 matrix,
%	maLabeled(i)=true if sample i is labeled, false else.
%
% param: Struct of hyper-parameters, please see the first cell of this
%	program ("default parameters") for details. You can set parameter p to
%	x by setting param.p = x. For parameters that are not set, default
%	values will be used.

% ftAllNew:	All samples in the learned subspace.
% transMdl:	A struct containing the model, transMdl.W is the projection
%	matrix.

% ref: Ke Yan, Lu Kou, and David Zhang, "Domain Adaptation via Maximum
%	Independence of Domain Features," http://arxiv.org/abs/1603.04535
% Copyright 2016 Ke YAN, Tsinghua Univ. http://yanke23.com , xjed09@gmail.com

%% default parameters
isRegress = 0; % 0 for classification problem, 1 for regression
kerName = 'lin'; % kernel name, see the next cell ("kernels")
kerSigma = 10; % kernel parameter, see the next cell ("kernels")
bSmida = true; % 0 for MIDA if no label information is considered for all
	% samples, 1 for semisupervised MIDA (SMIDA) if some labels are
	% considered.
doSample = false; % when there are too many unlabeled data, eigenvalue
	% decomposition can be slow. Setting this variable to true makes the
	% code to sample some unlabeled data.
nSmpRatio = 1; % if doSample=true, the number of unlabeled data to sample
	% will be ceil(ntr*nSmpRatio)
mu = 1; % the weight of the variance term, see the ref
m = 30; % the dimension of the subspace
gamma = .1; % the weight of the supervised term, only useful in SMIDA, see
	% the ref
ftAugType = 1; % feature augmentation, 0: no aug; 1: aug with domainFt;
	% 2: frustratingly easy aug, only for discrete domains, see the ref

defParam % set user-defined hyper-parameters

%% kernels
nm = @(X,p)repmat(sum(X.^2,2),1,p);
linKer = @(X1,X2)X1*X2';
rbfKer = @(X1,X2)exp(-(nm(X1,size(X2,1))+nm(X2,size(X1,1))'-2*X1*X2')/2/kerSigma^2);
lapKer = @(X1,X2)exp(-pdist2(X1,X2)/kerSigma);
polyKer = @(X1,X2)(1+kerSigma*X1*X2').^2;
if strcmpi(kerName,'lin'), kerFun = linKer; % linear kernel
elseif strcmpi(kerName,'poly'), kerFun = polyKer; % polynomial kernel
elseif strcmpi(kerName,'rbf'), kerFun = rbfKer;
elseif strcmpi(kerName,'lap'), kerFun = lapKer; % Laplacian kernel
else error('unknown kernel'); end

%% sort samples
ntr = nnz(maLabeled);
nAll = size(ftAll,1);
Df = double(domainFtAll);
% Df = zscore(domainFtAll); % bad
if doSample && nSmpRatio<inf % sample some unlabeled data
	rng(0)
	nUnlabeledUsed = min(nAll-ntr,ceil(ntr*nSmpRatio));
	id = randperm(nAll-ntr,nUnlabeledUsed); % rand sel
	idTest = find(~maLabeled);
	ftUsed = [ftAll(maLabeled,:);ftAll(idTest(id),:)];
	dfUsed = [Df(maLabeled,:);Df(idTest(id),:)];
	nUsed = ntr+nUnlabeledUsed;
else
	ftUsed = [ftAll(maLabeled,:);ftAll(~maLabeled,:)];
	dfUsed = [Df(maLabeled,:);Df(~maLabeled,:)];
	nUsed = nAll;
end
Kd = linKer(dfUsed,dfUsed);

%% feature augmentation
if ftAugType==0
	ftUsedAug = ftUsed;
elseif ftAugType==1
	ftUsedAug = [ftUsed,dfUsed];
elseif ftAugType==2
	nFt = size(ftUsed,2);
	nDomain = size(dfUsed,2);
	ftUsedAug = zeros(nUsed,nFt*(nDomain+1));
	ftUsedAug(:,1:nFt) = ftUsed;
	for p = 1:nDomain
		ftUsedAug(dfUsed(:,p)==1,nFt*p+1:nFt*(p+1)) = ftUsed(dfUsed(:,p)==1,:);
	end
end

%% compute
Kx = kerFun(ftUsedAug,ftUsedAug);
H = eye(nUsed)-ones(nUsed)/(nUsed);

if ~bSmida % MIDA

	A = Kx*(mu*H-H*Kd*H)*Kx;

else % SMIDA

	target = double(target);
	if isRegress==0 && all(target==floor(target)) % classification
		Ky = repmat(target,1,ntr)==repmat(target,1,ntr)';
% 		Ky = (Kyy-.5)*2; % no need
	else % regression
		target = zscore(target);
		Ky = target*target'; % linear ker
	end

	Kyy_tilde = zeros(nUsed);
	Kyy_tilde(1:ntr,1:ntr) = Ky;

	A = Kx*(H*(-Kd+gamma*Kyy_tilde)*H+mu*H)*Kx;
end

A = (A+A')/2; % to compensate float num error
[V,D] = eig(A);
[D,I] = sort(diag(real(D)),'descend');
transMdl.W = real(V(:,I(1:m)));

%% project the samples
if ftAugType==0
	ftAllAug = ftAll;
elseif ftAugType==1
	ftAllAug = [ftAll,Df];
elseif ftAugType==2
	ftAllAug = zeros(nUsed,nFt*(nDomain+1));
	ftAllAug(:,1:nFt) = ftAll;
	for p = 1:nDomain
		ftAllAug(Df(:,p)==1,nFt*p+1:nFt*(p+1)) = ftAll(Df(:,p)==1,:);
	end
end

KxAll = kerFun(ftAllAug,ftUsedAug);
ftAllNew = KxAll*transMdl.W;

end

author: viggin

repo: domain-adaptation-toolbox


text:
function triangulation_test25 ( )

%*****************************************************************************80
%
%% TEST25 tests TRIANGULATION_ORDER6_BOUNDARY_NODE, TRIANGULATION_ORDER6_PLOT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  file_name = 'triangulation_order6_plot.eps';
  dim_num = 2;
  node_show = 2;
  triangle_show = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST25\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER6_BOUNDARY_COUNT counts the boundary\n' );
  fprintf ( 1, '    edges in an arbitrary triangulation.\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER6_PLOT plots the triangulation.\n' );

  [ node_num, triangle_num, hole_num ] = ...
    triangulation_order6_example1_size ( );

  [ node_xy, triangle_node, triangle_neighbor ] = ...
    triangulation_order6_example1 ( node_num, triangle_num );
%
%  Make the plot.
%
  triangulation_order6_plot ( file_name, node_num, node_xy, triangle_num, ...
    triangle_node, node_show, triangle_show );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  An Encapsulated PostScript image of this\n' );
  fprintf ( 1, '  triangulation is in "%s".\n', file_name );

  node_boundary = triangulation_order6_boundary_node ( node_num, triangle_num, ...
    triangle_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node  BN?\n' );
  fprintf ( 1, '\n' );

  for i = 1 : node_num
    fprintf ( 1, '  %6d  %1d\n', i, node_boundary(i) );
  end

  return
end


text:
function vec=getPermIndices(rank,n,dim)
%%GETPERMINDICES Get a vector of indices that can be used to rearrange the
%                subvectors in a vector consisting of n subvectors of
%                dimensionality dim according to the permutation of the
%                given rank in lexicographic ordering.
%
%INPUTS: rank The order of the desired permutation of 0:(n-1) (or 1:n) in
%             lexicographic order. Note  that 0<=rank<=(n!-1).
%           n The number of stacked vectors that are to be permuted.
%         dim The dimensionality of the n subvectors being permutated.
%
%OUTPUTS: vec A vector of indices that can be used to reorder the
%             subvectors of a vector of n subvectors of dimensionality dim.
%
%Given an (n*dim)X1 vector v of n stacked subvectors of dimension dim,
%v(vec) rearranges the subvectors in v according to the given permutation.
%This rearrangement can also be performed in a less computationally
%efficient manner using a permutation matrix generated by the function
%rank2PermMatrix.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    perm=unrankPermutation(rank,n);
    vec=zeros(dim*n,1);

    for curPerm=1:n
        baseVec=(curPerm-1)*dim;
        basePerm=(perm(curPerm)-1)*dim;

        for off=1:dim
           vec(baseVec+off)=basePerm+off;
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function ClickToSelectOffense(patch, ~, game, player)

N = numel(game.board);

%% RESET ALL PATCHES TO NEUTRAL
for i = 1:N
    occupant = get(game.board(i).Patch, 'UserData');
    if occupant(1) == player
        set(game.board(i).Patch, 'UserData', [player, 0], ...
                                 'FaceColor', game.cmap(player, :));
        set(game.board(i).Text, 'Color', [0 0 0]);
    end
end

%% SET CURRRENT PATCH TO OFFENSE
idx = str2double(get(patch, 'Tag'));
set(game.board(idx).Patch, 'UserData', [player 1], ...
                           'FaceColor', game.cmap(7,:));
set(game.board(idx).Text, 'Color', game.cmap(8,:));


text:
% biosig2eeglabevent() - convert biosig events to EEGLAB event structure
%
% Usage:
%   >> eeglabevent = biosig2eeglabevent( biosigevent, interval )
%
% Inputs:
%   biosigevent    - BioSig event structure
%   interval       - Period to extract events for, in frames.
%                    Default [] is all.
%
% Outputs:
%   eeglabevent    - EEGLAB event structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2006-

% Copyright (C) 13 2006- Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function event = biosig2eeglabevent(EVENT, interval)

if nargin < 2
    interval = [];
end;

event = [];
disp('Importing data events...');

% If the interval variable is empty, import all events.
if isempty(interval)
    if isfield(EVENT, 'Teeg')
        event = EVENT.Teeg;
    end
    if isfield(EVENT, 'TYP')
        for index = 1:length( EVENT.TYP )
            event(index).type = EVENT.TYP(index);
        end
    end
    if isfield(EVENT, 'POS')
        for index = 1:length( EVENT.POS )
            event(index).latency = EVENT.POS(index);
        end
    end
    if isfield(EVENT, 'DUR')
        if any( [ EVENT.DUR ] )
            for index = 1:length( EVENT.DUR )
                event(index).duration = EVENT.DUR(index);
            end
        end
    end
    if isfield(EVENT, 'CHN')
        if any( [ EVENT.CHN ] )
            for index = 1:length( EVENT.CHN )
                event(index).chanindex = EVENT.CHN(index);
            end
        end
    end
% If a subinterval was specified, select only events that fall in that range, and
% edit duration field if it exceeds that range.
elseif isfield(EVENT,'POS')
    count = 1;
    for index = 1:length(EVENT.POS)
        pos_tmp = EVENT.POS(index) - interval(1) + 1;
        if pos_tmp > 0 & EVENT.POS(index) <= interval(2)
            event(count).latency = pos_tmp;
            if isfield(EVENT, 'TYP'), event(count).type = EVENT.TYP(index); end
            if isfield(EVENT, 'CHN'), event(count).chanindex = EVENT.CHN(index); end
            if isfield(EVENT, 'DUR')
                event(count).duration = min(EVENT.DUR(index), interval(2) - EVENT.POS(index));
            end
            count = count + 1;
        end
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function R = euler2rotMat(phi, theta, psi)
    R(1,1,:) = cos(psi).*cos(theta);
    R(1,2,:) = -sin(psi).*cos(phi) + cos(psi).*sin(theta).*sin(phi);
    R(1,3,:) = sin(psi).*sin(phi) + cos(psi).*sin(theta).*cos(phi);

    R(2,1,:) = sin(psi).*cos(theta);
    R(2,2,:) = cos(psi).*cos(phi) + sin(psi).*sin(theta).*sin(phi);
    R(2,3,:) = -cos(psi).*sin(phi) + sin(psi).*sin(theta).*cos(phi);

    R(3,1,:) = -sin(theta);
    R(3,2,:) = cos(theta).*sin(phi);
    R(3,3,:) = cos(theta).*cos(phi);
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Calibrating kl, ku, mu and nu for a sabr model for Kienitz extrapolation
% the prices are calculated from known sabr parameters

f = 0.03; t =1;                         % forward time
a = 0.25; b = 0.5; r = -.5; n = 0.2;    % sabr parameters

k = 0.0001:0.0001:1;                         % strike range

call = sprice(a, b, r, n, f, k, t,1);   % call prices (standard sabr)
put = sprice(a, b, r, n, f, k, t,0);    % put prices (standard sabr)
put(1) = 0;
call(1) = f;                            % assures forward is matched

Nparam = 4;                             % number of calibrated params

% objective function
of = @(x) of_sabr(a, b, r, n, f, k, t, x(3), x(4), x(1), x(2), call, put);

x0 = [.25*f; 25.5*f; 1.8; 2.4];                   % starting values                  %

A = zeros(Nparam,Nparam); bc= zeros(Nparam,1);
Aeq = A; beq = bc;
lb = [.25*f; f; 1; 3];                        % lower bound
ub = [f; 30*f; 1.5; 5];                     % upper bound
y = fmincon(of,x0,A,bc,Aeq,beq,lb,ub);  % optimization

% verification of results
xval = 0:0.001:.25;                     % x-values
[cl, bl, al, cu, bu, au] = psabr_param_3(a, b, r, n, f, t,y(3),y(4),y(1),y(2));
yval = psabr_5(a, b, r, n, f, xval, t, y(1), y(2), ...
           y(3), cl, bl,al, y(4), cu,bu,au);% y-values calculated

plot(xval,yval);                        % plot the results
Factor = 1000000;                       % used for plotting

yval_call = sprice_5(a, b, r, n, f, xval, t, y(1), y(2), ...
    y(3), cl, bl,al, y(4), cu,bu,au, 1);    % SABR Call prices
figure; hold on;
    plot(xval,Factor*yval_call,'r');
    plot(k,Factor*call,'g');
hold off;

yval_put = sprice_5(a, b, r, n, f, xval, t, y(1), y(2), ...
    y(3), cl, bl,al, y(4), cu,bu,au, 0);    % SABR Put prices
figure; hold on;
    plot(xval,Factor*yval_put,'r');
    plot(k,Factor*put,'g');
hold off;

fval = sprice_5(a, b, r, n, f, 0, t, y(1), y(2), ...
    y(3), cl, bl,al, y(4), cu,bu,au, 1);    % calculate forward value
y                                       % calibrated values
f - fval                                % display difference to forward


text:
function[b]=allall(x)
%ALLALL  ALLALL(X)=ALL(X(:))
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004 J.M. Lilly --- type 'help jlab_license' for details
b=all(x(:));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function varargout = VScope_start(varargin)
% VSCOPE_START M-file for VScope_start.fig
%      VSCOPE_START, by itself, creates a new VSCOPE_START or raises the existing
%      singleton*.
%
%      H = VSCOPE_START returns the handle to a new VSCOPE_START or the handle to
%      the existing singleton*.
%
%      VSCOPE_START('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in VSCOPE_START.M with the given input arguments.
%
%      VSCOPE_START('Property','Value',...) creates a new VSCOPE_START or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before VScope_start_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to VScope_start_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help VScope_start

% Last Modified by GUIDE v2.5 19-Jan-2008 23:51:34

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @VScope_start_OpeningFcn, ...
                   'gui_OutputFcn',  @VScope_start_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before VScope_start is made visible.
function VScope_start_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to VScope_start (see VARARGIN)

% Choose default command line output for VScope_start
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes VScope_start wait for user response (see UIRESUME)
% uiwait(handles.figure1);

im=imread('osc.jpg');
axis(handles.axes1);
imshow(im);
axis off;

% --- Outputs from this function are returned to the command line.
function varargout = VScope_start_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.text3,'string','Initializing.....Please Wait');
pause(2)
set(handles.text3,'string',' ');
query
close('VScope_start');

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

close('VScope_start');


text:
function out = DN_Moments(y,theMom)
% DN_Moments    A moment of the distribution of the input time series.
%
% Normalizes by the standard deviation
% Uses the moment function from Matlab's Statistics Toolbox
%
%---INPUTS:
% y, the input data vector
% theMom, the moment to calculate (a scalar)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

out = moment(y,theMom) / std(y); % normalized

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function m = p03_m ( )

%*****************************************************************************80
%
%% P03_M returns the spatial dimension for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Output, integer M, the spatial dimension.
%
  m = 4;

  return
end


text:
function [nlzs log_losses times_cache] = GaussVarApproxHessian_alg2(compute_loss,guessMean,guessPrec,gradFun,nrSteps,nSamples,stepSize,X,y,X_te,y_te,init_a1,init_a2)
%Code is adapted from https://github.com/TimSalimans/LinRegVB/blob/master/probit/GaussVarApproxHessian.m

nlzs=zeros(nrSteps+1,1);
log_losses=zeros(nrSteps+1,1);
times_cache=zeros(nrSteps+1,1);

% dimension
k=length(guessMean);
N = size(X,1);

% initial statistics for stochastic approximation
dm1 = init_a1*ones(N,1);
dm2 = init_a2*ones(N,1);
tm11 = X'*dm1;
tmp1 = bsxfun(@times,X,dm2);%diag(dm2)*X
tm22 = X'*tmp1;

gamma = diag(guessPrec);
sW = 1 ./ sqrt( gamma );
L = chol(eye(k)+sW*sW'.*tm22);
half_V = L'\(diag(sW)); %V=half_V'*half_V; %%V= ( tm2+diag(gamma) )^{-1}
m = half_V'*(half_V*tm11);%note that prior_mean = 0
V = half_V'*half_V;

P = V\eye(k);
z = m;
a = zeros(k,1);
% do stochastic approximation
for i=1:nrSteps
    tic;
    % cholesky factor of current inverse variance matrix
    cholP = chol(P);
    % sample
    xMean = cholP\(cholP'\a) + z;

    post_mean = xMean;
    post_cov = cholP\(cholP'\eye(k));
    tt1 = toc;

    if compute_loss==1
        post_dist.mean = post_mean;
        post_dist.covMat = post_cov;
        [pred, log_lik]=get_loss(i-1, post_dist, X, y, gamma, X_te, y_te);
        nlzs(i) = -log_lik;
        log_losses(i) = pred.log_loss;
    end

    tic;
    [g H b] =gradFun(xMean, cholP, nSamples);

    % update statistics
    P = (1-stepSize)*P - stepSize*H;
    a = (1-stepSize)*a + stepSize*g;
    z = (1-stepSize)*z + stepSize*b;
    tt2 = toc;
    times_cache(i+1) = tt1 + tt2;
end
times_cache = cumsum(times_cache);

if compute_loss==1
    approxPrec = P;
    cholP = chol(approxPrec);
    approxMean = cholP\(cholP'\a) + z;
    approxV = cholP\(cholP'\eye(k));

    i = nrSteps+1;
    post_dist.mean = approxMean;
    post_dist.covMat = approxV;
    [pred, log_lik]=get_loss(i-1, post_dist, X, y, gamma, X_te, y_te);

    nlzs(i) = -log_lik;
    log_losses(i) = pred.log_loss;
end

end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function [ K, columns ] = spark( Phi )
%SPARK Calculates the spark of a given matrix Phi
% Let us get the rank of matrix
R = rank(Phi);
% Number of rows and columns of Phi
[~,N] = size(Phi);
for K=1:R+1
    % Possible number of choices of K columns
    % from N columns
    if K > N
        % This happens when we have
        % an NxN full rank matrix
        break;
    end
    numChoices = nchoosek(N,K);
    if (numChoices > 160000)
        error('N=%d, K=%d too large!',N,K);
    end
    % All possible choices of K columns out of N
    choices = nchoosek(1:N, K);
    % We iterate over choices
    for c=1:numChoices
        choice = choices(c, :);
        % We choose K columns out of N columns
        Phik = Phi (:, choice);
        % We compute the rank of this matrix
        r = rank(Phik);
        % We check if the columns are linearly
        % dependent.
        if r < K
            % We have found a set of columns
            % which are linearly dependent
            columns = choice;
            return
        end
    end
end
% All columns up to rank K are
% linearly independent
K = R+1;
columns = 1:N;
end

author: indigits

repo: sparse-plex


text:
function [I,J,K,IND]=logic2subind(L)

[J,I,K]=meshgrid(1:1:size(L,2),1:1:size(L,1),1:1:size(L,3));
I=I(L); J=J(L); K=K(L);
IND = sub2ind(size(L),I,J,K);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function y = tapas_condhalluc_obs2_sim(r, infStates, p)
% Simulates responses according to the condhalluc_obs model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Get parameters
be = p(1);
nu = p(2);

% Prediction trajectory
mu1hat = infStates(:,1,1);

% Get true-positive rate corresponding to stimuli
tp = r.u(:,2);

% Update belief using precision-weighted prediction error
% with nu the generalized precision
x = mu1hat + 1/(1 + nu)*(tp - mu1hat);

% Apply the logistic sigmoid to the inferred beliefs
prob = tapas_sgm(be.*(2.*x-1),1);

% Initialize random number generator
if isnan(r.c_sim.seed)
    rng('shuffle');
else
    rng(r.c_sim.seed);
end

% Simulate
y = binornd(1, prob);

return;

author: translationalneuromodeling

repo: tapas


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

%*****************************************************************************80
%
%% NCOH_NESTED_RULE computes the nested Newton Cotes Open Half rule.
%
%  Discussion:
%
%    Nested version of Newton Cotes Open Half rule on [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), W(N), the points and weights.
%

%
%  Get N equally spaced points that do not include [ 0, 1];
%
  num = vdc_numerator_sequence ( n );
  num = 2 * num - 1;

  den = 2 ^ ( i4_log_2 ( n ) + 2 ) - 2;

  x = num / den;
%
%  Linear transformation from [0,1] to [-1,+1].
%
  x = 2.0 * x - 1.0;
%
%  Compute the weights.
%
  x_min = -1.0;
  x_max = +1.0;

  w = nc_compute ( n, x_min, x_max, x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = chi2inv(p,nu)
% Inverse-chi-squared distribution
% x = Gamma(nu/2,1/(2*p)) / Gammma(nu/2)

x = gammainc(1/(2*p),nu/2) / gamma(nu/2);

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function value = angle_rad ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_RAD returns the angle swept out between two rays.
%
%  Discussion:
%
%    Except for the zero angle case, it should be true that
%
%      ANGLE_RAD(P1,P2,P3) + ANGLE_RAD(P3,P2,P1) = 2 * PI
%
%        P1
%        /
%       /
%      /
%     /
%    P2--------->P3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), P3(2,1), define the rays
%    P1 - P2 and P3 - P2 which in turn define the angle.
%
%    Output, real VALUE, the angle swept out by the rays, measured
%    in radians.  0 <= VALUE < 2*PI.  If either ray has zero length,
%    then VALUE is set to 0.
%
  p(1,1) = ( p3(1,1) - p2(1,1) ) * ( p1(1,1) - p2(1,1) ) ...
         + ( p3(2,1) - p2(2,1) ) * ( p1(2,1) - p2(2,1) );

  p(2,1) = ( p3(1,1) - p2(1,1) ) * ( p1(2,1) - p2(2,1) ) ...
         - ( p3(2,1) - p2(2,1) ) * ( p1(1,1) - p2(1,1) );

  if ( p(1,1) == 0.0 & p(2,1) == 0.0 )
    value = 0.0;
    return
  end

  value = atan2 ( p(2,1), p(1,1) );

  if ( value < 0.0 )
    value = value + 2.0 * pi;
  end

  return
end


text:
function t = tvec_even3 ( nt )

%*****************************************************************************80
%
%% TVEC_EVEN3 computes an evenly spaced set of angles between 0 and 2*PI.
%
%  Discussion:
%
%    The angles begin with 0 and end with 2*PI.
%
%  Example:
%
%    NT = 4
%
%    T = ( 0, 2*PI/3, 4*PI/3 2*PI )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, the number of values to compute.
%
%    Output, real TVEC(NT), the evenly spaced angles, in radians.
%
  if ( nt == 1 )
    t(1) = pi;
  else
    t(1:nt) = ( 0:2:(2*nt-1) ) * pi / ( nt - 1 );
  end

  return
end


text:
classdef TestPlotting < testUnfitted

    properties (Access = protected)
        testName;
        meshType;
        meshIncludeBoxContour;
    end

    properties (Access = private)
        varAdim
    end

    methods (Access = public)

        function obj = TestPlotting(cParams)
            obj.init(cParams);
            obj.createTopOpt();
            obj.createMesh();
            obj.plot();
        end

        function plot(obj)
            figure();
            if isequal(obj.meshType,'BOUNDARY')
                obj.unfittedMesh.plotBoundary();
                obj.unfittedMesh.plotNormals();
            elseif isequal(obj.meshType,'INTERIOR')
                obj.unfittedMesh.plot();
            end
            view(obj.getViewAngle());
        end

        function passed = computePassed(obj)
            d = load(obj.testName);
            unfittedMesh = obj.unfittedMesh;
            loaded = d.unfittedMesh;
            passed = isequaln(unfittedMesh,loaded);
            if ~passed
                save(obj.testName, 'unfittedMesh', '-append')
            end
        end

        function overwriteResults(obj)
            save(obj.testName, 'unfittedMesh', '-append')
        end

    end

    methods (Access = private)

        function init(obj, cParams)
            obj.testName        = cParams.testName;
            obj.meshType        = cParams.meshType;
            obj.meshIncludeBoxContour = cParams.meshIncludeBoxContour;
        end

        function angle = getViewAngle(obj)
            if isprop(obj,'viewAngle')
                angle = obj.viewAngle;
            else
                angle = [0 0 1];
            end
        end

    end

    %% Heredat de testUnfitted, s'ha de declarar buit
    methods (Access = protected)
        function printTestNotPassed()
        end
        function printTestPassed()
        end
        function hasPassed()
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% LFDecodeLensletImageSimple - decodes a 2D lenslet image into a 4D light field, called by LFUtilDecodeLytroFolder
%
% Usage:
%
%   [LF, LFWeight, DecodeOptions, DebayerLensletImage, CorrectedLensletImage] = ...
%      LFDecodeLensletImageSimple( LensletImage, WhiteImage, LensletGridModel, DecodeOptions )
%
% This function follows the simple decode process described in:
% D. G. Dansereau, O. Pizarro, and S. B. Williams, "Decoding, calibration and rectification for
% lenslet-based plenoptic cameras," in Computer Vision and Pattern Recognition (CVPR), IEEE
% Conference on. IEEE, Jun 2013.
%
% This involves demosaicing, devignetting, transforming and slicing the input lenslet image to
% yield a 4D structure. More sophisticated approaches exist which combine steps into joint
% solutions, and they generally yield superior results, particularly near the edges of lenslets.
% The approach taken here was chosen for its simplicity and flexibility.
%
% Input LensletImage is a raw lenslet iamge as found within Lytro's LFP picture files. Though
% this function is written with Lytro imagery in mind, it should be possible to adapt it for use
% with other lenslet-based cameras. LensletImage is ideally of type uint16.
%
% Input WhiteImage is a white image as found within Lytro's calibration data. A white image is an
% image taken through a diffuser, and is useful for removing vignetting (darkening near edges of
% images) and for locating lenslet centers. The white image should be taken under the same zoom and
% focus settings as the light field. In the case of Lytro imagery, the helper functions
% LFUtilProcessWhiteImages and LFSelectFromDatabase are provided to assist in forming a list of
% available white images, and selecting the appropriate image based on zoom and focus settings.
%
% Input LensletGridModel is a model of the lenslet grid, as estimated, for example, using
% LFBuildLensletGridModel -- see that function for more on the required structure.
%
% Optional Input DecodeOptions is a structure containing:
%          [Optional] ResampMethod : 'fast'(default) or 'triangulation', the latter is slower
%          [Optional]   Precision : 'single'(default) or 'double'
%          [Optional] LevelLimits : a two-element vector defining the black and white levels
%
% Output LF is a 5D array of size [Nj,Ni,Nl,Nk,3]. See [1] and the documentation accompanying this
% toolbox for a brief description of the light field structure.
%
% Optional output LFWeight is of size [Nj,Ni,Nl,Nk]. LFWeight contains a confidence measure
% suitable for use in filtering applications that accept a weight parameter. This parameter is kept
% separate from LF rather than building in as a fourth channel in order to allow an optimization:
% when the parameter is not requested, it is not computed, saving significant processing time and
% memory.
%
% Optional output DebayerLensletImage is useful for inspecting intermediary results. The debayered
% lenslet image is the result of devignetting and debayering, with no further processing. Omitting
% this output variable saves memory.
%
% Optional output CorrectedLensletImage is useful for inspecting intermediary results. The
% corrected lenslet image has been rotated and scaled such that lenslet centers lie on straight lines, and
% every lenslet center lies on an integer pixel spacing. See [1] for more detail. Omitting
% this output variable saves memory.
%
% A more direct version of this function LFDecodeLensletImageDirect does not generate the
% intermediate corrected lenslet image.
%
% User guide: LFToolbox.pdf
% See also:  LFLytroDecodeImage, LFUtilDecodeLytroFolder, LFDecodeLensletImageDirect

% Copyright (c) 2013-2020 Donald G. Dansereau

function [LF, LFWeight, DecodeOptions, DebayerLensletImage, CorrectedLensletImage] = ...
    LFDecodeLensletImageSimple( LensletImage, WhiteImage, LensletGridModel, DecodeOptions )

%---Defaults---
 % LevelLimits defaults are a failsafe, should be passed in / come from the white image metadata
DecodeOptions = LFDefaultField( 'DecodeOptions', 'LevelLimits', [min(WhiteImage(:)), max(WhiteImage(:))] );
DecodeOptions = LFDefaultField( 'DecodeOptions', 'ResampMethod', 'fast' ); %'fast', 'triangulation'
DecodeOptions = LFDefaultField( 'DecodeOptions', 'Precision', 'single' );
DecodeOptions = LFDefaultField( 'DecodeOptions', 'DoDehex', true );
DecodeOptions = LFDefaultField( 'DecodeOptions', 'DoSquareST', true );

%---Rescale image values, remove black level---
DecodeOptions.LevelLimits = cast(DecodeOptions.LevelLimits, DecodeOptions.Precision);
BlackLevel = DecodeOptions.LevelLimits(1);
WhiteLevel = DecodeOptions.LevelLimits(2);
WhiteImage = cast(WhiteImage, DecodeOptions.Precision);
WhiteImage = (WhiteImage - BlackLevel) ./ (WhiteLevel - BlackLevel);

LensletImage = cast(LensletImage, DecodeOptions.Precision);
LensletImage = (LensletImage - BlackLevel) ./ (WhiteLevel - BlackLevel);
LensletImage = LensletImage ./ WhiteImage; % Devignette
% Clip -- this is aggressive and throws away bright areas; there is a potential for an HDR approach here
LensletImage = min(1, max(0, LensletImage));

if( nargout < 2 )
    clear WhiteImage
end

%---Demosaic---
% This uses Matlab's demosaic, which is "gradient compensated". This likely has implications near
% the edges of lenslet images, where the contrast is due to vignetting / aperture shape, and is not
% a desired part of the image
LensletImage = cast(LensletImage.*double(intmax('uint16')), 'uint16');
LensletImage = demosaic(LensletImage, DecodeOptions.DemosaicOrder);
LensletImage = cast(LensletImage, DecodeOptions.Precision);
LensletImage = LensletImage ./  double(intmax('uint16'));
DecodeOptions.NColChans = 3;

if( nargout >= 2 )
    DecodeOptions.NWeightChans = 1;
else
    DecodeOptions.NWeightChans = 0;
end

if( nargout > 3 )
    DebayerLensletImage = LensletImage;
end

%---Tranform to an integer-spaced grid---
fprintf('\nAligning image to lenslet array...');
InputSpacing = [LensletGridModel.HSpacing, LensletGridModel.VSpacing];
NewLensletSpacing = ceil(InputSpacing);
% Force even so hex shift is a whole pixel multiple
NewLensletSpacing = ceil(NewLensletSpacing/2)*2;
XformScale = NewLensletSpacing ./ InputSpacing;  % Notice the resized image will not be square

NewOffset = [LensletGridModel.HOffset, LensletGridModel.VOffset] .* XformScale;
RoundedOffset = round(NewOffset);
XformTrans =  RoundedOffset-NewOffset;

NewLensletGridModel = struct('HSpacing',NewLensletSpacing(1), 'VSpacing',NewLensletSpacing(2), ...
    'HOffset',RoundedOffset(1), 'VOffset',RoundedOffset(2), 'Rot',0, ...
    'UMax', LensletGridModel.UMax, 'VMax', LensletGridModel.VMax, 'Orientation', LensletGridModel.Orientation, ...
    'FirstPosShiftRow', LensletGridModel.FirstPosShiftRow);

%---Fix image rotation and scale---
RRot = LFRotz( LensletGridModel.Rot );

RScale = eye(3);
RScale(1,1) = XformScale(1);
RScale(2,2) = XformScale(2);
DecodeOptions.OutputScale(1:2) = XformScale;
DecodeOptions.OutputScale(3:4) = [1,2/sqrt(3)];  % hex sampling

RTrans = eye(3);
RTrans(end,1:2) = XformTrans;

% The following rotation can rotate parts of the lenslet image out of frame.
% todo[optimization]: attempt to keep these regions, offer greater user-control of what's kept
FixAll = maketform('affine', RRot*RScale*RTrans);
NewSize = size(LensletImage(:,:,1)) .* XformScale(2:-1:1);
LensletImage = imtransform( LensletImage, FixAll, 'YData',[1 NewSize(1)], 'XData',[1 NewSize(2)]);
if( nargout >= 2 )
    WhiteImage = imtransform( WhiteImage, FixAll, 'YData',[1 NewSize(1)], 'XData',[1 NewSize(2)]);
end
if( nargout >= 4 )
    CorrectedLensletImage = LensletImage;
end

LF = SliceXYImage( NewLensletGridModel, LensletImage, WhiteImage, DecodeOptions );
clear WhiteImage LensletImage

%---Correct for hex grid and resize to square u,v pixels---
LFSize = size(LF);
HexAspect = 2/sqrt(3);
switch( DecodeOptions.ResampMethod )
    case 'fast'
        fprintf('\nResampling (1D approximation) to square u,v pixels');
        NewUVec = 0:1/HexAspect:(size(LF,4)+1);  % overshoot then trim
        NewUVec = NewUVec(1:ceil(LFSize(4)*HexAspect));
        OrigUSize = size(LF,4);
        LFSize(4) = length(NewUVec);
        %---Allocate dest and copy orig LF into it (memory saving vs. keeping both separately)---
        LF2 = zeros(LFSize, DecodeOptions.Precision);
        LF2(:,:,:,1:OrigUSize,:) = LF;
        LF = LF2;
        clear LF2

        if( DecodeOptions.DoDehex )
            ShiftUVec = -0.5+NewUVec;
            fprintf(' and removing hex sampling...');
        else
            ShiftUVec = NewUVec;
            fprintf('...');
        end
        for( ColChan = 1:size(LF,5) )
            CurUVec = ShiftUVec;
            for( RowIter = 1:2 )
                RowIdx = mod(NewLensletGridModel.FirstPosShiftRow + RowIter, 2) + 1;
                ShiftRows = squeeze(LF(:,:,RowIdx:2:end,1:OrigUSize, ColChan));
                SliceSize = size(ShiftRows);
                SliceSize(4) = length(NewUVec);
                ShiftRows = reshape(ShiftRows, [size(ShiftRows,1)*size(ShiftRows,2)*size(ShiftRows,3), size(ShiftRows,4)]);
                ShiftRows = interp1( (0:size(ShiftRows,2)-1)', ShiftRows', CurUVec' )';
                ShiftRows(isnan(ShiftRows)) = 0;
                LF(:,:,RowIdx:2:end,:,ColChan) = reshape(ShiftRows,SliceSize);
                CurUVec = NewUVec;
            end
        end
        clear ShiftRows
        DecodeOptions.OutputScale(3) = DecodeOptions.OutputScale(3) * HexAspect;

    case 'triangulation'
        fprintf('\nResampling (triangulation) to square u,v pixels');
        OldVVec = (0:size(LF,3)-1);
        OldUVec = (0:size(LF,4)-1) * HexAspect;

        NewUVec = (0:ceil(LFSize(4)*HexAspect)-1);
        NewVVec = (0:LFSize(3)-1);
        LFSize(4) = length(NewUVec);
        LF2 = zeros(LFSize, DecodeOptions.Precision);

        [Oldvv,Olduu] = ndgrid(OldVVec,OldUVec);
        [Newvv,Newuu] = ndgrid(NewVVec,NewUVec);
        if( DecodeOptions.DoDehex )
            fprintf(' and removing hex sampling...');
            FirstShiftRow = NewLensletGridModel.FirstPosShiftRow;
            Olduu(FirstShiftRow:2:end,:) = Olduu(FirstShiftRow:2:end,:) + HexAspect/2;
        else
            fprintf('...');
        end

        DT = delaunayTriangulation( Olduu(:), Oldvv(:) );  % use DelaunayTri in older Matlab versions
        [ti,bc] = pointLocation(DT, Newuu(:), Newvv(:));
        ti(isnan(ti)) = 1;

        for( ColChan = 1:size(LF,5) )
            fprintf('.');
            for( tidx= 1:LFSize(1) )
                for( sidx= 1:LFSize(2) )
                    CurUVSlice = squeeze(LF(tidx,sidx,:,:,ColChan));
                    triVals = CurUVSlice(DT(ti,:));
                    CurUVSlice = dot(bc',triVals')';
                    CurUVSlice = reshape(CurUVSlice, [length(NewVVec),length(NewUVec)]);

                    CurUVSlice(isnan(CurUVSlice)) = 0;
                    LF2(tidx,sidx, :,:, ColChan) = CurUVSlice;
                end
            end
        end
        LF = LF2;
        clear LF2
        DecodeOptions.OutputScale(3) = DecodeOptions.OutputScale(3) * HexAspect;

    otherwise
        fprintf('\nNo valid dehex / resampling selected\n');
end

%---Resize to square s,t pixels---
% Assumes only a very slight resampling is required, resulting in an identically-sized output light field
if( DecodeOptions.DoSquareST )
    fprintf('\nResizing to square s,t pixels using 1D linear interp...');

    ResizeScale = DecodeOptions.OutputScale(1)/DecodeOptions.OutputScale(2);
    ResizeDim1 = 1;
    ResizeDim2 = 2;
    if( ResizeScale < 1 )
        ResizeScale = 1/ResizeScale;
        ResizeDim1 = 2;
        ResizeDim2 = 1;
    end

    OrigSize = size(LF, ResizeDim1);
    OrigVec = floor((-(OrigSize-1)/2):((OrigSize-1)/2));
    NewVec = OrigVec ./ ResizeScale;

    OrigDims = [1:ResizeDim1-1, ResizeDim1+1:5];

    UBlkSize = 32;
    USize = size(LF,4);
    LF = permute(LF,[ResizeDim1, OrigDims]);
    for( UStart = 1:UBlkSize:USize )
        UStop = UStart + UBlkSize - 1;
        UStop = min(UStop, USize);
        LF(:,:,:,UStart:UStop,:) = interp1(OrigVec, LF(:,:,:,UStart:UStop,:), NewVec);
        fprintf('.');
    end
    LF = ipermute(LF,[ResizeDim1, OrigDims]);
    LF(isnan(LF)) = 0;

    DecodeOptions.OutputScale(ResizeDim2) = DecodeOptions.OutputScale(ResizeDim2) * ResizeScale;
end

%---Trim s,t---
LF = LF(2:end-1,2:end-1, :,:, :);

%---Slice out LFWeight if it was requested---
if( nargout >= 2 )
    LFWeight = LF(:,:,:,:,end);
    LFWeight = LFWeight./max(LFWeight(:));
    LF = LF(:,:,:,:,1:end-1);
end

end

%------------------------------------------------------------------------------------------------------
function LF = SliceXYImage( LensletGridModel, LensletImage, WhiteImage, DecodeOptions )
% todo[optimization]: The SliceIdx and ValidIdx variables could be precomputed

fprintf('\nSlicing lenslets into LF...');

USize = LensletGridModel.UMax;
VSize = LensletGridModel.VMax;
MaxSpacing = max(LensletGridModel.HSpacing, LensletGridModel.VSpacing);  % Enforce square output in s,t
SSize = MaxSpacing + 1; % force odd for centered middle pixel -- H,VSpacing are even, so +1 is odd
TSize = MaxSpacing + 1;

LF = zeros(TSize, SSize, VSize, USize, DecodeOptions.NColChans + DecodeOptions.NWeightChans, DecodeOptions.Precision);

TVec = cast(floor((-(TSize-1)/2):((TSize-1)/2)), 'int16');
SVec = cast(floor((-(SSize-1)/2):((SSize-1)/2)), 'int16');
VVec = cast(0:VSize-1, 'int16');
UBlkSize = 32;
for( UStart = 0:UBlkSize:USize-1 )  % note zero-based indexing
    UStop = UStart + UBlkSize - 1;
    UStop = min(UStop, USize-1);
    UVec = cast(UStart:UStop, 'int16');

    [tt,ss,vv,uu] = ndgrid( TVec, SVec, VVec, UVec );

    %---Build indices into 2D image---
    LFSliceIdxX = LensletGridModel.HOffset + uu.*LensletGridModel.HSpacing + ss;
    LFSliceIdxY = LensletGridModel.VOffset + vv.*LensletGridModel.VSpacing + tt;

    HexShiftStart = LensletGridModel.FirstPosShiftRow;
    LFSliceIdxX(:,:,HexShiftStart:2:end,:) = LFSliceIdxX(:,:,HexShiftStart:2:end,:) + LensletGridModel.HSpacing/2;

    %---Lenslet mask in s,t and clip at image edges---
    CurSTAspect = DecodeOptions.OutputScale(1)/DecodeOptions.OutputScale(2);
    R = sqrt((cast(tt,DecodeOptions.Precision)*CurSTAspect).^2 + cast(ss,DecodeOptions.Precision).^2);
    ValidIdx = find(R < LensletGridModel.HSpacing/2 & ...
        LFSliceIdxX >= 1 & LFSliceIdxY >= 1 & LFSliceIdxX <= size(LensletImage,2) & LFSliceIdxY <= size(LensletImage,1) );

    %--clip -- the interp'd values get ignored via ValidIdx--
    LFSliceIdxX = max(1, min(size(LensletImage,2), LFSliceIdxX ));
    LFSliceIdxY = max(1, min(size(LensletImage,1), LFSliceIdxY ));

    %---
    LFSliceIdx = sub2ind(size(LensletImage), cast(LFSliceIdxY,'int32'), ...
        cast(LFSliceIdxX,'int32'), ones(size(LFSliceIdxX),'int32'));

    tt = tt - min(tt(:)) + 1;
    ss = ss - min(ss(:)) + 1;
    vv = vv - min(vv(:)) + 1;
    uu = uu - min(uu(:)) + 1 + UStart;
    LFOutSliceIdx = sub2ind(size(LF), cast(tt,'int32'), cast(ss,'int32'), ...
        cast(vv,'int32'),cast(uu,'int32'), ones(size(ss),'int32'));

    %---
    for( ColChan = 1:DecodeOptions.NColChans )
        LF(LFOutSliceIdx(ValidIdx) + numel(LF(:,:,:,:,1)).*(ColChan-1)) = ...
            LensletImage( LFSliceIdx(ValidIdx) + numel(LensletImage(:,:,1)).*(ColChan-1) );
    end
    if( DecodeOptions.NWeightChans ~= 0 )
        LF(LFOutSliceIdx(ValidIdx) + numel(LF(:,:,:,:,1)).*(DecodeOptions.NColChans)) = ...
            WhiteImage( LFSliceIdx(ValidIdx) );
    end
    fprintf('.');
end
end

repo: LFToolbox


text:
function n = ndims(T)
% overloads size

n = max(2,ndims(T.M) - T.rank);


text:
%% (Internal) Alineate wave positions (onset, peak and offsets commonly)
%
%     slpos_out = groupnlineup_waves(slpos)
%
%
% Arguments:
%
%   + slpos: delineation structure,
%
% Output:
%
%   + slpos_out: delineation structure with all heartbeats lined up with its
%   respective waves
%
% Example:
%
%
% See also alineate_positions
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate: 17/05/2016
% Last update: 17/05/2016
% Copyright 2008-2016
%
function slpos_out = groupnlineup_waves(slpos)

cant_leads = length(slpos);
[slpos_mat, slpos_names] = positions2matrix(slpos);

for this_lead = 1:cant_leads

    cant_QRS = length(slpos(this_lead).qrs);
    this_pos_mat = slpos_mat{this_lead};
    cant_anns = size(this_pos_mat,1);

    if( isempty(this_pos_mat) )
        continue
    end

    prev_QRS = 1;

    for ii = 1:cant_QRS

        if( ii == cant_QRS)
            next_QRS = slpos(this_lead).qrs(cant_QRS) + 1.5* (slpos(this_lead).qrs(cant_QRS)-slpos(this_lead).qrs(cant_QRS-1));
        else
            next_QRS = slpos(this_lead).qrs(ii+1);
        end

        slpos_out(this_lead).Pon(ii) = match_value(slpos(this_lead).Pon, slpos(this_lead).qrs(ii), prev_QRS, slpos(this_lead).qrs(ii));
        slpos_out(this_lead).P(ii) = match_value(slpos(this_lead).P, slpos(this_lead).qrs(ii), prev_QRS, slpos(this_lead).qrs(ii));
        slpos_out(this_lead).Poff(ii) = match_value(slpos(this_lead).Poff, slpos(this_lead).qrs(ii), prev_QRS, slpos(this_lead).qrs(ii));
        slpos_out(this_lead).Ptipo(ii) = nan;

        slpos_out(this_lead).QRSon(ii) = match_value(slpos(this_lead).QRSon, slpos(this_lead).qrs(ii), prev_QRS, slpos(this_lead).qrs(ii));
        slpos_out(this_lead).qrs(ii) = slpos(this_lead).qrs(ii);
        slpos_out(this_lead).Q(ii) = match_value(slpos(this_lead).Q, slpos(this_lead).qrs(ii), prev_QRS, next_QRS);
        slpos_out(this_lead).R(ii) = match_value(slpos(this_lead).R, slpos(this_lead).qrs(ii), prev_QRS, next_QRS);
        slpos_out(this_lead).S(ii) = match_value(slpos(this_lead).S, slpos(this_lead).qrs(ii), prev_QRS, next_QRS);
        slpos_out(this_lead).QRSoff(ii) = match_value(slpos(this_lead).QRSoff, slpos(this_lead).qrs(ii), slpos(this_lead).qrs(ii), next_QRS);

        slpos_out(this_lead).Ton(ii) = match_value(slpos(this_lead).Ton, slpos(this_lead).qrs(ii), slpos(this_lead).qrs(ii), next_QRS);
        slpos_out(this_lead).T(ii) = match_value(slpos(this_lead).T, slpos(this_lead).qrs(ii), slpos(this_lead).qrs(ii), next_QRS);
        slpos_out(this_lead).Tprima(ii) = match_value(slpos(this_lead).Tprima, slpos(this_lead).qrs(ii), slpos(this_lead).qrs(ii), next_QRS);
        slpos_out(this_lead).Toff(ii) = match_value(slpos(this_lead).Toff, slpos(this_lead).qrs(ii), slpos(this_lead).qrs(ii), next_QRS);
        slpos_out(this_lead).Ttipo(ii) = nan;

        prev_QRS = slpos(this_lead).qrs(ii);

    end

    if( cant_QRS == 0 )
        slposout_mat = [];
    else
        % group remaining waves if any
        slposout_mat = positions2matrix(slpos_out(this_lead));
    end

    cant_names = length(slpos_names);
    other_waves = cell(1,cant_names);

    for ii = 1:cant_names
        if( all(isnan(this_pos_mat(:,ii))) )
            other_waves{ii} = [];
        else
            if( isempty(slposout_mat) )
                aux_val = this_pos_mat(:,ii);
            else
                aux_val = setdiff(this_pos_mat(:,ii), slposout_mat(:,ii));
            end
            other_waves{ii} = aux_val(~isnan(aux_val));
        end
    end

    if( any(cellfun(@(a)(~isempty(a)), other_waves)) )
        % Group P waves and add as a new heartbeat
        [~, wave_idx] = intersect(slpos_names, {'Pon' 'P' 'Poff'});
        if( any(cellfun(@(a)(~isempty(a)), other_waves(wave_idx) )) )
            other_waves(wave_idx) = alineate_positions( other_waves(wave_idx) );
            aux_val = nan(length(other_waves{wave_idx(1)}),cant_names);
            aux_val(:,wave_idx) = cell2mat(other_waves(wave_idx));
            slposout_mat = [ slposout_mat; aux_val ];
        end

        % Group QRS waves and add as a new heartbeat
        [~, wave_idx] = intersect(slpos_names, {'QRSon' 'qrs' 'Q' 'R' 'S' 'QRSoff'});
        if( any(cellfun(@(a)(~isempty(a)), other_waves(wave_idx) )) )
            other_waves(wave_idx) = alineate_positions( other_waves(wave_idx) );
            aux_val = nan(length(other_waves{wave_idx(1)}),cant_names);
            aux_val(:,wave_idx) = cell2mat(other_waves(wave_idx));
            slposout_mat = [ slposout_mat; aux_val ];
        end

        % Group T waves and add as a new heartbeat
        [~, wave_idx] = intersect(slpos_names, {'Ton' 'T' 'Toff'});
        if( any(cellfun(@(a)(~isempty(a)), other_waves(wave_idx) )) )
            other_waves(wave_idx) = alineate_positions( other_waves(wave_idx) );
            aux_val = nan(length(other_waves{wave_idx(1)}),cant_names);
            aux_val(:,wave_idx) = cell2mat(other_waves(wave_idx));
            slposout_mat = [ slposout_mat; aux_val ];
        end

        slpos_out(this_lead) = matrix2positions(slposout_mat, slpos_names);
    end

    % sort heartbeats and waves and add as a new heartbeat
    [~, aux_idx] = sort((sum(~isnan(this_pos_mat))), 'descend');
    aux_val = slpos_out(this_lead).( slpos_names{aux_idx(1)} );

    [~, aux_idx] = sort(aux_val);

    for fname = slpos_names
        if( isfield(slpos_out(this_lead), fname{1}) )
            aux_val = slpos_out(this_lead).(fname{1});
            if( isempty(aux_val) )
                slpos_out(this_lead).(fname{1}) = nan(cant_anns,1);
            else
                slpos_out(this_lead).(fname{1}) = colvec(aux_val(aux_idx));
            end
        else
            slpos_out(this_lead).(fname{1}) = nan(cant_anns,1);
        end
    end

end

function close_val = get_closer(values, reference)

    close_val = values( min_index(abs(values-reference)) );

function vals_in_range = get_values_in_range(values, start_range, end_range)

    vals_in_range = values( values >= start_range & values <= end_range );

function matched_val = match_value(values, reference, start_range, end_range)

    aux_val = get_values_in_range(values, start_range, end_range);

    if( isempty(aux_val) )
        matched_val = nan;
    else
        matched_val = get_closer(aux_val, reference);
        if( isempty(matched_val) )
            matched_val = nan;
        end
    end

author: marianux

repo: ecg-kit

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


text:
function [ know, x ] = p10_sol ( n )

%*****************************************************************************80
%
%% P10_SOL returns the solution for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = [ 1.0E+06, 2.0E-06 ]';

  return
end


text:
function lines = assign_lineclass(lines, vp)

THRES_THETA = 10;
if length(vp)>=1 && ~isempty(vp{1})
    lineclass1 = line_belongto_vp(lines, vp{1}, THRES_THETA);
else
    lineclass1 = zeros(1, length(lines));
end
if length(vp)>=2 && ~isempty(vp{2})
    lineclass2 = line_belongto_vp(lines, vp{2}, THRES_THETA);
else
    lineclass2 = zeros(1, length(lines));
end
if length(vp)>=3 && ~isempty(vp{3})
    lineclass3 = line_belongto_vp(lines, vp{3}, THRES_THETA);
else
    lineclass3 = zeros(1, length(lines));
end

for i = 1:length(lines)
	lines(i).lineclass1 = lineclass1(i);
	lines(i).lineclass2 = lineclass2(i);
	lines(i).lineclass3 = lineclass3(i);

	if lineclass1(i) + lineclass2(i) + lineclass3(i) == 1
		lines(i).lineclass = 1*lineclass1(i) + 2*lineclass2(i) + 3*lineclass3(i);
	else
		lines(i).lineclass = 0;
	end
end

author: zouchuhang

repo: LayoutNet


text:
function [vec,istart,iend] = SHCreateVec(lmax)

% [vec,start,end] = SHCreateVec(lmax)
%
% For a given degree lmax, creates an array of zero-valued spherical
% harmonic coefficients. If lmax is a vector of length N, returns a
% concatenated vector of spherical harmonics of degrees lmax(1) .. lmax(N).
% The arrays 'start' and 'end' contain the first and the last indices
% for each of the N sections, so that vec(start(i):end(i)) exactly
% corresponds to the spherical harmonics degree lmax(i).

N = length(lmax);

nmax=zeros(1,N);
istart=zeros(1,N);
vec=[];

for k=1:N
    if lmax(k)<0
        error('invalid usage: lmax must be a non-negative integer');
    end
    nmax(k)=SHl2n(lmax(k));
    section=zeros(1,nmax(k))';
    vec = [vec;section];
end

istart(1)=1;
for k=2:N
    istart(k)=istart(k-1)+nmax(k-1);
end
iend = istart+nmax-1;

repo: mcbench-benchmarks


text:
function view = update3ViewLoc(view);
% view = update3ViewLoc(view);
%
% Updates the location of a volume 3-view window
% based on the settings of the slice view/orientation
% fields at the bottom of the window.

% 03/07/03 by ras

axiSlice=str2num(get(view.ui.sliceNumFields(1),'String'));
corSlice=str2num(get(view.ui.sliceNumFields(2),'String'));
sagSlice=str2num(get(view.ui.sliceNumFields(3),'String'));

loc = [axiSlice corSlice sagSlice];

view = refreshScreen(view,loc);

return

author: vistalab

repo: vistasoft


text:
function q = ch_qnormlz(q)
% 四元数归一化
    q = q/norm(q);
    if(q(1)<0)
        q(1) = -q(1);
        q(2) = -q(2);
        q(3) = -q(3);
        q(4) = -q(4);
    end

author: yandld


text:
function str=vec2str(vec)
    %       string=vec2str(vector)
    % convert a numeric vector to an equivalent string
    % so that eval(string)=vector                                rcobb 9/95
    %
    if min(size(vec)) > 1,error('error, expecting vector input'),end
    if isempty(vec)
        str = '';
    else
        str=['[',num2str(vec(1))];
        len=length(vec);
        indx = 2;
        while indx <= len
            sindx = indx;
            while vec(sindx-1)+1 == vec(sindx) & sindx <= len-1
                sindx = sindx + 1;
                if sindx == len && vec(sindx-1)+1 == vec(sindx)
                    str=[str,':',num2str(vec(sindx))];
                    str=[str,']'];
                    return
                end
            end
            if sindx == indx
                str=[str,',',num2str(vec(indx))];
                indx =sindx + 1;
            else
                str=[str,':',num2str(vec(sindx-1))];
                indx =sindx;
            end
        end
        str=[str,']'];
    end

author: CelsoReyes


text:
function value = p16_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P16_F evaluates the integrand for problem 16.
%
%  Discussion:
%
%    The integrand can be regarded as the L1 norm of X - Z.
%
%    It would be nice to allow the use to specify several
%    base points Z, to make the function more jagged more places%
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    The integrand can be regarded as the L1 norm of X - Z.
%
%    There is a basis point Z associated with the integrand.
%    Z(1:DIM_NUM) defaults to ( 0.5, 0.5, ..., 0.5 ).
%    The user can set, get, or randomize this value by calling
%    P16_R8VEC.
%
%  Integrand:
%
%    sum ( abs ( x(1:dim_num) - z(1:dim_num) ) )
%
%  Exact Integral:
%
%    The integral is separable into
%
%       Int ( A(1) <= X(1) <= B(1) ) abs ( X(1) - Z(1) )
%     * Product ( B(1:N)-A(1:N), skip index 1 )
%     + Int ( A(2) <= X(2) <= B(2) ) abs ( X(2) - Z(2) )
%     * Product ( B(1:N)-A(1:N), skip index 2 )
%     ...
%     + Int ( A(N) <= X(N) <= B(N) ) abs ( X(N) - Z(N) )
%     * Product ( B(1:N)-A(1:N), skip index N )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  z = [];
  z = p16_r8vec ( 'G', 'Z', dim_num, z );

  value(1:point_num) = 0.0;

  for point = 1 : point_num

    value(point) = sum ( abs ( x(1:dim_num,point) - z(1:dim_num)' ) );

  end

  p16_i4 ( 'I', '#', point_num );

  return
end


text:
function [maps,data,stats] = RippleStats(filtered,ripples,varargin)

%RippleStats - Compute descriptive stats for ripples (100~200Hz oscillations).
%
%  USAGE
%
%    [maps,data,stats] = RippleStats(filtered,ripples,<options>)
%
%    filtered       ripple-band filtered samples (one channel)
%    ripples        ripple timing information (obtained using FindRipples)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'frequency'   sampling rate (in Hz) (default = 1250Hz)
%     'durations'   durations before and after ripple peak (in s)
%                   (default = [-0.5 0.5])
%    =========================================================================
%
%  OUTPUT
%
%    maps.ripples               instantaneous amplitude (one ripple per row)
%    maps.frequency             instantaneous frequency
%    maps.phase                 instantaneous phase
%    maps.amplitude             enveloppe amplitude
%    data.peakFrequency         frequency at peak
%    data.peakAmplitude         amplitude at peak
%    data.duration              durations
%    stats.durationAmplitude    correlation between duration and amplitude (rho, p)
%    stats.durationFrequency    correlation between duration and frequency (rho, p)
%    stats.amplitudeFrequency   correlation between amplitude and frequency (rho, p)
%    stats.acg.data             autocorrelogram data
%    stats.acg.t                autocorrelogram time bins
%
%  SEE
%
%    See also FindRipples, SaveRippleEvents, PlotRippleStats.

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

% Default values
samplingRate = 1250;
durations = [-50 50]/1000;
nCorrBins = 50;
%  corrBinSize = 400;
corrDuration = 20;
corrBinSize = 0.01;

% Check number of parameters
if nargin < 2,
  error('Incorrect number of parameters (type ''help RippleStats'' for details).');
end

% Check parameter sizes
if size(filtered,2) ~= 2,
	error('Parameter ''filtered'' is not a Nx2 matrix (type ''help RippleStats'' for details).');
end
%  if size(ripples,2) ~= 4,
%  	error('Parameter ''ripples'' is not a Nx4 matrix (type ''help RippleStats'' for details).');
%  end
if size(ripples,2) < 3,
	error('Parameter ''ripples'' is not a Nx4 or Nx4 matrix (type ''help RippleStats'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help RippleStats'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'frequency',
			samplingRate = varargin{i+1};
			if ~isdscalar(samplingRate,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help RippleStats'' for details).');
			end
		case 'durations',
			durations = varargin{i+1};
			if ~isdvector(durations,'#2','<'),
				error('Incorrect value for property ''durations'' (type ''help RippleStats'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help RippleStats'' for details).']);
	end
end

nBins = floor(samplingRate*diff(durations)/2)*2+1; % must be odd
nHalfCenterBins = 3;
centerBin = ceil(nBins/2);
centerBins = centerBin-nHalfCenterBins:centerBin+nHalfCenterBins;

% Compute instantaneous phase and amplitude
h = hilbert(filtered(:,2));
phase(:,1) = filtered(:,1);
phase(:,2) = angle(h);
amplitude(:,1) = filtered(:,1);
amplitude(:,2) = abs(h);
unwrapped(:,1) = filtered(:,1);
unwrapped(:,2) = unwrap(phase(:,2));
% Compute instantaneous frequency
frequency = Diff(unwrapped,'smooth',0);
frequency(:,2) = frequency(:,2)/(2*pi);

% Compute ripple map
[r,i] = Sync(filtered,ripples(:,2),'durations',durations);
maps.ripples = SyncMap(r,i,'durations',durations,'nbins',nBins,'smooth',0);

% Compute frequency Map
[f,i] = Sync(frequency,ripples(:,2),'durations',durations);
maps.frequency = SyncMap(f,i,'durations',durations,'nbins',nBins,'smooth',0);

% Compute phase map
[p,i] = Sync(phase,ripples(:,2),'durations',durations);
maps.phase = SyncMap(p,i,'durations',durations,'nbins',nBins,'smooth',0);

% Compute amplitude map
[a,i] = Sync(amplitude,ripples(:,2),'durations',durations);
maps.amplitude = SyncMap(a,i,'durations',durations,'nbins',nBins,'smooth',0);

% Ripple frequency and amplitude at peak
data.peakFrequency = maps.frequency(:,centerBin);
data.peakAmplitude = maps.amplitude(:,centerBin);

% Ripple durations
data.duration = ripples(:,3)-ripples(:,1);

% Autocorrelogram and correlations
%  if nargin > 2,
%  	[stats.acg.data,stats.acg.t] = CCG(ripples(:,2),1,'binSize',corrBinSize,'halfBins',nCorrBins/2);
	[stats.acg.data,stats.acg.t] = CCG(ripples(:,2),1,'binSize',corrBinSize);
	[stats.amplitudeFrequency.rho,stats.amplitudeFrequency.p] = corrcoef(data.peakAmplitude,data.peakFrequency);
	[stats.durationFrequency.rho,stats.durationFrequency.p] = corrcoef(data.duration,data.peakFrequency);
	[stats.durationAmplitude.rho,stats.durationAmplitude.p] = corrcoef(data.duration,data.peakAmplitude);
%  end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Author: Relja Arandjelovic (relja@relja.info)

function Y= relja_l1normalize_col( X )
    Y= bsxfun(@rdivide, X, sum(abs(X),1) + 1e-12 );
end

author: HajimeTaira


text:
function figure_num = cube3d_grid_plot ( x1, x2, x3, filename, figure_num )

%*****************************************************************************80
%
%% CUBE3D_GRID_PLOT plots hypersphere gridpoints onto the surface of a 3D cube.
%
%  Discussion:
%
%    The X1, X2, X3 data has, presumably, been computed by CUBE_GRID for
%    a 3D cube.
%
%    This function plots the surface, and then the surface plus the
%    projected grid points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X1(2*N+1,N+1), X2(2*N+1,N+1), X3(2*N+1,N+1), the coordinates
%    of points on the cube surface, projected from the hypersphere.
%
%    Input, string FILENAME, the "first name" of the two files to be created:
%    filename_surface.png and filename_points.png.
%
%    Input/output, int FIGURE_NUM, the current figure index.
%
  if ( nargin < 4 )
    filename = 'cube';
  end
%
%  Draw the surface.
%
  figure_num = figure_num + 1;
  figure ( figure_num )
  mesh ( x3, x2, x1, 'FaceColor', 'interp' );
  axis equal
  grid on
  xlabel ( '<---X--->', 'FontSize', 16 );
  ylabel ( '<---Y--->', 'FontSize', 16 );
  zlabel ( '<---Z--->', 'FontSize', 16 );
  title ( 'Cube transition surface', 'FontSize', 24 )
  hold off
  filename1 = sprintf ( '%s_surface.png', filename );
  print ( '-dpng', filename1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created plotfile "%s".\n', filename1 );
%
%  Draw the surface and points.
%
%  Listing the points in the order X1, X2, X3, rational as it may seem,
%  is apparently NOT the way to go here.
%
  figure_num = figure_num + 1;
  figure ( figure_num )
  hold on
  mesh ( x3, x2, x1, 'FaceColor', 'interp' );
  plot3 ( x3, x2, x1, 'k.', 'MarkerSize', 5 );
  axis equal
  grid on
  xlabel ( '<---X--->', 'FontSize', 16 );
  ylabel ( '<---Y--->', 'FontSize', 16 );
  zlabel ( '<---Z--->', 'FontSize', 16 );
  title ( 'Grid points on transition surface', 'FontSize', 24 )
  hold off
  filename2 = sprintf ( '%s_points.png', filename );
  print ( '-dpng', filename2 );
  fprintf ( 1, '  Created plotfile "%s".\n', filename2 );

  return
end


text:
function e = calcError(pf,rec,varargin)
% RP and mean square error
%
% *calcError(pf,rec)* calculates reconstruction error between meassured
% intensities and the recalcuated ODF or between two meassured pole
% figures. It can be specified whether the RP
% error or the mean square error is calculated. The scaling coefficients
% are calculated by the function PoleFigure/calcNormalization
%
% Syntax
%   e = calcError(pf,pf2) % compares two different @PoleFigure with same @S2Grid
%   e = calcError(pf,rec) % compares @PoleFigure with the Recalculated @SO3Fun
%
% Input
%  pf,pf2 - @PoleFigure
%  rec    - @SO3Fun
%
% Output
%  e - error
%
% Flags
%  RP    - RP value (default) |pfmeas - pfcalc|./ pfcalc
%  l1    - l1 error           |pfmeas - pfcalc|
%  l2    - l2 error           |pfmeas - pfcalc|.^2
%
% See also
% SO3Fun/calcError PoleFigure/calcNormalization PoleFigure/scale

argin_check(rec,{'SO3Fun','PoleFigure'});

% calc difference PoleFigure
errorpf = calcErrorPF(pf,rec,varargin{:});

% calc error
e = zeros(1,pf.numPF);
for i = 1:pf.numPF

  e(i) = mean(errorpf.allI{i}(:)); % RP error

  if check_option(varargin,'l1')
    e(i) = e(i)/mean(abs(pf.allI{i}(:))); % L^1 error
  elseif check_option(varargin,'l2')
    e(i) = e(i)/mean((pf.allI{i}(:)).^2); % L^2 error
  end
end

% TODO: implement a nice default output
%if nargout == 0
%  disp('TODO')
%  clear e;
%end


text:
% Main camera calibration toolbox:
format compact

repo: ekfmonoslam