text:
function y = ArbitraryOctaveFilt(x, SPECT, FREQS, N, fs, octBandwidth)
% Filters a signal with any arbitrary spectrum smoothed with any fractional octave band average
%
% Syntax:	Y = ARBITRARYOCTAVEFILT(X, SPECT, FREQS, N, FS, OCTBANDWIDTH)
%
% Inputs:
% 	x - Input signal to filter as a vector
% 	SPECT - The spectrum to shape the input signal to
% 	FREQS - The frequencies of each SPECT element
% 	N - The length of the filter to usee
% 	fs - Description
% 	octBandwidth - Description
%
% Outputs:
% 	y - Description
%
% Example:
%   fs = 16000;
%   T = 10;
%   N = 1000;
%   f = linspace(0,fs/2,N);
%   s = 1./f;
%   x = wgn(T*fs,1,0);
%   y = ArbitraryOctaveFilt(x,s,f,N,fs,1/3);
%   pwelch([x y]);
%
% See also: fir2, filter

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2017
% Date: 06 June 2016
% Revision: 0.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 6
    octBandwidth = 1/6;
end

% Find nth-octave averages
[MAG,f]=Tools.octaveBandMean(SPECT,FREQS,octBandwidth);

% Force even length filter
if isempty(N), if mod(length(SPECT),2), N=length(SPECT)-1; else N=length(SPECT); end; end
% Design arbitrary magnitude (linear-phase) filter
b = fir2(N,f/(fs/2),MAG);
% Apply filter
y = filter(b,1,x);

end

author: jtkim-kaist

repo: Speech-enhancement


text:
%function [] = test_fig2u3d_contour_surf

clf

res = [30, 31];
q = domain2vec([-2, 2, -2, 2], res);

x = q(1, :);
y = q(2, :);
z = -sin(x) .*cos(y);

ax = gca;
hold(ax, 'on')
%vcontourf(ax, q, z, res);
vcontour(ax, q, z, res);
vsurf(ax, q, z, res);
plot_scalings(ax, 0)

%{
h = get(gca, 'children');
hp = get(h(2), 'Children');
fc = get(hp(1), 'Faces');

x = v(:, 1);
y = v(:, 2);
x(isnan(x) ) = [];
y(isnan(y) ) = [];

dt = DelaunayTri(x, y);
triplot(dt, 'Parent', newax)
%}
fig2u3d(ax, 'test')

axis(ax, 'equal')

author: Sable

repo: mcbench-benchmarks


text:
function [soln,eqn,info] = PoissonP2(node,elem,bdFlag,pde,option)
%% POISSONP2 Poisson equation: P2 quadratic element.
%
% u = PoissonP2(node,elem,bdFlag,pde,option) produces the quadratic
%   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
%
% [soln,eqn,info] = PoissonP2(node,elem,pde,bdFlag,option)
%
% The usage is the same as Poisson. Quadratic element on a triangle is
% summarized in PoissonP2femrate for detail.
%
%   Example
%     clear all
%     node = [0,0; 1,0; 1,1; 0,1];
%     elem = [2,3,1; 4,1,3];
%     for k = 1:3
%       [node,elem] = uniformrefine(node,elem);
%     end
%     % Homogenous Dirichlet boundary condition
%     pde.f = inline('ones(size(p,1),1)','p');
%     pde.g_D = 0;
%     u = PoissonP2(node,elem,[],pde);
%     figure(1);
%     showresult(node,elem,u);
%
% See also Poisson, Poisson3, Poisson3P2
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Preprocess
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end

%% Construct data structure
time = cputime;  % record assembling time
[elem2dof,edge,bdDof] = dofP2(elem);
% important constants
N = size(node,1);  NT = size(elem,1); NE = size(edge,1);
Ndof = N + NE;
elem2edge = elem2dof(:,4:6) - N;

%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix
% generate sparse pattern
ii = zeros(21*NT,1);
jj = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        index = index + NT;
    end
end
% quadrature points
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'quadorder')
    % diffusion is piecewise constant
    option.quadorder = 2;        % default order
    if ~isempty(pde.d) && isnumeric(pde.d) % numerical diffusion
        option.quadorder = 3;    % exact for linear diffusion coefficient
    end
end
[lambda, w] = quadpts(option.quadorder);
nQuad = size(lambda,1);
% compute non-zeros
sA = zeros(21*NT,nQuad);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    index = 0;
    for i = 1:6
        for j = i:6
            Aij = 0;
            if isempty(pde.d) || isnumeric(pde.d)
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
            else
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
            end
            if ~isempty(pde.d) && isnumeric(pde.d) % d is piecewise constant
                Aij = pde.d.*Aij;
            end
            Aij = Aij.*area;
            sA(index+1:index+NT,p) = Aij;
            index = index + NT;
        end
    end
end
sA = sum(sA,2);
% 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 right hand side by high order quadrature rule
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 4;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
    phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
    phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
    phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
    phi(:,4) = 4*lambda(:,2).*lambda(:,3);
    phi(:,5) = 4*lambda(:,3).*lambda(:,1);
    bt = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        if isfield(pde,'f') && isnumeric(pde.f)
            fp = pde.f;        % piecewise constant
        else
            fp = pde.f(pxy);   % function handle
        end
        for j = 1:6
            bt(:,j) = bt(:,j) + w(p)*phi(p,j)*fp;
        end
    end
    bt = bt.*repmat(area,1,6);
    b = accumarray(elem2dof(:),bt(:),[Ndof 1]);
end

%% Boundary Conditions
[AD,b,u,freeDof,isPureNeumann] = getbdP2(b);

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

%% Solve the system of linear equations
if isempty(freeDof), return; end
% Set up solver
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % Multigrid-type  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeDof) = AD(freeDof,freeDof)\b(freeDof);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        option.x0 = u;
        option.solver = 'CG';
        option.tol = 1e-9;
        [u,info] = mg(AD,b,elem,option,edge);
    case 'amg'
        option.solver = 'CG';
        option.tol = 1e-9;
        [u(freeDof),info] = amg(AD(freeDof,freeDof),b(freeDof),option);
end
% post-process for pure Neumann problem
if isPureNeumann
      patchArea = accumarray(elem2edge(:),[area;area;area]/3, [NE 1]);
             uc = sum(u(N+1:end).*patchArea)/sum(area);
              u = u - uc;   % normalization for pure Neumann problem
% Here we use 3 middle points rule to compute the integral which is exact
% for quadratic elements. We can also use the following formulae
%      intuh = sum(1/3*u(elem2dof(:,4:6)),2).*area; %Compute the integrable of uh
%         uc = sum(intuh)/sum(area);
%          u = u -uc;
end

%% Compute Du
Du = [];

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdP2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,b,u,freeDof,isPureNeumann] = getbdP2(b)
    %% Boundary conditions for Poisson equation: P2 quadratic FEM.
    %
    % The set up of boundary condition consists of two parts:
    %
    % 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(fixedDof,fixedDof)=I, AD(fixedDof,freeDof)=0, AD(freeDof,fixedDof)=0.
    %
    % 2) Modify the right hand side b. The Neumann boundary integral is added
    % to b. For Dirichlet boundary ndoes, b(fixedDof) 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
    % fixedDof.
    %
    % 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.

    u = zeros(Ndof,1);
    %% Set up boundary and basic parameter
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    idxR = (bdFlag(:) == 3);      % index of Robin edges in bdFlag
    if any(idxR)
        isRobin = false(NE,1);
        isRobin(elem2edge(idxR)) = true;
        Robin = edge(isRobin,:);  % Robin edges
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        if ~isfield(option,'gRquadorder')
            option.gRquadorder = 5;   % default order exact for linear gR
        end
        [lambdagR,weightgR] = quadpts1(option.gRquadorder);
        nQuadgR = size(lambdagR,1);
        % quadratic bases (1---3---2)
        bdphi = zeros(nQuadgR,3);
        bdphi(:,1) = (2*lambdagR(:,1)-1).*lambdagR(:,1);
        bdphi(:,2) = (2*lambdagR(:,2)-1).*lambdagR(:,2);
        bdphi(:,3) = 4*lambdagR(:,1).*lambdagR(:,2);
        % length of edge
        el = sqrt(sum((node(Robin(:,1),:) - node(Robin(:,2),:)).^2,2));
        NR = size(Robin,1);
        ss = zeros(NR,3,3);
        for pp = 1:nQuadgR
            ppxy = lambdagR(pp,1)*node(Robin(:,1),:) ...
                 + lambdagR(pp,2)*node(Robin(:,2),:);
            gRp = pde.g_R(ppxy);
            for iR = 1:3
                for jR = iR:3   % only compute half of the off-diagonal part
                    ss(:,iR,jR) = ss(:,iR,jR) + ...
                    weightgR(pp)*gRp*bdphi(pp,iR).*bdphi(pp,jR);
                end
            end
        end
        ss(:) = ss(:).*repmat(el,9,1);
        index = 0;
        Robin(:,3) = find(isRobin)+N;  % the third one maps to corresponding dof
        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 dof: fixedDof
    fixedDof = [];
    freeDof = [];
    isFixedDof = false(Ndof,1);
    if ~isempty(bdFlag)
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;
        isFixedDof(N + find(isDirichlet')) = true;
        fixedDof = find(isFixedDof);
        freeDof = find(~isFixedDof);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        fixedDof = bdDof;
        isFixedDof(fixedDof) = true;
        freeDof = find(~isFixedDof);
    end
    isPureNeumann = false;
    if isempty(fixedDof) && isempty(Robin)  % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedDof = 1;
        freeDof = 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(fixedDof,fixedDof)=I, AD(fixedDof,freeDof)=0,
    % AD(freeDof,fixedDof)=0|.
    if ~isempty(fixedDof)
        bdidx = zeros(Ndof,1);
        bdidx(fixedDof) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    else
        AD = A;
    end

    %% Part 2: Find boundary edges and modify the load b
    % Find boundary edges: Neumann and Robin
    Neumann = [];
    if ~isempty(bdFlag)
        idxN = (bdFlag(:) == 2);      % all Neumann edges in bdFlag
        Neumannidx = elem2edge(idxN | idxR); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann   = edge(Neumannidx,:);
    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
        Neumannidx = find(bdDof>N);
        Neumann = edge(Neumannidx,:);
    end

    % Neumann boundary condition
    if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 4;   % default order
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        % quadratic bases on an edge (1---3---2)
        bdphi = zeros(nQuadgN,3);
        bdphi(:,1) = (2*lambdagN(:,1)-1).*lambdagN(:,1);
        bdphi(:,2) = (2*lambdagN(:,2)-1).*lambdagN(:,2);
        bdphi(:,3) = 4*lambdagN(:,1).*lambdagN(:,2);
        % length of edge
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        ge = zeros(size(Neumann,1),3);
        for pp = 1:nQuadgN
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            ge(:,1) = ge(:,1) + weightgN(pp)*gNp*bdphi(pp,1);
            ge(:,2) = ge(:,2) + weightgN(pp)*gNp*bdphi(pp,2);
            ge(:,3) = ge(:,3) + weightgN(pp)*gNp*bdphi(pp,3); % interior bubble
        end
        % update RHS
        ge = ge.*repmat(el,1,3);
        b(1:N) = b(1:N) + accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
        b(N+Neumannidx) = b(N+Neumannidx) + ge(:,3);
    end

    % Dirichlet boundary conditions
    if ~isPureNeumann && ~isempty(fixedDof) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
        idx = (fixedDof > N);                         % index of edge nodes
        u(fixedDof(~idx)) = pde.g_D(node(fixedDof(~idx),:)); % bd value at vertex dofs
        bdEdgeIdx = fixedDof(idx) - N;
        bdEdgeMid = (node(edge(bdEdgeIdx,1),:) + node(edge(bdEdgeIdx,2),:))/2;
        u(fixedDof(idx)) = pde.g_D(bdEdgeMid);
        b = b - A*u;
    end
    if ~isPureNeumann % non-empty Dirichlet boundary condition
        b(fixedDof) = u(fixedDof);
    end

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

repo: ifem


text:
function r8sm_print ( m, n, a, u, v, title )

%*****************************************************************************80
%
%% R8SM_PRINT prints a R8SM matrix.
%
%  Discussion:
%
%    The R8SM storage format is used for an M by N Sherman Morrison matrix B,
%    which is defined by an M by N matrix A, an M vector U, and
%    an N vector V, by B = A - U * V'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(M,N), the R8SM matrix.
%
%    Input, real U(M), V(N), the R8SM vectors.
%
%    Input, string TITLE, a title to be printed.
%
  r8sm_print_some ( m, n, a, u, v, 1, 1, m, n, title );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef AoATermConditionOptimizationVariable < AbstractOptimizationVariable
    %AoATermConditionOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) AngleOfAttackTermCondition = AngleOfAttackTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = AoATermConditionOptimizationVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.targetAoA;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.targetAoA = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Angle Of Attack Termination Condition', evtNum)};
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            varsStoredInRad = true;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ net,res,opts ] = qrnn_bp( net,res,opts )
%QRNN_BP Summary of this function goes here
%   Detailed explanation goes here

    n_frames=opts.parameters.n_frames;
    n_hidden_nodes=opts.parameters.n_hidden_nodes;
    batch_size=opts.parameters.batch_size;

    %1: calculate the gradients of the data fitting transform
    opts.dzdy=res.Fit(numel(net{end}.layers)+1).dzdx;
    [net{3},res.Fit,opts] = net_bp(net{3},res.Fit,opts);

    %2: BPTT: calculate the gradient wrt the hidden nodes
    res.Fit(1).dzdx=reshape(res.Fit(1).dzdx,n_hidden_nodes,batch_size,n_frames);
    for f=n_frames-1:-1:1
        res.Fit(1).dzdx(:,:,f)=res.Fit(1).dzdx(:,:,f)+(1-res.Gates(end).x(:,:,f+1)).*res.Fit(1).dzdx(:,:,f+1);
    end

    res.Gates(end).x=reshape(res.Gates(end).x,n_hidden_nodes,[]);
    res.Input(end).x=reshape(res.Input(end).x,n_hidden_nodes,[]);
    res.Fit(1).dzdx=reshape(res.Fit(1).dzdx,n_hidden_nodes,[]);
    res.Hidden_diff=reshape(res.Hidden_diff,n_hidden_nodes,[]);
    %3: calculate the gradients of the input transform
    opts.dzdy=res.Fit(1).dzdx.*res.Gates(end).x;
    [ net{1},res.Input,opts ] = net_bp( net{1},res.Input,opts );

    %4: calculate the gradients of the gates
    opts.dzdy=res.Fit(1).dzdx.*res.Hidden_diff;
    [ net{1},res.Gates,opts ] = net_bp( net{1},res.Gates,opts );

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function t = q2tr(q)

    q = double(q);
    s = q(1);
    x = q(2);
    y = q(3);
    z = q(4);

    r = [   1-2*(y^2+z^2)   2*(x*y-s*z) 2*(x*z+s*y)
        2*(x*y+s*z) 1-2*(x^2+z^2)   2*(y*z-s*x)
        2*(x*z-s*y) 2*(y*z+s*x) 1-2*(x^2+y^2)   ];
    t = eye(4,4);
    t(1:3,1:3) = r;
    t(4,4) = 1;
endfunction

author: petercorke

repo: robotics-toolbox-matlab


text:
function [uids, patient_ids] = MimImageImporter(filename_or_root_directory, database, reporting)
    % MimImageImporter. Recrusively imports image files into a MimImageDatabase object
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    reporting.ShowProgress('Importing data');
    tags_to_get = DMDicomDictionary.GroupingDictionary;

    [import_folder, filename_only] = CoreDiskUtilities.GetFullFileParts(filename_or_root_directory);

    % Find out file type. 0=does not exist; 2=file; 7=directory
    exist_result = exist(filename_or_root_directory, 'file');

    if exist_result == 7
        % For a directory, we import every file
        [uids, patient_ids] = ImportDirectoryRecursive(database, import_folder, tags_to_get, reporting);

    elseif exist_result == 2
        % For a single file, we check if it is Dicom.
        single_image_metainfo = MimGetSingleImageInfo(import_folder, filename_only, [], tags_to_get, reporting);
        if single_image_metainfo.ImageFileFormat == MimImageFileFormat.Dicom
            % For a Dicom file we import the whole folder
            [uids, patient_ids] = ImportDirectoryRecursive(database, import_folder, tags_to_get, reporting);
        else
            % For a non-Dicom file we import only the primary file
            database.AddImage(single_image_metainfo);
            uids = {single_image_metainfo.SeriesUid};
            patient_ids = {single_image_metainfo.PatientId};
        end

    else
        reporting.Error('MimImageImporter:FileDoesNotExist', ['The file or directory ' filename_or_root_directory ' does not exist.']);
    end

    reporting.CompleteProgress;
end

function [uids, patient_ids] = ImportDirectoryRecursive(database, import_folder, tags_to_get, reporting)
    uids = [];
    patient_ids = [];
    directories_to_do = CoreStack(import_folder);

    while ~directories_to_do.IsEmpty
        current_dir = directories_to_do.Pop;

        reporting.UpdateProgressMessage(['Importing data from: ' current_dir]);

        if reporting.HasBeenCancelled
            reporting.Error(CoreReporting.CancelErrorId, 'User cancelled');
        end

        [new_uids, new_patient_ids] = ImportFilesInDirectory(database, current_dir, tags_to_get, reporting);
        uids = [uids, new_uids];
        patient_ids = [patient_ids, new_patient_ids];
        next_dirs = CoreDiskUtilities.GetListOfDirectories(current_dir);
        for dir_to_add = next_dirs
            new_dir = fullfile(current_dir, dir_to_add{1});
            directories_to_do.Push(new_dir);
        end
    end
    uids = unique(uids);
end

function [uids, patient_ids] = ImportFilesInDirectory(database, directory, tags_to_get, reporting)
    uids = [];
    patient_ids = [];
    all_filenames = CoreTextUtilities.SortFilenames(CoreDiskUtilities.GetDirectoryFileList(directory, '*'));
    for filename = all_filenames
        if ~strcmp(filename{1}, 'DICOMDIR')
            try
                single_image_metainfo = MimGetSingleImageInfo(directory, filename{1}, [], tags_to_get, reporting);
                if ~isempty(single_image_metainfo.ImageFileFormat)
                    database.AddImage(single_image_metainfo);
                    new_uid = single_image_metainfo.SeriesUid;
                    new_patient_id = single_image_metainfo.PatientId;
                    if ~ismember(new_uid, uids)
                        uids{end + 1} = new_uid;
                    end
                    if ~ismember(new_patient_id, patient_ids)
                        patient_ids{end + 1} = new_patient_id;
                    end
                end
            catch ex
                reporting.ShowWarning('MimImageImporter:ImportFileFailed', ['Failed to import file ' filename{1} ' due to error: ' ex.message]);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% 避免大数吃小数
clear;

% 计算 (10^9+10^-9-10^9)/10^-9
S = (10^9 + 10^-9 - 10^9) / 10 ^ -9;
fprintf('直接计算的结果为: %.20e\n',S);
% S1 = S;

% 改变顺序，再次计算
S = (10^9 - 10^9 + 10^-9) / 10 ^ -9;
fprintf('改变顺序计算的结果为: %.20e\n',S);

repo: NumericalAnalysis


text:
function MapFig = drawLmk(MapFig,Lmk,MapOpt)

% DRAWLMK  Draw 3D landmark.
%   DRAWLMK(MapFig,Lmk) draws 3D landmark Lmk into the map figure by
%   updating the handles in MapFig structure. It accepts different types of
%   landmark by checking the field Lmk.type. Supported landmark types
%   (June 2009) are:
%       'idpPnt' - inverse depth points
%       'eucPnt' - Euclidean points
%       'hmgPnt' - homogeneous points.
%
%   DRAWLMK(...) draws mean and covariances of landmarks:
%     - For punctual landmarks, the mean is a dot and the covariance is the
%       3/sigma ellipsoid.
%     - For line landmarks, the mean is a segment between the two
%       landmark's endpoints, and the covariance is a 3-sigma ellipsoind at
%       each one of these endpoints.
%
%   See also DRAWIDPPNT, DRAWEUCPNT, DRAWHMGPNT.

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

switch (Lmk.type)

    % landmark types
    % --------------
    case {'idpPnt'}
        colors = MapOpt.colors.othPnt;
        drawIdpPnt(MapFig, Lmk, colors, MapOpt);

    case {'eucPnt'}
        colors = MapOpt.colors.defPnt;
        drawEucPnt(MapFig, Lmk, colors, MapOpt);

    case {'hmgPnt'}
        colors = MapOpt.colors.othPnt;
        drawHmgPnt(MapFig, Lmk, colors, MapOpt);

    case {'ahmPnt'}
        colors = MapOpt.colors.othPnt;
        drawAhmPnt(MapFig, Lmk, colors, MapOpt);

    case {'plkLin'}
        colors = MapOpt.colors.othLin;
        drawPlkLin(MapFig, Lmk, colors, MapOpt);

    case {'aplLin'}
        colors = MapOpt.colors.othLin;
        drawAplLin(MapFig, Lmk, colors, MapOpt);

    case 'idpLin'
        colors = MapOpt.colors.othLin;
        drawIdpLin(MapFig, Lmk, colors, MapOpt);

    case {'hmgLin'}
        colors = MapOpt.colors.othLin;
        drawHmgLin(MapFig, Lmk, colors, MapOpt);

    case {'ahmLin'}
        colors = MapOpt.colors.othLin;
        drawAhmLin(MapFig, Lmk, colors, MapOpt);

    otherwise
        error('??? Unknown landmark type ''%s''.',Lmk.type);

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [err_N, sG] = ptv_nuclear_metric(volfix, voldef, pix_resolution, mask, singular_coefs, norm_type)
    volsz = [size(voldef, 1), size(voldef, 2), size(voldef, 3)];
    npix = prod(volsz);
    Nimgs = size(voldef, 5);

%     nuc_G = zeros(size(voldef), 'like', voldef);
    err_N = 0;

%     imagesc([voldef(:,:, round(end/2), 1, 1), volfix(:,:, round(end/2), 1, 1)]); pause(0.02);
%     imagesc([ mean(voldef(:,:, round(end/2), 1, :), 5)]); pause(0.02);

    if isempty(singular_coefs)
        singular_coefs = fl(sqrt(1:Nimgs));
        singular_coefs(1) = 0;
    end

    if isempty(volfix)
        X = reshape(voldef, [npix, Nimgs]);
    else
        X = reshape(cat(5, volfix, voldef), [npix, Nimgs+1]);
        if numel(singular_coefs > 1)
%             singular_coefs = [0, 1, singular_coefs(2:end)];
%             singular_coefs = [fl(singular_coefs(1:end)); singular_coefs(end)];
%             singular_coefs = fl(sqrt(1:(Nimgs+1)));
%             singular_coefs
        end
    end
    if ~isempty(mask)
        mask = logical(mask > 0.5);
        X = X(mask, :);
    end
    singular_coefs = singular_coefs(1 : min(numel(singular_coefs), size(X, 1)));
    [sU, sS, sV] = svdecon(X);
%     sS
    diagS = diag(sS);
    if false
        sGt = sU * ( (diag(singular_coefs)) * sV');
        err_N = sum(diagS(:) .* singular_coefs(:));
    else
        eps = 1e-4;
%         eps=0.1;
        G_sS = diagS ./ sqrt(diagS.^2 + eps);
        sGt = sU * ( (diag(singular_coefs(:) .* G_sS(:))) * sV');
        err_N = sum( sqrt(diagS.^2 + eps) .* singular_coefs(:));
    end

%     if true
%         sk1 = sqrt(1 : size(singular_coefs, 1));
%         sk2 = sk1;
%         sk2(1) = 0;
%
%         sk3 = ones(size(singular_coefs, 1), 1);
%         sk4 = sk3;
%         sk4(1) = 0;
%
%         size(sk1)
%         subplot(221);
%         itmp = reshape(sU * (diag(sk1) * sV'), [volsz, Nimgs]);
%         imagesc(itmp(:,:, 70)); colorbar;
%
%         subplot(222);
%         itmp = reshape(sU * (diag(sk2) * sV'), [volsz, Nimgs]);
%         imagesc(itmp(:,:, 70)); colorbar;
%
%         subplot(223);
%         itmp = reshape(sU * (diag(sk3) * sV'), [volsz, Nimgs]);
%         imagesc(itmp(:,:, 70)); colorbar;
%
%         subplot(224);
%         itmp = reshape(sU * (diag(sk4) * sV'), [volsz, Nimgs]);
%         imagesc(itmp(:,:, 70)); colorbar;
%
%         pause(0.05);
%         pause();
%     end

    if isempty(volfix)
        if ~isempty(mask)
            sG = zeros([npix, Nimgs], 'like', voldef);
            sG(mask, :) = reshape(sGt, [nnz(mask), Nimgs]);
            sG = reshape(sG, [volsz, Nimgs]);
        else
            sG = reshape(sGt, [volsz, Nimgs]);
        end
    else
        if ~isempty(mask)
            sG = zeros([npix, Nimgs+1], 'like', voldef);
            sG(mask, :) = reshape(sGt, [nnz(mask), Nimgs+1]);
            sG = reshape(sG, [volsz, Nimgs+1]);
        else
            sG = reshape(sGt, [volsz, Nimgs+1]);
        end
        sG = sG(:,:,:, 2:end);
    end

end

repo: pTVreg


text:
function [p,ellipse]=phantom3dAniso(varargin)

%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom
%   P = PHANTOM3D(DEF,N) generates a 3D head phantom that can
%   be used to test 3-D reconstruction algorithms.
%
%   DEF is a string that specifies the type of head phantom to generate.
%   Valid values are:
%
%      'Shepp-Logan'            A test image used widely by researchers in
%                               tomography
%      'Modified Shepp-Logan'   (default) A variant of the Shepp-Logan phantom
%                               in which the contrast is improved for better
%                               visual perception.
%      'yu-ye-wang'             Another version of the modified Shepp-Logan
%                               phantom from "Katsevich-Type Algorithms for
%                               Variable Radius Spiral Cone-BeamCT"
%
%   N specifies the 3D grid size of P
%   If N is a scalar, P will have isotropic size [N, N, N]
%   If N is a 3-vector, P will have size [N(1) N(2) N(3)]
%   If you omit the argument, N defaults to [64 64 64].
%
%   P = PHANTOM3D(E,N) generates a user-defined phantom, where each row
%   of the matrix E specifies an ellipsoid in the image.  E has ten columns,
%   with each column containing a different parameter for the ellipsoids:
%
%     Column 1:  A      the additive intensity value of the ellipsoid
%     Column 2:  a      the length of the x semi-axis of the ellipsoid
%     Column 3:  b      the length of the y semi-axis of the ellipsoid
%     Column 4:  c      the length of the z semi-axis of the ellipsoid
%     Column 5:  x0     the x-coordinate of the center of the ellipsoid
%     Column 6:  y0     the y-coordinate of the center of the ellipsoid
%     Column 7:  z0     the z-coordinate of the center of the ellipsoid
%     Column 8:  phi    phi Euler angle (in degrees) (rotation about z-axis)
%     Column 9:  theta  theta Euler angle (in degrees) (rotation about x-axis)
%     Column 10: psi    psi Euler angle (in degrees) (rotation about z-axis)
%
%   For purposes of generating the phantom, the domains for the x-, y-, and
%   z-axes span [-1,1].  Columns 2 through 7 must be specified in terms
%   of this range.
%
%   [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom.
%
%   Class Support
%   -------------
%   All inputs must be of class double.  All outputs are of class double.
%
%   Remarks
%   -------
%   For any given voxel in the output image, the voxel's value is equal to the
%   sum of the additive intensity values of all ellipsoids that the voxel is a
%   part of.  If a voxel is not part of any ellipsoid, its value is 0.
%
%   The additive intensity value A for an ellipsoid can be positive or negative;
%   if it is negative, the ellipsoid will be darker than the surrounding pixels.
%   Note that, depending on the values of A, some voxels may have values outside
%   the range [0,1].
%
%   Example
%   -------
%        ph = phantom3d(128);
%        figure, imshow(squeeze(ph(64,:,:)))
%
%   Copyright 2005 Matthias Christian Schabel (matthias @ stanfordalumni . org)
%   University of Utah Department of Radiology
%   Utah Center for Advanced Imaging Research
%   729 Arapeen Drive
%   Salt Lake City, UT 84108-1218
%
%   This code is released under the Gnu Public License (GPL). For more information,
%   see : http://www.gnu.org/copyleft/gpl.html
%
%   Portions of this code are based on phantom.m, copyrighted by the Mathworks
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Modification May 25, 2015, by Patrick J. Bolan, University of Minnesota
% Added support for anisotropic phantom sizes: the phantom size can now be
% a vector.
%

[ellipse,n] = parse_inputs(varargin{:});

nx = n(1); ny = n(2); nz = n(3);
p = zeros([nx ny nz]);

rngx =  ( (0:nx-1)-(nx-1)/2 ) / ((nx-1)/2);
rngy =  ( (0:ny-1)-(ny-1)/2 ) / ((ny-1)/2);
rngz =  ( (0:nz-1)-(nz-1)/2 ) / ((nz-1)/2);

% PJB: Note the swap of the x and y with meshgrid parameters.
%[x,y,z] = meshgrid(rngx,rngy,rngz);
[x,y,z] = meshgrid(rngy,rngx,rngz);
x=single(x);y=single(y);z=single(z);
coord = [flatten(single(x)); flatten(single(y)); flatten(single(z))];
clear x y z;
p = flatten(p);

for k = 1:size(ellipse,1)
   A = ellipse(k,1);            % Amplitude change for this ellipsoid
   asq = ellipse(k,2)^2;        % a^2
   bsq = ellipse(k,3)^2;        % b^2
   csq = ellipse(k,4)^2;        % c^2
   x0 = ellipse(k,5);           % x offset
   y0 = ellipse(k,6);           % y offset
   z0 = ellipse(k,7);           % z offset
   phi = ellipse(k,8)*pi/180;   % first Euler angle in radians
   theta = ellipse(k,9)*pi/180; % second Euler angle in radians
   psi = ellipse(k,10)*pi/180;  % third Euler angle in radians

   cphi = cos(phi);
   sphi = sin(phi);
   ctheta = cos(theta);
   stheta = sin(theta);
   cpsi = cos(psi);
   spsi = sin(psi);

   % Euler rotation matrix
   alpha = [cpsi*cphi-ctheta*sphi*spsi   cpsi*sphi+ctheta*cphi*spsi  spsi*stheta;
            -spsi*cphi-ctheta*sphi*cpsi  -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta;
            stheta*sphi                  -stheta*cphi                ctheta];

   % rotated ellipsoid coordinates
   coordp = alpha*coord;

   idx = find((coordp(1,:)-x0).^2./asq + (coordp(2,:)-y0).^2./bsq + (coordp(3,:)-z0).^2./csq <= 1);
   p(idx) = p(idx) + A;
end

%p = reshape(p,[nx ny nz]);
p = reshape(p, [nx ny nz]);

return;

function out = flatten(in)

out = reshape(in,[1 numel(in)]);

return;

function [e,n] = parse_inputs(varargin)
%  e is the m-by-10 array which defines ellipsoids
%  n is a 3-vector with the size of the phantom brain image, [nx ny nz]

n = [64 64 64];     % The default size
e = [];
defaults = {'shepp-logan', 'modified shepp-logan', 'yu-ye-wang'};

for i=1:nargin
   if ischar(varargin{i})         % Look for a default phantom
      def = varargin{i};
      idx = strcmpi(def, defaults);
      if isempty(idx)
         eid = sprintf('Images:%s:unknownPhantom',mfilename);
         msg = 'Unknown default phantom selected.';
         error(eid,'%s',msg);
      end
      switch defaults{idx}
      case 'shepp-logan'
         e = shepp_logan;
      case 'modified shepp-logan'
         e = modified_shepp_logan;
      case 'yu-ye-wang'
         e = yu_ye_wang;
      end
   elseif numel(varargin{i})==1
      n = [varargin{i} varargin{i} varargin{i}];   % a scalar is the image size
   elseif numel(varargin{i})==3
      siz = varargin{i};
      n = [siz(1) siz(2) siz(3)]; % 3 integers specify image dimensions
   elseif ndims(varargin{i})==2 && size(varargin{i},2)==10
      e = varargin{i};            % user specified phantom
   else
      eid = sprintf('Images:%s:invalidInputArgs',mfilename);
      msg = 'Invalid input arguments.';
      error(eid,'%s',msg);
   end
end

% ellipse is not yet defined
if isempty(e)
   e = modified_shepp_logan;
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Default head phantoms:   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function e = shepp_logan

e = modified_shepp_logan;
e(:,1) = [1 -.98 -.02 -.02 .01 .01 .01 .01 .01 .01];

return;

function e = modified_shepp_logan
%
%   This head phantom is the same as the Shepp-Logan except
%   the intensities are changed to yield higher contrast in
%   the image.  Taken from Toft, 199-200.
%
%         A      a     b     c     x0      y0      z0    phi  theta    psi
%        -----------------------------------------------------------------
e =    [  1  .6900  .920  .810      0       0       0      0      0      0
        -.8  .6624  .874  .780      0  -.0184       0      0      0      0
        -.2  .1100  .310  .220    .22       0       0    -18      0     10
        -.2  .1600  .410  .280   -.22       0       0     18      0     10
         .1  .2100  .250  .410      0     .35    -.15      0      0      0
         .1  .0460  .046  .050      0      .1     .25      0      0      0
         .1  .0460  .046  .050      0     -.1     .25      0      0      0
         .1  .0460  .023  .050   -.08   -.605       0      0      0      0
         .1  .0230  .023  .020      0   -.606       0      0      0      0
         .1  .0230  .046  .020    .06   -.605       0      0      0      0 ];

return;

function e = yu_ye_wang
%
%   Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius Spiral Cone-Beam CT
%
%         A      a     b     c     x0      y0      z0    phi  theta    psi
%        -----------------------------------------------------------------
e =    [  1  .6900  .920  .900      0       0       0      0      0      0
        -.8  .6624  .874  .880      0       0       0      0      0      0
        -.2  .4100  .160  .210   -.22       0    -.25    108      0      0
        -.2  .3100  .110  .220    .22       0    -.25     72      0      0
         .2  .2100  .250  .500      0     .35    -.25      0      0      0
         .2  .0460  .046  .046      0      .1    -.25      0      0      0
         .1  .0460  .023  .020   -.08    -.65    -.25      0      0      0
         .1  .0460  .023  .020    .06    -.65    -.25     90      0      0
         .2  .0560  .040  .100    .06   -.105    .625     90      0      0
        -.2  .0560  .056  .100      0    .100    .625      0      0      0 ];

return;

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function L= DoPyrDec(x,level)

%% use the raised-cosine function to get a smooth transition band.
smooth_func = @rcos;
Pyr_mode=1.5;
L = PyrNDDec_mm(x, 'S', level, Pyr_mode, smooth_func);

author: thomaskuestner


text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returns the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

for j = 1:K
    centroids(j, :) = mean(X(idx==j, :), 1);
end

% =============================================================

end

author: imLogM


text:
function mse = MeanSquareError(x,y)
% This function computes the mean square error
% between x and y

mse = sum( sum( (y - x).^2 ) );
mse = mse / prod(size(x));


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

%*****************************************************************************80
%
%% FEJER2_COMPUTE computes a Fejer type 2 quadrature rule.
%
%  Discussion:
%
%    This method uses a direct approach.  The paper by Waldvogel
%    exhibits a more efficient approach using Fourier transforms.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%    Walter Gautschi,
%    Numerical Quadrature in the Presence of a Singularity,
%    SIAM Journal on Numerical Analysis,
%    Volume 4, Number 3, 1967, pages 357-362.
%
%    Joerg Waldvogel,
%    Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules,
%    BIT Numerical Mathematics,
%    Volume 43, Number 1, 2003, pages 1-18.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  if ( n == 1 )
    x(1) = 0.0;
    w(1) = 2.0;
    return
  elseif ( n == 2 )
    x(1) = -0.5;
    x(2) =  0.5;
    w(1:2) = 1.0;
    return
  end

  theta(1:n) = ( n : -1 : 1 ) * pi / ( n + 1 );
  x(1:n) = cos ( theta(1:n) );

  for i = 1 : n

    w(i) = 1;

    for j = 1 : floor ( ( n - 1 ) / 2 )
      w(i) = w(i) - 2 * cos ( 2 * j * theta(i) ) / ( 4 * j * j - 1 );
    end

    if ( 2 < n )
      p = 2 * floor ( ( n + 1 ) / 2 ) - 1;
      w(i) = w(i) - cos ( ( p + 1 ) * theta(i) ) / p;
    end

  end

  w(1:n) = 2 * w(1:n) / ( n + 1 );

  return
end


text:
function alpha = dfaScalingExponent(x, minBoxSize, maxBoxSize, pflag)
%
% varargout = dfaScalingExponent(xminBoxSize, midBoxSize, maxBoxSize, pflag)
% calculates the detrended fluctuation analysis estimate of the scaling
% exponent alpha.
%
% INPUTS
%         x          : A Nx1 vector containing the series to be analyzed
%         minBoxSize : Smallest box width (default: 4)
%         maxBoxSize : Largest box width (default: N/4)
%         pflag      : (Optional) pflag=1 plot,  pflag=0
% OUTPUTS
%         alpha      : estimate of scaling exponent, +
%                      minBoxSize <= n <= maxBoxSize
%
% The raw time series x(i) is first integrated to give y(i); i=1,...,N.
% For each length scale, n, y(i) is divided into segments of equal length, n.
% In each segment, the data is detrended by subtracting the local linear least
% squares fit, yn(k).  The root-mean-square fluctuation of this integrated
% and detrended time series is given by
% F(n) = sqrt( (1/N) sum_{k=1}^N [y(k) - yn(k)]^2 )
% We calculate the average fluctuation F(n) for each segment n.
% If the scaling approximately given by F(n) = c n^alpha,
% we can estimate alpha by calculating the slope of log F(n) versus log n.
% Such a linear relationship on a log-log plot indicates the presence of
% power law (fractal) scaling.
% A log-log plot of F(n) against n is provided when pflag=1.  Default: plag=0.
% Peng C-K, Buldyrev SV, Havlin S, Simons M, Stanley HE, Goldberger AL.
% Mosaic organization of DNA nucleotides. Phys Rev E 1994;49:1685-1689.
%
%
% 09-20-2017 Modified by Giulia Da Poian (GDP) to be included in the Physionet
%            HRV Toolkit for Matlab. (Original function name: dfa)
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
% Copyright (c) 2005 Patrick E. McSharry (patrick@mcsharry.net)
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

if nargin < 2 || isempty(minBoxSize)
    minBoxSize = 4;
end
if nargin < 3 || isempty(maxBoxSize)
    maxBoxSize = length(x)/4;
end
if nargin < 4
   pflag = 0;
end

if size(x,1)<size(x,2)
    x=x';
end

N = length(x);
y = cumsum(x);

n1 = round(log2(minBoxSize)); % modified GDP, was 3
n2 = round(log2(maxBoxSize)); % modified GDP, was n2 = round(log2(N/2))
ns = (2.^(n1:n2))';           % modified GDP, was ns =[2.^[n1:n2] N]'

nn = length(ns);
F = zeros(nn,1);
for n=1:nn
   t = trend(y, ns(n));
   z = y - t;
   F(n) = sqrt(mean(z.^2));

end

lns = log10(ns);
lF = log10(F);
A = ones(nn,2);
A(:,2) = lns;
a = pinv(A)*lF;
alpha = a(2);
lFpred = A*a;

if pflag == 1
    figure;
    loglog(10.^lns, 10.^lF,'b.-','MarkerSize',16);
    hold on;
    loglog(10.^[lns(1) lns(nn)], 10.^[lFpred(1) lFpred(nn)],'k');
    xlabel('n');
    ylabel('F(n)');
    title(['F(n) ~ n^{\alpha} with \alpha = ' num2str(a(2)) ]);
end

end % dfaScalingExponent function

function t = trend(y, n)
    N = length(y);
    t = zeros(N,1);
    r = floor(N/n);
    for i=1:r
       v = y((i-1)*n+1:i*n);
       t((i-1)*n+1:i*n) = linfit(v);
    end
    v = y(r*n+1:N);
    t(r*n+1:N) = linfit(v);
end % trend function

function up = linfit(v)
    k = length(v);
    A = ones(k,2);
    u = [1:k]';
    A(:,2) = u;
    a = pinv(A)*v;
    up = A*a;
end % linfit function

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [modelLOD] = setQualitativeConstraints(model,cond_uptake,cond_uptake_LODs, cond_secretion,cond_secretion_LODs,cellConc,t,cellWeight,ambiguous_metabolites,basisMedium)
% This function sets qualitative constraints (by enforcing minimal uptake or secretion based on individual detection limits), e.g., based on the uptake and secretion
% profile of metabolites measured through mass-spectrometry. Uptake is only possible if the lower bound has been set to a value >0 using `setMediumConstraints`
% The minimal allowable uptake and secretion flux is defined by enforcing uptake at or above the limit of detection (mass-spectrometry). If these values are not available, a very small value, e.g., 1.0E-06 can be used. Note that this
% value has to be below the concentrations defined in the medium. Otherwise the model will be infeasible. Alternatively to the LODs a small value can
% be used (e.g., 0.00001). Note that the value has to be above the threshold later called zero.
%
% USAGE:
%
%    [modelLOD] = setQualitativeConstraints(model, cond_uptake, cond_uptake_LODs, cond_secretion, cond_secretion_LODs, cellConc, t, cellWeight, ambiguous_metabolites, basisMedium)
%
% INPUTS:
%    model:                   Metabolic model (Recon), with set constraints (output of `setMediumConstraints`)
%    cond_uptake:             Vector of exchanges of metabolites consumed by the cells in the experiment
%    cond_uptake_LODs:        Vector of detection limits (LOD in Mm) for the compounds and in the experiment
%    cond_secretion:          Vector of metabolite exchanges consumed by the cells in the experiment
%    cond_secretion_LODs:     Vector of detection limits (LOD in Mm) for the compounds and in the experiment
%    cellConc:                Cell concentration (cells per 1 ml)
%    t:                       Time in hours
%    cellWeight:              gDW per cell
%    ambiguous_metabolites:   Since all exchanges, except the ones specified in `uptake` and `secretion` are closed in `modelLOD`, this input
%                             variable allows to specify metabolite exchanges that should remain open.
%                             Thus, if these exchanges are open in the starting model they can be uptaken
%                             or secreted by the `modelLOD`, e.g., `ambiguous_metabolites`
%                             This can for example be the case if it is suspected that an uptake or secretion might have taken place in concentrations below the
%                             detection limit.
%    basisMedium:             Vector defining the metabolite exchanges of the basic medium, i.e., `ions`, `mediumCompounds`
%
% OUTPUT:
%    modelLOD:                Model that is constrained qualitatively to the condition-specific uptake and secretion profile
%
% .. Author: - Maike K. Aurich 27/05/15

A = strfind(model.rxns, 'EX_');
%Find all exchange reactions
for i = 1:length(A);
    B= A{i};
       if  isempty(B);
       A{i}=[0];
       end
end
cc = cell2mat(A);
idx_EXreactions= find(cc);
ex_Rxns = model.rxns(idx_EXreactions); %make exchange reaction list

%% Calculate flux rates (LODs)
for i = 1 : length(cond_uptake_LODs)
    [ub_up_flux(i,1)] = conc2Rate(cond_uptake_LODs(i),cellConc, t, cellWeight);
end

for i = 1 : length(cond_secretion_LODs) % conversion from conc.(minimal detection value total, in Mm) to flux value per cell mmol/gDW/hr
        [lb_secr_flux(i,1)] = conc2Rate(cond_secretion_LODs(i),cellConc, t, cellWeight);
end

%% Define which exchanges to close and close them
conditional_exchanges = [cond_uptake;cond_secretion;basisMedium;ambiguous_metabolites];

no_exchange = ex_Rxns(find(~ismember(ex_Rxns,conditional_exchanges)));%find reactions not associated with detected exhange metabolites, basic medium metabolites, low quantity metabolites

model = changeRxnBounds(model,no_exchange,0,'b'); %block all other Exchange reactions

%% Apply LOD constraints to enforce uptake and secretion

for k=1:length(cond_uptake)
    ub=ub_up_flux(k);
    rxns = cond_uptake(k);
model = changeRxnBounds(model,rxns,-ub,'u'); %enforce uptake of metabolites taken up by cells in the experiment
end
clear k
for k=1:length(cond_secretion)
    lb=lb_secr_flux(k);
model = changeRxnBounds(model,cond_secretion(k),lb,'l');% enforce secretion of metabolites secreted by the cells in the experiment
end

modelLOD = model;
end

author: opencobra

repo: cobratoolbox


text:
% Copyright (C) Daphne Koller, Stanford Univerity, 2012
%
% Author: Huayan Wang

function VisualizeDataset(Dataset)

f = figure;
for i=1:size(Dataset,1)
    img = ShowPose(reshape(Dataset(i,:,:), [10 3]));
    imshow(img);
    pause(0.3)
    if (~ishandle(f)) break; end; % quit loop when user closes the figure
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [assignment,cost] = munkres(costMat)
% MUNKRES   Munkres Assign Algorithm
%
% [ASSIGN,COST] = munkres(COSTMAT) returns the optimal assignment in ASSIGN
% with the minimum COST based on the assignment problem represented by the
% COSTMAT, where the (i,j)th element represents the cost to assign the jth
% job to the ith worker.
%

% This is vectorized implementation of the algorithm. It is the fastest
% among all Matlab implementations of the algorithm.

% Examples
% Example 1: a 5 x 5 example
%{
[assignment,cost] = munkres(magic(5));
[assignedrows,dum]=find(assignment);
disp(assignedrows'); % 3 2 1 5 4
disp(cost); %15
%}
% Example 2: 400 x 400 random data
%{
n=400;
A=rand(n);
tic
[a,b]=munkres(A);
toc                 % about 6 seconds
%}

% Reference:
% "Munkres' Assignment Algorithm, Modified for Rectangular Matrices",
% http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html

% version 1.0 by Yi Cao at Cranfield University on 17th June 2008

assignment = false(size(costMat));
cost = 0;

costMat(costMat~=costMat)=Inf;
validMat = costMat<Inf;
validCol = any(validMat);
validRow = any(validMat,2);

nRows = sum(validRow);
nCols = sum(validCol);
n = max(nRows,nCols);
if ~n
    return
end

dMat = zeros(n);
dMat(1:nRows,1:nCols) = costMat(validRow,validCol);

%*************************************************
% Munkres' Assignment Algorithm starts here
%*************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   STEP 1: Subtract the row minimum from each row.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 dMat = bsxfun(@minus, dMat, min(dMat,[],2));

%**************************************************************************
%   STEP 2: Find a zero of dMat. If there are no starred zeros in its
%           column or row start the zero. Repeat for each zero
%**************************************************************************
zP = ~dMat;
starZ = false(n);
while any(zP(:))
    [r,c]=find(zP,1);
    starZ(r,c)=true;
    zP(r,:)=false;
    zP(:,c)=false;
end

while 1
%**************************************************************************
%   STEP 3: Cover each column with a starred zero. If all the columns are
%           covered then the matching is maximum
%**************************************************************************
    primeZ = false(n);
    coverColumn = any(starZ);
    if ~any(~coverColumn)
        break
    end
    coverRow = false(n,1);
    while 1
        %**************************************************************************
        %   STEP 4: Find a noncovered zero and prime it.  If there is no starred
        %           zero in the row containing this primed zero, Go to Step 5.
        %           Otherwise, cover this row and uncover the column containing
        %           the starred zero. Continue in this manner until there are no
        %           uncovered zeros left. Save the smallest uncovered value and
        %           Go to Step 6.
        %**************************************************************************
        zP(:) = false;
        zP(~coverRow,~coverColumn) = ~dMat(~coverRow,~coverColumn);
        Step = 6;
        while any(any(zP(~coverRow,~coverColumn)))
            [uZr,uZc] = find(zP,1);
            primeZ(uZr,uZc) = true;
            stz = starZ(uZr,:);
            if ~any(stz)
                Step = 5;
                break;
            end
            coverRow(uZr) = true;
            coverColumn(stz) = false;
            zP(uZr,:) = false;
            zP(~coverRow,stz) = ~dMat(~coverRow,stz);
        end
        if Step == 6
            % *************************************************************************
            % STEP 6: Add the minimum uncovered value to every element of each covered
            %         row, and subtract it from every element of each uncovered column.
            %         Return to Step 4 without altering any stars, primes, or covered lines.
            %**************************************************************************
            M=dMat(~coverRow,~coverColumn);
            minval=min(min(M));
            if minval==inf
                return
            end
            dMat(coverRow,coverColumn)=dMat(coverRow,coverColumn)+minval;
            dMat(~coverRow,~coverColumn)=M-minval;
        else
            break
        end
    end
    %**************************************************************************
    % STEP 5:
    %  Construct a series of alternating primed and starred zeros as
    %  follows:
    %  Let Z0 represent the uncovered primed zero found in Step 4.
    %  Let Z1 denote the starred zero in the column of Z0 (if any).
    %  Let Z2 denote the primed zero in the row of Z1 (there will always
    %  be one).  Continue until the series terminates at a primed zero
    %  that has no starred zero in its column.  Unstar each starred
    %  zero of the series, star each primed zero of the series, erase
    %  all primes and uncover every line in the matrix.  Return to Step 3.
    %**************************************************************************
    rowZ1 = starZ(:,uZc);
    starZ(uZr,uZc)=true;
    while any(rowZ1)
        starZ(rowZ1,uZc)=false;
        uZc = primeZ(rowZ1,:);
        uZr = rowZ1;
        rowZ1 = starZ(:,uZc);
        starZ(uZr,uZc)=true;
    end
end

% Cost of assignment
assignment(validRow,validCol) = starZ(1:nRows,1:nCols);
cost = sum(costMat(assignment));

repo: BUPT-ICS-Courseware


text:
function [nS, iNonMatch] = merge_spikes0(uBase, nS, uS, crit)

if ~isempty(uBase)
    cdot = uBase * uS';

    baseNorms = sum(uBase.^2, 2)';
    newNorms  = sum(uS.^2, 2)';

    cNorms = 1e-10 + repmat(baseNorms', 1, numel(newNorms)) + repmat(newNorms, numel(baseNorms), 1);

    cdot = 1 - 2*cdot./cNorms;

    [cdotmin, imin] = min(cdot, [], 1);

    iMatch = cdotmin<crit;

    nSnew = hist(imin(iMatch), 1:1:size(uBase,1));
    nS = nS + nSnew';

    iNonMatch = find(cdotmin>crit);
else
   iNonMatch = 1:size(uS,2);
   nS = [];
end

author: cortex-lab

repo: KiloSort

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


text:
% MacromoleculeUtility test
%
% Author: Derek Macklin, macklin@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated 6/29/2011
classdef MacromoleculeUtil_Test < TestCase
    properties
        simulation
    end

    methods
        function this = MacromoleculeUtil_Test(name)
            this = this@TestCase(name);
        end

        function setUp(this)
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load();
            this.simulation = sim;
        end

        function testGetMacromoleculeIndexsIDsNames(this)
            import edu.stanford.covert.cell.sim.util.MacromoleculeUtil;

            s = this.simulation;

            % Gene that doesn't exist
            [idx, compIdx, id, name, type, ...
                geneIdxs, geneCompIdxs, geneIds, geneNames, ...
                nascentRnaIdxs, nascentRnaCompIdxs, nascentRnaIds, nascentRnaNames, ...
                matureRnaIdxs, matureRnaCompIdxs, matureRnaIds, matureRnaNames, ...
                monomerIdxs, monomerCompIdxs, monomerIds, monomerNames, ...
                complexIdxs, complexCompIdxs, complexIds, complexNames] = ...
                MacromoleculeUtil.getMacromoleculeIndexsIDsNames('DOES_NOT_EXIST', s);

            assert(isequal(idx,[]));
            assert(isequal(compIdx,[]));
            assert(isequal(id,[]));
            assert(isequal(name,[]));
            assert(isequal(type,[]));
            assert(isequal(geneIdxs,[]));
            assert(isequal(geneCompIdxs,[]));
            assert(isequal(geneIds,[]));
            assert(isequal(geneNames,[]));
            assert(isequal(nascentRnaIdxs,[]));
            assert(isequal(nascentRnaCompIdxs,[]));
            assert(isequal(nascentRnaIds,[]));
            assert(isequal(nascentRnaNames,[]));
            assert(isequal(matureRnaIdxs,[]));
            assert(isequal(matureRnaCompIdxs,[]));
            assert(isequal(matureRnaIds,[]));
            assert(isequal(matureRnaNames,[]));
            assert(isequal(monomerIdxs,[]));
            assert(isequal(monomerCompIdxs,[]));
            assert(isequal(monomerIds,[]));
            assert(isequal(monomerNames,[]));
            assert(isequal(complexIdxs,[]));
            assert(isequal(complexCompIdxs,[]));
            assert(isequal(complexIds,[]));
            assert(isequal(complexNames,[]));

            % Two species present in a protein complex
            [idx, compIdx, id, name, type, ...
                geneIdxs, geneCompIdxs, geneIds, geneNames, ...
                nascentRnaIdxs, nascentRnaCompIdxs, nascentRnaIds, nascentRnaNames, ...
                matureRnaIdxs, matureRnaCompIdxs, matureRnaIds, matureRnaNames, ...
                monomerIdxs, monomerCompIdxs, monomerIds, monomerNames, ...
                complexIdxs, complexCompIdxs, complexIds, complexNames] = ...
                MacromoleculeUtil.getMacromoleculeIndexsIDsNames('MG_008_379_TETRAMER',s);

            assert(isequal(id,'MG_008_379_TETRAMER'));

            assert(isequal(name,'tRNA uridine 5-carboxymethylaminomethyl modification enzyme'));

            assert(isequal(type,'complex'));

            assert(isequal(size(geneIdxs),[2 1]));

            assert(isequal(size(geneCompIdxs),[2 1]));

            assert(isequal(size(geneIds),[2 1]));
            assert(isequal(geneIds,{'MG_008';'MG_379'}));

            assert(isequal(size(geneNames),[2 1]));
            assert(isequal(geneNames,{'tRNA modification GTPase TrmE';'tRNA uridine 5-carboxymethylaminomethyl modification enzyme GidA'}));

            assert(isequal(size(nascentRnaIdxs),[1 2]));

            assert(isequal(size(nascentRnaCompIdxs),[1 2]));

            assert(isequal(size(nascentRnaIds),[2 1]));
            assert(isequal(nascentRnaIds,{'TU_003';'TU_275'}));

            assert(isequal(size(nascentRnaNames),[2 1]));
            assert(isequal(nascentRnaNames,{'MG_005 - MG_009';'MG_379 - MG_381'}));

            assert(isequal(size(matureRnaIdxs),[1 2]));

            assert(isequal(size(matureRnaCompIdxs),[1 2]));

            assert(isequal(size(matureRnaIds),[2 1]));
            assert(isequal(matureRnaIds,{'TU_003';'TU_275'}));

            assert(isequal(size(matureRnaNames),[2 1]));
            assert(isequal(matureRnaNames,{'MG_005 - MG_009';'MG_379 - MG_381'}));

            assert(isequal(size(monomerIdxs),[2 1]));
            assert(isequal(size(monomerCompIdxs),[2 1]));

            assert(isequal(size(monomerIds),[2 1]));
            assert(isequal(monomerIds,{'MG_008_MONOMER';'MG_379_MONOMER'}));

            assert(isequal(size(monomerNames),[2 1]));
            assert(isequal(monomerNames,{'tRNA modification GTPase TrmE';'tRNA uridine 5-carboxymethylaminomethyl modification enzyme GidA'}));

            assert(isequal(complexIds,{'MG_008_379_TETRAMER'}));
            assert(isequal(complexNames,{'tRNA uridine 5-carboxymethylaminomethyl modification enzyme'}));
        end

        function testGetMacromoleculeCounts(~)
            import edu.stanford.covert.cell.sim.util.MacromoleculeUtil;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            %select a simulation to analyze
            [simDir, ~, sim] = SimulationDiskUtil.getLatestSimulation();

            %load data
            [cnts, ...
                geneCnts, ...
                nascentRnaCnts, processedRnaCnts, matureRnaCnts, boundRnaCnts, misfoldedRnaCnts, damagedRnaCnts, aminoacylatedRnaCnts, ...
                nascentMonCnts, processedIMonCnts, translocatedMonCnts, processedIIMonCnts, foldedMonCnts, matureMonCnts, inactivatedMonCnts, boundMonCnts, misfoldedMonCnts, damagedMonCnts, ...
                nascentCpxCnts, matureCpxCnts, inactivatedCpxCnts, boundCpxCnts, misfoldedCpxCnts, damagedCpxCnts, ...
                rnaPolTUs, rnaPolPosStrnds, ribMRNAs, ribPos] = ...
                MacromoleculeUtil.getMacromoleculeCounts('MG_357', sim, simDir);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
classdef RWMOP20 < PROBLEM
% <multi> <real> <constrained>
% Hydro-static thrust bearing design problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 4;
            obj.lower    = [ 1, 1,  1e-6,1];
            obj.upper    = [16, 16, 16*1e-6,16];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x = varargin{1};
            R = x(:,1); Ro = x(:,2);  mu = x(:,3); Q = x(:,4);
            gamma = 0.0307; C = 0.5; n = -3.55; C1 = 10.04;
            Ws = 101000; Pmax = 1000; delTmax = 50; hmin = 0.001;
            gg = 386.4; N = 750;
            P    = (log10(log10(8.122*1e6.*mu+0.8))-C1)./n;
            delT = 2.*(10.^P-560);
            Ef   = 9336.*Q.*gamma.*C.*delT;
            h    = (2.*pi.*N./60).^2.*2.*pi.*mu./Ef.*(R.^4./4-Ro.^4./4)-1e-5;
            Po   = (6.*mu.*Q./(pi.*h.^3)).*log(R./Ro);
            W    = pi.*Po./2.*(R.^2-Ro.^2)./(log(R./Ro)-1e-5);
            % Objective function
            f(:,1) = (Q.*Po./0.7+Ef)./12;
            f(:,2) = gamma./(gg.*Po).*(Q./(2.*pi.*R.*h));
            % Constraints
            g(:,1) = Ws-W;
            g(:,2) = Po-Pmax;
            g(:,3) = delT-delTmax;
            g(:,4) = hmin-h;
            g(:,5) = Ro-R;
            g(:,6) = f(:,2)-0.001;
            g(:,7) = W./(pi.*(R.^2-Ro.^2)+1e-5)-5000;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [2.6725846e+02  -2.7672651e-05];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [f,nxe_train,w] = train_linear_calibration(tar,non,prior,obj_func,niters,quiet)
% A function that uses 'train_binary_classifier' to train a linear
% fusion function.
% Inputs:
%   tar: a vector of scores for target trials
%   non: a vector of scores for non-target trials
%   prior: the effective target prior
%   obj_func: the objective function for the training algorithm.  If []
%     then cllr objective is used.
%   niters: The maximum number of training iterations
%   quiet: A boolean indicating whether the training algorithm
%     should print information on its progress.
% Outputs:
%   f: A function handle to the trained calibration function.
%     f(scores) will return a vector of the same length as scores
%       with each score scaled and shifted according to the values
%       learned during training.
%   nxe_train: normalized multiclass cross-entropy of the solution.
%       The range is 0 (good) to 1 (useless).
%   w: The calibration weigths (scale factor and offset).

% check the inputs
assert(nargin==6)
assert(size(tar,1)==1)
assert(size(non,1)==1)
assert(length(tar)>0)
assert(length(non)>0)

% create function handle for function that must be trained
[fusion,params] = linear_fuser([],[tar,non]);
% get a starting point for the calibration weights: 'w0'
w0 = params.get_w0();
% let the trainer know which scores are target scores and which are
% non-target scores
classf = [ones(1,length(tar)),-ones(1,length(non))];
% do the training to get the calibration weights 'w'
[w,nxe_train] = train_binary_classifier(fusion,classf,w0,obj_func,prior,[],0,niters,[],[],quiet);

% create a function handle that will calibrate input scores using
% the trained weights 'w'
f = @(scores) linear_fuser(w,scores);

author: nesl


text:
function gp_save(file,data,header)
%GP_SAVE saves x,y as a text file in a Gnuplot compatible format
%
%   Usage: gp_save(file,data,[header])
%
%   Input parameters:
%       file    - filename of the data file
%       data    - data matrix
%       header  - header comment added before the data
%
%   See also: gp_save_matrix

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 3;
narginchk(nargmin,nargmax);
if nargin==nargmax-1
    header = '';
end
isargchar(file);
isargmatrix(data);
isargchar(header);

%% ===== Main ============================================================
% Write header to the file
fid = fopen(file,'w');
fprintf(fid,'# Data file generated by gp_save.m\n');
fprintf(fid,['# ',header,'\n']);
fclose(fid);
% Append the data to the file using tabulator as a delimiter between the data
if isoctave
    dlmwrite(file,data,'\t','-append');
else
    dlmwrite(file,data,'delimiter','\t','-append');
end

author: sfstoolbox

repo: sfs-matlab


text:
function Prepare_data_Multi_PIE_all()

    % This bit collects all of the multi-pie labels into a single structure for
    % easy access
    labels_root = 'C:\Users\tbaltrus\Dropbox\AAM\test data\MultiPI_AAM/';

    % The location of the Multi-PIE data folder
    multi_pie_root = 'D:\MultiPIE/Image_Data/';

    multi_pie_labels = CollectMultiPieLabels(labels_root, multi_pie_root);

    %%

    % Make sure same images are used across scales
    rng(0);
    ExtractTrainingMultiPIE(0.25, multi_pie_labels);

    rng(0);
    ExtractTrainingMultiPIE(0.35, multi_pie_labels);

    rng(0);
    ExtractTrainingMultiPIE(0.5, multi_pie_labels);

    rng(0);
    ExtractTrainingMultiPIE(1.0, multi_pie_labels);

end
% Now extract the relevant information

function [multi_pie_labels] = CollectMultiPieLabels(labels_root, multi_pie_root)

    multi_pie_labels = struct;
    currentLabel = 0;
    left_to_frontal_map = [1, 28;
                        2, 29;
                        3, 30;
                        4, 31;
                        5, 34;
                        6, 36;
                        7, 27;
                        8, 26;
                        9, 25;
                        10,24;
                        11,46;
                        12,45;
                        13,44;
                        14,48;
                        15,47;
                        16,52;
                        17,53;
                        18,54;
                        19,55;
                        20,56;
                        21,57;
                        22,58;
                        23,63;
                        24,64;
                        26,66;
                        27,67;
                        30,9;
                        31,10;
                        32,11;
                        33,12;
                        34,13;
                        35,14;
                        36,15;
                        37,16;
                        38,17;];

    right_to_frontal_map = [1, 28;
                            2, 29;
                            3, 30;
                            4, 31;
                            5, 34;
                            6, 32;
                            7 ,18;
                            8 ,19;
                            9 ,20;
                            10,21;
                            11,37;
                            12,38;
                            13,39;
                            14,41;
                            15,42;
                            16,52;
                            17,51;
                            18,50;
                            19,49;
                            20,60;
                            21,59;
                            22,58;
                            23,63;
                            24,62;
                            26,68;
                            27,67;
                            30, 9;
                            31, 8;
                            32, 7;
                            33, 6;
                            34, 5;
                            35, 4;
                            36, 3;
                            37, 2;
                            38, 1];

    for i=1:4

        labelsDir = sprintf('%sSession%d/', labels_root, i);

        labels = dir([labelsDir '/*.mat']);

        double_label = false;

        for j = 1:numel(labels)
            if(double_label)
                double_label = false;
                continue;
            end

            load([labelsDir labels(j).name]);

            if(size(pts,1) ~= 68 && size(pts,1) ~= 39)
                continue;
            end

            userID = labels(j).name(1:3);
            recID = labels(j).name(8:9);
            camID = labels(j).name(11:12);
            camID2 = labels(j).name(13);
            viewID = labels(j).name(15:16);

            % doubling labels seem to be odd/incorrect
            if(j < numel(labels))

                userIDN = labels(j+1).name(1:3);
                recIDN = labels(j+1).name(8:9);
                camIDN = labels(j+1).name(11:12);
                camID2N = labels(j+1).name(13);

                if(strcmp(userIDN, userID) && strcmp(recID, recIDN) && strcmp(camIDN, camID) && strcmp(camID2N, camID2))
                    double_label = true;
                    continue;
                end

            end

            % camera id 08 is from a very different perspective
            if(strcmp(camID, '08') && strcmp(camID2, '1'))
                continue;
            end

            currentLabel = currentLabel + 1;

            multi_pie_labels(currentLabel).label_file = [labelsDir labels(j).name];
            multi_pie_labels(currentLabel).user_ID = userID;
            multi_pie_labels(currentLabel).rec_ID = recID;
            multi_pie_labels(currentLabel).cam_ID = camID;
            multi_pie_labels(currentLabel).cam_ID2 = camID2;

            multiPieImageDir = sprintf('%s/session0%d/multiview/%s/%s/%s_%s/', multi_pie_root, i,userID,recID,camID,camID2);
            multiPieImgs = dir([multiPieImageDir '*.png']);

            multi_pie_labels(currentLabel).img_dir = multiPieImageDir;
    %         multi_pie_labels(currentLabel).img_loc = [labels(j+1).name(1:end-7), '.png'];
            actualLabel = dir([multiPieImageDir '/*' viewID '.png']);
            multi_pie_labels(currentLabel).img_locs = {multiPieImgs.name};
            multi_pie_labels(currentLabel).actual_img = actualLabel.name;
            landmark_labels = zeros(68,2);

            % This is a profile image
            if(size(pts,1) == 39)
                % Determine if left or right

                if(pts(4,1) < pts(39,1))
                    multi_pie_labels(currentLabel).type = 'profile_left';
                    landmark_labels(left_to_frontal_map(:,2),:) = pts(left_to_frontal_map(:,1),:);
                else
                    multi_pie_labels(currentLabel).type = 'profile_right';
                    landmark_labels(right_to_frontal_map(:,2),:) = pts(right_to_frontal_map(:,1),:);
                end
            else
                landmark_labels = pts;
            end

            multi_pie_labels(currentLabel).landmark_labels = landmark_labels;

        end
    end
end

function ExtractTrainingMultiPIE( training_scale, multi_pie_labels)
%PREPARETRAININGIMAGEMPIE This function collects a subset of Multi-PIE
%images at different expressions and lighting conditions to store it in a
%format suitable for training patch experts

%   Detailed explanation goes here

    img_size = [400, 400] * training_scale;

    %%
    addpath('PDM_helpers/');
    load 'PDM_helpers/pdm_68_aligned_wild.mat';

    output_location = '../prepared_data/mpie_';
    output_location = [output_location num2str(training_scale,3)];

    num_landmarks = 68;

    % Use mirror images to provide extra training data
    mirror_inds = [1,17;2,16;3,15;4,14;5,13;6,12;7,11;8,10;18,27;19,26;20,25;21,24;22,23;...
                  32,36;33,35;37,46;38,45;39,44;40,43;41,48;42,47;49,55;50,54;51,53;60,56;59,57;...
                  61,65;62,64;68,66];

    % The centres of views we want to extract
    centres_all = [   0,   0,  0
                      0, -20,  0
                      0, -45,  0
                      0, -70,  0
                      0,  20,  0
                      0,  45,  0
                      0,  70,  0
                      ];

    num_centers = size(centres_all, 1);

    counter_colour = zeros(num_centers,1);

    % read in all of the labels, together with names of images used

    landmark_labels = zeros(68,2,numel(multi_pie_labels));
    img_locations = cell(numel(multi_pie_labels),1);
    extra_locations = cell(numel(multi_pie_labels),1);
    img_dirs = cell(numel(multi_pie_labels),1);

    for i=1:numel(multi_pie_labels)

        img_locations{i} = [multi_pie_labels(i).img_dir, multi_pie_labels(i).actual_img];
        extra_locations{i} =  multi_pie_labels(i).img_locs;
        img_dirs{i} =  multi_pie_labels(i).img_dir;

        landmark_labels(:,:,i) = multi_pie_labels(i).landmark_labels;

    end

    % go through all images, see which centres they match and then add to
    % the appropriate bin
    num_imgs = size(landmark_labels,3);
    scales = zeros(num_imgs, 1);
    views = zeros(num_imgs, 1);
    for lbl=1:num_imgs

        Msm = M;
        Vsm = V;
        labels = landmark_labels(:,:,lbl);

        % Find the best PDM parameters given the 2D labels
        [ a,R ] = fit_PDM_ortho_proj_to_2D( Msm, E, Vsm, labels);

        eul = Rot2Euler(R);
        eul_orig = eul * 180 / pi;

        if(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '24_0'))
            eul = [0, -90, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '01_0'))
            eul = [0, -75, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '20_0'))
            eul = [0, -60, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '19_0'))
            eul = [0, -45, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '04_1'))
            eul = [0, -30, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '05_0'))
            eul = [0, -15, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '05_1'))
            eul = [0, 0, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '14_0'))
            eul = [0, 15, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '13_0'))
            eul = [0, 30, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '08_0'))
            eul = [0, 45, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '09_0'))
            eul = [0, 60, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '12_0'))
            eul = [0, 75, 0];
        elseif(strcmp([multi_pie_labels(lbl).cam_ID,'_', multi_pie_labels(lbl).cam_ID2], '11_0'))
            eul = [0, 90, 0];
        end

        % find the closest view (also mirror img?)
        [~, view] = min(sum(abs(centres_all - repmat(eul, num_centers, 1)),2));
        counter_colour(view) = counter_colour(view) + 1;
        scales(lbl) = a;
        views(lbl) = view;
    end

    % preallocate data
    allExamplesColourAllViews = cell(size(centres_all,1),1);
    landmarkLocationsAllViews = cell(size(centres_all,1),1);

    % if we don't have enough labelled original images add ones from diff
    % lighting conditions (that might not be perfectly labelled, but they
    % are better than fewer images)
    images_aim = 8000;
    extra_factors = ones(size(counter_colour));

    for r=1:size(centres_all,1)

         % see if extra is needed
        mirrorIdx = find(sum(abs(centres_all - repmat([centres_all(r,1), -centres_all(r,2), -centres_all(r,3)], size(centres_all,1),1)),2)==0);
        count = counter_colour(r) + counter_colour(mirrorIdx);

        if(count < images_aim)
            extra_factors(r) = images_aim / count;
        end

    end

    for r=1:size(centres_all,1)
        counter_colour(r) = round(counter_colour(r) * extra_factors(r));

        allExamplesColourAllViews{r} = uint8(zeros(counter_colour(r), img_size(1), img_size(2)));
        landmarkLocationsAllViews{r} = zeros(counter_colour(r), num_landmarks, 2);

        actual_imgs_used_all_views{r} = cell(counter_colour(r), 1);

    end

    counter_colour = zeros(num_centers,1);

    % The shape fitting is performed in the reference frame of the
    % patch training scale
    refGlobal = [training_scale, 0, 0, 0, 0, 0]';

    % go through all images and add to corresponding container
    for lbl=1:num_imgs

        labels = landmark_labels(:,:,lbl);
        occluded = labels(:,1) == 0;

        % Convert the labels to matlab format (we expect 1,1 to represent
        % the center of the top left pixel)
        labels = labels + 1;

        labels(occluded,:) = 0;

        imgCol = imread(img_locations{lbl});

        if(size(imgCol,3) == 3)
            imgCol = rgb2gray(imgCol);
        end

        % the reference shape
        [ ~, ~, ~,~, local_params] = fit_PDM_ortho_proj_to_2D( Msm, E, Vsm, labels);
        refShape = GetShapeOrtho(M, Vsm, local_params, refGlobal);

        % Create transform using a slightly modified version of Kabsch that
        % takes scaling into account as well, in essence we get a
        % similarity transform from current estimate to reference shape
        [A_img2ref, T_img2ref, ~, ~] = AlignShapesWithScale(labels(~occluded,:), refShape(~occluded,1:2));

        T_img2ref = T_img2ref + [img_size(1)/2, img_size(2)/2];

        % Create a transform, from shape in image to reference shape
        T = affine2d([A_img2ref';T_img2ref]);

        % transform the current shape to the reference one
        shape2D_in_ref = bsxfun(@plus, (A_img2ref * labels')', T_img2ref);

        % warp the image
        resizeColImage = imwarp(imgCol, T, 'linear', 'OutputView', imref2d(img_size));

        shape2D_in_ref(occluded,1) = 0;
        shape2D_in_ref(occluded,2) = 0;

        counter_colour(views(lbl)) = counter_colour(views(lbl)) + 1;

        allExamplesColourAllViews{views(lbl)}(counter_colour(views(lbl)),:,:) = resizeColImage;

        landmarkLocationsAllViews{views(lbl)}(counter_colour(views(lbl)),:,:) = shape2D_in_ref;
        actual_imgs_used_all_views{views(lbl)}{counter_colour(views(lbl))} = img_locations{lbl};

        % Here to add extra images if missing and not filled yet
        if(extra_factors(views(lbl)) > 1 && counter_colour(views(lbl)) < numel(actual_imgs_used_all_views{views(lbl)}))

            factor = extra_factors(views(lbl));

            % pick if to use this image or not
            if(randi(100) > (factor - floor(factor)) * 100)
                factor = floor(factor);
            else
                factor = ceil(factor);
            end

            while(factor > 1)

                lighting_id =  randi(20);
                [~, img_orig, ext] = fileparts(img_locations{lbl});
                img_orig = [img_orig, ext];
                if( strcmp(img_orig, extra_locations{lbl}{lighting_id}))
                    lighting_id = lighting_id + 1;
                    if(lighting_id > 20)
                        lighting_id = 1;
                    end
                end
                img_loc = [img_dirs{lbl}, extra_locations{lbl}{lighting_id}];
                imgCol = imread(img_loc);

                if(size(imgCol,3) == 3)
                    imgCol = rgb2gray(imgCol);
                end

                % resize the image to desired scale
                resizeColImage = imwarp(imgCol, T, 'linear', 'OutputView', imref2d(img_size));

                counter_colour(views(lbl)) = counter_colour(views(lbl)) + 1;

                allExamplesColourAllViews{views(lbl)}(counter_colour(views(lbl)),:,:) = resizeColImage;

                landmarkLocationsAllViews{views(lbl)}(counter_colour(views(lbl)),:,:) = shape2D_in_ref;
                actual_imgs_used_all_views{views(lbl)}{counter_colour(views(lbl))} = img_loc;
                factor = factor - 1;
            end

        end
        % Make sure same one not added as well

        if(mod(lbl, 100) == 0)
            fprintf('%d/%d done\n', lbl, num_imgs);
        end

    end

    % write out the training data for each view
    for r=1:size(centres_all,1)
        mirrorIdx = find(sum(abs(centres_all - repmat([centres_all(r,1), -centres_all(r,2), -centres_all(r,3)], size(centres_all,1),1)),2)==0);

        % if the mirrored view already added no need to do it again
        % so if (0,-20,0) is done no need to compute (0,20,0)
        if(mirrorIdx < r)
            continue
        end

        % Cap to two thousand images for space reasons, and because more
        % wouldn't actually be used
        max_images = 2000;

        mirrorImgs = allExamplesColourAllViews{mirrorIdx}(1:counter_colour(mirrorIdx),:,:);
        mirrorLbls = landmarkLocationsAllViews{mirrorIdx}(1:counter_colour(mirrorIdx),:,:);

        for i=1:size(mirrorImgs,1)

            flippedImg = fliplr(squeeze(mirrorImgs(i,:,:)));

            flippedLbls = squeeze(mirrorLbls(i,:,:));
            flippedLbls(:,1) = img_size(1) - flippedLbls(:,1) + 1;

            tmp1 = flippedLbls(mirror_inds(:,1),:);
            tmp2 = flippedLbls(mirror_inds(:,2),:);
            flippedLbls(mirror_inds(:,2),:) = tmp1;
            flippedLbls(mirror_inds(:,1),:) = tmp2;

            mirrorImgs(i,:,:) = flippedImg;
            mirrorLbls(i,:,:) = flippedLbls;

            % Ensure that after mirroring invalid (occluded) feature points
            % are set to 0 in both x and y dims
            mirrorLbls(i,flippedLbls(:,2)==0,1) = 0;

        end

        all_images = cat(1, allExamplesColourAllViews{r}(1:counter_colour(r),:,:), mirrorImgs);
        landmark_locations = cat(1, landmarkLocationsAllViews{r}(1:counter_colour(r),:,:), mirrorLbls);
        actual_imgs_used = cat(1, actual_imgs_used_all_views{r}(1:counter_colour(r)), actual_imgs_used_all_views{mirrorIdx}(1:counter_colour(mirrorIdx)));

        centres = centres_all(r,:);

        % identify the visibility of a point
        num_visible = sum(landmark_locations(:,:,1)~=0);
        visible_max = max(num_visible);

        visiIndex = ones(1,68);
        visiIndex(num_visible < 0.5*visible_max) = 0;

        if(size(all_images,1) > max_images)
            im_to_select = randperm(size(all_images,1));
            im_to_select = im_to_select(1:max_images);
            all_images = all_images(im_to_select,:,:);
            landmark_locations = landmark_locations(im_to_select,:,:);
            actual_imgs_used = actual_imgs_used(im_to_select);
        end
        save([output_location '_' num2str(r) '.mat'], 'all_images', 'landmark_locations', 'training_scale', 'centres', 'actual_imgs_used', 'visiIndex', '-v7.3');

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%    EPSON C8.
%
% The author of this script is:
%		  Silvia Moreno Serrano
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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= '6AXIS_C8';

%Path where everything is stored for this robot
robot.path = 'robots/EPSON/6AXIS_C8';

robot.DH.theta='[q(1)+pi/2 q(2)+pi/2 q(3) q(4) q(5) q(6)]';
robot.DH.d='[0.472 0 0 0.310 0 0.080]';
robot.DH.a='[0.100 0.300 0.03 0 0 0]';
robot.DH.alpha='[pi/2 0 pi/2 -pi/2 pi/2 0]';

robot.J=[];

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

%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 =[deg2rad(-240) deg2rad(240); %Axis 1, minimum, maximum
                deg2rad(-158) deg2rad(65); %Axis 2, minimum, maximum
                deg2rad(-61) deg2rad(202); %Axis 3
                deg2rad(-200) deg2rad(200); %Axis 4: Unlimited (400? default)
                deg2rad(-135) deg2rad(135); %Axis 5
                deg2rad(-360) deg2rad(360)]; %Axis 6: Really Unlimited to (800? default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(331); %Axis 1, rad/s
                deg2rad(332); %Axis 2, rad/s
                deg2rad(450); %Axis 3, rad/s
                deg2rad(450); %Axis 4, rad/s
                deg2rad(450); %Axis 5, rad/s
                deg2rad(720)];%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 = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; %%cambiar para dar la vuelta?

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

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [255 102 51]./255;
%for transparency
robot.graphical.draw_transparent=1;
%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.5 1 -0.5 1 0 1];
%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=[22 13 6 5 2 1];

%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:
function output = F_weighting(input, weight, bias)
[dim, N] = size(input);
nWeight = length(weight);
D = dim/nWeight;

input2 = reshape(input, D, nWeight, N);
input2 = permute(input2, [2 1 3]);

output = sum(bsxfun(@times, input2, weight));

if N==1
    output = output';
else
    output = permute(output, [2 3 1]);
end

output = bsxfun(@plus, output, bias);

end

author: singaxiong

repo: SignalGraph


text:
function h = box3(vol, I2X, bbox)
% Cut out a box from volumetric data and display it using texture mapping.
% h = box3(vol, I2X, bbox)
%
% vol, is a volume that is either N x M x K if it consists of scalar values
% or N x M x K x 3 if it is in RGB. The behaviour of the color handling is
% is similar to image3 and slice3.
% vol is the volume
% I2X is a transformation matrix, see slice3 for an explanation.
% bbox is a vector [i1 i2 j1 j2 k1 k2], where i1 is the lowest index of the
% box in the first dimension and i2 is the highest index of the box in the
% first data dimension, j1 is the lowest in the second data dimension and
% so on.
%
% SEE ALSO: slice3, image3

vol = vol(bbox(1):bbox(2),bbox(3):bbox(4),bbox(5):bbox(6),:);

I2Xp = I2X; I2Xp(:,4) = I2Xp(:,4) + I2X*[bbox(1)-0.5 bbox(3) bbox(5) 1]';
h1 = slice3(vol,I2Xp,1,1);
I2Xp = I2X; I2Xp(:,4) = I2Xp(:,4) + I2X*[bbox(1) bbox(3)-0.5 bbox(5) 1]';
h2 = slice3(vol,I2Xp,2,1);
I2Xp = I2X; I2Xp(:,4) = I2Xp(:,4) + I2X*[bbox(1) bbox(3) bbox(5)-0.5 1]';
h3 = slice3(vol,I2Xp,3,1);
I2Xp = I2X; I2Xp(:,4) = I2Xp(:,4) + I2X*[bbox(1)+0.5 bbox(3) bbox(5) 1]';
he1 = slice3(vol,I2Xp,1,size(vol,1));
I2Xp = I2X; I2Xp(:,4) = I2Xp(:,4) + I2X*[bbox(1) bbox(3)+0.5 bbox(5) 1]';
he2 = slice3(vol,I2Xp,2,size(vol,2));
I2Xp = I2X; I2Xp(:,4) = I2Xp(:,4) + I2X*[bbox(1) bbox(3) bbox(5)+0.5 1]';
he3 = slice3(vol,I2Xp,3,size(vol,3));

h = [h1,h2,h3,he1,he2,he3];

author: Sable

repo: mcbench-benchmarks