text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   KINOVA JACO.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   09/01/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

robot.name= 'KINOVA_JACO';
robot.DH.theta= '[q(1) q(2)+pi/2 q(3)+pi/2 q(4) q(5) q(6)]';
robot.DH.d='[0.2755 0 0.0098 0.2492 0.0838 0.2106]';
robot.DH.a='[0 0.410 0 0 0 0]';
robot.DH.alpha= '[pi/2 pi pi/2 (11*pi)/36 (-11*pi)/36 0]';
robot.J=[];

robot.inversekinematic_fn = 'inverse_kinematics_jacobian(robot, T, q)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

robot.parameters.stop_iterations = 1500;
robot.parameters.step_time = 0.01;
robot.parameters.epsilonQ = 0.01;

%SPECIAL PARAMETERS TO SOLVE THE INVERSE KINEMATICS
robot.parameters.step_time=0.01;
%Error in XYZ to stop inverse kinematics
robot.parameters.epsilonXYZ=0.001;
%Error in Quaternion to stop inverse kinematics.
robot.parameters.epsilonQ=0.001;
robot.parameters.stop_iterations=1500;

% 1: maximize manipulability.
% 0: standard.
% -1: minimize manipulability.
robot.maximize_manipulability = 0;

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[-pi pi; %Axis 1, minimum, maximum
                deg2rad(-100) deg2rad(100); %Axis 2, minimum, maximum
                deg2rad(-220) deg2rad(60); %Axis 3
                deg2rad(-200) deg2rad(200); %Axis 4: Unlimited (400� default)
                deg2rad(-120) deg2rad(120); %Axis 5
                deg2rad(-400) deg2rad(400)]; %Axis 6: Really Unlimited to (800� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(200); %Axis 1, rad/s
                deg2rad(200); %Axis 2, rad/s
                deg2rad(260); %Axis 3, rad/s
                deg2rad(360); %Axis 4, rad/s
                deg2rad(360); %Axis 5, rad/s
                deg2rad(450)];%Axis 6, rad/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

% end effectors maximum velocity
robot.linear_velmax = 1; %m/s

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

% GRAPHICS
robot.graphical.has_graphics=0;
robot.graphical.color = [255 102 51]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-0.75 0.75 -0.75 0.75 0 1.5];
%read graphics files
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DYNAMIC PARAMETERS
%   WARNING! These parameters do not correspond to the actual IRB 140
%   robot. They have been introduced to demonstrate the necessity of
%   simulating the robot and should be used only for educational purposes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[25 27 15 10 2.5 1.5];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[0       0          0; %(rx, ry, rz) link 1
                     -0.05	 0.006	 0.1; %(rx, ry, rz) link 2
                    -0.0203	-0.0141	 0.070;  %(rx, ry, rz) link 3
                     0       0.019       0;%(rx, ry, rz) link 4
                     0       0           0;%(rx, ry, rz) link 5
                     0       0         0.032];%(rx, ry, rz) link 6

%Inertia matrices of each link with respect to its D-H reference system.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, for each row
robot.dynamics.Inertia=[0      0.35	0   	0	0	0;
    .13     .524	.539	0	0	0;
    .066	.086	.0125	0	0	0;
    1.8e-3	1.3e-3	1.8e-3	0	0	0;
    .3e-3	.4e-3	.3e-3	0	0	0;
    .15e-3	.15e-3	.04e-3	0	0	0];

robot.motors=load_motors([5 5 5 4 4 4]);
%Speed reductor at each joint
robot.motors.G=[300 300 300 300 300 300];

repo: ARTE


text:
classdef (Abstract) mm_shared_pfcpf_acp < mp.mm_shared_pfcpf_ac

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

%     properties
%     end

    methods
        function [vx_, z_, x_] = convert_x_m2n(obj, mmx, nm, only_v)
            %% x = obj.pf_convert(mmx, nm)
            %% [v, z] = obj.pf_convert(mmx, nm)
            %% [v, z, x] = obj.pf_convert(mmx, nm)
            %% ... = obj.pf_convert(mmx, nm, only_v)

            %% update v_, z_ from mmx
            nm_vars = obj.update_nm_vars(mmx, nm);
            vx_ = nm_vars.vm .* exp(1j * nm_vars.va);
            z_ = nm_vars.zr + 1j * nm_vars.zi;

            %% update z, if requested
            if nargin < 4 || ~only_v
                z_ = obj.update_z(nm, vx_, z_, obj.aux_data);
            end

            if nargout < 2
                vx_ = [vx_; z_];
            elseif nargout > 2
                x_ = [vx_; z_];
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [u,Du,eqn,info] = Poisson3T1(node,elem,pde,bdFlag, option)
%% POISSON3T1 Poisson equation: trilinear element.
%
%   u = Poisson3T1(node,elem,pde,bdFlag) produces the trilinear finite element
%   approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
% Author: Huayi Wei < huayiwei1984@gmail.com>.
%
% See also: PoissonQ1
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end
[N,Dim] = size(node);
[NT,NV] = size(elem);
Ndof = N;

tic;

%% Stiffness matrix and mass matrix in each direction
hx = node(elem(:,2),1) - node(elem(:,1),1);
hy = node(elem(:,4),2) - node(elem(:,1),2);
hz = node(elem(:,5),3) - node(elem(:,1),3);
volume = abs(hx.*hy.*hz);
Ax = 1./hx*[1 -1];
Mx = hx*[1/3 1/6];
Ay = 1./hy*[1 -1];
My = hy*[1/3 1/6];
Az = 1./hz*[1 -1];
Mz = hz*[1/3 1/6];

%% Assemble stiffness matrix
% generate sparse pattern
ii = zeros(36*NT,1); jj = zeros(36*NT,1);
index = 0;
for i = 1:NV
    for j = i:NV
        ii(index+1:index+NT) = double(elem(:,i));
        jj(index+1:index+NT) = double(elem(:,j));
        index = index + NT;
    end
end
% index map
itoixiyiz = [[1 2 2 1 1 2 2 1]' [1 1 2 2 1 1 2 2]' [1 1 1 1 2 2 2 2]'];

if ~isfield(pde,'d'), pde.d = []; end
sA = zeros(36*NT,1);
index = 0;
for i = 1:NV
    for j = i:NV
        ix = itoixiyiz(i,1); jx = itoixiyiz(j,1);
        iy = itoixiyiz(i,2); jy = itoixiyiz(j,2);
        iz = itoixiyiz(i,3); jz = itoixiyiz(j,3);
        Aij = Ax(:,2-(ix==jx)).*My(:,2-(iy==jy)).*Mz(:,2-(iz==jz)) ...
            + Mx(:,2-(ix==jx)).*Ay(:,2-(iy==jy)).*Mz(:,2-(iz==jz)) ...
            + Mx(:,2-(ix==jx)).*My(:,2-(iy==jy)).*Az(:,2-(iz==jz));
        if ~isempty(pde.d) && isnumeric(pde.d) % d is piecewise constant
            Aij = pde.d.*Aij;
        end
        sA(index+1:index+NT) = Aij;
        index = index + NT;
    end
end
% assemble the matrix
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
A = A + AU + AU';
clear Aij ii jj sA

%% Assemble the right hand side
b = zeros(N,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end

if ~isempty(pde.f)
    [GaussPt, weight] = quadptscube(option.fquadorder);
    lambdax = GaussPt(:,1);
    lambday = GaussPt(:,2);
    lambdaz = GaussPt(:,3);
    nQuad = size(GaussPt,1);
    bt = zeros(NT,NV);
    for p = 1:nQuad
		% quadrature points in the x-y-z coordinate
        phi2d = kron([lambday(p) 1-lambday(p)],[lambdax(p); 1-lambdax(p)]);
        phi = transpose(kron([lambdaz(p) 1-lambdaz(p)],phi2d([1 2 4 3])));
        J = volume;
%         [phi, tempvar, J] = hexbasis(node,elem, pts(p,:));
        pxyz = zeros(NT, Dim);
        for k = 1:Dim
            xk = node(:,k);
            pxyz(:,k) = xk(elem)*phi(:);
        end
		fp = pde.f(pxyz);
        bt = bt + (weight(p)*fp.*J)*phi';
    end
    b = accumarray(elem(:),bt(:),[Ndof 1]);
end
clear pxyz bt

%% Set up boundary conditions
if ~exist('bdFlag','var'), bdFlag = []; end
[AD,b,u,freeNode,isPureNeumann] = getbd3(b);

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

%% Solve the system of linear equations
if isempty(freeNode), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeNode) = AD(freeNode,freeNode)\b(freeNode);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
%         option.x0 = u;
        option.solver = 'CG';
%         [u,info] = mg(AD,b,elem,option);
        [u(freeNode),info] = amg(AD(freeNode,freeNode),b(freeNode),option);
    case 'amg'
        option.solver = 'CG';
        [u(freeNode),info] = amg(AD(freeNode,freeNode),b(freeNode),option);
end
% post-process for pure Neumann problem
if isPureNeumann
    patchArea = accumarray(elem(:),[J;J;J;J;J;J;J;J]/8, [N 1]); % TODO: Here is J?
    uc = sum(u.*patchArea)/sum(J);
    u = u - uc;   % int u = 0
end

%% Output information
eqn = struct('A',AD,'b',b,'freeNode',freeNode);
info.assembleTime = assembleTime;

%% Compute Du
Du = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbd3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,b,u,freeNode,isPureNeumann] = getbd3(b)
    %% Set up of boundary conditions.
    %
    % 1) Modify the matrix for Dirichlet boundary nodes, which are not degree
    % of freedom. Values at these nodes are evaluatation of pde.g_D. The
    % original stiffness matrix A is turn into the matrix AD by enforcing
    % AD(fixedNode,fixedNode)=I, AD(fixedNode,freeNode)=0, AD(freeNode,fixedNode)=0.
    %
    % 2) Modify the right hand side b. The Neumann boundary integral is added
    % to b. For Dirichlet boundary ndoes, b(fixedNode) is the evaluation of
    % pde.g_D.
    %
    % Special attentation should be given for the pure Neumann boundary
    % condition. To enforce the compatible condition, the vector b should have
    % mean value zero. To avoid a singular matrix, the 1st node is chosen as
    % fixedNode.
    %
    % The order of assigning Neumann and Dirichlet boundary condition is
    % important to get the right setting at the intersection nodes of Dirichlet
    % and Neumann boundary edges.
    %
    % Reference: Long Chen. Finite Element Methods and its Programming. Lecture
    % Notes.

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

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    isRobin = (bdFlag(:) == 3);
    if any(isRobin)
        allFace = [elem(:,[1 4 3 2]); elem(:,[1 2 6 5]); elem(:,[5 6 7 8]);...
                   elem(:,[8 7 3 4]); elem(:,[4 1 5 8]); elem(:,[2 3 7 6])];
        Robin = allFace(isRobin,:);
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        %% To do: this is for tet mesh. Modify for hex mesh later.
        v12 = node(Robin(:,2),:)-node(Robin(:,1),:);
        v13 = node(Robin(:,3),:)-node(Robin(:,1),:);
        area = 0.5*sqrt(abs(sum(mycross(v12,v13,2).^2,2)));
        if ~isfield(option,'gRquadorder')
            option.gRquadorder = 3;   % default order exact for linear gR
        end
        [lambdagR,weightgR] = quadpts(option.gRquadorder);
        nQuadgR = size(lambdagR,1);
        bdphi = lambdagR;
        NR = size(Robin,1);
        ss = zeros(NR,3,3);
        % int g_R phi_i phi_j
        for pp = 1:nQuadgR
            ppxyz = lambdagR(pp,1)*node(Robin(:,1),:) ...
                  + lambdagR(pp,2)*node(Robin(:,2),:) ...
                  + lambdagR(pp,3)*node(Robin(:,3),:);
            gRp = pde.g_R(ppxyz);
            for iR = 1:3
                for jR = iR:3
                    ss(:,iR,jR) = ss(:,iR,jR) + ...
                    weightgR(pp)*gRp*bdphi(pp,iR).*bdphi(pp,jR);
                end
            end
        end
        ss(:) = ss(:).*repmat(area,9,1);
        index = 0;
        for iR = 1:3
            for jR = 1:3
                iiR(index+1:index+NR) = double(Robin(:,iR));
                jjR(index+1:index+NR) = double(Robin(:,jR));
                if jR>=iR
                    ssR(index+1:index+NR) = ss(:,iR,jR);
                else
                    ssR(index+1:index+NR) = ss(:,jR,iR);
                end
                index = index + NR;
            end
        end
        A = A + sparse(iiR,jjR,ssR,Ndof,Ndof);
    end

    % Find Dirichlet boundary nodes: fixedNode
    fixedNode = []; freeNode = [];
    if ~isempty(bdFlag) % bdFlag specifies different bd conditions
        [fixedNode,bdFace,isBdNode] = findboundary3(elem,bdFlag);
        freeNode = find(~isBdNode);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given
        [fixedNode,bdFace,isBdNode] = findboundary3(elem);
        freeNode = find(~isBdNode);
    end
    isPureNeumann = false;
    if isempty(fixedNode) && isempty(Robin) % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedNode = 1;
        freeNode = 2:Ndof;    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedNode,fixedNode)=I, AD(fixedNode,freeNode)=0, AD(freeNode,fixedNode)=0.
    bdidx = zeros(Ndof,1);
    bdidx(fixedNode) = 1;
    Tbd = spdiags(bdidx,0,Ndof,Ndof);
    T = spdiags(1-bdidx,0,Ndof,Ndof);
    AD = T*A*T + Tbd;

    %% Part 2: Find boundary faces and modify the load b
    % Find boundary faces: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        Neumann = bdFace;
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        [tempvar,Neumann] = findboundary3(elem); %#ok<ASGLU>
    end

    % Neumann boundary condition
    if ~isempty(Neumann) && ~isempty(pde.g_N) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        v12 = node(Neumann(:,2),:)-node(Neumann(:,1),:);
        v14 = node(Neumann(:,4),:)-node(Neumann(:,1),:);
        area = sqrt(abs(sum(mycross(v12,v14,2).^2,2)));
        % three middle points rule
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 3;   % default order exact for linear gN
        end
        [GaussPt,weightgN] = quadptsquad(option.gNquadorder);
        xxi = GaussPt(:,1);
        eta = GaussPt(:,2);
        nQuadgN = size(GaussPt,1);
        % barycentrical coordinate for four basis
        lambdagN = zeros(nQuadgN,4);
        lambdagN(:,1) = (1 - xxi).*(1 - eta);
        lambdagN(:,2) = xxi.*(1 - eta);
        lambdagN(:,3) = xxi.*eta;
        lambdagN(:,4) = (1 - xxi).*eta;
        ge = zeros(size(Neumann,1),4);
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxyz = node(Neumann(:,1),:)*lambdagN(pp,1) + node(Neumann(:,2),:)*lambdagN(pp,2)...
                  + node(Neumann(:,3),:)*lambdagN(pp,3) + node(Neumann(:,4),:)*lambdagN(pp,4);
            gNp = pde.g_N(ppxyz);
            for iN = 1:4
                ge(:,iN) = ge(:,iN) + weightgN(pp)*lambdagN(pp,iN)*gNp;
            end
        end
        ge = ge.*repmat(area,1,4);
        b = b + accumarray(Neumann(:),ge(:),[Ndof,1]);
    end
    % The case with non-empty Neumann edges but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann edges and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if ~isPureNeumann && ~isempty(fixedNode) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && all(pde.g_D == 0))    % nonzero g_D
        if isnumeric(pde.g_D)
            u(fixedNode) = pde.g_D(fixedNode);
        else
            u(fixedNode) = pde.g_D(node(fixedNode,:));
        end
        b = b - A*u;
    end
    if ~isPureNeumann % non-empty Dirichlet boundary condition
        b(fixedNode) = u(fixedNode);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

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

repo: ifem


text:
classdef TestDualTVL1OpticalFlow
    %TestDualTVL1OpticalFlow

    methods (Static)
        function test_1
            im1 = 255*uint8([...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 1 0 1 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
            ]);
            im2 = circshift(im1, [0 1]);
            alg = cv.DualTVL1OpticalFlow();
            flow = alg.calc(im1, im2);
            validateattributes(flow, {'single'}, ...
                {'3d', 'size',[size(im1,1) size(im1,2) 2]});
        end

        function test_2
            prevImg = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale1.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            nextImg = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale2.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            alg = cv.DualTVL1OpticalFlow();
            alg.ScaleStep = 0.3;
            alg.MedianFiltering = 3;
            flow = alg.calc(prevImg, nextImg);
            validateattributes(flow, {'single'}, ...
                {'3d', 'size',[size(prevImg,1) size(prevImg,2) 2]});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Author : FUAT COGUN
% Date   : 16.07.2009
%
%
% Gets the video information of the specified video file.
%
%
% Inputs  : -Name of the video file                             (String)
% ======
%
% Outputs : -Number of frames in video file                     (Integer)
% =======   -Frames per Second                                  (Integer)
%           -Width of the video                                 (Integer)
%           -Height of the video                                (Integer)
%           -Width-to-Height ratio                              (Double)

function [numOfFrames fps videoWidth videoHeight ratio] = getVideoInfo(videoFile)

    info = aviinfo(videoFile);
    numOfFrames = info.NumFrames;
    fps = info.FramesPerSecond;
    videoWidth = info.Width;
    videoHeight = info.Height;
    ratio = videoWidth/videoHeight;

author: Sable

repo: mcbench-benchmarks


text:
function [edges,idx,edgemap]=uniqedges(elem)
%
% [edges,idx,edgemap]=uniqedges(elem)
%
% return the unique edge list from a surface or tetrahedral mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%     elem: a list of elements, each row is a list of nodes for an element.
%           elem can have 2, 3 or 4 columns
%
% output:
%     edge: unique edges in the mesh, denoted by a pair of node indices
%     idx:  index of the output in the raw edge list (returned by meshedge)
%     edgemap: index of the raw edges in the output list (for triangular mesh)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(elem)==2)
   edges=elem;
elseif(size(elem)>=3)
   edges=meshedge(elem);
else
   error('invalid input');
end

[uedges,idx,jdx]=unique(sort(edges,2),'rows');
edges=edges(idx,:);
if(nargout>2)
    edgemap=reshape(jdx,[size(elem,1) nchoosek(size(elem,2),2)]);
end

author: fangq


text:
function [xEst,PEst]=batchLSLinMeasLinDyn(z,H,F,R,kD,Q,numMeas)
%%BATCHLSLINMEASLINDYN Perform batch least squares state estimation under a
%                linear measurement model and a linear dynamic model with
%                optional process noise.
%
%INPUTS: z The zDim X N matrix of measurements for the whole batch. It is
%          assumed that the measurements have the same dimensionality over
%          the batch. If an empty matrix is passed, then it is assumed that
%          the user only wants the covariance matrix and not xEst.
%        H The zDim X xDim X N hypermatrix of measurement matrices such
%          that H(:,:,k)*x+w is the measurement at time k, where x is the
%          state and w is zero-mean Gaussian noise with covariance matrix
%          R(:,:,k). Alternatively, if all of the measurement matrices are
%          the same, one can just pass a single zDim X xDim matrix.
%        F An xDim X xDim X (N-1) hypermatrix of  matrices. The state at
%          discrete-time k+1 is modeled as F(:,:,k) times the state at time
%          k plus zero-mean Gaussian process noise with covariance matrix
%          Q(:,:,k). Alternatively, if all of the state transition matrices
%          are the same, one can just pass a single xDim X xDim matrix.
%          Note that all of the F matrices must be invertible if kD is not
%          equal to one.
%        R The zDim X zDim X N hypermatrix of measurement covariance
%          matrices. Alternatively, if all of the measurement covariance
%          matrices are the same, one can just pass a single zDimXzDim
%          matrix.
%       kD The discrete time-step for which the covariance matrix of an
%          estimate is desired, where z(:,1) is at discrete time-step 1
%          (not 0).
%        Q The xDim X xDim X (N-1) hypermatrix of process noise covariance
%          matrices. Alternatively, if all of the process noise covariance
%          matrices are the same, one can just pass a single xDimXxDim
%          matrix. If an empty matrix is passed for Q, then the estimation
%          is performed assuming there is no process noise.
%  numMeas If an empty matrix is passed for z, then the numMeas parameter
%          must be passed indicating the number of measurements in the
%          batch. Otherwise, this parameter is
%                ignored.
%
%OUTPUTS: xEst The batch state estimate at step kD, unless an empty matrix
%              was passed for z, in which case xEst is empty.
%         PEst A covariance matrix estimate that goes with the state
%              estimate.
%
%The algorithm is an implementation of the method of Section 3.3.2 of [1]
%for a linear dynamic model. Note that the cases in equation 28 to 32 of
%the paper did not cover all of the possibilities, so the correct
%generalization had to be derived from E[epsilon_{p,k}*epsilon_{q,k}'].
%
%REFERENCES:
%[1] A. B. Poore, B. J. Slocumb, B. J. Suchomel, F. H. Obermeyer, S. M.
%    Herman, and S. M. Gadaleta, "Batch maximum likelihood (ML) and maximum
%    a posteriori (MAP) estimation with process noise for tracking
%    applications," in Proceedings of SPIE: Signal and Data Processing of
%    Small Targets, vol. 5204, San Diego, CA, 3 Aug. 2003, pp. 188-199.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(~isempty(z))
    numMeas=size(z,2);
end

xDim=size(F,1);
zDim=size(R,1);

if(size(H,3)==1)
    HMats=repmat(H,[1,1,numMeas]);
else
    HMats=H;
end

if(size(F,3)==1)
    F=repmat(F,[1,1,numMeas-1]);
end

if(size(R,3)==1)
    R=repmat(R,[1,1,numMeas]);
end

if(isempty(Q))
    Q=zeros(xDim,xDim,numMeas);
elseif(size(Q,3)==1)
    Q=repmat(Q,[1,1,numMeas-1]);
end

transMats=getTransMats(F);

WMat=zeros(zDim*numMeas,zDim*numMeas);
for p=1:numMeas
    pIdxMin=(p-1)*zDim+1;
    pIdxMax=p*zDim;
    pSpan=pIdxMin:pIdxMax;
    for q=1:numMeas
        qIdxMin=(q-1)*zDim+1;
        qIdxMax=q*zDim;
        qSpan=qIdxMin:qIdxMax;
        if(p<kD&&q<kD)
            CumMat=zeros(xDim,xDim);

            for i=(max(p,q)+1):kD
                CumMat=CumMat+transMats(:,:,p,i)*Q(:,:,i-1)*transMats(:,:,q,i)';
            end

            WMat(pSpan,qSpan)=HMats(:,:,p)*CumMat*HMats(:,:,q)'+KDelta(p-q)*R(:,:,p);
        elseif(p>kD&&q>kD)
            CumMat=zeros(xDim,xDim);

            for i=(kD+1):min(p,q)
                CumMat=CumMat+transMats(:,:,p,i)*Q(:,:,i-1)*transMats(:,:,q,i)';
            end

            WMat(pSpan,qSpan)=KDelta(p-q)*R(:,:,p)+HMats(:,:,p)*CumMat*HMats(:,:,q)';
        elseif(p==kD&&q==kD)
            WMat(pSpan,qSpan)=R(:,:,p);
        end
    end
end

%The propagated measurement matrices.
bigH=zeros(zDim*numMeas,xDim);
for curMeas=1:numMeas
    idxMin=(curMeas-1)*zDim+1;
    idxMax=curMeas*zDim;

    bigH(idxMin:idxMax,:)=HMats(:,:,curMeas)*transMats(:,:,curMeas,kD);
end

WInv=inv(WMat);
PEstInv=bigH'*WInv*bigH;
PEst=inv(PEstInv);
if(isempty(z))
    xEst=[];
else
    xEst=PEstInv\bigH'*(WMat\z(:));
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function title = p12_title ( title )

%*****************************************************************************80
%
%% P12_TITLE returns the title of problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Flat Stanley (ALL derivatives are zero at the root.)';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function handles = plotDataStarCoordinates(hObject, handles)

plotDimensions = get(handles.lstPlotDimensions, 'Value');

sFigureTitle = 'Data';
sFigureTitle = [sFigureTitle ' - Star Coordinates'];
[handles, handles.figData] = openPlotFigure(hObject, handles, ...
    'Data Plot (Star Coordinates)', sFigureTitle);

workData = normalizeData(handles.Data);
workData = workData(plotDimensions, :);

d = size(workData, 1);
unitRoot = exp(2 * pi * 1i / d);

starAxes = zeros(1, d);
for ii = 1:d
    starAxes(ii) = unitRoot^ii;
end

starAxes = repmat(starAxes', 1, size(workData, 2));
plotPoints = sum(workData .* starAxes, 1);
scatter(real(plotPoints), imag(plotPoints), ...
    getPlotMarkerSize(), ['k' getPlotMarkerStyle()]);

hold off;

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
classdef CarlaEnvironment < matlab.System & matlab.system.mixin.Propagates
    % Untitled Add summary here
    %
    % This template includes the minimum set of functions required
    % to define a System object with discrete state.

    % Public, tunable properties
    properties
        % Other actors
        npc_to_spawn = 20;
    end

    properties(DiscreteState)

    end

    % Pre-computed constants
    properties(Access = private)
        % Sensors resolution
        width = 1920;
        height = 1080;

        tesla;
        npc_list;

        keyboard_handle;

        rgb;
        module_rgb;
        time_last_update = 0;
    end

    methods(Access = protected)
        function setupImpl(obj)
            % Perform one-time calculations, such as computing constants
            port = int16(2000);
            client = py.carla.Client('localhost', port);
            client.set_timeout(2.0);
            world = client.get_world();

            % Spawn Vehicle
            blueprint_library = world.get_blueprint_library();
            car_list = py.list(blueprint_library.filter("model3"));
            car_bp = car_list{1};
            spawn_point = py.random.choice(world.get_map().get_spawn_points());
            obj.tesla = world.spawn_actor(car_bp, spawn_point);
            obj.tesla.set_autopilot(true);

            % Spawn NPC's
            npc_bps = blueprint_library.filter("vehicle");

            i = 1;
            while i <= obj.npc_to_spawn
                try
                    npc_bp = py.random.choice(npc_bps);
                    spawn_point = py.random.choice(world.get_map().get_spawn_points());
                    obj.npc_list{i} = world.spawn_actor(npc_bp, spawn_point);
                    obj.npc_list{i}.set_autopilot(true);
                catch
                    % In case spawing fails due to collision, try again
                    i = i - 1;
                end
                i = i + 1;
            end

            % RGB
            blueprint = world.get_blueprint_library().find('sensor.camera.rgb');
            blueprint.set_attribute('image_size_x', num2str(obj.width))
            blueprint.set_attribute('image_size_y', num2str(obj.height))

            transform = py.carla.Transform(py.carla.Location(pyargs('x',0.8, 'z',1.7)));
            obj.rgb = world.spawn_actor(blueprint, transform, pyargs('attach_to',obj.tesla));

            obj.module_rgb = sensorBind(obj.rgb, "rgb", "rgb", "array");

            % Handle keyboard
            callstr = 'set(gcbf, ''Userdata'', get(gcbf, ''Currentkey'')) ; uiresume ' ;
            obj.keyboard_handle = figure( 'name', 'Press a key', ...
                                           'keypressfcn', callstr, ...
                                           'windowstyle', 'modal', ...
                                           'numbertitle', 'off', ...
                                           'position', [0 0 1 1], ...
                                           'userdata', 'timeout');
            set(obj.keyboard_handle,'Visible','on');
        end

        function [RGB] = stepImpl(obj)
            % RGB
            RGB = uint8(py.getattr(obj.module_rgb, 'array'));

            % Handle keyboard
            dt = cputime - obj.time_last_update;

            try
                keyPressed = get(obj.keyboard_handle,'Userdata');
                processKeyboard(obj.tesla, keyPressed, dt);

                % Reset the pressed key
                set(obj.keyboard_handle,'Userdata', char(0));
            catch
                obj.tesla.set_autopilot(true);
            end

            % Update time
            obj.time_last_update = obj.time_last_update + dt;
            pause(0.025);
        end

        function [RGB] = isOutputComplexImpl(~)
            RGB = false;
        end

        function [RGB] = getOutputSizeImpl(obj)
            RGB = [obj.height obj.width 3];
        end

        function [RGB] = getOutputDataTypeImpl(~)
            RGB = 'uint8';
        end

        function [RGB] = isOutputFixedSizeImpl(~)
            RGB = true;
        end

        function resetImpl(~)
            % Initialize / reset discrete-state properties
        end
    end

    methods(Access= public)
        function delete(obj)
            % Delete the car from the Carla world
            if ~isempty(obj.tesla)
                obj.tesla.destroy();
            end

            % Delete NPCs
            if ~isempty(obj.npc_list)
                for i=1:obj.npc_to_spawn
                    obj.npc_list{i}.destroy();
                end
            end

            % RGB camera
            if ~isempty(obj.rgb)
                obj.rgb.destroy();
            end

            % Close the figure if it still exists
            if ~isempty(obj.keyboard_handle)
                close(obj.keyboard_handle);
            end

        end
    end
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function [obj] = modelReadTable(obj,varargin)
	% Read in table, decides whether to do a single or multiple tables. If multiple tables, should have the same column names.
	% Biafra Ahanonu
	% branched from controllerAnalysis: 2014.08.01 [16:09:16]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% which table to read in
	options.table = 'discreteStimulusTable';
	% whether to add file information to the table
	options.addFileInfoToTable = 0;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	if isempty(obj.delimiter)
		delimiter = 'tab';
	end
	pathClass = class(obj.(options.table));
	% pathClass
	switch pathClass
		case 'char'
			display(['loading: ' obj.(options.table)])
			outputTable = readtable(obj.(options.table),'Delimiter',obj.delimiter,'FileType','text');
		case 'cell'
			nPaths = length(obj.(options.table));
			for i=1:nPaths
				thisTablePath = obj.(options.table){i};
				if(exist(thisTablePath, 'file'))
					display(['loading: ' thisTablePath])
					if ~exist('outputTable')
						outputTable = readtable(char(obj.(options.table)(i)),'Delimiter',obj.delimiter,'FileType','text','TreatAsEmpty',{'NA','N/A'});
						if options.addFileInfoToTable==1
							[outputTable] = addFileInfoToTable(outputTable,thisTablePath);
						end
					else
						tmpTable = readtable(char(obj.(options.table)(i)),'Delimiter',obj.delimiter,'FileType','text','TreatAsEmpty',{'NA','N/A'});
						if options.addFileInfoToTable==1
							[tmpTable] = addFileInfoToTable(tmpTable,thisTablePath);
						end
						outputTable = [outputTable;tmpTable];
					end
				end
			end
		otherwise
	end

	obj.(options.table) = outputTable;

	function [inputTable] = addFileInfoToTable(inputTable,inputFileName)
		fileInfo = getFileInfo(inputFileName);
		repeatSize = [size(inputTable,1) 1];
		inputTable.subject = repmat(fileInfo.subjectNum,repeatSize);
		% inputTable.assay = repmat({fileInfo.assay},repeatSize);
		inputTable.trial = repmat({fileInfo.assay},repeatSize);
		inputTable.date = repmat({fileInfo.date},repeatSize);
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function C = spm_atranspa(A)
% Multiplies the transpose of a matrix by itself
% FORMAT C = spm_atranspa(A)
% A - real matrix
% C - real symmetric matrix resulting from A'A
%_______________________________________________________________________
%
% This compiled routine was written to save both memory and CPU time but
% is now deprecated. Use A'*A directly instead.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_atranspa.m 4418 2011-08-03 12:00:13Z guillaume $

persistent runonce
if isempty(runonce)
    warning('spm_atranspa is deprecated. Use A''*A instead.');
    runonce = 1;
end

C = A'*A;

author: spm


text:
function datapt = lvd_EngineTasks(stateLogEntry, subTask, engine)
%lvd_EngineTasks Summary of this function goes here
%   Detailed explanation goes here

    switch subTask
        case 'active'
            engineStates = stateLogEntry.getAllEngineStates();
            engineState = engineStates([engineStates.engine] == engine);

            if(not(isempty(engineState)))
                engineState = engineState(1);
                datapt = double(engineState.active);
            else
                datapt = -1;
            end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Reconstruction of 2D Cartesian Pulseq data
% provides an example on how data reordering can be detected from the MR
% sequence with almost no additional prior knowledge
%
% it loads Matlab .mat files with the rawdata in the format
%     adclen x channels x readouts
% if Matlab .mat file not available it attempt to load Siemens .dat (which needs mapVBVD in the path)
% but first it seeks the accompanying .seq file with the same name to interpret
%     the data

%% Load the latest file from the specified directory
path='../IceNIH_RawSend/'; % directory to be scanned for data files
%path='/data/Dropbox/ismrm2021pulseq_liveDemo/dataLive/Vienna_7T_Siemens'; % directory to be scanned for data files

if path(end)~=filesep, path=[path filesep]; end

pattern='*.seq';
D=dir([path pattern]);
[~,I]=sort([D(:).datenum]);
seq_file_path=[path D(I(end-0)).name]; % use end-1 to reconstruct the second-last data set, etc...
                                                % or replace I(end-0) with I(1) to process the first dataset, I(2) for the second, etc...
%seq_file_path='../interpreters/siemens/data_example/gre_example.seq'

% keep basic filename without the extension
[p,n,e] = fileparts(seq_file_path);
basic_file_path=fullfile(p,n);

% try loading Matlab data
data_file_path=[basic_file_path '.mat'];
if isfile(data_file_path)
    fprintf(['loading `' data_file_path '´ ...\n']);
    data_unsorted = load(data_file_path);
    if isstruct(data_unsorted)
        fn=fieldnames(data_unsorted);
        assert(length(fn)==1); % we only expect a single variable
        data_unsorted=data_unsorted.(fn{1});
    end
else
    % revert to Siemens .dat file
    data_file_path=[basic_file_path '.dat'];
    fprintf(['loading `' data_file_path '´ ...\n']);
    twix_obj = mapVBVD(data_file_path);
    if iscell(twix_obj)
        data_unsorted = twix_obj{end}.image.unsorted();
        seqHash_twix=twix_obj{end}.hdr.Dicom.tSequenceVariant;
    else
        data_unsorted = twix_obj.image.unsorted();
        seqHash_twix=twix_obj.hdr.Dicom.tSequenceVariant;
    end

    if length(seqHash_twix)==32
        fprintf(['raw data contain pulseq-file signature ' seqHash_twix '\n']);
    end

end
[adc_len,channels,readouts]=size(data_unsorted);

%% Load sequence from file
fprintf(['loading `' seq_file_path '´ ...\n']);
seq = mr.Sequence();              % Create a new sequence object
seq.read(seq_file_path,'detectRFuse');
seqName=seq.getDefinition('Name');
if ~isempty(seqName), fprintf('sequence name: %s\n',seqName); end
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();
figure; plot(ktraj(1,:),ktraj(2,:),'b',...
             ktraj_adc(1,:),ktraj_adc(2,:),'r.'); % a 2D plot
axis('equal'); title('2D k-space trajectory');

%% Analyze the trajectory data (ktraj_adc)
fprintf('analyzing the k-space trajectory ...\n');
k_extent=max(abs(ktraj_adc),[],2);
k_scale=max(k_extent);
k_threshold=k_scale/5000;

% detect unused dimensions and delete them
if any(k_extent<k_threshold)
    ktraj_adc(k_extent<k_threshold,:)=[]; % delete rows
    k_extent(k_extent<k_threshold)=[];
end

% detect dK, k-space reordering and repetitions (or slices, etc)
kt_sorted=sort(ktraj_adc,2);
dk_all=kt_sorted(:,2:end)-kt_sorted(:,1:(end-1));
dk_all(dk_all<k_threshold)=NaN;
dk_min=min(dk_all,[],2);
dk_max=max(dk_all,[],2);
dk_all(dk_all-dk_min(:,ones(1,size(dk_all,2)))>k_threshold)=NaN;
dk_all_cnt=sum(isfinite(dk_all),2);
dk_all(~isfinite(dk_all))=0;
dk=sum(dk_all,2)./dk_all_cnt;
dk(~isfinite(dk))=0;
[~,k0_ind]=min(sum(ktraj_adc.^2,1));
kindex=round((ktraj_adc-ktraj_adc(:,k0_ind*ones(1,size(ktraj_adc,2))))./dk(:,ones(1,size(ktraj_adc,2))));
kindex(~isfinite(kindex))=0;
kindex_min=min(kindex,[],2);
kindex_mat=kindex-kindex_min(:,ones(1,size(ktraj_adc,2)))+1;
kindex_end=max(kindex_mat,[],2);
sampler=zeros(kindex_end');
repeat=zeros(1,size(ktraj_adc,2));
for i=1:size(kindex_mat,2)
    if (size(kindex_mat,1)==3)
        ind=sub2ind(kindex_end,kindex_mat(1,i),kindex_mat(2,i),kindex_mat(3,i));
    else
        ind=sub2ind(kindex_end,kindex_mat(1,i),kindex_mat(2,i));
    end
    repeat(i)=sampler(ind);
    sampler(ind)=repeat(i)+1;
end
if (max(repeat(:))>0)
    kindex=[kindex;(repeat+1)];
    kindex_mat=[kindex_mat;(repeat+1)];
    kindex_end=max(kindex_mat,[],2);
end
%figure; plot(kindex(1,:),kindex(2,:),'.-');

%% sort the k-space data into the data matrix
% the incoming data order is [kx coils acquisitions]
data_coils_last = permute(data_unsorted, [1, 3, 2]);
data_coils_last = reshape(data_coils_last, [adc_len*readouts, channels]);

data=zeros([kindex_end' channels]);
if (size(kindex,1)==3)
    for i=1:size(kindex,2)
        data(kindex_mat(1,i),kindex_mat(2,i),kindex_mat(3,i),:)=data_coils_last(i,:);
    end
else
    for i=1:size(kindex,2)
        data(kindex_mat(1,i),kindex_mat(2,i),:)=data_coils_last(i,:);
    end
end

if size(kindex,1)==3
    nImages=size(data,3);
else
    nImages=1;
    data=reshape(data, [size(data,1) size(data,2) 1 size(data,3)]); % we need a dummy images/slices dimension
end

% account for the matlab's strange convention with data dimensions
data=flip(data,1); % left-right orientation, works on TRIO
data=flip(data,2);

%figure; imab(log(abs(data))); title('k-space data');

%% Reconstruct coil images

images = zeros(size(data));
%figure;

for ii = 1:channels
    images(:,:,:,ii) = ifftshift(ifft2(ifftshift(data(:,:,:,ii)))); % 1.4.0 does not need read inversion
end

% Phase images (possibly channel-by-channel and echo-by-echo)
%figure;imab(angle(images));colormap('jet');
%figure;imab(abs(images));colormap('gray');

%% Image display with optional sum of squares combination
figure;
if channels>1
    sos=abs(sum(images.*conj(images),ndims(images))).^(1/2);
    sos=sos./max(sos(:));
    imab(sos); title('reconstructed image(s), sum-of-squares');
    %imwrite(sos, ['img_combined.png']
else
    imab(abs(images)); title('reconstructed image(s)');
end
colormap('gray');
saveas(gcf,[basic_file_path '_image_2dfft'],'png');

%% reconstruct field map (optional)

if size(images,3)>=2
    cmplx_diff=images(:,:,2,:).*conj(images(:,:,1,:));
    phase_diff_image=angle(sum(cmplx_diff,4));
    figure;
    imab(phase_diff_image);colormap('jet');
    title('phase difference(s) echo2-echo1');
end

%% gif movie export

scale=1.2;
filename='fftReconMovie';
fps=5;

if size(images,3)>4

    clm=gray(256);

    if channels>1
        imex=sos;
    else
        imex=abs(images);
        imex=imex/max(imex(:));
    end

    imex=permute(imex(:,end:-1:1,:),[2,1,3]);

    imind=uint8(scale*imex*(size(clm,1)-1)+1);
    imwrite(imind(:,:,1),clm, [filename, '.gif'],'DelayTime',1/fps,'Loopcount',inf);
    for i=2:size(imind,3),
        imwrite(imind(:,:,i),clm, [filename, '.gif'],'DelayTime',1/fps, 'WriteMode','append');
    end

end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
%SETLABELS Reset labels of dataset or mapping
%
%   A = SETLABELS(A,LABELS,J)
%   W = SETLABELS(W,LABELS,J)
%
% INPUT
%   A          Input dataset, size [M,K]
%   W          Input mapping (classifier), size [K,C]
%   LABELS     Desired labels. M or length(J) labels for a dataset.
%              C or length(J) labels for a mapping.
% OUTPUT
%   A          Dataset with reset labels
%   W          Mapping with reset labels
%
% DESCRIPTION
% The labels of the dataset A are reset by LABELS. If supplied, the
% index vector J defines the objects for wich LABELS applies. If in
% LABELS just a single label is given all the objects defined by J
% are given that label. If LABELS is empty ([]) or NaN all the objects
% defined by J are marked as unlabeled.
%
% If A has soft labels (label type is 'soft') or has no labels but
% targets (label type is 'targets'), these soft labels or targets are
% replaced by LABELS, provided it has the right size.
%
% For soft labels and targets supplied to relabel a dataset, LABELS may be
% supplied as a dataset of which the data are used for the soft labels or
% targets and the feature labels are used to set LABLIST of A.
%
% The labels stored in a classifier mapping W are assigned as feature
% labels of the resulting dataset D in case a dataset B is applied to W:
% D = A*W.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, MULTI_LABELING

author: marianux

repo: ecg-kit

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


text:
% QUATERNION DIVISION
function [qDiv] = qDivide(q,r)
% The function compute the division of quaternion q, by
% quaterion v.
% Associated block:
% "Quaternion Division"

qDiv = zeros(4,1);
% Normalise the second vector
qDiv(1) =  q(1)*r(1) + q(2)*r(2) + q(3)*r(3) + q(4)*r(4);
qDiv(2) = -q(1)*r(2) + q(2)*r(1) + q(3)*r(4) - q(4)*r(3);
qDiv(3) =  q(1)*r(3) + q(2)*r(4) + q(3)*r(1) + q(4)*r(2);
qDiv(4) = -q(1)*r(4) + q(2)*r(3) - q(3)*r(2) + q(4)*r(1);
qDiv = qDiv./(norm(r));
end

repo: OpenMAS


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rfAmp,rfPhase,rfFreq,rfCoil,GzAmp,GyAmp,GxAmp,ADC,Ext,uts,ts,flags]=PSD_SPGR3DGM
global VCtl
global VVar
CV1=2e-3;
CV10=0;
CV11=0;
CV12=0;
CV13=0;
CV14=0;
CV2=3e-3;
CV3=abs(rem(VVar.TRCount,2)*2-1);
CV4=100e-3;
CV5=4000;
CV6=-4.0*42.576*7;
CV7=0;
CV8=0;
CV9=0;
rfAmpAll=[];
rfPhaseAll=[];
rfFreqAll=[];
rfCoilAll=[];
GzAmpAll=[];
GyAmpAll=[];
GxAmpAll=[];
ADCAll=[];
ExtAll=[];
rfTimeAll=[];
GzTimeAll=[];
GyTimeAll=[];
GxTimeAll=[];
ADCTimeAll=[];
ExtTimeAll=[];
SEtAll=[];
uts=[];
ts=[];
flags=[];
if VCtl.PlotSeq == 1
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
ts = [ts tS tE];
end
ts = [0 max(ts)-min(ts)];
return;
end
%==============Pulses 1==============
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
rfTime=[];
GzTime=[];
GyTime=[];
GxTime=[];
ADCTime=[];
ExtTime=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if isempty(tS) | isempty(tE) | (tS>=tE)
error('SE setting is incorrect for Pulses 1!');
end
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{1};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=CV4+0.6e-3;
p.tStart=CV4+0.1e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp1,rfPhase1,rfFreq1,rfCoil1,rfTime1]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil1(1)
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
else
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=0;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=0.2e-3;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=VCtl.TR;
p.tStart=CV4+0.7e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp2,rfPhase2,rfFreq2,rfCoil2,rfTime2]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil2(1)
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
else
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=CV5;
p.Notes='Gaussian rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.dt=100e-6;
p.rfFreq=CV6;
p.rfPhase=0;
p.tEnd=CV4;
p.tStart=0;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp3,rfPhase3,rfFreq3,rfCoil3,rfTime3]=rfGaussian(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil3(1)
rfAmp=[rfAmp rfAmp3];
rfPhase=[rfPhase rfPhase3];
rfFreq=[rfFreq rfFreq3];
rfCoil=[rfCoil rfCoil3];
rfTime=[rfTime rfTime3];
end
else
rfAmp=[rfAmp rfAmp3];
rfPhase=[rfPhase rfPhase3];
rfFreq=[rfFreq rfFreq3];
rfCoil=[rfCoil rfCoil3];
rfTime=[rfTime rfTime3];
end
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gz1Sign=1;
p.Gz2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV4+CV2;
p.t1Start=CV4+CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GzAmp1,GzTime1]=GzCartesian(p);
GzAmp=[GzAmp GzAmp1];
GzTime=[GzTime GzTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gy1Sign=1;
p.Gy2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV4+CV2;
p.t1Start=CV4+CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GyAmp1,GyTime1]=GyCartesian(p);
GyAmp=[GyAmp GyAmp1];
GyTime=[GyTime GyTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gx1Sign=-1;
p.Gx2Sign=1;
p.Gx3Sign=0;
p.Notes='cartesian frequency';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1Start=CV4+CV1;
p.t2Middle=VCtl.TE;
p.t3Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GxAmp1,GxTime1]=GxCartesian(p);
GxAmp=[GxAmp GxAmp1];
GxTime=[GxTime GxTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Notes='cartesian readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tMiddle=VCtl.TE;
if strcmp(p.Switch,'on')
[ADC1,ADCTime1]=ADCCartesian(p);
ADC=[ADC ADC1];
ADCTime=[ADCTime ADCTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=5;
p.Notes='calculate remaining scan time';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=0;
if strcmp(p.Switch,'on')
[Ext1,ExtTime1]=ExtBit(p);
Ext=[Ext Ext1];
ExtTime=[ExtTime ExtTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=1;
p.Notes='reset K space location';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV4+0.6e-3;
if strcmp(p.Switch,'on')
[Ext2,ExtTime2]=ExtBit(p);
Ext=[Ext Ext2];
ExtTime=[ExtTime ExtTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=VCtl.TR*(99/100);
if strcmp(p.Switch,'on')
[Ext3,ExtTime3]=ExtBit(p);
Ext=[Ext Ext3];
ExtTime=[ExtTime ExtTime3];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV4+0.05e-3;
if strcmp(p.Switch,'on')
[Ext4,ExtTime4]=ExtBit(p);
Ext=[Ext Ext4];
ExtTime=[ExtTime ExtTime4];
end
p=[];
%--------------------
SEt=[tS tE];
rfAmp(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfPhase(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfFreq(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfCoil(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzAmp(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyAmp(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxAmp(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADC(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
Ext(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfTime(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzTime(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyTime(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxTime(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADCTime(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
ExtTime(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfAmp(abs(rfAmp)<eps) = 0;
rfTime = rfTime + SEt(1);
GzTime = GzTime + SEt(1);
GyTime = GyTime + SEt(1);
GxTime = GxTime + SEt(1);
ADCTime = ADCTime + SEt(1);
ExtTime = ExtTime + SEt(1);
rfAmpAll=[rfAmpAll rfAmp];
rfPhaseAll=[rfPhaseAll rfPhase];
rfFreqAll=[rfFreqAll rfFreq];
rfCoilAll=[rfCoilAll rfCoil];
GzAmpAll=[GzAmpAll GzAmp];
GyAmpAll=[GyAmpAll GyAmp];
GxAmpAll=[GxAmpAll GxAmp];
ADCAll=[ADCAll ADC];
ExtAll=[ExtAll Ext];
rfTimeAll=[rfTimeAll rfTime];
GzTimeAll=[GzTimeAll GzTime];
GyTimeAll=[GyTimeAll GyTime];
GxTimeAll=[GxTimeAll GxTime];
ADCTimeAll=[ADCTimeAll ADCTime];
ExtTimeAll=[ExtTimeAll ExtTime];
SEtAll=[SEtAll SEt];
end
%====================================
if isempty(rfTimeAll)
error('rf sequence line can not be empty! Master Tx coil element must be used.');
end
if isempty(GzTimeAll)
error('GzSS sequence line can not be empty!');
end
if isempty(GyTimeAll)
error('GyPE sequence line can not be empty!');
end
if isempty(GxTimeAll)
error('GxR sequence line can not be empty!');
end
if isempty(ADCTimeAll)
error('ADC sequence line can not be empty!');
end
if isempty(ExtTimeAll)
error('Ext sequence line can not be empty!');
end
SEflag=repmat([0 0 0 0 0 0]',[1 2]);
rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTimeAll))]);
Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTimeAll))]);
Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTimeAll))]);
Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTimeAll))]);
ADCflag=repmat([0 0 0 0 1 0]',[1 max(size(ADCTimeAll))]);
Extflag=repmat([0 0 0 0 0 1]',[1 max(size(ExtTimeAll))]);
ts=[[min(SEtAll) max(SEtAll)] rfTimeAll GzTimeAll GyTimeAll GxTimeAll ADCTimeAll ExtTimeAll]-min(SEtAll);
flags=[SEflag rfflag Gzflag Gyflag Gxflag ADCflag Extflag];
[ts,ind]=sort(ts);
uts=unique(ts);
flags=flags(:,ind);
[rfTime,ind]=sort(rfTimeAll-min(SEtAll));
rfAmp=rfAmpAll(:,ind);
rfPhase=rfPhaseAll(:,ind);
rfFreq=rfFreqAll(:,ind);
rfCoil=rfCoilAll(:,ind);
[GzTime,ind]=sort(GzTimeAll-min(SEtAll));
GzAmp=GzAmpAll(:,ind);
[GyTime,ind]=sort(GyTimeAll-min(SEtAll));
GyAmp=GyAmpAll(:,ind);
[GxTime,ind]=sort(GxTimeAll-min(SEtAll));
GxAmp=GxAmpAll(:,ind);
[ADCTime,ind]=sort(ADCTimeAll-min(SEtAll));
ADC=ADCAll(:,ind);
[ExtTime,ind]=sort(ExtTimeAll-min(SEtAll));
Ext=ExtAll(:,ind);
rfAmp(1) = 0;
rfPhase(1) = 0;
rfFreq(1) = 0;
GzAmp(1) = 0;
GyAmp(1) = 0;
GxAmp(1) = 0;
ADC(1) = 0;
Ext(1) = 0;
rfAmp(end) = 0;
rfPhase(end) = 0;
rfFreq(end) = 0;
GzAmp(end) = 0;
GyAmp(end) = 0;
GxAmp(end) = 0;
ADC(end) = 0;
Ext(end) = 0;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [m,other] = esvm_update_dfun(m)
%% Perform Distance Function learning for a single exemplar model.  We
% assume that the exemplar has a set of detections loaded in
% m.model.svxs and m.model.svbbs.
% Returns: model [m] with updated classifier
% If no arguments are given returns the suffix '-dfun' and
% classifier type 'dfun'
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

other = 'dfun';
%if no inputs are specified, just return the suffix of current method
if nargin==0
  m = '-dfun';
  return;
end

%If no mask is specified, then we assume that every single bin
%should be used.
if ~isfield(m.model,'mask') | length(m.model.mask)==0
  m.model.mask = logical(ones(numel(m.model.w),1));
end

%If the mask is only MxNx1, where the feature are MxNxF, then we
%repmat the mask to match the features.  This is important because
%it is sometimes smarter to save a 2D mask is 1/F the size of the
%full mask.
if length(m.model.mask(:)) ~= numel(m.model.w)
  m.model.mask = repmat(m.model.mask,[1 1 features]);
  m.model.mask = logical(m.model.mask(:));
end

mining_params = m.mining_params;
xs = m.model.svxs;
bbs = m.model.svbbs;

% A trick in case we have too many detections
MAXSIZE = 2000;
if size(xs,2) >= MAXSIZE
  HALFSIZE = MAXSIZE/2;
  %NOTE: random is better than top 5000
  r = m.model.w(:)'*xs;
  [tmp,r] = sort(r,'descend');
  r1 = r(1:HALFSIZE);

  r = HALFSIZE+randperm(length(r((HALFSIZE+1):end)));
  r = r(1:HALFSIZE);
  r = [r1 r];
  xs = xs(:,r);
  bbs = bbs(r,:);
end

newx = cat(2,m.model.x,xs);
newx = bsxfun(@minus,newx,m.model.x(:,1)).^2;

newy = cat(1,ones(size(m.model.x,2),1),-1*ones(size(xs,2),1));

number_positives = sum(newy==1);
number_negatives = sum(newy==-1);

wpos = mining_params.train_positives_constant;
wneg = 1;

% if mining_params.BALANCE_POSITIVES == 1
%   fprintf(1,'balancing positives\n');
%   wpos = 1/number_positives;
%   wneg = 1/number_negatives;
%   wpos = wpos / wneg;
%   wneg = wneg / wneg;
% end

newx = newx(logical(m.model.mask),:);

fprintf(1,' -----\nStarting SVM dim=%d... s+=%d, s-=%d ',...
        size(newx,1), number_positives, number_negatives);
starttime = tic;

svm_model = libsvmtrain(newy, newx',sprintf(['-s 0 -t 0 -c' ...
                    ' %f -w1 %.9f -q'], mining_params.train_svm_c, wpos));

if length(svm_model.sv_coef) == 0
  %learning had no negatives
  wex = m.model.w;
  b = m.model.b;
  fprintf(1,'reverting to old model...\n');

else

  %convert support vectors to decision boundary
  svm_weights = full(sum(svm_model.SVs .* ...
                         repmat(svm_model.sv_coef,1, ...
                                size(svm_model.SVs,2)),1));

  wex = svm_weights';
  b = svm_model.rho;

  %With libsvm, the first datapoint is interpreted as class +1, so if
  %we started with a negative data point, we have to flip everything.
  if newy(1) == -1
    wex = wex*-1;
    b = b*-1;
  end

  wex2 = zeros(size(newx,1),1);
  wex2(m.model.mask) = wex;

  wex = wex2;

  %% issue a warning if the norm is very small
  if norm(wex) < .00001
    fprintf(1,'learning broke down!\n');
  end
end

maxpos = max(wex'*newx - b);
fprintf(1,' --- Max positive is %.3f\n',maxpos);

fprintf(1,'took %.3f sec\n',toc(starttime));

m.model.w = reshape(wex, size(m.model.w));
m.model.b = b;

%Take top train_max_negatives detections
r = m.model.w(:)'*bsxfun(@minus,m.model.svxs,m.model.x(:,1)).^2 - m.model.b;

[alpha,beta] = sort(r,'descend');
svs = beta(1:min(length(beta),mining_params.train_max_negatives_in_cache));
m.model.svxs = m.model.svxs(:,svs);
m.model.svbbs = m.model.svbbs(svs,:);

s = length(svs);
fprintf(1,' ---length of kept svs in dfun is: %d\n',s);

author: quantombone

repo: exemplarsvm


text:
function SO3F = conj(SO3F)
% Construct the complex conjugate function $\overline{f}$ of an SO3Fun $f$.
%
% Syntax
%   SO3F = conj(F)
%
% Input
%  F - @SO3Fun
%
% Output
%  SO3F - @SO3FunHarmonic
%

SO3F = SO3FunHandle(@(rot) conj(SO3F.eval(rot)),SO3F.SRight,SO3F.SLeft);
% SO3F = SO3FunHarmonic(SO3F);
% SO3F = conj(SO3F);

end

author: mtex-toolbox

repo: mtex

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


text:
function out = makeLabel(type, label, value)
%makeLabel Create a ADC Label.
%   label=makeLabel() Create a Label event for ADC line
%         Possible values for type are: 'SET','INC'.
%         Label may be a counter of a flag and should be one of
%         (counters) 'SLC','SEG','REP','AVG','SET','ECO','PHS','LIN','PAR',
%         (flags)    'NAV','REV','SMS'.
%         Value: numeric value of the parameter (increment may be negative)
%                or true/false for a flag
%
%   See also  Sequence.addBlock, mr.getSupportedLabels

supported_labels=mr.getSupportedLabels();
if nargin~=3
    error('makeLabel:invalidArguments','Must supply exactly 3 parameters');
end
if ~any(ismember(supported_labels,label))
    error('makeLabel:invalidArguments','Must supply a valid label');
end
if ~any(ismember({'SET','INC'},type))
    error('makeLabel:invalidArguments','Must supply a valid type');
end
if ~isnumeric(value) && ~islogical(value)
    error('makeLabel:invalidArguments','Must supply a valid numeric or logical value');
end

switch(type)
    case 'SET'
        out.type = 'labelset';
    case 'INC'
        out.type = 'labelinc';
    otherwise
        error('unknown label type');
end
out.label=label;
out.value=value;
end


text:
function [S,f,R,Serr]=mtspectrumpt(data,params,fscorr,t)
% Multi-taper spectrum - point process times
%
% Usage:
%
% [S,f,R,Serr]=mtspectrumpt(data,params,fscorr,t)
% Input:
%       data        (structure array of spike times with dimension channels/trials;
%                   also accepts 1d array of spike times) -- required
%       params: structure with fields tapers, pad, Fs, fpass, err, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%           trialave (average over channels/trials when 1, don't average when 0) - optional. Default 0
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
%                1 (use finite size corrections) - optional
%                (available only for spikes). Defaults 0.
%       t        (time grid over which the tapers are to be calculated:
%                      this argument is useful when calling the spectrum
%                      calculation routine from a moving window spectrogram
%                      calculation routine). If left empty, the spike times
%                      are used to define the grid.
% Output:
%       S       (spectrum with dimensions frequency x channels/trials if trialave=0;
%               dimension frequency if trialave=1)
%       f       (frequencies)
%       R       (rate)
%       Serr    (error bars) - only if err(1)>=1
if nargin < 1; error('Need data'); end;
if nargin < 2; params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear params
data=change_row_to_column(data);
if nargout > 3 && err(1)==0; error('cannot compute error bars with err(1)=0; change params and run again'); end;
if nargin < 3 || isempty(fscorr); fscorr=0;end;
if nargin < 4 || isempty(t);
   [mintime,maxtime]=minmaxsptimes(data);
   dt=1/Fs; % sampling time
   t=mintime-dt:dt:maxtime+dt; % time grid for prolates
end;
N=length(t); % number of points in grid for dpss
nfft=max(2^(nextpow2(N)+pad),N); % number of points in fft of prolates
[f,findx]=getfgrid(Fs,nfft,fpass); % get frequency grid for evaluation
tapers=dpsschk(tapers,N,Fs); % check tapers
[J,Msp,Nsp]=mtfftpt(data,tapers,nfft,t,f,findx); % mt fft for point process times
S=squeeze(mean(conj(J).*J,2));
if trialave; S=squeeze(mean(S,2));Msp=mean(Msp);end;
R=Msp*Fs;
if nargout==4;
   if fscorr==1;
      Serr=specerr(S,J,err,trialave,Nsp);
   else
      Serr=specerr(S,J,err,trialave);
   end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function hf_out = lhs_operation(hf, samplesf, reg_filter, sample_weights, feature_reg)

% This is the left-hand-side operation in Conjugate Gradient

% size of the padding
num_features = length(hf);
output_sz = [size(hf{1},1), 2*size(hf{1},2)-1];
pad_sz = cellfun(@(hf) (output_sz - [size(hf,1), 2*size(hf,2)-1]) / 2, hf, 'uniformoutput',false);

% Compute the operation corresponding to the data term in the optimization
% (blockwise matrix multiplications)
%implements: A' diag(sample_weights) A f

% sum over all features in each block
sh_cell = cellfun(@(hf,samplesf) mtimesx(samplesf, permute(hf, [3 4 1 2]), 'speed'), hf, samplesf, 'uniformoutput', false);

% sum over all feature blocks
sh = sh_cell{1};    % assumes the feature with the highest resolution is first
for k = 2:num_features
    sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) = ...
        sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) + sh_cell{k};
end

% weight all the samples
sh = bsxfun(@times,sample_weights,sh);

% multiply with the transpose
hf_out = cellfun(@(samplesf,pad_sz) permute(conj(mtimesx(sh(:,1,1+pad_sz(1):end-pad_sz(1), 1+pad_sz(2):end), 'C', samplesf, 'speed')), [3 4 2 1]), ...
    samplesf, pad_sz, 'uniformoutput', false);

% compute the operation corresponding to the regularization term (convolve
% each feature dimension with the DFT of w, and the tramsposed operation)
% add the regularization part

reg_pad = cellfun(@(hf, reg_filter) min(size(reg_filter,2)-1, size(hf,2)-1), hf, reg_filter, 'uniformoutput', false);

% add part needed for convolution
hf_conv = cellfun(@(hf,reg_pad) cat(2, hf, conj(rot90(hf(:, end-reg_pad:end-1, :), 2))), hf, reg_pad, 'uniformoutput', false);

% do first convolution
hf_conv = cellfun(@(hf_conv, reg_filter) convn(hf_conv, reg_filter), hf_conv, reg_filter, 'uniformoutput', false);

% do final convolution and put toghether result
hf_out = cellfun(@(hf_conv, hf_data, reg_filter, hf, reg_pad, feature_reg) hf_data + convn(hf_conv(:,1:end-reg_pad,:), reg_filter, 'valid') + feature_reg * hf, ...
    hf_conv, hf_out, reg_filter, hf, reg_pad, feature_reg, 'uniformoutput', false);

end

author: martin-danelljan

repo: Continuous-ConvOp


text:
function FPDF = TabulatedPDF (x, p)
% Return function handles to routines to calculate the area, mean, and
% second moment of a unit-variance, zero-mean, uniform probability
% density function.
%                 b
%   Farea(a,b) = Int p(x) dx
%                 a
%                 b
%   Fmean(a,b) = Int x p(x) dx
%                 a
%                b
%   Fvar(a,b) = Int x^2 p(x) dx
%                a
% where p(x) is a tabulated function. The PDF is assumed to be zero outside
% of the limits, and defined by linear interpolation between points. The
% tabulated function need not be initially normalized to unit area.
%
% x - Abscissa values (in increasing order)
% p - Vector of probability values at the corresponding abscissa values.
%     The PDF is assumed to be zero outside of [x(1), x(end)] and to be
%     linearly interpolated between abscissa points.

% global xpdf pdf

% Test for increasing x
if (any(diff(x) < 0))
  error('TabulatedPDF - Non-increasing abscissa values');
end
if (any(p < 0))
  error('TabulatedPDF - Negative probability value');
end
if (length(x) ~= length(p) || length(x) <= 1)
  error('TabulatedPDF - Invalid table length');
end

% Normalize the tabulated data
xpdf = x;
pdf = p;
A = Tarea(-Inf, Inf);
pdf = pdf / A;

FPDF = {@Tarea, @Tmean, @Tvar};

return

%----- ----- begin nested functions
function v = Tarea (a, b)

% global xpdf pdf

v = TabulatedPDFInt(@TLarea, a, b, xpdf, pdf);

return
end

% ----- ------
function v = Tmean (a, b)

% global xpdf pdf

v = TabulatedPDFInt(@TLmean, a, b, xpdf, pdf);

return
end

% ----- ------
function v = Tvar (a, b)

% global xpdf pdf

v = TabulatedPDFInt(@TLvar, a, b, xpdf, pdf);

return
end

% ----- -----
function v = TabulatedPDFInt (Fn, a, b, x, p)
% The abscissa values are assumed to be in increasing order

% Limit the evaluation interval and flip the limits if a > b
br = min(x(end), max(a, b));
ar = max(x(1), min(a, b));

if (br <= ar)
  v = 0;
  return
end

% Search for the intervals included in the integral
iL = max(find(x <= ar));
iU = min(find(x >= br));
if (iL == iU)
  error('TabulatedPDFInt - iL == iU');
end
x = x(iL:iU);       % Keep only intervals of interest
p = p(iL:iU);

% Lop off the end intervals
p(1) = Alin(ar, x(1:2), p(1:2));
x(1) = ar;
p(end) = Alin(br, x(end-1:end), p(end-1:end));
x(end) = br;

% Integrate
v = feval(Fn, x, p);

% Negate the integral if a > b
if (a > b)
  v = -v;
end

return
end

% ----- -----
function v = TLarea (x, p)

xl = x(1:end-1);
xu = x(2:end);
pl = p(1:end-1);
pu = p(2:end);
% syms p pl pu xl xu x
% p = pl+(x-xl)*(pu-pl)/(xu-xl)
% factor(int(p,xl,xu)
% ans = -1/2*(pl+pu)*(-xu+xl)
A = (xu-xl) .* (pl+pu);

v = 0.5 * sum(A);

return
end

% ----- ------
function v = TLmean (x, p)

xl = x(1:end-1);
xu = x(2:end);
pl = p(1:end-1);
pu = p(2:end);
% syms p pl pu xl xu x
% p = pl+(x-xl)*(pu-pl)/(xu-xl)
% factor(int(x*p,xl,xu)
% ans = -1/6*(-xu+xl)*(2*pl*xl+pu*xl+pl*xu+2*pu*xu)
A = (xu-xl) .* ((xl+xu) .* (pl+pu) + pl.*xl + pu.*xu);

v = sum(A) / 6;

return
end

% ----- -----
function v = TLvar (x, p)

xl = x(1:end-1);
xu = x(2:end);
pl = p(1:end-1);
pu = p(2:end);
% syms p pl pu xl xu x
% p = pl+(x-xl)*(pu-pl)/(xu-xl)
% factor(int(x^2*p,xl,xu)
% ans = -1/12*(-xu+xl)
%          *(3*pl*xl^2+pu*xl^2+2*xl*pu*xu+2*pl*xu*xl+pl*xu^2+3*pu*xu^2)
A =  (xu-xl) .* ((pu+pl) .* (xu+xl).^2 + 2*(pu.*xu.^2 + pl.*xl.^2));

v = sum(A) / 12;

return
end

% ----- -----
function pa = Alin(a, x, p)

slope = (p(2) - p(1)) / (x(2) - x(1));
pa = (a - x(1)) * slope + p(1);

return
end

% ---- end of the nested functions
end


text:
%FEATSETC Set classifier

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

function [out1,out2] = featsetc(a,objclassf,fsetindex,fsetcombc,fsetclassf,fsetlab)

error('featsetc has been replaced by bagc')


text:
function [data] = fif_setup_raw(sFile, fid, allow_maxshield)
% FIF_SETUP_RAW:  Read information about raw data file
%
% USAGE:  [data] = fif_setup_raw(sFile, fid, allow_maxshield)
%
% INPUT:
%    - sFile           : Brainstorm structure for importing files
%    - fid             : Pointer to an open file
%    - allow_maxshield : Accept unprocessed MaxShield data

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
%   Copyright 2006
%
%   Matti Hamalainen
%   Athinoula A. Martinos Center for Biomedical Imaging
%   Massachusetts General Hospital
%   Charlestown, MA, USA
%
%   No part of this program may be photocopied, reproduced,
%   or translated to another program language without the
%   prior written consent of the author.
%
%   Adaptations for Brainstorm by Francois Tadel, 2009-2019
%   Section about file links copied from MNE-Python 0.19.1

%% ===== PARSE INPUT =====
% FIF Constants
global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end
% Other constants
FIFFT_SHORT=2;
FIFFT_INT=3;
FIFFT_FLOAT=4;
FIFFT_DAU_PACK16=16;
me = 'MNE-BST:fif_setup_raw';
% Arguments
if (nargin < 3)
    allow_maxshield = 1;
end
% Get file id
info = sFile.header.info;
meas = sFile.header.meas;

%% ===== LOCATE DATA OF INTEREST =====
raw = fiff_dir_tree_find(meas,FIFF.FIFFB_RAW_DATA);
if isempty(raw) && allow_maxshield
    raw = fiff_dir_tree_find(meas,FIFF.FIFFB_IAS_RAW_DATA);
    %warning('Reading FIFFB_IAS_RAW_DATA.');
end
if isempty(raw)
    raw = fiff_dir_tree_find(meas,FIFF.FIFFB_CONTINUOUS_DATA);
    disp([10 '--------' 10 ...
          'WARNING: This file contains raw Internal Active Shielding data. It may be distorted.' 10 ...
          'Elekta recommends it be run run through MaxFilter to produce reliable results.' 10 ...
          'Consider closing the file and running MaxFilter on the data.' 10 ...
          '--------' 10]);
end
if isempty(raw)
    error(me,'No raw data in file');
end

%
%   Set up the output structure
%
data.first_samp = 0;
data.last_samp  = 0;
%
%   Process the directory
%
dir          = raw.dir;
nent         = raw.nent;
first        = 1;
first_samp   = 0;
first_skip   = 0;
%
%  Get first sample tag if it is there
%
if dir(first).kind == FIFF.FIFF_FIRST_SAMPLE
    tag = fiff_read_tag(fid,dir(first).pos);
    first_samp = tag.data;
    first = first + 1;
end
%
%  Omit initial skip
%
if dir(first).kind == FIFF.FIFF_DATA_SKIP
    %
    %  This first skip can be applied only after we know the buffer size
    %
    tag = fiff_read_tag(fid,dir(first).pos);
    first_skip = tag.data;
    first = first + 1;
end
data.first_samp = first_samp;
%
%   Go through the remaining tags in the directory
%
rawdir = struct('ent',{},'first',{},'last',{},'nsamp',{});
nskip = 0;
ndir  = 0;
for k = first:nent
    ent = dir(k);
    if ent.kind == FIFF.FIFF_DATA_SKIP
        tag = fiff_read_tag(fid,ent.pos);
        nskip = tag.data;
    elseif ent.kind == FIFF.FIFF_DATA_BUFFER
        %
        %   Figure out the number of samples in this buffer
        %
        switch ent.type
            case FIFFT_DAU_PACK16
                nsamp = ent.size/(2*info.nchan);
            case FIFFT_SHORT
                nsamp = ent.size/(2*info.nchan);
            case FIFFT_FLOAT
                nsamp = ent.size/(4*info.nchan);
            case FIFFT_INT
                nsamp = ent.size/(4*info.nchan);
            otherwise
                fclose(fid);
                error(me,'Cannot handle data buffers of type %d',ent.type);
        end
        %
        %  Do we have an initial skip pending?
        %
        if first_skip > 0
            first_samp = first_samp + nsamp*first_skip;
            data.first_samp = first_samp;
            first_skip = 0;
        end
        %
        %  Do we have a skip pending?
        %
        if nskip > 0
            ndir        = ndir+1;
            rawdir(ndir).ent   = [];
            rawdir(ndir).first = first_samp;
            rawdir(ndir).last  = first_samp + nskip*nsamp - 1;
            rawdir(ndir).nsamp = nskip*nsamp;
            first_samp = first_samp + nskip*nsamp;
            nskip = 0;
        end
        %
        %  Add a data buffer
        %
        ndir               = ndir+1;
        rawdir(ndir).ent   = ent;
        rawdir(ndir).first = first_samp;
        rawdir(ndir).last  = first_samp + nsamp - 1;
        rawdir(ndir).nsamp = nsamp;
        first_samp = first_samp + nsamp;
    end
end
data.last_samp  = first_samp - 1;

data.rawdir = rawdir;
%
% fprintf(1,'\tRange : %d ... %d  =  %9.3f ... %9.3f secs\n',...
%     data.first_samp,data.last_samp,...
%     double(data.first_samp)/info.sfreq,...
%     double(data.last_samp)/info.sfreq);

%% ===== GET LINK TO NEXT FILE =====
% Initialize references
data.next_fname = [];
data.next_num = [];
data.next_id = [];
% Find the reference block
ref = fiff_dir_tree_find(meas,FIFF.FIFFB_REF);
if ~isempty(ref)
    for iBlock = 1:length(ref)
        for k = 1:length(ref(iBlock).dir)
            ent = ref(iBlock).dir(k);
            switch (ent.kind)
                case FIFF.FIFF_REF_ROLE
                    % Check the role of the reference: accept only "next file"
                    tag = fiff_read_tag(fid, ent.pos);
                    if (tag.data ~= FIFF.FIFFV_ROLE_NEXT_FILE)
                        break;
                    end
                case FIFF.FIFF_REF_FILE_NAME
                    % Get filename of the next linked file
                    tag = fiff_read_tag(fid, ent.pos);
                    data.next_fname = tag.data;
                case FIFF.FIFF_REF_FILE_NUM
                    % Some files don't have the name, just the number. So we construct the name from the current name.
                    tag = fiff_read_tag(fid, ent.pos);
                    data.next_num = tag.data;
                case FIFF.FIFF_REF_FILE_ID
                    tag = fiff_read_tag(fid, ent.pos);
                    data.next_id = tag.data;
            end
        end
    end
end

author: brainstorm-tools


text:
function [sp, dcell] = deconvolution_standalone2(ops, ca, neu, gt)
% takes as input the calcium and (optionally) neuropil traces,
% both NT by NN (number of neurons).
% outputs a cell array dcell containing spike times (dcell.st) and amplitudes
% (dcell.c). dcell.B(3) is the neuropil contamination coefficient.
% dcell.B(2) is an estimate of the baseline.

% this version also estimates a single timescale

% specify in ops the following options, or leave empty for defaults
%       fs = sampling rate
%       recomputeKernel = whether to estimate kernel from data
%       sensorTau  = timescale of sensor, if recomputeKernel = 0
% additional options can be specified (mostly for linking with Suite2p, see below).

ops.imageRate    = getOr(ops, {'imageRate'}, 30); % total image rate (over all planes)
ops.nplanes      = getOr(ops, {'nplanes'}, 1); % how many planes at this total imaging rate
ops.sensorTau    = getOr(ops, {'sensorTau'}, 2); % approximate timescale in seconds
ops.sameKernel   = getOr(ops, {'sameKernel'}, 1); % 1 for same kernel per plane, 0 for individual kernels (does not work right now)
ops.maxNeurop    = getOr(ops, {'maxNeurop'}, Inf); % maximum allowed neuropil contamination coefficient.
ops.recomputeKernel = getOr(ops, {'recomputeKernel'}, 1); % whether to estimate kernel from data
lam = getOr(ops, 'lam', 3);

% the kernel should depend on timescale of sensor and imaging rate
ops.fs           = getOr(ops, 'fs', ops.imageRate/ops.nplanes);
mtau             = ops.fs * ops.sensorTau;

if nargin<3 || isempty(neu)
    neu = zeros(size(ca));
end

Params = [1 lam 1 2e4]; %parameters of deconvolution

% f0 = (mtau/2); % resample the initialization of the kernel to the right number of samples
kernel = exp(-[1:ceil(5*mtau)]'/mtau) ;
%
npad        = 250;
[NT, NN]    = size(ca);
coefNeu     = .8 * ones(1,NN); % initialize neuropil subtraction coef with 0.8

caCorrected = ca - bsxfun(@times, neu, coefNeu);

if ops.recomputeKernel
    tlag                     = 1;
    [kernel, mtau]           = estimateKernel(ops, caCorrected, tlag);

%     [kernel, mtau, coefNeu]  = estimateKernel(ops, ca - coefNeu * neu, tlag);

    fprintf('Timescale determined is %4.4f samples \n', mtau);
end
kernel = normc(kernel(:));
%%
kernelS     = repmat(kernel, 1, NN);
dcell       = cell(NN,1);

% running minimum subtraction here???

Fsort       = my_conv2(caCorrected, ceil(ops.fs), 1);

Fsort       = sort(Fsort, 1, 'ascend');
baselines   = Fsort(ceil(NT/20), :);

% determine and subtract the neuropil
F1 = caCorrected - bsxfun(@times, ones(NT,1), baselines);

% normalize signal
% sd   = 1/2 * std(F1 - my_conv2(F1, max(2, ops.fs/4), 1), [], 1);

sd   = 1/2 * 1/sqrt(2) * std(F1(2:end, :) - F1(1:end-1, :), [], 1);

F1   = bsxfun(@rdivide, F1 , 1e-12 + sd);

if 1
    % get new kernel
    kerns = exp(-bsxfun(@rdivide, [1:ceil(5*mtau)]',  mtau * [.05 .125 .25 .5 1 2]));
    spfilt = ones([size(kerns,2) size(F1)]);
    for j = 1:size(kerns,2)
        spfilt(j, :,:) = filter(kerns(:,j), 1, gt);
    end

    sts = zeros(size(spfilt,1));
    stF = zeros(size(spfilt,1), 1);
    for k = 1:size(spfilt,3)
        sts = sts + spfilt(:,:,k) * spfilt(:,:,k)';
        stF = stF + spfilt(:,:,k) * F1(:,k);
    end
    coefs = sts \ stF;
    kernel = normc(kerns * coefs);
    kernelS = repmat(kernel, 1, NN);

    plot(kernel)
    drawnow
end

sp = zeros(size(F1));

% run the deconvolution to get fs etc\
parfor icell = 1:size(ca,2)
    [sp(:,icell),dcell{icell}] = ...
        single_step_single_cell(F1(:,icell), Params, kernelS(:,icell), NT, npad,dcell{icell});
end

% rescale baseline contribution
for icell = 1:size(ca,2)
    dcell{icell}.c                      = dcell{icell}.c * sd(icell);
    dcell{icell}.baseline               = baselines(icell);
    dcell{icell}.neuropil_coefficient   = coefNeu(icell);
end

%%

author: cortex-lab


text:
function [Jh,dJh] = calcJacobian(obj, x, dx)
    % calculate the holonomic constraints, Jacobian and its time
    % derivatives
    %
    % Parameters:
    % x: the states @type colvec
    % dx: the first order derivatives @type colvec
    %
    % Return values:
    %  Jh: the first-order Jacobian of the constraint @type matrix
    %  dJh: the time derivative of the Jacobian matrix @type matrix

    switch obj.Model.Type
        case 'FirstOrder'
            if obj.DerivativeOrder == 1
                Jh = feval(obj.Jh_name, x);
                dJh = [];
            else
                Jh = feval(obj.Jh_name, x);
                dJh = feval(obj.dJh_name, x);
            end
        case 'SecondOrder'
            Jh = feval(obj.Jh_name, x);
            dJh = feval(obj.dJh_name, x, dx);
    end

end

author: ayonga

repo: frost-dev

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


text:
function aIdx = getPaddingIndices22(aSize,padSize,method,direction)
%getPaddingIndices is used by padarray and blockproc.
%   Computes padding indices of input image.  This is function is used to
%   handle padding of in-memory images (via padarray) as well as
%   arbitrarily large images (via blockproc).
%
%   aSize : result of size(I) where I is the image to be padded
%   padSize : padding amount in each dimension.
%             numel(padSize) can be greater than numel(aSize)
%   method : X or a 'string' padding method
%   direction : pre, post, or both.
%
%   See the help for padarray for additional information.

% Copyright 2010 The MathWorks, Inc.

% make sure we have enough image dims for the requested padding
if numel(padSize) > numel(aSize)
    singleton_dims = numel(padSize) - numel(aSize);
    aSize = [aSize ones(1,singleton_dims)];
end

switch method
    case 'circular'
        aIdx = CircularPad(aSize, padSize, direction);
    case 'symmetric'
        aIdx = SymmetricPad(aSize, padSize, direction);
    case 'replicate'
        aIdx = ReplicatePad(aSize, padSize, direction);
end

%%%
%%% CircularPad
%%%
function idx = CircularPad(aSize, padSize, direction)

numDims = numel(padSize);

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
for k = 1:numDims
    M = aSize(k);
    dimNums = uint32(1:M);
    p = padSize(k);

    switch direction
        case 'pre'
            idx{k}   = dimNums(mod(-p:M-1, M) + 1);

        case 'post'
            idx{k}   = dimNums(mod(0:M+p-1, M) + 1);

        case 'both'
            idx{k}   = dimNums(mod(-p:M+p-1, M) + 1);

    end
end

%%%
%%% SymmetricPad
%%%
function idx = SymmetricPad(aSize, padSize, direction)

numDims = numel(padSize);

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
for k = 1:numDims
    M = aSize(k);
    dimNums = uint32([1:M M:-1:1]);
    p = padSize(k);

    switch direction
        case 'pre'
            idx{k}   = dimNums(mod(-p:M-1, 2*M) + 1);

        case 'post'
            idx{k}   = dimNums(mod(0:M+p-1, 2*M) + 1);

        case 'both'
            idx{k}   = dimNums(mod(-p:M+p-1, 2*M) + 1);
    end
end

%%%
%%% ReplicatePad
%%%
function idx = ReplicatePad(aSize, padSize, direction)

numDims = numel(padSize);

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
for k = 1:numDims
    M = aSize(k);
    p = padSize(k);
    onesVector = uint32(ones(1,p));

    switch direction
        case 'pre'
            idx{k}   = [onesVector 1:M];

        case 'post'
            idx{k}   = [1:M M*onesVector];

        case 'both'
            idx{k}   = [onesVector 1:M M*onesVector];
    end
end

repo: Deep-ADMM-Net


text:
%% Bonus Examples
% Tucker McClure
% Copyright 2013 The MathWorks, Inc.

%% Bonus 1: Please
% If it doesn't work, try saying "please".

please = @(f) feval(curly({fprintf('Since you asked nicely...\n'), f}, 2));

% Example:
please(@why)

%% Bonus 2: Anonymous Do-While
% In the examples, we didn't show how to create a do-while loop
% anonymously. Here it is.

dowhile = @(x0, cont, fcn, cleanup) recur(... % Make the recursion.
    @(f, x) use(fcn(x), @(z) ...              %   Run the fcn, save as z.
                iif(cont(z), @() f(f, z), ... %   Continue?
                    true,    cleanup(z))),... % Stop, clean up.
    x0);                                      % Initial state

% Create a function to draw sets of numbers until it draws n positive
% numbers in a row.
draw_until_all_positive = @(n) dowhile([], ...               % No x0
                                       @(x) any(x < 0), ...  % Until all >0
                                       @(~) randn(1, n), ... %   Draw set
                                       @(x) x);              % Return set

fprintf('Draw numbers until we get 5 positive numbers in a row.\n');
draw_until_all_positive(5)

%% Bonus 3: A Large Program
% This is probably not the world's best example of good programming
% practices, but it does serve to illustrate that entire programs *can* be
% built with anonymous functions (whether or not they *should* be built
% this way or whether this is the best way to do it are very different
% questions).
%
% This example uses most of the techniques developed in the
% functional_programming_examples (though it uses the .m versions instead
% of the anonymous versions), including |forloop|, |iif|, |last|, |use|,
% |recur|, |tern|, and |void|. Creating this example also served to prompt
% me to make functions that simplify programming this way.

% Set up the figure.
figure(1);
clf();
set(1, 'Color', 'k');
axis([-20 20 0 25], 'off');
set(gca(), 'Position', [0 0 1 1]);
hold('on');

% Create some constants.
nf = 6;       % Number of fireworks
np = 100;     % Number of particles per firework
n  = nf * np; % Total number of particles
dt = 1/24;    % Time step

% Create some fireworks.
create_fireworks = @() forloop([], nf, @(p,~) ...
    use(2*randn, rand, @(x, y) ...
    [p repmat([2*x, -5*y, 5*x, 20*y+30]', 1, np)]));

% Plot the fireworks.
random_color = @() 1./(1+exp(-20*(rand(1, 3)-0.4))); % It's a sigmoid fcn.
create_plot = @(x0) forloop([], nf, @(h,~) ...
    [h plot(x0(1,:), x0(2,:), '.', 'Color', random_color())]);
limit = @(x) min(max(x, 0), 1);
update_plot = @(x, h, k) forloop([], nf, @(~,z) ...
    void(@() set(h(z), 'XData', x(1, np*(z-1)+1:np*z), ...
                       'YData', x(2, np*(z-1)+1:np*z), ...
                       'Color', limit(1.03-1e-3*k) * get(h(z),'Color'))));

% Propagate the fireworks' particles. At k == 25, make them explode.
explode  = @(x) x + [zeros(2, n); 5*randn(2, n)];
Phi      = expm(dt*[0 0 1 0; 0 0 0 1; 0 0 -1 0; 0 0 0 -1]);
g        = repmat([0 -dt*9.81 0 -0.5*9.81*dt^2]', 1, n);
update_x = @(k, x) Phi * tern(k==25, @() explode(x), @() x) + g;

% Functions to change callback and show text.
set_callback = @(cb) void(@() set(gcf, 'WindowButtonDownFcn', cb));
show_text    = @() text(0, 1, 'Click for fireworks.', ...
                        'HorizontalAlignment', 'center', ...
                        'Color', 'w');

% Create the anonymous program.
fireworks = @() recur(@(self) {...                 % To reference itself.
    set_callback([]), ...                          % Clear the callback.
    void(@() delete(findobj('Type', 'text'))), ... % Delete the text.
    use(create_fireworks(), @(x0) ...              % Create positions, x0.
    use(create_plot(x0),    @(h) ...               % Plot, keep handle, h.
    {forloop(x0, 100, @(x,k) ...                   % Loop from k=1:100.
             last(update_plot(x,h,k), ...          %   Update plot.
                  void(@() pause(dt)), ...         %   Wait.
                  {update_x(k,x)})), ...           %   Update state.
	 void(@() delete(h))} ...                      % Delete fireworks.
    )), ...                                        % (End of x0, h scope.)
	show_text(), ...                               % Create text again.
    set_callback(@(~, ~) self(self))});            % Restore callback.

% Execute the anonymous program, which will also attach 'fireworks' to the
% button click callback.
fireworks();

author: Sable

repo: mcbench-benchmarks


text:
function out = simplifyMap(map, scalling)
% rescale map

mapB = map;
[a,b] = size(mapB);
out = zeros(a/scalling,b/scalling);
[c,d] = size(out);

s = scalling;
a = floor((scalling-1)/2);
b = scalling - a;

for x=1:1:c-1
    for y=1:1:d-1
        tmp = sum( sum(  mapB( x*s-a:x*s+b, y*s-a:y*s+b  )   ));
        if tmp<12
            out(x,y) = 0;
        else
            out(x,y) = 1;
        end
    end
end

author: LazyFalcon


text:
% GIfTI Class for MATLAB
%
% Geometry format under the Neuroimaging Informatics Technology Initiative
% (NIfTI):
%                 http://www.nitrc.org/projects/gifti/
%                      http://nifti.nimh.nih.gov/
%
% This MATLAB class is part of SPM:
%                  http://www.fil.ion.ucl.ac.uk/spm/
%
% It relies on external libraries:
%   Base64, by Peter J. Acklam:
%      http://home.online.no/~pjacklam/
%   miniz, by Rich Geldreich:
%      http://code.google.com/p/miniz/
%   XMLTree, by Guillaume Flandin:
%      http://www.artefact.tk/software/matlab/xml/
%__________________________________________________________________________
% Copyright (C) 2008-2015 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: Contents.m 6404 2015-04-13 14:29:53Z guillaume $

% GIfTI file format for MATLAB (The Mathworks, Inc.).
% Copyright (C) 2008-2015 Wellcome Trust Centre for Neuroimaging
%
% 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 any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation Inc, 59 Temple Pl. - Suite 330, Boston, MA 02111-1307, USA.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
row=112;
column=92;
NN=row*column;

Class_Train_NUM=5;
Class_Sample_NUM=10; % total
Class_Test_NUM=Class_Sample_NUM-Class_Train_NUM;

Class_NUM=40;
Train_NUM=Class_NUM*Class_Train_NUM; %
Test_NUM=Class_NUM*(Class_Sample_NUM-Class_Train_NUM); %

Eigen_NUM=10;
Disc_NUM=Eigen_NUM;

Train_DAT=zeros(NN,Train_NUM);
s=1;
for r=1:Class_NUM
   for t=1:Class_Train_NUM
     string=['E:\ORL_face\orlnumtotal\s' int2str(r) '_' int2str(t)];
     A=imread(string,'bmp');
     B=im2double(A);
     Train_DAT(:,s)=B(:);
     s=s+1;
   end
end

Test_DAT=zeros(NN,Test_NUM);
s=1;
for r=1:Class_NUM
   for t=Class_Train_NUM+1:Class_Sample_NUM
     string=['E:\ORL_face\orlnumtotal\s' int2str(r) '_' int2str(t)];
     A=imread(string,'bmp');
     B=im2double(A);
     Test_DAT(:,s)=B(:);
     s=s+1;
   end
end

% to center the each training sample and testing sample
% !!! Note that: Centralization have great effection when
%  Cos distance is used, but it has no impact when L2 or L1 distance is used
Mean_Image=mean(Train_DAT,2);
Train_DAT=Train_DAT-Mean_Image*ones(1,Train_NUM);
Test_DAT=Test_DAT-Mean_Image*ones(1,Test_NUM);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Projection,disc_value]=Eigenface_f(Train_DAT, Disc_NUM);
Eigen_face=zeros(NN,10);

for k=1:1:10
    Eigen_face(:,k)=Projection(:,k);
    Eigen_face(:,k)=mat2gray(Eigen_face(:,k));
end

Eigen_face=reshape(Eigen_face,[row,column,10]);

for count=1:1:10
    subplot(1,10,count);imshow(Eigen_face(:,:,count))
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [data,age,ndiscarded] = onl_lslreceive(inlet,timeout,maxage)
if nargin<2
    timeout = -1;
end
if nargin<3
    maxage = Inf;
end
ndiscarded = 0;

[tmpsample,timestamp] = inlet.pull_sample(timeout,true);

if isempty(timestamp) || timestamp == 0
    data = [];
    age = NaN;
    return;
end

% tic;
% calculate the age of the sample
age = lsl_local_clock(inlet.get_libhandle()) - (timestamp+inlet.time_correction());
% fprintf('age: %0.5g\n',age);
while age > maxage
    % discard samples that have expired
    [tmpsample,timestamp] = inlet.pull_sample(timeout,true);
    age = lsl_local_clock(inlet.get_libhandle()) - (timestamp+inlet.time_correction());
    ndiscarded = ndiscarded + 1;
end
% toc

if isempty(tmpsample)
    data = [];
    age = NaN;
    return;
end

% now we have a sufficiently current sample
% ...decode the data
% tic;
data = hlp_deserialize(tmpsample{1});
% toc

repo: REST


text:
%FbaLPWriter test
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updaetd 2/3/2011
classdef FbaLPWriter_Test < TestCase
    methods
        function this = FbaLPWriter_Test(name)
            this = this@TestCase(name);
        end
    end

    methods
        %checks for run time errors
        function testWriter(~)
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], true);
            edu.stanford.covert.cell.sim.util.FbaLPWriter.write(sim, 'output/runSmallTests');
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell