text:
function [funs, ends] = constructor(op, dom, data, pref)
%CONSTRUCTOR   CHEBFUN constructor.
%   FUNS = CONSTRUCTOR(OP, DOM) constructs the piecewise components (known as
%   "FUNS") used by a CHEBFUN object to represent the function OP on the
%   interval DOM. OP must be a function_handle, string, numerical vector, or a
%   cell array containing a combination of these first three data types. In the
%   later case, the number of elements in the array must be one less than the
%   length of the DOM vector.
%
%   It is not expected that CHEBFUN.CONSTRUCTOR() be called directly,
%
%   If OP is a function_handle or a string, it should be vectorised in that it
%   accepts a column vector of length N and return a matrix of size N x M. If M
%   ~= 1, we say the resulting CHEBFUN is "array-valued".
%
%   CONSTRUCTOR(OP, DOM, DATA, PREF), where DATA is a MATLAB structure and PREF
%   is a CHEBFUNPREF object, allows construction data and alternative
%   construction preferences to be passed to the constructor.  See CHEBFUNPREF
%   for more details on preferences.
%
%   In particular, if PREF.SPLITTING = TRUE and OP is a function_handle or a
%   string, then the constructor adaptively introduces additional breakpoints
%   into the domain so as to better represent the function. These are returned
%   as the second output argument in [FUNS, END] = CONSTRUCTOR(OP, DOM).
%
%   The DATA structure input contains information which needs to be passed to
%   the lower layers about parameters which may affect the construction process.
%   Presently, the only fields CONSTRUCTOR expects DATA to have on input are
%   DATA.EXPONENTS and DATA.SINGTYPE, which convey information about endpoint
%   singularities. These fields are populated by CHEBFUN.PARSEINPUTS as need be.
%   Before calling the FUN constructor, DATA will be augmented to include
%   information about the construction domain as well as the horizontal and
%   vertical scales involved in the construction procedure.
%
% See also CHEBFUN, CHEBFUNPREF.

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

% Initial setup:
numIntervals = numel(dom) - 1;

% Initialise hscale and vscale:
data.hscale = norm(dom, inf);
if ( isinf(data.hscale) )
    data.hscale = 1;
end
if ( isempty(data.vscale) )
    data.vscale = 0;
end

% Sanity check:
if ( iscell(op) && (numel(op) ~= numIntervals) )
    error('CHEBFUN:CHEBFUN:constructor:cellInput', ...
        ['Number of cell elements in OP must match the number of ', ...
         'intervals in DOMAIN.'])
end

% Construct the FUNs.
if ( pref.splitting )
    [funs, ends] = constructorSplit(op, dom, data, pref);
else
    [funs, ends] = constructorNoSplit(op, dom, data, pref);
end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%  SPLITTING OFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [funs, ends] = ...
    constructorNoSplit(op, domain, data, pref)
% In 'OFF' mode, seek only one piece with length < maxLength.

% Initial setup:
numIntervals = numel(domain) - 1;
ends = domain;

% Initialise the FUN array:
funs = cell(1, numIntervals);

% We only want to throw the warning 'CHEBFUN:constructor:notResolved'once:
warningThrown = false;

singDetect = pref.blowup;
exps = data.exponents;
singTypes = data.singType;

% Loop over the intervals:
for k = 1:numIntervals
    endsk = ends(k:k+1);
    % Unwrap if OP is a cell:
    if ( iscell(op) )
        opk = op{k};
    else
        opk = op;
    end

    if ( singDetect )
        % Extract the exponents for this interval:
        if ( ~isempty(exps) )
            data.exponents = exps(2*k-1:2*k);
        end
        % Replace the information for the singularity type in the preference:
        if ( ~isempty(singTypes) && ~ischar(singTypes) )
            data.singType = singTypes(2*k-1:2*k);
        end
    end

    % Call GETFUN() (which calls FUN.CONSTRUCTOR()):
    [funs{k}, ishappy, data.vscale] = getFun(opk, endsk, data, pref);

    % Warn if unhappy (as we're unable to split the domain to improve):
    if ( ~ishappy && ~warningThrown )
        if ( strcmpi(func2str(pref.tech), 'trigtech') )
            str = 'a non-trig representation';
        else
            str = '''splitting on''';
        end

        warning('CHEBFUN:CHEBFUN:constructor:notResolved', ...
            ['Function not resolved using %d pts.', ...
            ' Have you tried ' str, '?'], pref.techPrefs.maxLength);
        warningThrown = true;
    end
end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%  SPLITTING ON %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [funs, ends] = constructorSplit(op, dom, data, pref)
% In 'ON' mode, seek only many pieces with total length < maxlength.

% Initial setup:
numIntervals = numel(dom) - 1;
ends = dom;

% Set the maximum length (i.e., number of sample points for CHEBTECH):
pref.techPrefs.maxLength = pref.splitPrefs.splitLength;

% We extrapolate when splitting so that we can construct functions like
% chebfun(@sign,[-1 1]), which otherwise would not be happy at x = 0.
pref.techPrefs.extrapolate = true;

% Initialise the FUN array:
funs = cell(1, numIntervals);
% Initialise happiness:
ishappy = ones(1, numel(ends) - 1);

singDetect = pref.blowup;
exps = data.exponents;
singTypes = data.singType;

% Try to get one smooth piece for the entire domain before splitting:
for k = 1:numIntervals
    % Unwrap if OP is a cell:
    if ( iscell(op) )
        opk = op{k};
    else
        opk = op;
    end

    if ( singDetect )
        % Extract the singularity information for this interval:
        data = getSingInfo(exps, singTypes, 2*k-1:2*k, data);
    end

    [funs{k}, ishappy(k), data.vscale] = ...
        getFun(opk, ends(k:k+1), data, pref);

    % For the case where vscale is Inf due to blowup in the interior of the
    % domain:
    if ( isinf(data.vscale) )
        % An infinite vscale doesn't help us at all, but ruin the consequential
        % constructions at the lower levels:
        data.vscale = 0;
    end

end
sad = ~ishappy;

% MAIN LOOP. If the above didn't work, enter main loop and start splitting.
% (Note that at least one new breakpoint will be introduced).
while ( any(sad) )
    % If a FUN is sad in a subinterval, split this subinterval.

    % Choose a subinterval to improve:

%     % Old choice = the first sad interval:
%     k = find(sad, 1, 'first');

    % New choice = the largest sad interval:
    diffEnds = diff(ends);
    diffEnds(~sad) = 0;
    [ignored, k] = max(diffEnds);

    % Ends of this subinterval:
    a = ends(k);
    b = ends(k+1);

    % Unwrap if OP is a cell:
    if ( iscell(op) )
        opk = op{k};
    else
        opk = op;
    end

    % Look for an edge:
    edge = fun.detectEdge(funs{k}, opk, data.hscale, data.vscale, pref);

    if ( singDetect )
        % Update singularity info:
        if ( ~isempty(exps) )
            exps = [exps(1:2*k-1), zeros(1,2), exps(2*k:end)];
        end
        if ( ~isempty(singTypes) )
            singTypes = [singTypes(1:2*k-1), singTypes(2*k-1), ...
                singTypes(2*k) singTypes(2*k:end)];
        end
    end

    if ( singDetect )
        % Extract the singularity information for this interval:
        data = getSingInfo(exps, singTypes, 2*k-1:2*k, data);
    end

    % Try to obtain happy child FUN objects on each new subinterval:
    [childLeft, happyLeft, data.vscale] = getFun(opk, [a, edge], data, pref);

    if ( singDetect )
        % Extract the singularity information for this interval:
        data = getSingInfo(exps, singTypes, 2*k+1:2*k+2, data);
    end

    [childRight, happyRight, data.vscale] = getFun(opk, [edge, b], data, pref);

    % Check for happiness/sadness:
    sad = [sad(1:k-1), ~happyLeft, ~happyRight, sad(k+1:end)];

    % Insert new pieces in to existing funs:
    funs = [funs(1:k-1), {childLeft, childRight}, funs(k+1:end)];
    ends = [ends(1:k), edge, ends(k+1:end)];

    % If a cell was given, we must store pieces on new intervals:
    if ( iscell(op) )
        op = [op(1:k), {opk}, op(k+1:end)];
    end

    % Fail if too many points are required:
    len = sum(cellfun(@length, funs));
    if ( len > pref.splitPrefs.splitMaxLength )
        warning('CHEBFUN:CHEBFUN:constructor:funNotResolved', ...
            'Function not resolved using %d pts.', ...
            sum(cellfun(@length, funs)));
        return
    end

end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GETFUN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [g, ishappy, vscale] = getFun(op, dom, data, pref)
%GETFUN   Call the FUN constructor.

% If the interval is very small then skip adaptation and treat OP as a constant:
if ( diff(dom) < 4*1e-14*data.hscale && ~isnumeric(op) )
    op = op(mean(dom));
end

% Bolt domain information onto the data structure.
data.domain = dom;

% Call the FUN constructor:
g = fun.constructor(op, data, pref);

% See if the construction was happy:
ishappy = get(g, 'ishappy');

% Update the vertical scale:
if ( ishappy )
    vscale = max([data.vscale, get(g, 'vscale')]);
else
    vscale = data.vscale;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GETSINGINFO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function data = getSingInfo(exps, singTypes, kk, data)
% Place information about the singularity type in the data structure.
if ( ~isempty(exps) )
    data.exponents = exps(kk);
end
if ( ~isempty(singTypes) )
    data.singType = singTypes(kk);
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function img = showTemplates(img, A, A_mean, A_norm, tsize, nt)
if(size(img,3) == 1)
    for i=1:3
        img_color(:,:,i) = img;
    end
    img = img_color;
end

for i=1:nt
    tt = A(:,i);
    tt = tt*A_norm(i)+A_mean(:,i);
    tt = reshape(tt, [tsize, 3]);
    for j=1:3
        img(end-tsize(1)+1:end,(i-1)*tsize(2)+1:i*tsize(2),j) = tt(:,:,j);
    end
end
%save templates templates;

author: yuxng


text:
function mr=v_rotqr2mr(qr)
%V_ROTQR2MR converts a matrix of real quaternion vectors to quaternion matrices
% Inputs:
%
%     QR(4m,n,...)   mxn matrix of real quaternion vectors (each 4x1)
%
% Outputs:
%
%     MR(4m,4n,...)   mxn matrix of real quaternion matrices (each 4x4)
%
% In matrix form, quaternions can be multiplied and added using normal matrix
% arithmetic. Each element of an mxn matrix of quaternions is itself a 4x4 block
% so the total dimension of MR is 4m x 4n.

%
%      Copyright (C) Mike Brookes 2000-2018
%      Version: $Id: v_rotqr2mr.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent a b c
if isempty(a)
    a=[1 2 3 3 1 2];    % destination row of +ve entries (from 0)
    b=[1 2 3 2 3 1];    % destination col of +ve entries (from 0)
    c=[0 0 0 1 2 3];    % source row of +ve entries (from 0)
end
s=size(qr);
m=s(1);
mr=repmat(reshape(qr,s(1),[]),4,1);
n=size(mr,2);
mn=s(1)*n;
j=repmat(4*m*(0:n-1),m/4,1);
i=repmat((1:4:m)',n,1)+j(:);
ni=length(i);
i6=repmat(i,1,6);
mr(i6+repmat(a+m*b,ni,1))=mr(i6+repmat(c,ni,1));
mr(i6+repmat(c+m*b,ni,1))=-mr(i6+repmat(a,ni,1));
s(2)=4*s(2); % output array size
mr=reshape(mr,s);
if ~nargout
    qr=qr(1:4);
    v_rotqr2ro(qr(:)); % plot a rotated cube
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [desired_state] = circle(t, qn)
% CIRCLE trajectory generator for a circle

% =================== Your code goes here ===================
% You have to set the pos, vel, acc, yaw and yawdot variables
time_tol = 12;
radius = 5;
dt = 0.0001;

    function pos = pos_from_angle(a)
        pos = [radius*cos(a); radius*sin(a); 2.5*a/(2*pi)];
    end

    function vel = get_vel(t)
        angle1 = tj_from_line(0, 2*pi, time_tol, t);
        pos1 = pos_from_angle(angle1);
        angle2 = tj_from_line(0, 2*pi, time_tol, t+dt);
        vel = (pos_from_angle(angle2) - pos1)/dt;
    end

if t > time_tol
    pos = [radius; 0; 2.5];
    vel = [0;0;0];
    acc = [0;0;0];
else
    angle = tj_from_line(0, 2*pi, time_tol, t);
    pos = pos_from_angle(angle);
    vel = get_vel(t);
    acc = (get_vel(t+dt) - get_vel(t))/dt;
end

yaw = 0;
yawdot = 0;
% =================== Your code ends here ===================

desired_state.pos = pos(:);
desired_state.vel = vel(:);
desired_state.acc = acc(:);
desired_state.yaw = yaw;
desired_state.yawdot = yawdot;

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function [pnew,timing,x_min,upper] = diagonalize_quadratic_program(p,timing,x_min,upper)
pnew = p;
pnew.V = [];
pnew.diagonalized = 0;

if ~p.options.bmibnb.diagonalize
    return
end

% Preprocessing has moved nonliner constraints
% to bounds, so don't be fooled...
if p.originallyNonlinearConstraints
    return
end

% No quadratic terms
if all(p.variabletype == 0)
    return
end

% Any polynomial terms or simple linear by some reason
if any(p.variabletype > 2) & ~all(p.variabletype == 0) | p.options.bmibnb.diagonalize==0
    return
end

if ~isempty(p.evalVariables)
    return
end

if ~isempty(p.binary_variables) | ~isempty(p.integer_variables)
    return
end

if isempty(p.F_struc)
    % Bound-constrained QP is faster in original
    % form as we can exploit concavity
    return
end

nonlinear = find(p.variabletype > 0);
linear = find(p.variabletype == 0);
if ~isempty(p.F_struc)
    % Nonlinear terms in constraints
    if nnz(p.F_struc(:,1 + nonlinear))> 0
        return
    end
end

% Find quadratic and linear terms
used_in_c = find(p.c);
quadraticterms = used_in_c(find(ismember(used_in_c,nonlinear)));
Q = zeros(length(p.c),length(p.c));
if ~isempty(quadraticterms)
    usedinquadratic = zeros(1,length(p.c));
    for i = 1:length(quadraticterms)
        Qij = p.c(quadraticterms(i));
        power_index = find(p.monomtable(quadraticterms(i),:));
        if length(power_index) == 1
            Q(power_index,power_index) = Qij;
        else
            Q(power_index(1),power_index(2)) = Qij/2;
            Q(power_index(2),power_index(1)) = Qij/2;
        end
    end
end
Qlin = Q(linear,linear);
clin = p.c(linear);

if nnz(Qlin-diag(diag(Qlin)))==0
    % already diagonal
    return
end

% Decompose Q
[V,D] = eig(full(Qlin));
V = real(V);
D = real(D);
V(abs(V)<1e-11) = 0;
D(abs(D)<1e-11) = 0;
lb = p.lb(linear);
ub = p.ub(linear);
Z = V';
newub = sum([Z>0].*Z.*repmat(ub',length(Z),1),2)+sum([Z<0].*Z.*repmat(lb',length(Z),1),2);
newlb = sum([Z>0].*Z.*repmat(lb',length(Z),1),2)+sum([Z<0].*Z.*repmat(ub',length(Z),1),2);
newub(isnan(newub)) = inf;
newlb(isnan(newlb)) = -inf;

% Create new problem
clin = V'*clin;

n = length(linear);
pnew.original_linear = linear;
pnew.original_n = length(p.c);
pnew.V = V;
pnew.c = [clin;diag(D)];
pnew.Q = spalloc(2*n,2*n,0);

% find constraint polytope
if size(p.F_struc,1)>0
    A = -p.F_struc(:,1 + linear);
    b = p.F_struc(:,1);
    pnew.F_struc = [b -A*V zeros(length(b),n)];
    Abounds = [V;-V];
    bbounds = [ub;-lb];
    keep = find(~isinf(bbounds));
    if ~isempty(keep)
        pnew.F_struc = [pnew.F_struc;bbounds(keep) -Abounds(keep,:) zeros(length(keep),n)];
        pnew.K.l = pnew.K.l + length(keep);
    end
end

pnew.variabletype = [zeros(1,n) ones(1,n)*2];
pnew.monomtable = [eye(n);2*eye(n)];
pnew.monomtable(2*n,2*n) = 0;
pnew.lb =-inf(2*n,1);
pnew.ub = inf(2*n,1);
pnew.lb(1:n) = newlb;
pnew.ub(1:n) = newub;
pnew.lb(n+1:end) = 0;
pnew.ub(n+1:end) = max(newlb.^2,newub.^2);
if length(p.x0)>0
    pnew.x0 = [(V'*p.x0(p.linears));(V'*p.x0(p.linears)).^2];
end
x_min = [V'*x_min(p.linears);(V'*x_min(p.linears)).^2];
pnew.diagonalized = 1;

% Ok ,we now got a new diagonalized model
% Sometimes it performs much worse though
% Make trial calls to judge it
if p.options.bmibnb.diagonalize == -1
    tstart = tic;
    % If the model lacks bounds, these relaxations will
    % fail at this point, as bound propagation hasn't been performed yet.
    % Add fake bounds and see which performs best on fake model
    p.socpcuts.F_struc = [];
    pnew.socpcuts.F_struc = [];
    p.delayedconvex = [];
    pnew.delayedconvex = [];
    p.upper=inf;
    pnew.upper=inf;
    [output1,cost1,~,timing] = solvelower(fakeBounds(defineQuadratics(fakeLowerModel(pnew))),p.options,p.solver.lowersolver.call,[],[],timing);
    [output2,cost2,~,timing] = solvelower(fakeBounds(fakeLowerModel(p)),p.options,p.solver.lowersolver.call,[],[],timing);
    [upper_cost1,x_u1] = quickCost(pnew,output1.Primal);
    [upper_cost2,x_u2] = quickCost(p,output2.Primal);
    timing.lowersolve = timing.lowersolve + toc(tstart);
    if output1.problem && output2.problem
        % both had issues. Let's assume diagonal is better
        if p.options.verbose>=0;display('* -Diagonalized QP');end
    elseif output1.problem && ~output2.problem
        % new model has issues. use old
        pnew = p;
        pnew.V = [];
        pnew.diagonalized = 0;
        upper = upper_cost2;x_min = x_u2;
        if p.options.verbose>=0;display('* -Diagonalized QP but switched back to original.');end
    elseif ~output1.problem && output2.problem
        % old has issue, so keep new
        upper = upper_cost1;x_min = x_u1;
        if p.options.verbose>=0;display('* -Diagonalized QP.');end
    else
        % Both solved, so which had the best lower bound
        if cost1 < cost2
            % New model weaker, keep old.
            % Take best upper though!
            x_min = x_u2;
            upper = upper_cost2;
            if upper_cost1 < upper
                upper = upper_cost1;
                x_min(p.linears)=V*x_u1(pnew.linears);
            end
            pnew = p;
            pnew.V = [];
            pnew.diagonalized = 0;
            if p.options.verbose>=0;display('* -Diagonalized QP but switched back to original.');end
        else
            % New model better
            % Take best upper though!
            x_min = x_u1;
            upper = upper_cost1;
            if upper_cost2 < upper
                upper = upper_cost2;
                x_min(pnew.linears)=V'*x_u2(p.linears);
            end
            if p.options.verbose>=0;display('* -Diagonalized QP.');end
        end
    end
else
	if p.options.verbose>=0;display('* -Diagonalized QP.');end
end

function pnew = fakeLowerModel(p)
pnew = p;
pnew.complementary=[];
pnew.shiftedQP=[];
pnew.originalModel.variabletype=[];
pnew.concavityEqualities=[];
pnew.EqualityConstraintState=[];
pnew.EqualityConstraintState=[];
pnew.InequalityConstraintState=[];

function [q,z] = quickCost(p,x)
q = 0;
z=[];
for i = 1:size(p.monomtable,1)
    j = find(p.monomtable(i,:));
    v = x(j);
    w = p.monomtable(i,j);
   	z =[z;prod(v(:).^w(:))];
    q = q + p.c(i)*z(i);
end

function p = defineQuadratics(p)
n = length(p.c)/2;
p.Quadratics = n+1:2*n;
p.nonlinears = n+1:2*n;
p.QuadraticsList = [zeros(n,2);(1:n)' (1:n)'];
p.bilinears = [(n+1:2*n)' (1:n)' (1:n)'];
%p.solver.lowersolver.constraint.inequalities.secondordercone.linear=0;

function p = fakeBounds(p)
U = 1e3;
p.ub(p.linears) = min(p.ub(p.linears),U);
p.lb(p.linears) = max(p.lb(p.linears),-U);
p.ub(p.Quadratics) = min(p.ub(p.Quadratics),U^2);
p.lb(p.Quadratics) = max(p.lb(p.Quadratics),-U^2);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function volume = simplex_unit_volume_nd ( dim_num )

%*****************************************************************************80
%
%% SIMPLEX_UNIT_VOLUME_ND computes the volume of the unit simplex in ND.
%
%  Discussion:
%
%    The formula is simple: volume = 1/N!.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Output, real VOLUME, the volume of the cone.
%
  volume = 1.0;
  for i = 1 : dim_num
    volume = volume / i;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef RMMEDA_F9 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing RM-MEDA

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, and Y. Jin, RM-MEDA: A regularity model-based
% multiobjective estimation of distribution algorithm, IEEE Transactions on
% Evolutionary Computation, 2008, 12(1): 41-63.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = [1,zeros(1,obj.D-1)+10];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,2:end).^2 - repmat(X(:,1),1,size(X,2)-1);
            g = sum(t.^2/4000,2) - prod(cos(t./repmat(sqrt(1:size(X,2)-1),size(X,1),1)),2) + 2;
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-sqrt(PopObj(:,1)./g));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function dlY = model(dlX,parameters,hyperparameters,doTraining)

numBlocks = hyperparameters.NumBlocks;
dropoutFactor = hyperparameters.DropoutFactor;

dlY = dlX;

% Residual blocks.
for k = 1:numBlocks
    dilationFactor = 2^(k-1);
    parametersBlock = parameters.("Block"+k);

    dlY = residualBlock(dlY,dilationFactor,dropoutFactor,parametersBlock,doTraining);
end

% Fully connect
weights = parameters.FC.Weights;
bias = parameters.FC.Bias;
dlY = fullyconnect(dlY,weights,bias,'DataFormat','CBT');

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [varargout] = plus(varargin)

% PLUS compile the missing mex file on the fly

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

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);
  mex(mexsrc);
  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was just created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a=ighmap(b)
%IGHMAP Single-level inverse discrete 2-D multiwavelet transform.
%   IGHMAP performs a single-level 2-D multiwavelet reconstruction
%   using GHM multiwavelet with four multi-filters
%
%   Y = IGHMAP(X) computes the original matrix from the approximation
%   coefficients matrix LL and details coefficients matrices
%   LH, HL, HH, obtained by a multiwavelet decomposition of the
%   original input matrix and puts the result in Y.
%
%   The size of Y is the same as that of X which should be
%   a square matrix of size NxN where N is power of 2 since
%   X is de-vectorized by critical sampling preprocessing.
%   The postprocessing filter is of order 2 degree 1.
%   LL, LH, HL, and HH should have the size N/2xN/2.
%
%   X should be arranged as [LL,LH;HL,HH].
%
%   See also GHM, IGHM, GHMAP, GHMAP2, IGHMAP2, WAVEDEC2, WAVEINFO.

%   Auth: Dr. Bessam Z. Hassan
%   Last Revision: 27-Feb-2004.
%   Copyright 1995-2002 The MathWorks, Inc.
% $Revision: 1.0 $

if nargin == 0,
	error('Not enough input arguments.');
end
if isempty(b)
   a = [];
   return
end
H0=[3/(5*sqrt(2)),4/5;-1/20,-3/(10*sqrt(2))];
H1=[3/(5*sqrt(2)),0;9/20,1/sqrt(2)];
H2=[0,0;9/20,-3/(10*sqrt(2))];
H3=[0,0;-1/20,0];
G0=[-1/20,-3/(10*sqrt(2));1/(10*sqrt(2)),3/10];
G1=[9/20,-1/sqrt(2);-9/(10*sqrt(2)),0];
G2=[9/20,-3/(10*sqrt(2));9/(10*sqrt(2)),-3/10];
G3=[-1/20,0;-1/(10*sqrt(2)),0];
[N,M]=size(b);
if N~=2^round(log(N)/log(2))
    error('size of the input should be power of 2');
end
if M~=N
    error('the input matrix must be square');
end
w=[H0,H1,H2,H3;G0,G1,G2,G3];
for i=1:N/4-1
    W(4*(i-1)+1:4*i,4*i-3:4*i+4)=w;
end
W=[W;[[H2,H3;G2,G3],zeros(4,N-8),[H0,H1;G0,G1]]];
p=[];X=[];

% shuffle

N=N/2;
b1=b(1:N,1:N);b2=b(1:N,N+1:2*N);b3=b(N+1:2*N,1:N);b4=b(N+1:2*N,N+1:2*N);
b1=b1';b2=b2';b3=b3';b4=b4';
T(1:2:N,:)=b1(1:N/2,:);T(2:2:N,:)=b1(N/2+1:N,:);T=T';b1(1:2:N,:)=T(1:N/2,:);b1(2:2:N,:)=T(N/2+1:N,:);
T(1:2:N,:)=b2(1:N/2,:);T(2:2:N,:)=b2(N/2+1:N,:);T=T';b2(1:2:N,:)=T(1:N/2,:);b2(2:2:N,:)=T(N/2+1:N,:);
T(1:2:N,:)=b3(1:N/2,:);T(2:2:N,:)=b3(N/2+1:N,:);T=T';b3(1:2:N,:)=T(1:N/2,:);b3(2:2:N,:)=T(N/2+1:N,:);
T(1:2:N,:)=b4(1:N/2,:);T(2:2:N,:)=b4(N/2+1:N,:);T=T';b4(1:2:N,:)=T(1:N/2,:);b4(2:2:N,:)=T(N/2+1:N,:);
b=[b1,b2;b3,b4];
p=b';

%column vector permutation

ii=0:4:2*N-1;jj=sort([ii+1,ii+2]);kk=sort([ii+3,ii+4]);
z(jj,:)=p(1:N,:);z(kk,:)=p(N+1:2*N,:);
N=N*2;

%inverse column transform

X=W'*z;

%postprocess columns

aa(2:2:N,:)=X(2:2:N,:)/(sqrt(2)-1);
aa(1:2:N,:)=(X(1:2:N,:)-0.11086198019724*aa(2:2:N,:)-0.11086198019724*[zeros(1,N);aa(2:2:N-1,:)])/0.37361535735427;
p=aa';N=N/2;

%row vector permutation

ii=0:4:2*N-1;jj=sort([ii+1,ii+2]);kk=sort([ii+3,ii+4]);z=[];
z(jj,:)=p(1:N,:);z(kk,:)=p(N+1:2*N,:);

%inverse row transform

X=W'*z;N=N*2;

%postprocess rows

a(2:2:N,:)=X(2:2:N,:)/(sqrt(2)-1);
a(1:2:N,:)=(X(1:2:N,:)-0.11086198019724*a(2:2:N,:)-0.11086198019724*[zeros(1,N);a(2:2:N-1,:)])/0.37361535735427;

author: Sable

repo: mcbench-benchmarks


text:
function area = sphericalCapsAreaC6(varargin)
%SPHERICALCAPSAREAC6 compute area of spherical caps on the unit sphere
%
%   AREA = sphericalCapsAreaC6
%   Compute the area of the spherical caps associated to the voronoi
%   diagram on the unit sphere, when the germs correspond to the 26
%   discrete directions in the unit cube.
%   Result is a 1x3 array, containing fraction of area of each type of
%   spherical cap.
%   area(1) concerns the direction [1 0 0]
%   area(2) concerns the direction [0 1 0]
%   area(3) concerns the direction [0 0 1]
%   Result is formatted such that 2*(area(1)+area(2)+area(3))=1.
%   For homogeneous lattices, all areas are equal.
%
%   AREA = sphericalCapsAreaC6(DELTA)
%   where DELTA = [DELTA1 DELTA2 DELTA3]  specifies the resolution of unit
%   voxel in each direction.
%
%
%   Algorithm :
%   - separate 3 basic case, corresponding on the type of directions
%   - for each cases, specify manually germ and neighbours
%   - compute great circle between each couple of germ
%   - compute intersection points of these circles
%   - sort intersection points around the central germ, giving a spherical
%      polygon
%   - compute triangulation of the resulting polygon
%   - compute spherical area of each triangle
%   - and area of each polygon
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY
%   27/04/2007 extends to non uniform grid spacing

%% Initializations

% grid resolution
delta = [1 1 1];
if ~isempty(varargin)
    delta = varargin{1};
end

% poins in the 26 discrete directions
%pt000 = normalize([-1 -1 -1].*delta);
%pt100 = normalize([ 0 -1 -1].*delta);
%pt200 = normalize([+1 -1 -1].*delta);
%pt010 = normalize([-1  0 -1].*delta);
pt110 = normalize([ 0  0 -1].*delta);
% pt210 = normalize([+1  0 -1].*delta);
% pt020 = normalize([-1 +1 -1].*delta);
% pt120 = normalize([ 0 +1 -1].*delta);
% pt220 = normalize([+1 +1 -1].*delta);

%pt001 = normalize([-1 -1  0].*delta);
pt101 = normalize([ 0 -1  0].*delta);
% pt201 = normalize([+1 -1  0].*delta);
pt011 = normalize([-1  0  0].*delta);
%pt111 = [ 0  0  0];             % origin point
pt211 = normalize([+1  0  0].*delta);
% pt021 = normalize([-1 +1  0].*delta);
pt121 = normalize([ 0 +1  0].*delta);
% pt221 = normalize([+1 +1  0].*delta);

% pt002 = normalize([-1 -1 +1].*delta);
% pt102 = normalize([ 0 -1 +1].*delta);
% pt202 = normalize([+1 -1 +1].*delta);
% pt012 = normalize([-1  0 +1].*delta);
pt112 = normalize([ 0  0 +1].*delta);
% pt212 = normalize([+1  0 +1].*delta);
% pt022 = normalize([-1 +1 +1].*delta);
% pt122 = normalize([ 0 +1 +1].*delta);
% pt222 = normalize([+1 +1 +1].*delta);

% basic unit sphere
sphere = [0 0 0  1];

%% Spherical cap type 1 direction [1 0 0]
% ----------------------------------------------------------------------

% Compute area of voronoi cell for a point on the Ox axis, i.e. a point
% in the 6-neighboorhoud of the center.
centre1 = pt211;
% neighbours of chosen point, sorted by angle
%voisins1 = [pt200;pt201;pt202;pt212;pt222;pt221;pt220;pt210];
voisins1 = [pt112;pt121;pt110;pt101];
n1 = size(voisins1, 1);

% compute separating circles
planes1 = zeros(n1, 9);
for i=1:n1
    planes1(i,1:9) = normalizePlane(medianPlane(centre1, voisins1(i,:)));
%    circles1(i,1:7) = intersectPlaneSphere(planes1(i,:), sphere);
end

% compute circle intersections
lines1 = zeros(n1, 6);
for i=1:n1
    lines1(i,1:6) = intersectPlanes(planes1(i,:), ...
        planes1(mod(i,n1)+1,:));
    points1(2*i-1:2*i,1:3) = intersectLineSphere(lines1(i,:), sphere);
end

% keep only points with x>0
ind = dot(points1, repmat(centre1, [2*n1 1]), 2)>0;
points1 = points1(ind,:);
n1 = size(points1, 1);

% compute spherical area of each triangle [center  pt[i+1]%4   pt[i] ]
angles1 = zeros(n1, 1);
for i=1:n1
    pt1 = points1(i, :);
    pt2 = points1(mod(i  , n1)+1, :);
    pt3 = points1(mod(i+1, n1)+1, :);

    angles1(i) = sphericalAngle(pt1, pt2, pt3);
end

% compute area of spherical polygon
area1 = sum(angles1)-pi*(n1-2);

%% Spherical cap type 1 direction [0 1 0]
% ----------------------------------------------------------------------

% Compute area of voronoi cell for a point on the Oy axis, i.e. a point
% in the 6-neighboorhoud of the center.
centre1 = pt121;
% neighbours of chosen point, sorted by angle
%voisins1 = [pt200;pt201;pt202;pt212;pt222;pt221;pt220;pt210];
%voisins1 = [pt221;pt222;pt122;pt022;pt021;pt020;pt120;pt220];
voisins1 = [pt211;pt112;pt011;pt110];
n1 = size(voisins1, 1);

% compute separating circles
planes1 = zeros(n1, 9);
for i=1:n1
    planes1(i,1:9) = normalizePlane(medianPlane(centre1, voisins1(i,:)));
%    circles1(i,1:7) = intersectPlaneSphere(planes1(i,:), sphere);
end

% compute circle intersections
lines1 = zeros(n1, 6);
for i=1:n1
    lines1(i,1:6) = intersectPlanes(planes1(i,:), ...
        planes1(mod(i,n1)+1,:));
    points1(2*i-1:2*i,1:3) = intersectLineSphere(lines1(i,:), sphere);
end

% keep only points with x>0
ind = dot(points1, repmat(centre1, [2*n1 1]), 2)>0;
points1 = points1(ind,:);
n1 = size(points1, 1);

% compute spherical area of each triangle [center  pt[i+1]%4   pt[i] ]
angles1 = zeros(n1, 1);
for i=1:n1
    pt1 = points1(i, :);
    pt2 = points1(mod(i  , n1)+1, :);
    pt3 = points1(mod(i+1, n1)+1, :);

    angles1(i) = sphericalAngle(pt1, pt2, pt3);
end

% compute area of spherical polygon
area2 = sum(angles1)-pi*(n1-2);

%% Spherical cap type 1 direction [0 0 1]
% ----------------------------------------------------------------------

% Compute area of voronoi cell for a point on the Oz axis, i.e. a point
% in the 6-neighboorhoud of the center.
centre1 = pt112;
% neighbours of chosen point, sorted by angle
%voisins1 = [pt200;pt201;pt202;pt212;pt222;pt221;pt220;pt210];
%voisins1 = [pt212;pt222;pt122;pt022;pt012;pt002;pt102;pt202];
voisins1 = [pt121;pt011;pt101;pt211];

n1 = size(voisins1, 1);

% compute separating circles
planes1 = zeros(n1, 9);
for i=1:n1
    planes1(i,1:9) = normalizePlane(medianPlane(centre1, voisins1(i,:)));
%    circles1(i,1:7) = intersectPlaneSphere(planes1(i,:), sphere);
end

% compute circle intersections
lines1 = zeros(n1, 6);
for i=1:n1
    lines1(i,1:6) = intersectPlanes(planes1(i,:), ...
        planes1(mod(i,n1)+1,:));
    points1(2*i-1:2*i,1:3) = intersectLineSphere(lines1(i,:), sphere);
end

% keep only points with z>0
ind = dot(points1, repmat(centre1, [2*n1 1]), 2)>0;
points1 = points1(ind,:);
n1 = size(points1, 1);

% compute spherical area of each triangle [center  pt[i+1]%4   pt[i] ]
angles1 = zeros(n1, 1);
for i=1:n1
    pt1 = points1(i, :);
    pt2 = points1(mod(i  , n1)+1, :);
    pt3 = points1(mod(i+1, n1)+1, :);

    angles1(i) = sphericalAngle(pt1, pt2, pt3);
    angles1(i) = min(angles1(i), 2*pi-angles1(i));
end

% compute area of spherical polygon
area3 = sum(angles1)-pi*(n1-2);

%% Results
% -------------------------------------------------------------------

% display some results
disp('results : ')
pattern = 'area of cell in direction %d: %12.10f sr, = %12.10f %%';
disp(sprintf(pattern, 1, area1, area1*100/4/pi));
disp(sprintf(pattern, 2, area2, area2*100/4/pi));
disp(sprintf(pattern, 3, area3, area3*100/4/pi));

% return computed areas, formatted as fraction of sphere surface
area = [area1 area2 area3]/4/pi;

%% Internal functions
% -------------------------------------------------------------------
% This functions are part of a more general geometric library, but are
% included here for avoiding dependencies

function theta = angle3Points(varargin)
%ANGLE3POINTS return oriented angle made by 3 points
%
%   ALPHA = ANGLE3POINTS(P1, P2, P3).
%   Pi are either [1*2] arrays, or [N*2] arrays, in this case ALPHA is a
%   [N*1] array. The angle computed is the directed angle between line
%   (P2P1) and line (P2P3).
%   Result is always given in radians, between 0 and 2*pi.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 23/02/2004.
%

%   HISTORY :
%   25/09/2005 : enable single parameter

if length(varargin)==3
    p1 = varargin{1};
    p2 = varargin{2};
    p3 = varargin{3};
elseif length(varargin)==1
    var = varargin{1};
    p1 = var(1,:);
    p2 = var(2,:);
    p3 = var(3,:);
end

% angle line (P2 P1)
theta = lineAngle(createLine(p2, p1), createLine(p2, p3));

function line = createLine(varargin)
%CREATELINE create a line with various inputs.
%
%   Line is represented in a parametric form : [x0 y0 dx dy]
%   x = x0 + t*dx
%   y = y0 + t*dy;
%
%
%   l = CREATELINE(p1, p2) return the line going through the two given
%   points.
%
%   l = CREATELINE(x0, y0, dx, dy) the line going through point (x0, y0)
%   and with direction vector(dx, dy).
%
%   l = CREATELINE(param) where param is an array of 4 values, create the
%   line oing through the point (param(1) param(2)), and with direction
%   vector (param(3) param(4)).
%
%   l = CREATELINE(theta) create a line originated at (0,0) and
%   with angle theta.
%
%   l = CREATELINE(rho, theta) create a line with normal theta, and with
%   min distance to origin equal to rho. rho can be negative, in this case,
%   the line is the same as with CREATELINE(-rho, theta+pi), but the
%   orientation is different.
%
%
%   Note : in all cases, parameters can be vertical arrays of the same
%   dimension. The result is then an array of lines, of dimensions [N*4].
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

if length(varargin)==1
    % Only one input parameter. It can be :
    % - line angle
    % - array of four parameters
    var = varargin{1};

    if size(var, 2)==4
        % 4 parameters of the line in a single array.
        line = var;
    elseif size(var, 2)==1
        % 1 parameter : angle of the line, going through origin.
        line = [zeros(size(var)) zeros(size(var)) cos(var) sin(var)];
    else
        error('wrong number of dimension for arg1 : can be 1 or 4');
    end

elseif length(varargin)==2
    % 2 input parameters. They can be :
    % - line angle and signed distance to origin.
    % - 2 points, then 2 arrays of 1*2 double.
    v1 = varargin{1};
    v2 = varargin{2};
    if size(v1, 2)==1
        % first param is angle of line, and second param is signed distance
        % to origin.
        line = [v1.*cos(v2) v1.*sin(v2) -sin(v2) cos(v2)];
    else
        % first input parameter is first point, and second input is the
        % second point.
        line = [v1(:,1), v1(:,2), v2(:,1)-v1(:,1), v2(:,2)-v1(:,2)];
    end

elseif length(varargin)==3
    % 3 input parameters :
    % first one is a point belonging to the line,
    % second and third ones are direction vector of the line (dx and dy).
    p = varargin{1};
    line = [p(:,1) p(:,2) varargin{2} varargin{3}];

elseif length(varargin)==4
    % 4 input parameters :
    % they are x0, y0 (point belongng to line) and dx, dy (direction vector
    % of the line).
    % All parameters should have the same size.
    line = [varargin{1} varargin{2} varargin{3} varargin{4}];
else
    error('Wrong number of arguments in ''createLine'' ');
end

function plane = createPlane(varargin)
%CREATEPLANE create a plane in parametrized form
%
%   Create a plane in the following format :
%   PLANE = [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where :
%   - (X0, Y0, Z0) is a point belonging to the plane
%   - (DX1, DY1, DZ1) is a first direction vector
%   - (DX2, DY2, DZ2) is a second direction vector
%
%
%
%   PLANE = createPlane(P1, P2, P3)
%   create a plane containing the 3 points
%
%   PLANE = createPlane(PTS)
%   The 3 points are packed into a single 3x3 array.
%
%   PLANE = createPlane(P0, N);
%   create a plane from a point and from a normal to the plane.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY :
%   24/11/2005 : add possibility to pack points for plane creation

if length(varargin)==1
    var = varargin{1};

    if iscell(var)
        plane = zeros(length(var), 9);
        for i=1:length(var)
            plane(i,:) = createPlane(var{i});
        end
    elseif size(var, 1)==3
        % 3 points in a single array
        p1 = var(1,:);
        p2 = var(2,:);
        p3 = var(3,:);

        % create direction vectors
        v1 = p2-p1;
        v2 = p3-p1;

        % create plane
        plane = [p1 v1 v2];
        return;
    end

elseif length(varargin)==2

    p0 = varargin{1};

    var = varargin{2};
    if size(var, 2)==2
        n = sph2cart2([var repmat(1, [size(var, 1) 1])]);
    elseif size(var, 2)==3
        n  = normalize3d(var);
    else
        error ('wrong number of parameters in createPlane');
    end

    % find a vector not colinear to the normal
    v0 = repmat([1 0 0], [size(p0, 1) 1]);
    if abs(cross(n, v0, 2))<1e-14
        v0 = repmat([0 1 0], [size(p0, 1) 1]);
    end

    % create direction vectors
    v1 = normalize3d(cross(n, v0, 2));
    v2 = -normalize3d(cross(v1, n, 2));

    plane = [p0 v1 v2];
    return;

elseif length(varargin)==3
    p1 = varargin{1};
    p2 = varargin{2};
    p3 = varargin{3};

    % create direction vectors
    v1 = p2-p1;
    v2 = p3-p1;

    plane = [p1 v1 v2];
    return;

else
    error('wrong number of arguments in "createPlane".');
end

function point = intersectLineSphere(line, sphere)
%INTERSECTLINESPHERE return intersection between a line and a sphere
%
%   GC = intersectLineSphere(LINE, SPHERE) return the two points which are
%   the intersection of the given line and sphere.
%   LINE   : [x0 y0 z0  dx dy dz]
%   SPHERE : [xc yc zc  R]
%   GC     : [x1 y1 z1 ; x2 y2 z2]
%
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY

% difference between centers
dc = line(1:3)-sphere(1:3);

a = sum(line(:, 4:6).*line(:, 4:6), 2);
b = 2*sum(dc.*line(4:6), 2);
c = sum(dc.*dc, 2) - sphere(:,4).*sphere(:,4);

delta = b.*b -4*a.*c;

if delta>1e-14
    % find two roots of second order equation
    u1 = (-b -sqrt(delta))/2/a;
    u2 = (-b +sqrt(delta))/2/a;

    % convert into 3D coordinate
    point = [line(1:3)+u1*line(4:6) ; line(1:3)+u2*line(4:6)];

elseif abs(delta) > 1e-14
    % find unique root, and convert to 3D coord.
    u = -b/2./a;
    point = line(1:3) + u*line(4:6);

else
    point = zeros(0, 3);
    return;
end

function point = intersectPlaneLine(plane, line)
%INTERSECTPLANELINE return intersection between a plane and a line
%
%   PT = intersectPlaneSphere(PLANE, LINE) return the intersection point of
%   the given line and the given plane.
%   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   LINE :  [x0 y0 z0 dx dy dz]
%   PT :    [XI YI ZI]
%
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/02/2005.
%

%   HISTORY
%   24/11/2005 add support for multiple input

% unify sizes of data
if size(plane, 1)~=size(line, 1)
    if size(plane, 1)==1
        plane = repmat(plane, [size(line, 1) 1]);
    elseif size(line,1)==1
        line = repmat(line, [size(plane, 1) 1]);
    else
        error('line and plane do not have the same dimension');
    end
end

% plane normal
n = cross(plane(:,4:6), plane(:, 7:9), 2);

% test if line and plane are parallel
if abs(dot(n, line(:,4:6), 2))<1e-14
    point = [NaN NaN NaN];
    return;
end

% difference between origins of plane and line
dp = plane(:,1:3) - line(:,1:3);

% relative position of intersection on line
t = dot(n, dp, 2)/dot(n, line(:,4:6), 2);

% compute coord of intersection point
point = line(:,1:3) + t*line(:,4:6);

function line = intersectPlanes(plane1, plane2)
%INTERSECTPLANES return intersection between 2 planes in space
%
%   PT = intersectPlanes(PLANE1, PLANE2) return the straight line belonging
%   to both planes
%   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   LINE :  [x0 y0 z0 dx dy dz]
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/02/2005.
%

%   HISTORY

% plane normal
n1 = normalize3d(cross(plane1(:,4:6), plane1(:, 7:9), 2));
n2 = normalize3d(cross(plane2(:,4:6), plane2(:, 7:9), 2));

% test if planes are parallel
if abs(cross(n1, n2, 2))<1e-14
    line = [NaN NaN NaN NaN NaN NaN];
    return;
end

% Uses hessian form, ie : N.p = d
% I this case, d can be found as : -N.p0, when N is normalized
d1 = dot(n1, plane1(:,1:3), 2);
d2 = dot(n2, plane2(:,1:3), 2);

% compute dot products
dot1 = dot(n1, n1, 2);
dot2 = dot(n2, n2, 2);
dot12 = dot(n1, n2, 2);

det = dot1*dot2 - dot12*dot12;
c1 = (d1*dot2 - d2*dot12)./det;
c2 = (d2*dot1 - d1*dot12)./det;

p0 = c1*n1 + c2*n2;
dp = cross(n1, n2, 2);

line = [p0 dp];

function theta = lineAngle(varargin)
%LINEANGLE return angle between lines
%
%   a = LINEANGLE(line) return the angle between horizontal, right-axis
%   and the given line. Angle is fiven in radians, between 0 and 2*pi,
%   in counter-clockwise direction.
%
%   a = LINEANGLE(line1, line2) return the directed angle between the
%   two lines. Angle is given in radians between 0 and 2*pi.
%
%   see createLine for more details on line representation.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY :
%   19/02/2004 : added support for multiple lines.

nargs = length(varargin);
if nargs == 1
    % one line
    line = varargin{1};
    theta = mod(atan2(line(:,4), line(:,3)) + 2*pi, 2*pi);
elseif nargs==2
    % two lines
    theta1 = lineAngle(varargin{1});
    theta2 = lineAngle(varargin{2});
    theta = mod(theta2-theta1+2*pi, 2*pi);
end

function plane = medianPlane(p1, p2)
%MEDIANPLANE create a plane in the middle of 2 points
%
%   plane = medianPlane(P1, P2)
%   plane is perpendicular to line (P1 P2) and contains the midpoint of p1
%   and p2.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

p0 = (p1 + p2)/2;
n = p2-p1;

plane = createPlane(p0, n);

function vn = normalize(v)
%NORMALIZE normalize a vector
%
%   V2 = normalize(V);
%   return the normalization of vector V, such that ||V|| = 1. V can be
%   either a row or a column vector.
%
%   When V is a MxN array, normalization is performed for each row of the
%   array.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 29/11/2004.
%

dim = size(v);
if dim(1)==1 || dim(2)==1
    vn = v/sqrt(sum(v.*v));
else
    vn = v./repmat(sqrt(sum(v.*v, 2)), [1 dim(2)]);
end

function vn = normalize3d(v)
%NORMALIZE3D normalize a 3D vector
%
%   V2 = normalize3d(V);
%   return the normalization of vector V, such that ||V|| = 1. Vector V is
%   given in vertical form.
%
%   When V is a Nx3 array, normalization is performed for each row of the
%   array.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 29/11/2004.
%

%   HISTORY
%   30/11/2005  : correct a bug

n = sqrt(v(:,1).*v(:,1) + v(:,2).*v(:,2) + v(:,3).*v(:,3));
vn = v./[n n n];

function plane2 = normalizePlane(plane1)
%NORMALIZEPLANE normalize parametric form of a plane
%
%   plane2 = normalizePlane(plane1);
%   transform the plane PANE1 in the following format :
%   [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where :
%   - (X0, Y0, Z0) is a point belonging to the plane
%   - (DX1, DY1, DZ1) is a first direction vector
%   - (DX2, DY2, DZ2) is a second direction vector
%   into another plane, with the same format, but with :
%   - (x0 y0 z0) is the closest point of plane to origin
%   - (DX1 DY1 DZ1) has norm equal to 1
%   - (DX2 DY2 DZ2) has norm equal to 1 and is orthogonal to (DX1 DY1 DZ1)
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY :

% compute origin point of the plane
p0 = projPointOnPlane([0 0 0], plane1);

% compute first direction vector
d1 = normalize3d(plane1(:,4:6));

% compute second direction vector
n = normalize3d(planeNormal(plane1));
d2 = -normalize3d(cross(d1, n));

% create the resulting plane
plane2 = [p0 d1 d2];

function n = planeNormal(plane)
%PLANENORMAL compute the normal to a plane
%
%   N = planeNormal(PLANE)
%   compute the normal of the given plane
%   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   N : [dx dy dz]
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/02/2005.
%

%   HISTORY

% plane normal
n = cross(plane(:,4:6), plane(:, 7:9), 2);

function pos = planePosition(point, plane)
%PLANEPOSITION compute position of a point on a plane
%
%   PT2 = PLANEPOSITION(POINT, PLANE)
%   POINT has format [X Y Z], and plane has format
%   [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where :
%   - (X0, Y0, Z0) is a point belonging to the plane
%   - (DX1, DY1, DZ1) is a first direction vector
%   - (DX2, DY2, DZ2) is a second direction vector
%
%   Result PT2 has the form [XP YP], with [XP YP] coordinate of the point
%   in the coordinate system of the plane.
%
%
%   CAUTION :
%   WORKS ONLY FOR PLANES WITH ORTHOGONAL DIRECTION VECTORS
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY :
%   24/11/2005 add support for multiple input

% unify size of data
if size(point, 1)~=size(plane, 1)
    if size(point, 1)==1
        point = repmat(point, [size(plane, 1) 1]);
    elseif size(plane, 1)==1
        plane = repmat(plane, [size(point, 1) 1]);
    else
        error('point and plane do not have the same dimension');
    end
end

p0 = plane(:,1:3);
d1 = plane(:,4:6);
d2 = plane(:,7:9);

s = dot(point-p0, d1, 2)./vecnorm3d(d1);
t = dot(point-p0, d2, 2)./vecnorm3d(d2);

pos = [s t];

function point = projPointOnPlane(point, plane)
%PROJPOINTONPLANE return the projection of a point on a plane
%
%   PT2 = PROJECTEDPOINT(PT1, PLANE).
%   Compute the (orthogonal) projection of point PT1 onto the line PLANE.
%
%   Function works also for multiple points and planes. In this case, it
%   returns multiple points.
%   Point PT1 is a [N*3] array, and PLANE is a [N*9] array (see createPlane
%   for details). Result PT2 is a [N*3] array, containing coordinates of
%   orthogonal projections of PT1 onto planes PLANE.
%
%   See also planePosition
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

n = planeNormal(plane);
line = [point repmat(n, [size(point, 1) 1])];
point = intersectPlaneLine(plane, line);

function varargout = sph2cart2(varargin)
%SPH2CART2 convert spherical coordinate to cartesian coordinate
%
%   usage :
%   C = SPH2CART2(S)
%   C = SPH2CART2(PHI, THETA)       (assume rho = 1)
%   C = SPH2CART2(PHI, THETA, RHO)
%   [X, Y, Z] = SPH2CART2(PHI, THETA, RHO);
%
%   S = [phi theta rho] (sphercial coordiante).
%   C = [X Y Z]  (cartesian coordinate)
%
%   Math convention is used : theta is angle with vertical, 0 for north
%   pole, +pi for south pole, pi/2 for points with z=0.
%   phi is the same as matlab cart2sph : angle from Ox axis, counted
%   counter-clockwise.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY
%   22/03/2005 : make test for 2 args, and add radius if not specified for
%       1 arg.

if length(varargin)==1
    var = varargin{1};
    if size(var, 2)==2
        var = [var ones(size(var, 1), 1)];
    end
elseif length(varargin)==2
    var = [varargin{1} varargin{2} ones(size(varargin{1}))];
elseif length(varargin)==3
    var = [varargin{1} varargin{2} varargin{3}];
end

[x y z] = sph2cart(var(:,1), pi/2-var(:,2), var(:,3));

if nargout == 1 || nargout == 0
    varargout{1} = [x, y, z];
else
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end

function alpha = sphericalAngle(p1, p2, p3)
%SPHERICALANGLE compute angle on the sphere
%
%   ALPHA = sphericalAngle(P1, P2, P3)
%   compute angle (P1, P2, P2), in radians, between 0 and 2*PI.
%
%   Points are given either as [x y z] (there will be normalized to lie on
%   the unit sphere), or as [phi theta], with phi being the longitude in [0
%   2*PI] and theta being the elevation on horizontal [-pi/2 pi/2].
%
%
%   NOTE :
%   this is an 'oriented' version of the angle computation, that is, the
%   result of sphericalAngle(P1, P2, P3) equals
%   2*pi-sphericalAngle(P3,P2,P1). To have the more classical relation
%   (with results given betwen 0 and PI), it suffices to take the minimum
%   of angle and 2*pi-angle.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY

% test if points are given as matlab spherical coordinate
if size(p1, 2) ==2
    [x y z] = sph2cart(p1(:,1), p1(:,2));
    p1 = [x y z];
    [x y z] = sph2cart(p2(:,1), p2(:,2));
    p2 = [x y z];
    [x y z] = sph2cart(p3(:,1), p3(:,2));
    p3 = [x y z];
end

% normalize points
p1 = normalize3d(p1);
p2 = normalize3d(p2);
p3 = normalize3d(p3);

% create the plane tangent to the unit sphere and containing central point
plane = createPlane(p2, p2);

% project the two other points on the plane
pi1 = planePosition(intersectPlaneLine(plane, [0 0 0 p1]), plane);
pi3 = planePosition(intersectPlaneLine(plane, [0 0 0 p3]), plane);

% compute angle on the tangent plane
alpha = angle3Points(pi1, [0 0], pi3);

function n = vecnorm3d(v)
%VECNORM3D compute euclidean norm of vector or of set of 3D vectors
%
%   n = vecnorm(V);
%   return euclidean norm of vector V.
%
%   When V is a Nx3 array, compute norm for each vector of the array.
%   Vector are given as rows. Result is then a [N*1] array.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY

n = sqrt(sum(v.*v, 2));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [L] = lap2DPeriodic(N, h);
%
% [L] = lap2DPeriodic(N, h)
%
%  Constructs the 2D Laplacian for a periodic square mesh using
%     the standard 5-point stencil.
%
%  Returns:
%     L = discrete Laplacian
%
%  Input:
%     N = number of mesh points in each direction
%     h = mesh width
%
%
%
%  License: This code is free to use for any purposes, provided
%           any publications resulting from the use of this code
%           reference the original code/author.
%
%  Author:  Samuel Isaacson (isaacson@math.utah.edu)
%  Date:    11/2007
%
%  Please notify the author of any bugs, and contribute any
%  modifications or bug fixes back to the original author.
%
%  Disclaimer:
%   This code is provided as is. The author takes no responsibility
%   for its results or effects.

M = N * N;

% Periodic 2D Laplace Operator on Square:
isPeriodic     = 1;
e              = ones(M,1);
ed             = -4*e;
if( ~isPeriodic )
  ed(N:N:M)      = -3;
  ed(1:N:M)      = -3;
  ed(1)          = -2;
  ed(N)          = -2;
  ed((N-1)*N+1)  = -2;
  ed(N*N)        = -2;
end
eu1            = ones(M,1);
eu1((N+1):N:M) = 0;
ed1            = ones(M,1);
ed1(N:N:M)     = 0;

L = spdiags([e ed1 ed eu1 e], [-N -1:1 N], M, M);

if( isPeriodic )
  for(i = 1:N)
    L( i, N*N-N+i )     = 1;
    L( N*N-N+i, i )     = 1;
    L( (i-1)*N+1, i*N ) = 1;
    L( i*N, (i-1)*N+1 ) = 1;
  end
end

L = L ./ (h*h);

repo: Aero-matlab


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function ind = selectPosesIdx(len,sub,d)

if nargin<3
    d = 2;
end
assert(mod(len,d)==0);
gap = len/d;
ind = [];
sub = sub(:)';
for i = 1:d
    ind = [ind sub + (i-1)*gap];
end

end

author: zhusz


text:
%ISINT Test on number(s) on integer >= 0
% $Id: isint.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function n = isint(m)

	if all(all(m == round(m) & m >= 0.5))
		n = 1;
	else
		n = 0;
	end

	return

author: marianux

repo: ecg-kit

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


text:
function s = resampdet(p,m,n);
%RESAMPDET Deterministic resampling
%
%   Description
%   S = RESAMPDET(P) returns a new set of indices according to the
%   probabilities P. P is array of probabilities, which are not
%   necessarily normalized, though they must be non-negative, and
%   not all zero. The size of S is the size of P.
%
%   S = RESAMPDET(P,M,N) returns M by N matrix.
%
%   S = RESAMPDET(P,M) returns M by M matrix.
%
%   Default is to use no-sort resampling. For sorted resampling use
%    [PS,PI]=SORT(P);
%    S=PI(RESAMPDET(PS));
%   Sorted re-sampling is slower but has smaller variance. Note
%   that deterministic resampling is not unbiased. Stratified
%   resampling (RESAMPSTR) is unbiased, almost as fast as
%   deterministic resampling, and has only slightly larger
%   variance.
%
%   In deterministic resampling indices are sampled using
%   deterministic numbers u_j~(j-a)/n, for fixed a in [0,1) and
%   n is length of P. Compare this to simple random resampling
%   where u_j~U[0,n]. See, Kitagawa, G., Monte Carlo Filter and
%   Smoother for Non-Gaussian Nonlinear State Space Models,
%   Journal of Computational and Graphical Statistics, 5(1):1-25,
%   1996.
%
%   See also RESAMPSIM, RESAMPRES, RESAMPSTR
%
% Copyright (c) 2003-2004 Aki Vehtari

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

if nargin<2
    [m,n]=size(p);
elseif nargin==2
    n=m;
end
mn=m.*n;
pn=p./sum(p(:)).*mn;
fpn=floor(pn);
s=zeros(m,n);
k=0;
c=0.5;
for i=1:numel(p)
  if pn(i)>=1
    a=fpn(i);
    pn(i)=pn(i)-a;
    s(k+[1:a])=i;
    k=k+a;
  end
  c=c+pn(i);
  if c>=1
    k=k+1;
    s(k)=i;
    c=c-1;
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
% POP_INTERP - interpolate data channels
%
% Usage: EEGOUT = pop_interp(EEG, badchans, method, t_range);
%
% Inputs:
%     EEG      - EEGLAB dataset
%     badchans - [integer array] indices of channels to interpolate.
%                For instance, these channels might be bad.
%                [chanlocs structure] channel location structure containing
%                either locations of channels to interpolate or a full
%                channel structure (missing channels in the current
%                dataset are interpolated).
%     method   - [string] method used for interpolation (default is 'spherical').
%                'invdist'/'v4' uses inverse distance on the scalp
%                'spherical' uses superfast spherical interpolation.
%                'spacetime' uses griddata3 to interpolate both in space
%                and time (very slow and cannot be interrupted).
%     t_range  - [integer array with just two elements] time interval of the
%                badchans which should be interpolated. First element is
%                the start time and the second element is the end time.
% Output:
%     EEGOUT   - data set with bad electrode data replaced by
%                interpolated data
%
% Author: Arnaud Delorme, CERCO, CNRS, 2009-

% Copyright (C) Arnaud Delorme, CERCO, 2009, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG com] = pop_interp(EEG, bad_elec, method, t_range)

    com = '';
    t_range = '';
    if nargin < 1
        help pop_interp;
        return;
    end
    if nargin > 1 && nargin < 3
        method = 'spherical';
    end

    if nargin < 2
        disp('Warning: interpolation can be done on the fly in studies');
        disp('         this function will actually create channels in the dataset');
        disp('Warning: do not interpolate channels before running ICA');
        disp('You may define channel location to interpolate in the channel');
        disp('editor and declare such channels as non-data channels');

        enablenondat = 'off';
        if isfield(EEG.chaninfo, 'removedchans')
            if ~isempty(EEG.chaninfo.removedchans)
                enablenondat = 'on';
            end
        end
        if isempty(EEG.epoch)
            uilist = { { 'Style' 'text' 'string' 'What channel(s) do you want to interpolate' 'fontweight' 'bold' } ...
                       { 'style' 'text' 'string' 'none selected' 'tag' 'chanlist' } ...
                       { 'style' 'pushbutton' 'string' 'Select from removed channels' 'callback' 'pop_interp(''nondatchan'',gcbf);' 'enable' enablenondat } ...
                       { 'style' 'pushbutton' 'string' 'Select from data channels'    'callback' 'pop_interp(''datchan'',gcbf);' } ...
                       { 'style' 'pushbutton' 'string' 'Use specific channels of other dataset' 'callback' 'pop_interp(''selectchan'',gcbf);'} ...
                       { 'style' 'pushbutton' 'string' 'Use all channels from other dataset' 'callback' 'pop_interp(''uselist'',gcbf);'} ...
                       { } ...
                       { 'style' 'text'  'string' 'Interpolation method'} ...
                       { 'style' 'popupmenu'  'string' 'Spherical|Planar (slow)'  'tag' 'method' } ...
                       { } ...
                       { 'Style', 'text', 'string', 'Time range [min max] (s)', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') } ...
                       { 'Style', 'edit', 'string', '', 'enable', fastif(length(EEG)>1, 'off', 'on') } ...
                       {} { 'Style' 'text' 'string' 'Note: for group level analysis, interpolate in STUDY' } ...
                       };

            geom     = { 1 1 1 1 1 1 1 [1.1 1] 1 [1.1 1] 1   1 };
            geomvert = [ 1 1 1 1 1 1 1 1       0.5 1      0.5  1 ];
        else
            uilist = { { 'Style' 'text' 'string' 'What channel(s) do you want to interpolate' 'fontweight' 'bold' } ...
                   { 'style' 'text' 'string' 'none selected' 'tag' 'chanlist' } ...
                   { 'style' 'pushbutton' 'string' 'Select from removed channels' 'callback' 'pop_interp(''removedchans'',gcbf);' 'enable' enablenondat } ...
                   { 'style' 'pushbutton' 'string' 'Select from data channels'    'callback' 'pop_interp(''datchan'',gcbf);' } ...
                   { 'style' 'pushbutton' 'string' 'Use specific channels of other dataset' 'callback' 'pop_interp(''selectchan'',gcbf);'} ...
                   { 'style' 'pushbutton' 'string' 'Use all channels from other dataset' 'callback' 'pop_interp(''uselist'',gcbf);'} ...
                   { } ...
                   { 'style' 'text'  'string' 'Interpolation method'} ...
                   { 'style' 'popupmenu'  'string' 'Spherical|Planar (slow)'  'tag' 'method' } ...
                   {} { 'Style' 'text' 'string' 'Note: for group level analysis, interpolate in STUDY' } ...
                   };

            geom     = { 1 1 1 1 1 1 1 [1.1 1] 1   1 };
            geomvert = [ 1 1 1 1 1 1 1 1       0.5 1 ];
        end
        [res, userdata, ~, restag ] = inputgui( 'uilist', uilist, 'title', 'Interpolate channel(s) -- pop_interp()', 'geometry', geom, 'geomvert', geomvert, 'helpcom', 'pophelp(''pop_interp'')');
        if isempty(res) || isempty(userdata), return; end

        if restag.method == 1
             method = 'spherical';
        else method = 'invdist';
        end
        bad_elec = userdata.chans;
        if nargin < 4
            if numel(res) > 1
                t_range = res{2};
            else
                t_range = '';
            end
        end
        if isempty(t_range)
            t_range = [EEG.xmin EEG.xmax];
        else
            t_range = eval( [ '[' t_range ']' ] );
        end
        if size(t_range,2) ~= 2
            error('Time/point range must contain 2 columns exactly');
        end
        if floor(max(t_range)) > EEG.xmax
            error('Time/point range exceed upper data limits');
        end
        if min(t_range) < EEG.xmin
            error('Time/point range exceed lower data limits');
        end

        com = sprintf('EEG = pop_interp(EEG, %s, ''%s'');', userdata.chanstr, method);
        if ~isempty(findstr('removedchans', userdata.chanstr))
            eval( [ userdata.chanstr '=[];' ] );
        end

    elseif ischar(EEG)
        command = EEG;
        clear EEG;
        fig = bad_elec;
        userdata = get(fig, 'userdata');

        if strcmpi(command, 'removedchans')
            global EEG;
            tmpchaninfo = EEG.chaninfo;
            [chanlisttmp, chanliststr] = pop_chansel( { tmpchaninfo.removedchans.labels } );
            if ~isempty(chanlisttmp)
                userdata.chans   = EEG.chaninfo.removedchans(chanlisttmp);
                userdata.chanstr = [ 'EEG.chaninfo.removedchans([' num2str(chanlisttmp) '])' ];
                set(fig, 'userdata', userdata);
                set(findobj(fig, 'tag', 'chanlist'), 'string', chanliststr);
            end
        elseif strcmpi(command, 'datchan')
            global EEG;
            tmpchaninfo = EEG.chanlocs;
            [chanlisttmp, chanliststr] = pop_chansel( { tmpchaninfo.labels } );
            if ~isempty(chanlisttmp)
                userdata.chans   = chanlisttmp;
                userdata.chanstr = [ '[' num2str(chanlisttmp) ']' ];
                set(fig, 'userdata', userdata);
                set(findobj(fig, 'tag', 'chanlist'), 'string', chanliststr);
            end
        else
            global ALLEEG EEG;
            tmpanswer = inputdlg2({ 'Dataset index' }, 'Choose dataset', 1, { '' });
            if ~isempty(tmpanswer)
                tmpanswernum = round(str2num(tmpanswer{1}));
                if ~isempty(tmpanswernum)
                    if tmpanswernum > 0 && tmpanswernum <= length(ALLEEG)
                        TMPEEG = ALLEEG(tmpanswernum);

                        tmpchans1 = TMPEEG.chanlocs;
                        if strcmpi(command, 'selectchan')
                            chanlist = pop_chansel( { tmpchans1.labels } );
                        else
                            chanlist = 1:length(TMPEEG.chanlocs); % use all channels
                        end

                        % look at what new channels are selected
                        tmpchans2 = EEG.chanlocs;
                        [tmpchanlist, chaninds] = setdiff_bc( { tmpchans1(chanlist).labels }, { tmpchans2.labels } );
                        if ~isempty(tmpchanlist)
                            if length(chanlist) == length(TMPEEG.chanlocs)
                                userdata.chans   = TMPEEG.chanlocs;
                                userdata.chanstr = [ 'ALLEEG(' tmpanswer{1} ').chanlocs' ];
                            else
                                userdata.chans   = TMPEEG.chanlocs(chanlist(sort(chaninds)));
                                userdata.chanstr = [ 'ALLEEG(' tmpanswer{1} ').chanlocs([' num2str(chanlist(sort(chaninds))) '])' ];
                            end
                            set(fig, 'userdata', userdata);
                            tmpchanlist(2,:) = { ' ' };
                            set(findobj(gcbf, 'tag', 'chanlist'), 'string', [ tmpchanlist{:} ]);
                        else
                            warndlg2('No new channels selected');
                        end
                    else
                        warndlg2('Wrong index');
                    end
                end
            end
        end
        return;
    end

    % remove from removedchans if interpolated
    if isfield(EEG.chaninfo, 'removedchans') && ~isempty(EEG.chaninfo.removedchans) && isstruct(bad_elec)
        for iChan = 1:length(bad_elec)
            ind = strmatch(lower(bad_elec(iChan).labels), lower({EEG.chaninfo.removedchans.labels}), 'exact');
            if ~isempty(ind)
                EEG.chaninfo.removedchans(ind) = [];
            end
        end
    end

    EEG = eeg_interp(EEG, bad_elec, method, t_range);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [xout,yout]=rejectsample(f, n)
% [xout,yout]=rejectsample(f, n)
% rejection sampling from the 2D discrete distribution f
% n: number of samples
% the proposal distriburion is a uniform distribution on the range of f
sizef=size(f);
maxf=max(f(:));
nsofar=0;
while nsofar<n
    xrand=sizef(1)*rand;
    yrand=sizef(2)*rand;
    u=maxf*rand;
    if u<f(ceil(xrand), ceil(yrand))
        nsofar=nsofar+1;
        x(nsofar)=xrand;
        y(nsofar)=yrand;
    end
end

%to the middle of the pixels...
xout=x-0.5;
yout=y-0.5;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ p, l, u ] = maxij_plu ( n )

%*****************************************************************************80
%
%% MAXIJ_PLU returns the PLU factors of the MAXIJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real P(N,N), L(N,N), U(N,N), the PLU factors.
%
  p = zeros ( n, n );

  for j = 1 : n
    for i = 1 : n
      if ( i4_wrap ( j - i, 1, n ) == 1 )
        p(i,j) = 1.0;
      end
    end
  end

  l = zeros ( n, n );

  i = 1;
  j = 1;
  l(i,j) = 1.0;

  j = 1;
  for i = 2 : n
    l(i,j) = ( i - 1 ) / n;
  end

  for j = 2 : n
    l(j,j) = 1.0;
  end

  u = zeros ( n, n );

  i = 1;
  for j = 1 : n
    u(i,j) = n;
  end

  for i = 2 : n
    for j = i : n
      u(i,j) = j + 1 - i;
    end
  end

  return
end


text:
function varargout = help1(varargin)
% HELP1 M-file for help1.fig
%      HELP1, by itself, creates a new HELP1 or raises the existing
%      singleton*.
%
%      H = HELP1 returns the handle to a new HELP1 or the handle to
%      the existing singleton*.
%
%      HELP1('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in HELP1.M with the given input arguments.
%
%      HELP1('Property','Value',...) creates a new HELP1 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before help1_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to help1_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

% Copyright 2002-2003 The MathWorks, Inc.

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

% Last Modified by GUIDE v2.5 06-Feb-2008 22:05:36

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @help1_OpeningFcn, ...
                   'gui_OutputFcn',  @help1_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 help1 is made visible.
function help1_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 help1 (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = help1_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)


text:
function test_quadratic_constrained
%% Tests the solvers on a simple constrained quadratic function
%{
    Solve:

    minimize_x  c'x + x'Dx/2
    subject to  ||x||_1 <= 10

    as an example of using TFOCS without the "SCD" interface
    (since the objective is smooth and we can project, there's
     no need to smooth and solve the dual)

    It's also an example of using 2 objective functions

%}

%% Now, add in constraints

randn( 'state', sum('quadratic test') );
N = 100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CONSTRUCT THE TEST PROBLEM %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c = randn(N,1);
D = randn(N,N);
D = D * D' + .5*eye(N);
s = svd(D);
% x_star = - D \ c;
% L = max(s);
% mu = min(s);
% if mu < eps, disp('WARNING: may have 0 eigenvalues'); end
% f_star = 0.5 * c' * x_star;
% n_star = norm( x_star );
x0 = zeros(N,1);

% f       = @(x) c'*x + x'*D*x/2;
% grad_f  = @(x) c + D*x;
% smoothF = @(x) wrapper_objective( f, grad_f, x );
%   -- or --
% smoothF = smooth_quad(D,c);
%   -- or (demonstrating how to use 2 objective functions) --
f1      = @(x) c'*x;
f2      = @(x) x'*D*x/2;
g1      = @(x) c;
g2      = @(x) D*x;
smoothF = { @(x) wrapper_objective(f1,g1,x); @(x) wrapper_objective(f2,g2,x) };
linearF = { 1 ; 1 };

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET UP THE TEST PARAMETERS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opts = [];
opts.tol        = 1e-16;
opts.maxits     = 3000;
opts.restart    = 100;

%%%%%%%%%%%%%%%%%%%%%
% ADD IN CONSTRAINT %
%%%%%%%%%%%%%%%%%%%%%
% projectorF      = proj_simplex(10); % x >=0, sum(x) = 10
projectorF      = proj_l1(10);  % l1 ball, radius 10 (e.g. norm(x,1) <= 10)

%%%%%%%%%%%%%%%%%
% RUN THE TESTS %
%%%%%%%%%%%%%%%%%
[ x, out, optsOut ] = tfocs( smoothF, linearF, projectorF,x0, opts );
% Check that we are within allowable bounds

function [v,gr] = wrapper_objective(f,g,x)
v = f(x);
if nargout > 1
    gr = g(x);
end
% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [score]=crps(ygibbs,yo)

% compute first the dimension of the gibbs sampler draws (normally, 'It' iterations minus 'Bu' discarded burn iterations)
It_Bu=size(ygibbs,1);

% compute the first summation term
sum1=sum(abs(ygibbs-repmat(yo,It_Bu,1)));

% compute the second summation term
 temp=abs(repmat(ygibbs,1,It_Bu)-repmat(ygibbs',It_Bu,1));
 sum2=sum(temp(:));

% eventually compute the score
score=(1/It_Bu)*sum1-(1/(2*It_Bu^2))*sum2;

author: european-central-bank

repo: BEAR-toolbox


text:
function data = prepare_batch_feat_input_data(pooler, feature_paths, windows, data_param)
num_windows   = size(windows,1);
feat_dim      = data_param.feat_dim;
data          = zeros(feat_dim, num_windows, 'single');
num_threads   = 6;

if isfield(data_param, 'num_threads') && data_param.num_threads > 0
    num_threads = data_param.num_threads;
end

select_img_idx  = unique(windows(:,1));
num_sel_img     = length(select_img_idx);
feature_paths   = feature_paths(select_img_idx);
window_this_img = cell(num_sel_img, 1);
img_window_map  = cell(num_sel_img, 1);
for i = 1:num_sel_img
    indices = find(windows(:,1) == select_img_idx(i));
    window_this_img{i} = windows(indices, :);
    img_window_map{i}  = indices(:);
end
assert(sum(cellfun(@(x) size(x, 1), window_this_img, 'UniformOutput', true)) == num_windows);
assert(sum(cellfun(@(x) size(x, 1), img_window_map, 'UniformOutput', true)) == num_windows);

random_scale = data_param.random_scale;

num_blocks   = 3;
block_size   = ceil(num_sel_img / num_blocks);
for b = 1:num_blocks
    start_idx       = (b-1) * block_size + 1;
    stop_idx        = min(b*block_size,num_sel_img);
    this_block_size = stop_idx - start_idx + 1;

    block_sel_img_idx     = select_img_idx(start_idx:stop_idx);
    block_window_this_img = window_this_img(start_idx:stop_idx);
    block_feature_paths   = feature_paths(start_idx:stop_idx);

    block_feats    = cell(1,this_block_size);
    block_data_map = cell2mat(img_window_map(start_idx:stop_idx));

    assert(all(block_data_map > 0 & block_data_map <= num_windows));
    assert(length(unique(block_data_map)) == length(block_data_map));
    if num_threads == 1
        for i = 1:this_block_size
            assert(all(block_window_this_img{i}(:,1) == block_sel_img_idx(i)));
            %feat_cache = load(block_feature_paths{i}, 'feat');
            feat_cache  = read_feat_conv_data( block_feature_paths{i}, true );
            block_feats{i} = cell2mat(convFeat_to_poolFeat_multi_region(pooler, ...
                feat_cache.feat, block_window_this_img{i}(:, 4:7), random_scale));
        end
    else
        parfor (i = 1:this_block_size,num_threads)
            assert(all(block_window_this_img{i}(:,1) == block_sel_img_idx(i)));
            feat_cache  = read_feat_conv_data( block_feature_paths{i}, true );
            block_feats{i} = cell2mat(convFeat_to_poolFeat_multi_region(pooler, ...
                feat_cache.feat, block_window_this_img{i}(:, 4:7), random_scale));
        end
    end
    data(:,block_data_map) = cell2mat(block_feats);
    for i = 1:this_block_size, block_feats{i} = []; end
end
end

author: gidariss

repo: mrcnn-object-detection


text:
function [p] = config_hog3x3(c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

p.grid_spacing = 6; % distance between grid centers
p.patch_size = 24;
p.Mw = 2; % number of spatial scales for spatial pyramid histogram
p.descriptor = 'hog';
p.w = p.patch_size; % boundary for HOG

p.pyramid_levels = 2;
p.llcknn = 3;
p.maxsize = 400;

% dictionary parameters
p.dictionary_size = 256;
p.num_images = 500;
p.descPerImage = 2000;
p.num_desc = 500000;

p.train_file = '%s/train_hog3x3_%d.mat';
p.test_file = '%s/test_hog3x3_%d.mat';
p.dictionary_file = '%s/dictionary_hog3x3_%d.mat';

author: adikhosla

repo: feature-extraction


text:
function [posterior, out] = VBA_multisession_factor(posterior,out)
% [posterior, out] = VBA_multisession_factor(posterior,out)
% take a multisession model, extract within session parameters and
% trajectories (in posterior.perSession), and clean up the out structure
% by undoing the multisession wrapping performed by VBA_multisession_expand.

%% skip if no multisession
if ~(isfield(out.options,'multisession') && isfield(out.options.multisession,'expanded') && out.options.multisession.expanded)
    return;
end

%% extract session_wise posteriors
Ts = [0 cumsum(out.options.multisession.split)];

multisession = out.options.inF{2};
out.options.inF = out.options.inF{1:end-1};
out.options.inG = out.options.inG{1:end-1};

out.options.f_fname = multisession.f_fname;
out.options.g_fname = multisession.g_fname;

for i=1:numel(out.options.multisession.split)

    idx_X0 = multisession.indices.X0(:,i);
    idx_theta = multisession.indices.theta(:,i);
    idx_phi = multisession.indices.phi(:,i);

    posteriors(i).muX0 = posterior.muX0(idx_X0);
    posteriors(i).SigmaX0 = posterior.SigmaX0(idx_X0,idx_X0);

    posteriors(i).muTheta = posterior.muTheta(idx_theta);
    posteriors(i).SigmaTheta = posterior.SigmaTheta(idx_theta,idx_theta);

    posteriors(i).muPhi = posterior.muPhi(idx_phi);
    posteriors(i).SigmaPhi = posterior.SigmaPhi(idx_phi,idx_phi);

    if isfield(posterior,'a_sigma')
        posteriors(i).a_sigma =  posterior.a_sigma ;
        posteriors(i).b_sigma =  posterior.b_sigma ;
    end
    if isfield(posterior,'a_alpha')
        posteriors(i).a_alpha =  posterior.a_alpha ;
        posteriors(i).b_alpha =  posterior.b_alpha ;
    end

    idx_t = (Ts(i)+1):Ts(i+1);

    posteriors(i).iQy = posterior.iQy(idx_t,:);

    try
    posteriors(i).iQx =  posterior.iQx(idx_t);%cellfun(@(Q) Q(idx_X0,idx_X0),posterior.iQx(idx_t),'UniformOutput',false);
    posteriors(i).muX = posterior.muX(idx_X0,idx_t);
    posteriors(i).SigmaX.current =  cellfun(@(Q) Q(idx_X0,idx_X0),posterior.SigmaX.current(idx_t,:),'UniformOutput',false);
    posteriors(i).SigmaX.inter =  cellfun(@(Q) Q(idx_X0,idx_X0),posterior.SigmaX.inter(:,idx_t),'UniformOutput',false);
    end
end

posterior.perSession = posteriors;

out.options.multisession.expanded = false ;

end

author: MBB-team

repo: VBA-toolbox


text:
function r8mat_det_5d_test ( )

%*****************************************************************************80
%
%% R8MAT_DET_5D_TEST tests R8MAT_DET_5D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  x = [ 1.0, 10.0, 4.0, 2.0, 3.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_DET_5D_TEST\n' );
  fprintf ( 1, '  R8MAT_DET_5D: determinant of 5 by 5 matrix.\n' );

  a = r8mat_vand2 ( n, x );
  det = r8mat_det_5d ( a );

  r8mat_print ( n, n, a, '  Matrix:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R8MAT_DET_5D computes determinant: %f\n', det );
%
%  Special formula for the determinant of a Vandermonde matrix:
%
  det = 1.0;
  for i = 1 : n
    for j = 1 : i-1
      det = det * ( x(i) - x(j) );
    end
  end

  fprintf ( 1, '  Exact determinant is %f\n', det );

  return
end


text:
function f = allenAtlasBrowser_original(templateVolume, annotationVolume, structureTree)
% Browser for the allen atlas ccf data in matlab.
%
% Inputs templateVolume, annotationVolume, and structureTree are the data describing the atlas.
% The annotation volume should be the "by_index" version
%
% Features:
% - Use scroll wheel to go between slices
% - Apply 1mm-spacing grid lines ('g' to enable/disable)
% - View overlay of the region hovered over ('o' to enable/disable)
% - Switch to viewing the atlas annotations ('v' to toggle)
% - Click to register points in 3D ('p' to toggle)
%   - use:
%     >> ud = get(f, 'UserData'); myPoints = ud.pointList;
%     to retrieve these points

% todo/future features:
% - switch to sagittal/horizontal slice types

fprintf(1, 'Controls: \n');
fprintf(1, '--------- \n');
fprintf(1, 'scroll: move between slices \n');
fprintf(1, 'g: add/remove Gridlines \n');
fprintf(1, 'o: add/remove Overlay of current region extent \n');
fprintf(1, 'v: switch to Viewing annotations (or switch back) \n');
fprintf(1, 'p: enable/disable mode where clicks are logged as Points \n');

f = figure;

ud.bregma = allenCCFbregma;

ud.currentSlice = min(ud.bregma(1), size(templateVolume, 1));
ud.oldContour = [];
ud.showContour = false;
ud.showOverlay = false; ud.overlayAx = [];
ud.getPoint = false; ud.pointList = zeros(0,3); ud.pointHands = [];
ud.showAtlas = false;

ud.atlasAx = axes('Position', [0.05 0.05 0.9 0.9]);

ud.im = plotTVslice(squeeze(templateVolume(ud.currentSlice,:,:)));

set(ud.im, 'ButtonDownFcn', @(f,k)atlasClickCallback(f, k));

ud.bregmaText = annotation('textbox', [0 0.95 0.4 0.05], ...
    'String', '[coords]', 'EdgeColor', 'none', 'Color', 'k');

allData.tv = templateVolume;
allData.av = annotationVolume;
allData.st = structureTree;

hold(ud.atlasAx, 'on');
set(ud.atlasAx, 'HitTest', 'off');

set(f, 'UserData', ud);

set(f, 'KeyPressFcn', @(f,k)hotkeyFcn(f, k, allData));
set(f, 'WindowScrollWheelFcn', @(src,evt)updateSlice(f, evt, allData))
set(f, 'WindowButtonMotionFcn',@(f,k)fh_wbmfcn(f, allData)); % Set the motion detector.

function hotkeyFcn(f, keydata, allData)

ud = get(f, 'UserData');
switch lower(keydata.Key)
    case 'g' % toggle showing Gridlines
        if ~isfield(ud, 'gridlines') || isempty(ud.gridlines)
            axes(ud.atlasAx); hold on;
            gridY = 100:100:700; % assuming the size of the atlas for this for now
            gridX = 70:100:1140;
            xl = xlim(); yl = ylim();
            gx = arrayfun(@(x)plot(x*[1 1], yl, 'w'), gridX, 'uni', false);
            gy = arrayfun(@(y)plot(xl, y*[1 1], 'w'), gridY, 'uni', false);
            ud.gridlines = [gx gy];
        elseif strcmp(get(ud.gridlines{1}, 'Visible'), 'on');
            cellfun(@(x)set(x, 'Visible', 'off'), ud.gridlines);
        elseif strcmp(get(ud.gridlines{1}, 'Visible'), 'off');
            cellfun(@(x)set(x, 'Visible', 'on'), ud.gridlines);
        end
    case 'c' % toggle Contour
        % contours are disabled for now - they are slow and buggy...
%         if ud.showContour
%             delete(ud.oldContour); ud.oldContour = [];
%         end
%         ud.showContour = ~ud.showContour;
    case 'o' % toggle Overlay
        ud.showOverlay = ~ud.showOverlay;
        if ~ud.showOverlay
            % weird bug means that point clicking stops working when the
            % overlayAx exists, even though I turned all the hittests off??
            % hack solution is to just delete it every time... :/
%             ovIm = get(ud.overlayAx, 'Children');
%             set(ovIm, 'AlphaData', zeros(size(get(ovIm, 'CData'))));
            delete(ud.overlayAx); ud.overlayAx = [];
        end
    case 'p' % toggle mode to register clicks as Points
        ud.getPoint = ~ud.getPoint;
    case 'v' % toggle View between template/annotation
        ud.showAtlas = ~ud.showAtlas;
        if ud.showAtlas
            set(ud.im, 'CData', squeeze(allData.av(ud.currentSlice,:,:)));
            cmap = allen_ccf_colormap('2017');
            colormap(ud.atlasAx, cmap); caxis(ud.atlasAx, [1 size(cmap,1)]);
        else
            set(ud.im, 'CData', squeeze(allData.tv(ud.currentSlice,:,:)));
            colormap(ud.atlasAx, 'gray'); caxis(ud.atlasAx, [0 400]);
        end
        set(f, 'UserData', ud);

end
set(f, 'UserData', ud);

function updateSlice(f, evt, allData)

ud = get(f, 'UserData');
ud.currentSlice = ud.currentSlice+evt.VerticalScrollCount*3;
if ud.currentSlice>size(allData.tv,1); ud.currentSlice = 1; end %wrap around
if ud.currentSlice<1; ud.currentSlice = size(allData.tv,1); end %wrap around
set(f, 'UserData', ud);

if ud.showAtlas
    set(ud.im, 'CData', squeeze(allData.av(ud.currentSlice,:,:)));
else
    set(ud.im, 'CData', squeeze(allData.tv(ud.currentSlice,:,:)));
end

if ~isempty(ud.pointHands)
    % any points that are at this AP position, make them visible
    arrayfun(@(x)pointVis(ud.pointHands(x), ud.currentSlice, ud.pointList(x,:)), 1:numel(ud.pointHands));
end

pixel = getPixel(ud.atlasAx);
updateStereotaxCoords(ud.currentSlice, pixel, ud.bregma, ud.bregmaText);

[name, acr, ann] = getPixelAnnotation(allData, pixel, ud.currentSlice);

updateTitle(ud.atlasAx, name, acr)

if ud.showOverlay
    updateOverlay(f, allData, ann);
end

function pointVis(h, sliceInd, coord)
    if sliceInd==coord(3)
        set(h, 'Visible', 'on');
    else
        set(h, 'Visible', 'off');
    end

function atlasClickCallback(im, keydata)
f = get(get(im, 'Parent'), 'Parent');
ud = get(f, 'UserData');

if ud.getPoint
    clickX = round(keydata.IntersectionPoint(1));
    clickY = round(keydata.IntersectionPoint(2));
    clickZ = ud.currentSlice;

    ud.pointList(end+1, :) = [clickX, clickY, clickZ];
    ud.pointHands(end+1) = plot(ud.atlasAx, clickX, clickY, 'ro');
end
set(f, 'UserData', ud);

function fh_wbmfcn(f, allData)
% WindowButtonMotionFcn for the figure.

ud = get(f, 'UserData');

ax = ud.atlasAx;

pixel = getPixel(ax);

% show bregma coords
updateStereotaxCoords(ud.currentSlice, pixel, ud.bregma, ud.bregmaText);

% get annotation for this pixel
[name, acr, ann] = getPixelAnnotation(allData, pixel, ud.currentSlice);

updateTitle(ax, name, acr);

if ~isempty(name)
    if ud.showContour
        if ~isempty(ud.oldContour)
            delete(ud.oldContour);
        end
        [~,ch] = contour(squeeze(allData.av(ud.currentSlice,:,:)==ann), 1, 'r');
        ud.oldContour = ch;
        set(f, 'UserData', ud);
    end

    if ud.showOverlay
        updateOverlay(f, allData, ann)
    end
end

function updateStereotaxCoords(currentSlice, pixel, bregma, bregmaText)
atlasRes = 0.010; % mm
ap = -(currentSlice-bregma(1))*atlasRes;
dv = (pixel(1)-bregma(2))*atlasRes;
ml = (pixel(2)-bregma(3))*atlasRes;
set(bregmaText, 'String', sprintf('%.2f AP, %.2f DV, %.2f ML', ap, dv, ml));

function pixel = getPixel(ax)

currPoint = get(ax,'currentpoint');  % The current point w.r.t the axis.

Cx = currPoint(1,1); Cy = currPoint(1,2);
pixel = round([Cy Cx]);

function updateOverlay(f, allData, ann)
ud = get(f, 'UserData');
if isempty(ud.overlayAx) % first time
    avo = plotAVoverlay(fliplr(squeeze(allData.av(ud.currentSlice,:,:))'), ann, ud.atlasAx);
    ud.overlayAx = avo;
    set(ud.overlayAx, 'HitTest', 'off');
    set(f, 'UserData', ud);
else
    ovIm = get(ud.overlayAx, 'Children');
    set(ovIm, 'HitTest', 'off');
    thisSlice = squeeze(allData.av(ud.currentSlice,:,:));
    if ud.showAtlas
        % if showing the atlas annotations, make the cdata 1302 so you get
        % a black overlay (it's the number for the parafloccular sulcus, which
        % is unused) in current annotations so I stole its entry in the
        % colormap to be black.
        adat = ones(size(get(ovIm, 'CData')));
        set(ovIm, 'CData', 1302*adat);
        adat(thisSlice~=ann) = 0;
        set(ovIm, 'AlphaData', flipud(100*adat));
    else
        % otherwise show in the color of the atlas
        set(ovIm, 'CData', flipud(thisSlice));
        plotAVoverlay(fliplr(thisSlice'), ann, ud.atlasAx, ud.overlayAx);
    end

end

function [name, acr, ann] = getPixelAnnotation(allData, pixel, currentSlice)
if pixel(1)>0&&pixel(1)<size(allData.av,2) && pixel(2)>0&&pixel(2)<=size(allData.av,3)
    ann = allData.av(currentSlice,pixel(1),pixel(2));
    name = allData.st.safe_name(ann);
    acr = allData.st.acronym(ann);
else
    ann = []; name = []; acr = [];
end

function updateTitle(ax, name, acr)
if ~isempty(name)
    title(ax, [name{1} ' (' acr{1} ')']);
else
    title(ax, 'not found');
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function obj = segmentation(obj, r, varargin)
% SEGMENTATION Segmentation of point cloud in plane elements.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% * This method adds an attribute 'segId' for each point, which contains the
%   segment id. Two segment ids have a special meaning:
%   * segId = 0 -> for isolated points (no neighbours within radius r)
%   * segId = 1 -> for points from small segments (where the no. of points is
%                  smaller than minNoPoints)
% * This segmentation algorithm is modified from [1].
% * This algorithm needs the normal vector for each point (attributes nx, ny,
%   nz). The normal vectors can be calculated with the method 'normals' (for
%   help run 'help pointCloud.normals').
% ------------------------------------------------------------------------------
% INPUT
% 1 [r]
%   Radius for segment growing.
%
% 2 ['dAngleMax', dAngleMax]
%   A point is only added to an existing segment, if the difference between its
%   normal and the segment normal is smaller than dAngleMax.
%
% 3 ['MinNoPoints', minNoPoints]
%   Minimum number of points in a segments. To segments with a lower number of
%   points the segmentation id 0 is assigned.
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Segment point cloud.
%   pc = pointCloud('Lion.xyz');
%   pc.segmentation(3, 5);
% ------------------------------------------------------------------------------
% REFERENCES
% [1] Rabbani T., 2006: Automatic Reconstruction of Industrial Installations
%     Using Point Clouds and Images
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired( 'r'          , @(x) isscalar(x) && x>0);
p.addParameter('dAngleMax'  , 10, @(x) isscalar(x) && x>0);
p.addParameter('MinNoPoints', 100, @(x) isscalar(x) && x>0);
p.parse(r, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear r

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'SEGMENTATION'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));

% Preparations -----------------------------------------------------------------

% Indices of activated points
idxAct = find(obj.act);

X = obj.X(idxAct,:);
N = [obj.A.nx(idxAct) obj.A.ny(idxAct) obj.A.nz(idxAct)];
segId = NaN(numel(idxAct),1);

% Search neighbors of each point -----------------------------------------------

msg('S', {procHierarchy{:} 'NNSEARCH'});

% NN of each point
[idxNN, dist] = knnsearch(X, X, 'K', 2);

idxNN = idxNN(:,2);
dist  = dist(:,2);

msg('E', {procHierarchy{:} 'NNSEARCH'});

% Find segments ----------------------------------------------------------------

segId = findSegments(X, N, segId, idxNN, dist, p);

% Report results ---------------------------------------------------------------

for i = 0:max(obj.A.segId)
    switch i
        case 0
            addInfo = '(isolated points)';
        case 1
            addInfo = '(too small segments)';
        otherwise
            addInfo = '';
    end

    msg('V', sum(obj.A.segId == i), sprintf('number of points in segment with segId=%d %s', i, addInfo), 'Prec', 0);
end

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

end

function percent = percentSegmentedPoints(obj)

    percent = sum(~isnan(obj.A.segId)) / obj.noPoints * 100;

end

author: pglira


text:
% 机动车刹车距离
%% 插值数据处理
clear; clc;
v = 20 : 10 : 150;  % 速度，插值节点，单位：千米/小时
d2=[3.15,7.08,12.59,19.68,28.34,38.57,50.4,63.75,...
    78.71,95.22,113.29,132.93,154.12,176.87];  % 制动距离，插值节点
vi = 20 : 1 : 150;  % 速度，待插值点
d2i = interp1(v, d2, vi, 'spline');  % 样条插值得到的制动距离

%% 根据某驾驶员的实际视力和视觉习惯，其驾驶时的有效视距为120m，则其在该路面行车时，时速最高不能超过多少(结果取整)?
dl = 120;  % 某路段的有效视距，单位：米
time = 10;  % 驾驶员的平均反应时间，单位：秒
vs = vi * 1000 / 3600; % 速度，插值节点，单位：米/秒
d1 = vs * time; % 驾驶员反应距离
d3 = 10;  % 预留安全距离，单位：米
di = d1 + d2i + d3;  % 有效视距 = 驾驶员反应距离 + 制动距离 + 预留安全距离
x = abs(di - dl);   % 有效视距差的绝对值
[y,i]=sort(x);    % 排序，i是排序后元素的序号
fprintf('有效视距%.4f米某路段的最高时速: %.f千米/小时\n',dl,vi(i(1)));
plot(vi,di,vi(i(1)),di(i(1)),'rp');

%% 若以表中数据为参考，设计一条最高时速为125km/h的高速公路，则设计人员应该保 证驾驶者在公路上任一点的可视距离为多少米?
v_max = 125; % 最高限速，单位：千米/小时
i = find(vi - 125 == 0);
fprintf('最高时速%.4f千米/小时的路段需保证有效视距: %.4f米\n', v_max, di(i));

repo: NumericalAnalysis


text:
% SQRT_COV  -  Factorizes a positive semi-definite symmetric matrix.
%
% V = SQRT_COV(A)
%
% The matrix A must be positive semi-definite and symmetric. The resulting
% matrix V is such that V*V'=A.

% Last modified 2011-10-20
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@aalto.fi)

function V = sqrt_cov(A)

[L,D] = ldl(A);
D(D<0) = 0;
V = L * sqrt(D);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function styleLabelAsError(hLabel, errorStr)
%styleLabelAsError Summary of this function goes here
%   Detailed explanation goes here
    lblStr = strWrapForLabel(hLabel, errorStr);

    set(hLabel,'Visible','on');
    set(hLabel,'BackgroundColor',[0.847,0.161,0]);
    set(hLabel,'ForegroundColor',[1,1,1]);
    set(hLabel,'String',lblStr);
    set(hLabel,'TooltipString','');
    set(hLabel,'FontWeight','bold');

    setappdata(hLabel,'tooltipstring',errorStr);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef TestImgHash
    %TestImgHash

    methods (Static)
        function test_1
            img1 = imread(fullfile(mexopencv.root(), 'test', 'tsukuba_l.png'));
            img2 = imread(fullfile(mexopencv.root(), 'test', 'tsukuba_r.png'));

            algs = {'AverageHash', 'BlockMeanHash', 'ColorMomentHash',
                'MarrHildrethHash', 'PHash', 'RadialVarianceHash'};
            for i=1:numel(algs)
                obj = cv.ImgHash(algs{i});
                hash1 = obj.compute(img1);
                hash2 = obj.compute(img2);
                validateattributes(hash1, {'numeric'}, {'vector'});
                validateattributes(hash1, {'numeric'}, {'vector'});
                val = obj.compare(hash1, hash2);
                validateattributes(val, {'numeric'}, {'scalar'});
            end
        end

        function test_2
            img = imread(fullfile(mexopencv.root(), 'test', 'tsukuba.png'));

            hash = cv.ImgHash.averageHash(img);
            validateattributes(hash, {'uint8'}, {'vector'}); % 'numel',8

            hash = cv.ImgHash.blockMeanHash(img, 'Mode','Mode1');
            validateattributes(hash, {'uint8'}, {'vector'}); % 'numel',32|121 (Mode0|Mode1)

            hash = cv.ImgHash.colorMomentHash(img);
            validateattributes(hash, {'double'}, {'vector'}); % 'numel',42

            hash = cv.ImgHash.marrHildrethHash(img, 'Alpha',2);
            validateattributes(hash, {'uint8'}, {'vector'}); % 'numel',72

            hash = cv.ImgHash.pHash(img);
            validateattributes(hash, {'uint8'}, {'vector'}); % 'numel',8

            hash = cv.ImgHash.radialVarianceHash(img, 'Sigma',1);
            validateattributes(hash, {'uint8'}, {'vector'}); % 'numel',40
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function P = cubic_eval(C,t);
  % CUBIC_EVAL Evaluate a cubic Bezier curve.
  %
  % P = cubic_eval(C,t);
  %
  % Inputs:
  %   C  4 by dim list of control points
  %   t  #t list of evaluation parameters
  % Outputs:
  %   P  #t by dim list of evaluated points
  %
  % See also: readSVG_cubics, cubic_split
  %
  t = reshape(t,numel(t),1);
  %B = bernstein_eval([0 1 2 3],3,t);
  %P = B*C;
  P =  ...
    1*(1-t).^3.*t.^0.*C(1,:) +  ...
    3*(1-t).^2.*t.^1.*C(2,:) +  ...
    3*(1-t).^1.*t.^2.*C(3,:) +  ...
    1*(1-t).^0.*t.^3.*C(4,:);
end