text:
function [x,info] = amgMaxwellinterface(A,b,node,edge,option)
%% AMGMAXWELL algebraic multigrid solver for Maxwell equations.
%
% x = amgMaxwell(A,b,node,edge) attempts to solve the system of
% linear equations Ax = b using multigrid type solver. The linear system is
% obtained by either the first or second family of linear edge element
% discretization of the Maxwell equation; See doc Maxwell.
%
% amgMaxwell is more algebraic than mgMaxwell but still requires geometric
% information node and edge. Grapha Laplacian of vertices are used as
% auxiliary Poisson operator and amg is used as Poisson solver.
%
% Input
%   -  A: curl(alpha curl) + beta I
%   -  b: right hand side
%   -  node,edge: mesh information
%   -  options: extra structures
%
% By default, the HX preconditioned PCG is used which works well for
% symmetric positive definite matrices (e.g. arising from eddy current
% simulation). For symmetric indefinite matrices (e.g. arising from time
% harmonic Maxwell equation), set option.solver = 'minres' or 'bicg' or
% 'gmres' to try other Krylov method with HX preconditioner.
%
% See also mg, mgMaxwell, mgMaxwellsaddle
%
% Reference:
% R. Hiptmair and J. Xu, Nodal Auxiliary Space Preconditioning in
% H(curl) and H(div) Spaces. SIAM J. Numer. Anal., 45(6):2483-2509, 2007.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

t = cputime;
%% Initial check
Ndof = length(b);                  % number of dof
N = size(node,1);                  % number of nodes
NE = size(edge,1);                 % number of edge;
dim = size(node,2);
% Assign default values to unspecified parameters
if ~exist('option','var')
    option = [];
end
option = mgoptions(option,length(b));    % parameters
x0 = option.x0;
% N0 = option.N0;
tol = option.tol;
%tol = 5*10^(-7);
% mu = option.smoothingstep; preconditioner = option.preconditioner; coarsegridsolver = option.coarsegridsolver;
printlevel = option.printlevel; %setupflag = option.setupflag;
maxIt = 2000; %400;   % increase the default step (200) for Maxwell equations
% tol = 1e-7;    % reset the tol
% Check for zero right hand side
if (norm(b) == 0)                  % if rhs vector is all zeros
    x = zeros(Ndof,1);             % then solution is all zeros
    flag = 0;
    itStep = 0;
    err = 0;
    time = cputime - t;
    info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));
    return
end

%% Transfer operators from nodal element to edge element
if isfield(option,'isBdEdge')
    isBdEdge = option.isBdEdge;
else
    deg = sum(spones(A(1:NE,1:NE)),2);
    isBdEdge = (deg == 1);
end
if Ndof == NE        % lowest order edge element
    II = node2edgematrix(node,edge,isBdEdge);
elseif Ndof >= 2*NE  % first or second order edge element
    II = node2edgematrix1(node,edge,isBdEdge);
end
IIt = II';
[grad,isBdNode] = gradmatrix(edge,isBdEdge);
gradt = grad';

%% Block smoother
%if option.blklevel == 0 || ~isfield(option,'blklevel')

interfaceEdgeID = option.blkId+1:NE;
interfaceEdge = false(NE, 1);
interfaceEdge(option.blkId+1:NE) = true;
regularEdge = false(NE, 1);
regularEdge(1:option.blkId) = true;

if isfield(option,'blklevel')

    ExtLevel = 0;
    while (ExtLevel<option.blklevel)

        V2Emat= sparse([1:NE,1:NE]',[edge(:,1);edge(:,2)],1,NE,N);
        interfaceNodeID = union(edge(interfaceEdgeID,1),edge(interfaceEdgeID,2));
        interfaceNode = false(N, 1); interfaceNode(interfaceNodeID) = true;
        interfaceEdge = (V2Emat*interfaceNode>0);
        regularEdge = ~interfaceEdge;
        interfaceEdgeID = find(interfaceEdge==1);
        ExtLevel = ExtLevel+1;

    end

end

Ainterface = A(interfaceEdge,interfaceEdge);
Aregular = A(regularEdge,regularEdge);
B = A(interfaceEdge, regularEdge);
D = diag(Aregular);

tic
if isfield(option,'fact') == 1
    switch option.fact
        case 'lu'
            [L,U,P,Q] = lu(Ainterface);
        case 'chol'
            [R,flag,P] = chol(Ainterface);
    end
else
    [L,U,P,Q] = lu(Ainterface);
end
toc

%% Auxiliary Poisson matrix
%   -  A: curl(alpha curl) + beta I
%   - AP: - div(alpha grad) + beta I
%   - BP: - div(beta grad)

if isfield(option,'AP')
    AP = option.AP;
   edgeVec = node(edge(:,2),:) - node(edge(:,1),:);
%     edgeLength = sqrt(sum(edgeVec.^2,2));
%     if isfield(option,'beta') % resacle by the dielectric coefficients
%         if isreal(option.beta) && (length(option.beta) == NE)
%            beta = option.beta;
%         else % option.beta is a function
%            edgeMiddle = (node(edge(:,2),:) + node(edge(:,1),:))/2;
%            beta = option.beta(edgeMiddle);
%         end
%     end
%     M = accumarray(edge(:),repmat((edgeLength.^3).*beta,2,1),[N 1]);
%     AP = AP + spdiags(M,0,N,N);
else
    % build graph Laplacian to approximate AP
    edgeVec = node(edge(:,2),:) - node(edge(:,1),:);
    edgeLength = sqrt(sum(edgeVec.^2,2));
    if isfield(option,'alpha') % resacle by the magnetic coefficients
        if isreal(option.alpha) && (length(option.alpha) == NE)
           alpha = option.alpha;
        else % option.alpha is a function
           edgeMiddle = (node(edge(:,2),:) + node(edge(:,1),:))/2;
           alpha = option.alpha(edgeMiddle);
        end
    end
    % edge weight: h*alpha
    % AP = gradt*spdiags(edgeLength.*alpha,0,NE,NE)*grad;
    i1 = (1:NE)'; j1 = double(edge(:,1)); s1 = sqrt(edgeLength.*alpha);
    i2 = (1:NE)'; j2 = double(edge(:,2)); s2 = -s1;
    isFreeEdge = ~isBdEdge;
    G = sparse([i1(isFreeEdge);i2(isFreeEdge)],...
        [j1(isFreeEdge);j2(isFreeEdge)],...
        [s1(isFreeEdge);s2(isFreeEdge)],NE,N);
    %i1 = (1:NE)'; j1 = double(edge(:,1)); s1 = ones(size(alpha)).*sqrt(alpha);
%     i2 = (1:NE)'; j2 = double(edge(:,2)); s2 = -s1;
%     isFreeEdge = ~isBdEdge;
%     G = sparse([i1(:);i2(:)],...
%         [j1(:);j2(:)],...
%         [s1(:);s2(:)],NE,N);
    AP = G'*G;
    %E = G'*A*G;
    % lumped mass matrix: h^3*beta
    if isfield(option,'beta') % resacle by the dielectric coefficients
        if isreal(option.beta) && (length(option.beta) == NE)
           beta = option.beta;
        else % option.beta is a function
           edgeMiddle = (node(edge(:,2),:) + node(edge(:,1),:))/2;
           beta = option.beta(edgeMiddle);
        end
    end
    M = accumarray(edge(:),repmat((edgeLength.^3).*beta,2,1),[N 1]);
    AP = AP + spdiags(M,0,N,N);
end
% BP is Galerkin projection to the free node space
% boundary nodes
bdidx = zeros(N,1);
bdidx(isBdNode) = 1;
Tbd = spdiags(bdidx,0,N,N);
BP = gradt*A(1:NE,1:NE)*grad + Tbd;
%BP = gradt*spdiags(edgeLength.*beta,0,NE,NE)*grad + Tbd;
% if strcmp(option.outsolver,'minres')
%     BP = gradt*option.BPP(1:NE,1:NE)*grad + Tbd;
% end

setupOption.solver = 'NO';
[x,info,APi,Ri,RRi,ResAP,ProAP,clA] = amg(AP,ones(N,1),setupOption); %#ok<ASGLU>
[x,info,BPi,Si,SSi,ResBP,ProBP,clB] = amg(BP,ones(N,1),setupOption); %#ok<ASGLU>

%% Krylov iterative methods with HX preconditioner
k = 1;
err = 1;
switch upper(option.outsolver)
    case 'CG'
        if printlevel>=1
            fprintf('Conjugate Gradient Method using HX preconditioner \n');
        end
        x = x0;
        r = b - A*x;
        nb = norm(b);
        err = zeros(maxIt,2);
        err(1,:) = norm(r)/nb;
        while (max(err(k,:)) > tol) && (k <= maxIt)
            % compute Br by HX preconditioner
            Br = HXpreconditioner(r);
            % update tau, beta, and p
            rho = r'*Br;  % r'*Br = e'*ABA*e approximates e'*A*e
            if k==1
                p = Br;
            else
                beta = rho/rho_old;
                p = Br + beta*p;
            end
            % update alpha, x, and r
            Ap = A*p;
            alpha = rho/(Ap'*p);
            r = r - alpha*Ap;
            x = x + alpha*p;
            rho_old = rho;
            % compute err for the stopping criterion
            k = k + 1;
            err(k,1) = sqrt(abs(rho/(x'*b))); % approximate relative error in energy norm
            err(k,2) = norm(r)/nb; % relative error of the residual in L2-norm
            if printlevel >= 2
                fprintf('#dof: %8.0u, HXCG iter: %2.0u, err = %12.8g\n',...
                         Ndof, k, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
        if k > maxIt || (max(err(end,:))>tol)
            flag = 1;
        else
            flag = 0;
        end
    case 'MINRES'
        fprintf('Minimum Residual Method with HX preconditioner \n')
        [x,flag,err,itStep] = minres(A,b,tol,maxIt,@HXpreconditioner,[],x0);
%         x = minres(A,b,tol,maxIt,@HXpreconditioner,[],x0);
    case 'GMRES'
        fprintf('General Minimum Residual Method with HX preconditioner \n')
        tic
        [x,flag,err,itStep] = gmres(A,b,10,tol,maxIt,@HXpreconditioner,[],x0);
        toc
        itStep = (itStep(1) -1)*10 + itStep(2);
end

%% Output
time = cputime - t;
if printlevel >= 1
    fprintf('#dof: %8.0u,   #nnz: %8.0u,   iter: %2.0u,   err = %8.4e,   time = %4.2g s\n',...
                 Ndof, nnz(A), itStep, max(err(end,:)), time)
end
if  (flag == 1) && (printlevel>0)
    fprintf('NOTE: the iterative method does not converge');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions HXpreconditioner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function Br = HXpreconditioner(r)
    %% 1. Smoothing in the finest grid of the original system
    if strcmp(option.smoother,'BD')
        if isfield(option,'fact') == 1
            switch option.fact
                case 'lu'
                    eh = zeros(NE, 1);
                    eh(regularEdge,:) = 0.75*r(regularEdge,:)./D;
                    %        eh(interfaceEdge,:) = Ainterface\r(interfaceEdge,:);
                    rtmp = P*r(interfaceEdge,:);
                    rtmp = L\rtmp;
                    rtmp = U\rtmp;
                    rtmp = Q*rtmp;
                    eh(interfaceEdge,:) = rtmp;
                case 'chol'
                    eh = zeros(NE, 1);
                    eh(regularEdge,:) = 0.75*r(regularEdge,:)./D;
                    %        eh(interfaceEdge,:) = Ainterface\r(interfaceEdge,:);
                    rtmp = P'*r(interfaceEdge,:);
                    rtmp = R'\rtmp;
                    rtmp = R\rtmp;
                    rtmp = P*rtmp;
                    eh(interfaceEdge,:) = rtmp;
            end
        else
            eh = zeros(NE, 1);
            eh(regularEdge,:) = 0.75*r(regularEdge,:)./D;
            %        eh(interfaceEdge,:) = Ainterface\r(interfaceEdge,:);
            rtmp = P*r(interfaceEdge,:);
            rtmp = L\rtmp;
            rtmp = U\rtmp;
            rtmp = Q*rtmp;
            eh(interfaceEdge,:) = rtmp;
        end
    elseif strcmp(option.smoother,'GS')
        %         eh = triu(A)\(D.*(tril(A)\r)); # does not work if interface edge is
        %         not treated differently
        rN = r(regularEdge); rI = r(interfaceEdge);
        rN = rN./D; rI = Ainterface\rI;
        ehN = rN + B'*(Ainterface\(B*rN))./D - B'*rI./D;
        ehI = -Ainterface\(B*rN) + rI;
        eh = [ehN; ehI];
    else
        eh = 0.75*r./D;  % Jacobi method. less computational time
    end
    %% 2. amg solver for auxiliary operators
    amgoption.solvermaxit = 3;  % 3 for SPD matrix
    amgoption.printlevel = 0;
    rc = reshape(IIt*r(1:size(IIt,2)),N,dim);   % transfer to the nodal linear element space
    %eaux = II*reshape(AP\rc,dim*N,1);
    %eaux = II*reshape(amg(AP,rc,amgoption),dim*N,1);
    rAP = amgInterface(AP,rc,APi,Ri,RRi,ResAP,ProAP,clA,amgoption);
    eaux = II*reshape(rAP,dim*N,1);
    eh = eh + eaux;
    rb = gradt*r(1:NE);
    %eauxb = grad*(BP\rb);
    %eauxb = grad*(amg(BP,rb,amgoption));
    rBP = amgInterface(BP,rb,BPi,Si,SSi,ResBP,ProBP,clB,amgoption);
    eauxb = grad*rBP;
    Br = eh + eauxb;
    end

    function z = amgInterface(M,r0,Ai,Bi,BBi,Res,Pro,cl,amgoption)

        maxItz = amgoption.solvermaxit;
        Nb = size(r0,2);
        z = repmat(zeros(size(r0,1),1),1,Nb);
        kz = 1;
        rz = r0 - M*z;
        nbz = max(sqrt(sum(r0.^2,1)));
        errz = zeros(maxItz,2);
        errz(1,:) = max(sqrt(sum(rz.^2,1)))/nbz;

        amgprintlevel = amgoption.printlevel;
        level = max(min(ceil(log2(N)/2-4),8),2);
        prefunc = @wcycle;
        if amgprintlevel >= 1
            fprintf('Conjugate Gradient Method\n')
        end
        if isfield(amgoption,'smoothingstep')  % smoothing steps
            mu = amgoption.smoothingstep;
        else
            mu = 2;
        end

        while (max(errz(kz,:)) > tol) && (kz <= maxItz)
            % compute Br by MG
            Brz = prefunc(r0);
            % update tau, beta, and p
            rhoz = dot(Brz,r0);  % e'*ABA*e approximates e'*A*e
            if kz == 1
                pz = Brz;
            else
                betaz = rhoz./rho_oldz;
                pz = Brz + betaz.*pz;
            end
            % update alpha, x, and r
            Apz = M*pz;
            alpha = rhoz./dot(Apz,pz);
            r0 = r0 - alpha.*Apz;
            z = z + alpha.*pz;
            rho_oldz = rhoz;
            kz = kz + 1;
            % compute err for the stopping criterion
        %     err(k,1) = alpha*sqrt(p'*Ap/(x'*A*x)); % increamental error in energy norm
            errz(kz,1) = max(sqrt(abs(rhoz./dot(z,r0)))); % approximate relative error in energy norm
            errz(kz,2) = max(sqrt(sum(r0.^2,1)))/nbz; % relative error of the residual in L2-norm
            if amgprintlevel >= 2
                fprintf('#dof: %8.0u, MGCG iter: %2.0u, err = %8.4e\n',...
                         N, kz-1, max(errz(kz,:)));
            end
        end
        errz = errz(1:kz,:);
        itStepz = kz-1;

        function e = wcycle(r,J)        % solve equations Ae = r in each level
            if nargin<=1
                J = level;
            end
            if J == cl
                e = Ai{cl}\r;   % exact solver in the coaresest grid
                return
            end
            % fine grid pre-smoothing
            e = Bi{J}\r;   % pre-smoothing
            for s = 1:mu           % extra mu steps smoothing
                e = e + Bi{J}\(r-Ai{J}*e);
            end
            rc = Res{J}*(r - Ai{J}*e);
            % coarse grid correction twice
            ec = wcycle(rc,J-1);
            ec = ec + wcycle(rc - Ai{J-1}*ec,J-1);
            % fine grid post-smoothing
            e = e + Pro{J-1}*ec;
            e = e + BBi{J}\(r-Ai{J}*e);
            for s = 1:mu
                e = e + BBi{J}\(r-Ai{J}*e); % post-smoothing
            end
        end

    end

end

repo: ifem


text:
classdef FileLoadSaveClass < matlab.mixin.Copyable

    properties (Access = private)
        filename;
        fileformat;
        supportedFomats;
        err;
        dataStorageScheme;
    end

    methods

        % ----------------------------------------------------------------------------------
        function obj = FileLoadSaveClass()
            obj.filename = '';
            obj.fileformat = '';
            obj.supportedFomats = struct( ...
                'matlab', {{'.mat','matlab','mat'}}, ...
                'hdf5', {{'hdf','.hdf','hdf5','.hdf5','hf5','.hf5','h5','.h5'}} ...
                );
            obj.err = 0;
            obj.dataStorageScheme = 'memory';
        end

        % ---------------------------------------------------------
        function Load(obj, filename, params, format)
            if ~exist('filename','var')
                filename = obj.filename;
            end
            if ~exist('format','var')
                format = obj.fileformat;
            elseif obj.Supported(format)
                obj.fileformat = format;
            end
            if ~exist('params','var')
                params = [];
            end

            switch(lower(format))
                case obj.supportedFomats.matlab
                    if ismethod(obj, 'LoadMat')
                        obj.err = obj.LoadMat(filename, params);
                    end
                case obj.supportedFomats.hdf5
                    if ismethod(obj, 'LoadHdf5')
                        obj.err = obj.LoadHdf5(filename, params);
                    end
            end
        end

        % ---------------------------------------------------------
        function Save(obj, filename, params, format)
            if ~exist('filename','var')
                filename = obj.filename;
            end
            if ~exist('params','var')
                params = [];
            end
            if ~exist('format','var')
                format = obj.fileformat;
            end

            switch(lower(format))
                case obj.supportedFomats.matlab
                    if ismethod(obj, 'SaveMat')
                        obj.SaveMat(filename, params);
                    end
                case obj.supportedFomats.hdf5
                    if ismethod(obj, 'SaveHdf5')
                        obj.SaveHdf5(filename, params);
                    end
            end
        end

        % ---------------------------------------------------------
        function b = Supported(obj, format)
            b = true;
            switch(lower(format))
                case obj.supportedFomats.matlab
                    return;
                case obj.supportedFomats.hdf5
                    return;
            end
            b = false;
        end

        % -------------------------------------------------------
        function B = ne(obj, obj2)
            if obj==obj2
                B = false;
            else
                B = true;
            end
        end

        % -------------------------------------------------------
        function SetFileFormat(obj, fmt)
            obj.fileformat = fmt;
        end

        % -------------------------------------------------------
        function fmt = GetFileFormat(obj)
            fmt = obj.fileformat;
        end

        % -------------------------------------------------------
        function SetFilename(obj, fname)
            obj.filename = fname;
        end

        % -------------------------------------------------------
        function fname = GetFilename(obj)
            fname = obj.filename;
        end

        % -------------------------------------------------------
        function SetDataStorageScheme(obj, scheme)
            obj.dataStorageScheme = scheme;
        end

        % -------------------------------------------------------
        function scheme = GetDataStorageScheme(obj)
            scheme = obj.dataStorageScheme;
        end

        % -------------------------------------------------------
        function SetError(obj, err)
            obj.err = err;
        end

        % -------------------------------------------------------
        function err = GetError(obj)
            err = [];
            if isempty(obj)
                return
            end
            err = obj.err;
        end

        % -------------------------------------------------------
        function supportedFomats = GetSupportedFormats(obj)
            supportedFomats = obj.supportedFomats;
        end

    end

end

author: BUNPC


text:
function [pOutIdxs,chOutIdxs] = treeWpBFrange(wt)
%TREEWPBFRANGE Wavelet packet tree output ranges in BF order
%   Usage: [pOutIdxs,chOutIdxs] = treeBFranges(wt);
%
%   Input parameters:
%         wt       : Filterbank tree struct.
%   Output parameters:
%         pOutIdxs  : Array of parent nodes in BF order
%         chOutIdxs : Cell array of children nodes in BF order
%
%   `[pOutIdxs,chOutIdxs] = treeBFranges(wt)` is a helper function
%   determining direct relationship between nodes in tree *wt*.
%   Elements in both returned arrays are ordered according to the BF order.
%   `pOutIdxs` is array of indices in the subbands
%

treePath = nodeBForder(0,wt);
trLen = numel(treePath);
pOutIdxs = zeros(1,trLen);
chOutIdxs = cell(1,trLen);
pRunIdx = [0];
chRunIdx = 1;
% do trough tree and look for nodeNo and its parent
for ii=1:trLen
    tmpfiltNo = length(wt.nodes{treePath(ii)}.g);
    locRange = nodesLocOutRange(treePath(ii),wt);
    diffRange = 1:tmpfiltNo;
    diffRange(locRange{1})=[];
    chOutIdxs{ii} = chRunIdx:chRunIdx+tmpfiltNo-1;
    chRunIdx = chRunIdx + tmpfiltNo;
    pOutIdxs(ii) = pRunIdx(1);
    pRunIdx = [pRunIdx(2:end),chOutIdxs{ii}(diffRange)];
end

pOutIdxs = pOutIdxs(end:-1:1);
chOutIdxs = chOutIdxs(end:-1:1);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function varargout = drawGrid3d(varargin)
%DRAWGRID3D Draw a 3D grid on the current axis.
%
%   drawGrid3d
%   draws a 3D square grid, with origin (0,0,0) and spacing 1 in each
%   direction, with bounds corresponding to the bounds of current axis.
%
%   drawGrid3d(SPACING)
%   where spacing is either a scalar or a [1x3] matrix, specifies the size
%   of the unit cell.
%
%   drawGrid3d(ORIGIN, SPACING)
%   Also specify origin of grid. ORIGIN is a [1x3] array.
%
%   drawGrid3d(..., EDGE)
%   specifies whether function should draw edges touching edges of axis.
%   EDGE is a characheter string, which can be :
%   - 'OPEN' : each line start from one face of window to the opposite
%   face. This results in a 'spiky' grid.
%   - 'CLOSED' (default value) : each line stops at the last visible point
%   of the grid for this line. The result looks like a box (no free spikes
%   around the grid).
%
%   H = drawGrid3d(...);
%   return a vector of handles for each LINE object which was crated.
%

%   ------
%   Author: David Legland
%   e-mail: david.legland@grignon.inra.fr
%   Created: 2005-11-17
%   Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%% initialize variables -----

% default values
closed = true;
origin = [0 0 0];
spacing = [1 1 1];

% check if grid is open or not
str = '';
if ~isempty(varargin)
    str = varargin{end};
end
if ischar(str)
    if strncmpi(str, 'open', 4)
        closed = false;
    end
    varargin = varargin(1:end-1);
end

% check origin and grid spacing
if length(varargin)==1
    spacing = varargin{1};
elseif length(varargin)==2
    origin = varargin{1};
    spacing = varargin{2};
end

%% Compute internam data -----

% get axis limits
ax = axis;
x0 = ax(1); x1 = ax(2);
y0 = ax(3); y1 = ax(4);
z0 = ax(5); z1 = ax(6);

% get first and last coordinates of the grid in each direction
dx = spacing(1); dy = spacing(2); dz = spacing(3);
xe = x0 + mod(origin(1) - x0, dx);
xf = x1 - mod(x1 - origin(1), dx);
ye = y0 + mod(origin(2) - y0, dy);
yf = y1 - mod(y1 - origin(2), dy);
ze = z0 + mod(origin(1) - z0, dz);
zf = z1 - mod(z1 - origin(1), dz);

% update first and last coordinate if grid is 'closed'
if closed
    x0 = xe; x1 = xf;
    y0 = ye; y1 = yf;
    z0 = ze; z1 = zf;
end

%% Draw the grid -----

h = [];
%TODO: rewrite code, avoiding loops

% draw lines parallel to x axis
for y = ye:dy:yf
    for z = ze:dz:zf
        h = [h; drawEdge3d([x0 y z x1 y z])]; %#ok<AGROW>
    end
end

% draw lines parallel to y axis
for x = xe:dx:xf
    for z = ze:dz:zf
        h = [h; drawEdge3d([x y0 z x y1 z])]; %#ok<AGROW>
    end
end

% draw lines parallel to z axis
for x = xe:dx:xf
    for y = ye:dy:yf
        h = [h; drawEdge3d([x y z0 x y z1])]; %#ok<AGROW>
    end
end

%% Check output arguments -----

if nargout>0
    varargout{1} = h;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function F = normalizePivots(F)
%NORMALIZEPIVOTS   Scale rows and cols of a SEPARABLEAPPROX so that all pivots are 1.
%
% Additionally, the norm of the kth row and column will be the same.

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

% TODO: Document
% TODO: is this useful?

F = normalizeRowsAndCols(F);

d = F.pivotValues(:).';
s = sign(d);
sqrtp = sqrt(abs(d));
F.cols = F.cols./(s.*sqrtp);
F.rows = F.rows./sqrtp;
F.pivotValues = ones(1, length(d));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Y=trapz(X,I,dummy)
%TRAPZ (overloaded)

try
    n = X.dim(1);
    m = X.dim(2);
    if nargin==1
        Y = X;
        if n==1 || m==1
            % Spezialized code...
            Y.basis = trapz(Y.basis,1);
            Y.dim(1) = 1;
            Y.dim(2) = 1;
            temp = 1;
        else
            % Standard case

            temp = trapz(reshape(X.basis(:,1),n,m));
            Y.basis =  kron(speye(m),ones(1,n))*X.basis;

        end
    else
        Y = X;
        temp = trapz(reshape(X.basis(:,1),n,m),I);
        Y.basis = temp(:);
        for i = 1:length(Y.lmi_variables)
            temp = trapz(reshape(X.basis(:,i+1),n,m),I);
            Y.basis(:,i+1) = temp(:);
        end
    end
catch
    error(lasterr)
end
Y.dim(1) = size(temp,1);
Y.dim(2) = size(temp,2);
% Reset info about conic terms
Y.conicinfo = [0 0];
Y.extra.opname = '';
Y = clean(Y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function fieldnames = listAllFieldNames(planC)
% List all fieldnames in the plan.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
cellNames = fields(indexS);
for i=1:length(cellNames)
    cellVals(i) = indexS.(cellNames{i});
end

fieldnames = [];
%Call listFields on each cell in planC.  Manually generate the index alias.
for i=1:length(cellNames)
    fieldnames = [fieldnames;listFields(planC{cellVals(i)}, ['planC{indexS.' cellNames{i} '}'])];
end
fieldnames = unique(fieldnames);

function fieldnames = listFields(structure, myName)
% list fields in struct given the name of the structure above.
fieldnames = {myName};
dataType = class(structure);
switch dataType
    case 'struct'
        for i=1:length(structure(:))
            fieldNames = fields(structure);
            for j=1:length(fieldNames)
                childName = [myName '.' fieldNames{j}];
                childFields = listFields(structure(i).(fieldNames{j}), childName);
                fieldnames = [fieldnames;childFields];
            end
        end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [EToE,EToF]= tiConnect2D(EToV)

% function [EToE,EToF]= tiConnect2D(EToV)
% Purpose: triangle face connect algorithm due to Toby Isaac

Nfaces=3;
K = size(EToV,1);
Nnodes = max(max(EToV));

% create list of all faces 1, then 2, & 3
fnodes = [EToV(:,[1,2]);EToV(:,[2,3]);EToV(:,[3,1])];
fnodes = sort(fnodes,2)-1;

% set up default element to element and Element to faces connectivity
EToE= (1:K)'*ones(1,Nfaces); EToF= ones(K,1)*(1:Nfaces);

% uniquely number each set of three faces by their node numbers
id = fnodes(:,1)*Nnodes + fnodes(:,2)+1;
spNodeToNode=[id, (1:Nfaces*K)', EToE(:), EToF(:)];

% Now we sort by global face number.
sorted=sortrows(spNodeToNode,1);

% find matches in the sorted face list
[indices,dummy]=find( sorted(1:(end-1),1)==sorted(2:end,1) );

% make links reflexive
matchL = [sorted(indices,:)   ;sorted(indices+1,:)];
matchR = [sorted(indices+1,:) ;sorted(indices,:)];

% insert matches
EToE(matchL(:,2)) = matchR(:,3); EToF(matchL(:,2)) = matchR(:,4);
return;

repo: Aero-matlab


text:
classdef TestMinEnclosingCircle
    %TestMinEnclosingCircle

    properties (Constant)
        p = [0 0; 1 0; 2 2; 3 3; 3 4];
    end

    methods (Static)
        function test_1
            [c,r] = cv.minEnclosingCircle(TestMinEnclosingCircle.p);
            validateattributes(c, {'numeric'}, {'vector', 'numel',2});
            validateattributes(r, {'numeric'}, {'scalar'});
        end

        function test_2
            [c,r] = cv.minEnclosingCircle(num2cell(TestMinEnclosingCircle.p,2));
            validateattributes(c, {'numeric'}, {'vector', 'numel',2});
            validateattributes(r, {'numeric'}, {'scalar'});
        end

        function test_error_argnum
            try
                cv.minEnclosingCircle();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function predicted=cosmo_classify_libsvm(samples_train, targets_train, samples_test, opt)
% libsvm-based SVM classifier
%
% predicted=cosmo_classify_libsvm(samples_train, targets_train, samples_test, opt)
%
% Inputs
%   samples_train      PxR training data for P samples and R features
%   targets_train      Px1 training data classes
%   samples_test       QxR test data
%   opt                (optional) struct with options for svmtrain
%     .autoscale       If true (default), z-scoring is done on the training
%                      set; the test set is z-scored using the mean and std
%                      estimates from the training set.
%     ?                any option supported by libsvm's svmtrain.
%
% Output
%   predicted          Qx1 predicted data classes for samples_test
%
% Notes:
%  - this function requires libsvm version 3.18 or later:
%    https://github.com/cjlin1/libsvm
%  - by default a linear kernel is used ('-t 0')
%  - this function uses LIBSVM's svmtrain function, which has the same
%    name as matlab's builtin version. Use of this function is not
%    supported when matlab's svmtrain precedes in the matlab path; in
%    that case, adjust the path or use cosmo_classify_matlabsvm instead.
%  - for a guide on svm classification, see
%      http://www.csie.ntu.edu.tw/~cjlin/papers/guide/guide.pdf
%  - By default this function performs z-scoring of the data. To switch
%    this off, set 'autoscale' to false
%  - cosmo_crossvalidate and cosmo_crossvalidation_measure
%    provide an option 'normalization' to perform data scaling
%
%
% See also svmtrain, svmclassify, cosmo_classify_svm,
%          cosmo_classify_matlabsvm
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<4
        opt=struct();
    end

    % support repeated testing on different data after training every time
    % on the same data. This is achieved by caching the training data
    % and associated model
    persistent cached_targets_train;
    persistent cached_samples_train;
    persistent cached_opt;
    persistent cached_model;

    cache_limit=1e5; % avoid caching huge training sets

    if isequal(cached_targets_train, targets_train) && ...
            isequal(cached_opt, opt) && ...
            numel(samples_train) < cache_limit && ...
            isequal(cached_samples_train, samples_train)
        % use cache
        model=cached_model;
    else
        model=train(samples_train,targets_train,opt);

                % store model
        cached_targets_train=targets_train;
        cached_samples_train=samples_train;
        cached_opt=opt;
        cached_model=model;
    end

    predicted=test(model, samples_test);

function model=train(samples_train,targets_train,opt)
    [ntrain, nfeatures]=size(samples_train);

    model=struct();
    model.nfeatures=nfeatures;
    model.normalize=[]; % off by default

    % check input size
    ntrain_=numel(targets_train);
    if ntrain_~=ntrain
        error('illegal input size');
    end

    % construct options string for svmtrain
    opt_str=libsvm_opt2str(opt);

    % auto-scale is the default
    autoscale= ~isstruct(opt) || ...
                ~isfield(opt,'autoscale') || ...
                isempty(opt.autoscale) || ...
                opt.autoscale;

    % perform autoscale if necessary
    if autoscale
        [samples_train, params]=cosmo_normalize(samples_train, ...
                                                    'zscore', 1);
        model.normalize=params;
    end
    train_func=@()svmtrain(targets_train(:), ...
                            samples_train, ...
                            opt_str);

    model.libsvm_model=eval_with_check_external(train_func);

function output=eval_with_check_external(func)
% Evaluates func() in try-catch block. If it fails, this may be due
% to missing libsvm and/or libsvm conflicting with Matlab's svm.
% Therefore first cosmo_check_external is called, which will give an
% informative error message if that is the case. Otherwise the original
% message is shown.

    try
        output=func();
    catch
        cosmo_check_external('libsvm');
        rethrow(lasterror());
    end

function predicted=test(model, samples_test)
    [ntest, nfeatures]=size(samples_test);
    if nfeatures~=model.nfeatures
        error(['Number of features in train set (%d) and '...
                    'test set (%d) do not match'],...
                    model.nfeatures,nfeatures);
    end

    if ~isempty(model.normalize)
        samples_test=cosmo_normalize(samples_test, model.normalize);
    end

    test_opt_str='-q'; % quiet (no output)
    test_func=@()svmpredict(NaN(ntest,1), samples_test, ...
                            model.libsvm_model, test_opt_str);
    predicted=eval_with_check_external(test_func);

function opt_str=libsvm_opt2str(opt)
    persistent cached_opt
    persistent cached_opt_str

    if ~isequal(opt, cached_opt)
        default_opt={'t','q';...
                        '0',''};
        n_default=size(default_opt,2);

        libsvm_opt_keys={'s','t','d','g','r','c','n','p',...
                         'm','e','h','n','wi','v'};
        opt_struct=cosmo_structjoin(opt);

        keys=intersect(fieldnames(opt_struct),libsvm_opt_keys);
        n_keys=numel(keys);

        use_defaults=true(1,n_default);

        libsvm_opt=cell(2,n_keys);
        for k=1:n_keys
            key=keys{k};
            default_pos=find(cosmo_match(default_opt(1,:),(key)),1);

            if ~isempty(default_pos)
                use_defaults(default_pos)=false;
            end

            value=opt_struct.(key);
            if isnumeric(value)
                value=sprintf('%d',value);
            end

            libsvm_opt(:,k)={key;value};
        end

        all_opt=[libsvm_opt default_opt(:, use_defaults)];

        cached_opt_str=sprintf('-%s %s ', all_opt{:});
        cached_opt=opt;
    end

    opt_str=cached_opt_str;

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [mResult, mResult2]=calc_pwtime(a,fMinTime, fMaxTime, fTstep, timef,bootloops,maepi, nMod)
% function [mResult] = calc_pwtime(a,fMinTime, fMaxTime, fTstep, timef,bootloops,maepi, nMod);
% -------------------------------------------------------------------------
% Calculate and plot p,c,k-value evolution with time
%
% Input variables:
% a         : earthquake catalog
% fMinTime  : Starting time after mainshock to calculate Omori-parameters
%             in days
% fMaxTime  : End time after mainshock to calculate Omori-parameters
%             in days
% fTstep    : Time step length fMinTime:fTstep:fMaxTime [days]
% timef     : dummy for forecast period, usually =1
% bootloops : Number of bootstraps
% maepi     : mainshock
% nMod      : Models to fit aftershock sequence
%
% j.woessner@sed.ethz.ch
% last update: 29.09.2004

% Surpress warnings from fmincon
warning off;

% Initialize
mResult = [];
mResult2 = [];

% Time loop
for time = fMinTime:fTstep:fMaxTime
    %[m_main, main] = max(a.Magnitude);
    date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
    date_main = datenum(floor(maepi(3)),maepi(4),maepi(5),maepi(8),maepi(9),0);
    time_aftershock = date_matlab-date_main;
    % Select biggest aftershock earliest in time, but more than 1 day after mainshock
    fDay = 1;
    ft_c=fDay/365; % Time not considered to find biggest aftershock
    vSel = (a.Date > maepi(:,3)+ft_c & a.Date<= maepi(:,3)+time/365);
    mCat = a.subset(vSel);
    vSel = mCat(:,6) == max(mCat(:,6));
    vBigAf = mCat(vSel,:);
    if length(mCat(:,1)) > 1
        vSel = vBigAf(:,3) == min(vBigAf(:,3));
        vBigAf = vBigAf(vSel,:);
    end

    date_biga = datenum(floor(vBigAf(3)),vBigAf(4),vBigAf(5),vBigAf(8),vBigAf(9),0);
    fT1 = date_biga - date_main; % Time of big aftershock

    % Aftershock times
    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = a.subset(l);

    % time_as: Learning period
    l = tas <= time;
    time_as=tas(l);
    if length(time_as) > 50

        % Times up to the forecast time
        lf = tas <= time+timef ;
        time_asf= [tas(lf) ];
        time_asf=sort(time_asf);

        % Calculate uncertainty and mean values of p,c,and k
        [mMedModF, mStdL, loopout] = brutebootloglike_a2(time_as, time_asf, bootloops,fT1,nMod);
        pmed1 = mMedModF(1,1);
        pstd1 = mMedModF(1,2);
        pmed2 = mMedModF(1,3);
        pstd2 = mMedModF(1,4);
        cmed1 = mMedModF(1,5);
        cstd1 = mMedModF(1,6);
        cmed2 = mMedModF(1,7);
        cstd2 = mMedModF(1,8);
        kmed1 = mMedModF(1,9);
        kstd1 = mMedModF(1,10);
        kmed2 = mMedModF(1,11);
        kstd2 = mMedModF(1,12);

        % Compute model according to model choice
        if nMod == 1
            [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
        elseif nMod == 2
            [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
        elseif nMod == 3
            [pval1, pval2, cval1, cval2, kval1, kval2 , fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
        else
            [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
        end

        % Calculate aftershock sequence ...
        cumnrf = (1:length(time_as))';
        cumnr_modelf = [];
        if nMod == 1
            for i=1:length(time_as)
                if pval1 ~= 1
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1));
                else
                    cm = kval1*log(time_as(i)/cval1+1);
                end
                cumnr_modelf = [cumnr_modelf; cm];
            end % END of FOR on length(time_as)
        else
            for i=1:length(time_as)
                if time_as(i) <= fT1
                    if pval1 ~= 1
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1));
                    else
                        cm = kval1*log(time_as(i)/cval1+1);
                    end
                    cumnr_modelf = [cumnr_modelf; cm];
                else
                    if (pval1 ~= 1 & pval2 ~= 1)
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_as(i)-fT1+cval2)^(1-pval2));
                    else
                        cm = kval1*log(time_as(i)/cval1+1) + kval2*log((time_as(i)-fT1)/cval2+1);
                    end
                    cumnr_modelf = [cumnr_modelf; cm];
                end; %END of IF on fT1
            end; % End of FOR length(time_as)
        end; % End of if on nMod
        time_as=sort(time_as);
        cumnr_modelf=sort(cumnr_modelf);

        % Calculate KSTEST2 as a measure of the goodness of fit
        [H,P,KSSTAT] = kstest2(cumnr_modelf,cumnrf);

        % Calculate RMS
        i=(1:1:length(time_as))';
        fRMS = (sum((i-cumnr_modelf).^2)/length(i))^0.5;

        % Result matrix
        mResult = [mResult; time pval1 pval2 cval1 cval2 kval1 kval2 fAIC fL H P KSSTAT fRMS];
        mResult2 = [mResult2; pmed1 pstd1 pmed2 pstd2 cmed1 cstd1 cmed2 cstd2 kmed1 kstd1 kmed2 kstd2];
    else
        mResult = [mResult; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
        mResult2 = [mResult2; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ];
    end % ENDIF on length(time_as)
end % FOR time

%save Result_pwtime.mat mResult mResult2;

% figure_w_normalized_uicontrolunits('tag','p-value time series', 'visible','on')
% hPplot4 = plot(mResult(:,1),mResult2(:,1));
% hold on
% hPplot=plot(mResult(:,1),mResult(:,2),'-k','Linewidth',2)
% hPplot2 = plot(mResult(:,1),mResult2(:,1)+mResult2(:,2));
% hPplot3 = plot(mResult(:,1),mResult2(:,1)-mResult2(:,2));
% set(hPplot2,'Linestyle','--','Linewidth',2,'Color',[0.5 0.5 0.5]);
% set(hPplot3,'Linestyle','--','Linewidth',2,'Color',[0.5 0.5 0.5]);
% set(hPplot4,'Linestyle','--','Linewidth',2,'Color',[0.5 0.5 0.5]);
% set(hPplot,'Linewidth',2)
% xlabel('Time after mainshock [days]','Fontweight','bold','FontSize',12)
% ylabel('p-value','Fontweight','bold','FontSize',12)
% set(gca,'Fontweight','bold','FontSize',12,'Box','on','Tickdir','out')

% Smoothing the plot
nWindowSize = 5;
figure_w_normalized_uicontrolunits('tag','p-value time series smooth', 'visible','on')
mPmean = filter(ones(1,nWindowSize)/nWindowSize,1,mResult2(:,1));
mPmean(1:nWindowSize,1)=mResult2(1:nWindowSize,1);
mPstd1 = filter(ones(1,nWindowSize)/nWindowSize,1,mResult2(:,2));
mPstd1(1:nWindowSize,1)=mResult2(1:nWindowSize,2);

mP = filter(ones(1,nWindowSize)/nWindowSize,1,mResult(:,2));
mP(1:nWindowSize,1)=mResult(1:nWindowSize,2);
% Plotting
hp1=plot(mResult(:,1),mPmean,'--','Linewidth',2,'Color',[0.3 0.3 0.3]);
hold on;
hp2=plot(mResult(:,1),mPmean-mPstd1,'-.','Linewidth',2,'Color',[0.5 0.5 0.5]);
plot(mResult(:,1),mPmean+mPstd1,'-.','Linewidth',2,'Color',[0.5 0.5 0.5]);
hPplot=plot(mResult(:,1),mResult(:,2),'-k','Linewidth',2)
ylabel('p-value','Fontweight','bold','FontSize',12)
% xlim([(min(mResult(:,1))) (max(mResult(:,1)))])
% ylim([floor(min(mP(:,1))) ceil(max(mP(:,1)))])
xlabel('Time / [days after main shock]','Fontweight','bold','FontSize',12)
l1=legend([hp1 hp2 hPplot],'p-value (mean)','\sigma p','p-value');
set(l1,'Fontweight','bold')
set(gca,'Fontweight','bold','FontSize',10,'Linewidth',2,'Tickdir','out')

author: CelsoReyes


text:
function F = acscd(F, varargin)
%ACSCD   Inverse cosecant of a CHEBFUN, result in degrees.
%   ACSCD(F) computes the inverse cosecant (in degrees) of the CHEBFUN F.
%
%   ACSCD(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also CSCD, ACSC.

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

% Call the compose method:
F = compose(F, @acscd, varargin{:});

end


text:
function [success,msg] = ConvertMotr2JAABA(varargin)

success = false;
msg = {};  %#ok<NASGU>

% Parse inputs
[inmoviefile,seqindexfile,...
  intrxfile,...
  expdir,moviefilestr,trxfilestr,perframedirstr,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  pxpermm,...
  dofliplr,doflipud,dotransposeimage,...
  dosoftlink,frameinterval,...
  sex,roi2] = myparse(varargin,...
  'inmoviefile','','seqindexfile','',...
  'intrxfile','',...
  'expdir','','moviefilestr','movie.ufmf','trxfilestr','trx.mat','perframedirstr','perframe',...
  'arenatype','None','arenacenterx',0,'arenacentery',0,...
  'arenaradius',123,'arenawidth',123,'arenaheight',123,...
  'pxpermm',1,...
  'fliplr',false,'flipud',false,'dotransposeimage',false,...
  'dosoftlink',false,...
  'frameinterval',[],...
  'sex',{},'roi2',[]);  %#ok<NASGU,ASGLU>

% Check that required inputs are given
if isempty(intrxfile),
  msg = 'Input trx mat file is empty';
  return
end

% Figure out whether it's a raw Motr _tracks.mat file or a Catalytic .ctc file by looking
% for the critical fields
td = load('-mat',intrxfile);
if isfield(td,'astrctTrackers') ,
    [success,msg] = ConvertRawMotr2JAABA(varargin{:});
elseif isfield(td,'trx') ,
    [success,msg] = ConvertMotrCtc2JAABA(varargin{:});
else
    msg='Input trx file doesn''t have astrctTrackers or trx field';
    return
end

end

author: kristinbranson

repo: JAABA


text:
function oeprint1(mu, oev)

% print six classical orbital elements
% and orbital period in minutes

% input

%  mu     = gravitational constant (km**3/sec**2)
%  oev(1) = semimajor axis (kilometers)
%  oev(2) = orbital eccentricity (non-dimensional)
%           (0 <= eccentricity < 1)
%  oev(3) = orbital inclination (radians)
%           (0 <= inclination <= pi)
%  oev(4) = argument of perigee (radians)
%           (0 <= argument of perigee <= 2 pi)
%  oev(5) = right ascension of ascending node (radians)
%           (0 <= raan <= 2 pi)
%  oev(6) = true anomaly (radians)
%           (0 <= true anomaly <= 2 pi)

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rtd = 180 / pi;

% unload orbital elements array

sma = oev(1);
ecc = oev(2);
inc = oev(3);
argper = oev(4);
raan = oev(5);
tanom = oev(6);

arglat = mod(tanom + argper, 2.0 * pi);

if (sma > 0.0)
    period = 2.0d0 * pi * sma * sqrt(sma / mu);
else
    period = 99999.9;
end

% print orbital elements

fprintf ('\n        sma (km)              eccentricity          inclination (deg)         argper (deg)');

fprintf ('\n %+16.14e  %+16.14e  %+16.14e  %+16.14e \n', sma, ecc, inc * rtd, argper * rtd);

fprintf ('\n       raan (deg)          true anomaly (deg)         arglat (deg)            period (min)');

fprintf ('\n %+16.14e  %+16.14e  %+16.14e  %+16.14e \n', raan * rtd, tanom * rtd, arglat * rtd, period / 60);

author: Sable

repo: mcbench-benchmarks


text:
function [D,B]=distance_wei(L)
%DISTANCE_WEI       Distance matrix
%
%   D = distance_wei(L);
%   [D,B] = distance_wei(L);
%
%   The distance matrix contains lengths of shortest paths between all
%   pairs of nodes. An entry (u,v) represents the length of shortest path
%   from node u to node v. The average shortest path length is the
%   characteristic path length of the network.
%
%   Input:      L,      Directed/undirected connection-length matrix.
%   *** NB: The length matrix L isn't the weights matrix W (see below) ***
%
%   Output:     D,      distance (shortest weighted path) matrix
%               B,      number of edges in shortest weighted path matrix
%
%   Notes:
%       The input matrix must be a connection-length matrix, typically
%   obtained via a mapping from weight to length. For instance, in a
%   weighted correlation network higher correlations are more naturally
%   interpreted as shorter distances and the input matrix should
%   consequently be some inverse of the connectivity matrix.
%       The number of edges in shortest weighted paths may in general
%   exceed the number of edges in shortest binary paths (i.e. shortest
%   paths computed on the binarized connectivity matrix), because shortest
%   weighted paths have the minimal weighted distance, but not necessarily
%   the minimal number of edges.
%       Lengths between disconnected nodes are set to Inf.
%       Lengths on the main diagonal are set to 0.
%
%   Algorithm: Dijkstra's algorithm.
%
%
%   Mika Rubinov, UNSW/U Cambridge, 2007-2012.
%   Rick Betzel and Andrea Avena, IU, 2012

%Modification history
%2007: original (MR)
%2009-08-04: min() function vectorized (MR)
%2012: added number of edges in shortest path as additional output (RB/AA)
%2013: variable names changed for consistency with other functions (MR)

n=length(L);
D=inf(n);
D(1:n+1:end)=0;                             %distance matrix
B=zeros(n);                                 %number of edges matrix

for u=1:n
    S=true(1,n);                            %distance permanence (true is temporary)
    L1=L;
    V=u;
    while 1
        S(V)=0;                             %distance u->V is now permanent
        L1(:,V)=0;                          %no in-edges as already shortest
        for v=V
            T=find(L1(v,:));                %neighbours of shortest nodes
            [d,wi]=min([D(u,T);D(u,v)+L1(v,T)]);
            D(u,T)=d;                       %smallest of old/new path lengths
            ind=T(wi==2);                   %indices of lengthened paths
            B(u,ind)=B(u,v)+1;              %increment no. of edges in lengthened paths
        end

        minD=min(D(u,S));
        if isempty(minD)||isinf(minD),      %isempty: all nodes reached;
            break,                          %isinf: some nodes cannot be reached
        end;

        V=find(D(u,:)==minD);
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function idx = selectRef(I, winSize, exposureThres)

if ~exist('winSize', 'var')
   winSize = 3;
end

if ~exist('exposureThres', 'var')
   exposureThres = 0.01;
end

I = double(I);
I = reorderByLum(I);
[~, ~, s3, s4] = size(I);

if s4 == 3
    idx = 2;
else
    window = ones(winSize, winSize, s3);
    window = window / sum(window(:));
    p = zeros(s4,1);
    for i = 1 : s4
        mu = convn(I(:, :, :, i), window, 'valid');
        p(i) = sum( sum( mu < exposureThres | mu > 1 - exposureThres) );
    end
    [~, idx] = min(p);
end

author: thfylsty


text:
function [imc] = undistort_image(im, KK, kc, nx, ny, KK_new)

fc = [ KK(1, 1); KK(2, 2) ];
alpha_c = 0;
cc = [ KK(1,3); KK(2, 3) ];

dist_amount = 1; %(1+kc(1)*r2_extreme + kc(2)*r2_extreme^2);
fc_new = dist_amount * fc;

if ~exist('KK_new','var')
    KK_new = [fc_new(1) alpha_c*fc_new(1) cc(1);0 fc_new(2) cc(2) ; 0 0 1];
end

I = im;

if (size(I,1)>ny)|(size(I,2)>nx),
    I = I(1:ny,1:nx);
end;

%% SHOW THE ORIGINAL IMAGE:
%{
if size(I,3) == 3
    figure(2);
    image(uint8(I));
    title('Original image (with distortion) - Stored in array I')
    axis equal;
    axis tight;
    drawnow;
else
    figure(2);
    image(I);
    colormap(gray(256));
    axis equal;
    axis tight;
    title('Original image (with distortion) - Stored in array I');
    drawnow;
end
%}

%% UNDISTORT THE IMAGE:

fprintf(1,'Computing the undistorted image...')

if size(I,3) == 3
    [Ipart_1] = rect(I(:,:,1),eye(3),fc,cc,kc,alpha_c,KK_new);
    [Ipart_2] = rect(I(:,:,2),eye(3),fc,cc,kc,alpha_c,KK_new);
    [Ipart_3] = rect(I(:,:,3),eye(3),fc,cc,kc,alpha_c,KK_new);

    I2 = ones(ny, nx,3);
    I2(:,:,1) = Ipart_1;
    I2(:,:,2) = Ipart_2;
    I2(:,:,3) = Ipart_3;

    fprintf(1,'done\n')

    %{
    figure(3);
    image(uint8(I2));
    axis equal;
    axis tight;
    %}
else
    [I2] = rect(I,eye(3),fc,cc,kc,alpha_c,KK_new);

    fprintf(1,'done\n');

    %{
    figure(3);
    image(I2);
    colormap(gray(256));
    axis equal;
    axis tight;
    %}
end;

%title('Undistorted image - Stored in array I2')
%drawnow;

imc = I2;

%

author: jianxiongxiao

repo: ProfXkit


text:
function bwn = bw_filter(bw, keepnum)
if nargin < 2
    keepnum = 15;
end
[L, num] = bwlabel(bw, 8);
Ln = zeros(1, num);
stats = regionprops(L, 'Area');
Ln = cat(1, stats.Area);
[Ln, ind] = sort(Ln);
if num>keepnum || num==keepnum
    for i = 1 : num-keepnum
        bw(L == ind(i)) = 0;
    end
end
bwn = bw;

author: decouples


text:
function [ w, x ] = line_o05 ( )

%*****************************************************************************80
%
%% LINE_O05 returns a 5 point quadrature rule for the unit line.
%
%  Discussion:
%
%    The integration region is:
%
%    - 1.0 <= X <= 1.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(5), the weights.
%
%    Output, real X(5), the abscissas.
%
  w(1:5) = [ ...
    0.118463442528095, ...
    0.239314335249683, ...
    0.284444444444444, ...
    0.239314335249683, ...
    0.118463442528095 ];

  x(1:5) = [ ...
    -0.90617984593866399280, ...
    -0.53846931010568309104, ...
     0.00000000000000000000, ...
     0.53846931010568309104, ...
     0.90617984593866399280 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gammaAB = diffCoherence(k, r_A, r_B, a_nm, b_nm, G_mtx)
%DIFFCOHERENCE Computes diffuse field coherence of two directional patterns
%
%   Does that. For more details see
%
%       Politis, A., 2016.
%       Diffuse-field coherence of sensors with arbitrary directional responses. arXiv preprint arXiv:1608.07713.
%
%   Inputs:
%       k:  vector of wavenumbers to compute coherence at
%       r_A:    [x y z] vector of location of first sensor/beamformer
%       r_B:    [x y z] vector of location of second sensor/beamformer
%       a_nm:   SHD beamweights of the first directional pattern
%       b_nm:   SHD beamweights of the second directional pattern
%       G_mtx:  (optional) pre-computed matrix of Gaunt coefficients used
%           in the derivation. If not passed then they are computed here,
%           but that can be very slow for high orders.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DIFFCOHERENCE.M - 7/06/2015
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get the vector connecting the two patterns
r_AB = r_B-r_A;
[phi_r, elev_r, d] = cart2sph(r_AB(1),r_AB(2),r_AB(3));
theta_r = pi/2-elev_r;

orderA = sqrt(length(a_nm)) -1;
orderB = sqrt(length(b_nm)) -1;
orderC = orderA+orderB;
minOrder = min(orderA, orderB);

% convert second pattern to its conjugate for coherence between complex
% patterns (doesn't do anything if the patterns are real)
b_nm = conjCoeffs(b_nm);

% denominator of the coherence function
denAB = sqrt((a_nm'*a_nm) * (b_nm'*b_nm));

if d==0
    % numerator of the coherence function
    numAB = dot(a_nm(1:(minOrder+1)^2), b_nm(1:(minOrder+1)^2));
    numAB = numAB*ones(size(k));

else
    % get coefficients of the product function
    if nargin<6
        c_nm = sphMultiplication(a_nm, b_nm);
    else
        c_nm = sphMultiplication(a_nm, b_nm, G_mtx);
    end

    % numerator of the coherence function
    y_nm_r = getSH(orderC, [phi_r theta_r], 'complex')';
    kd = k*d;
    imag_n = 1i.^(0:orderC).';
    bessel_n = zeros(orderC+1,length(kd));
    for no=0:orderC, bessel_n(no+1,:) = sph_besselj(no, kd).'; end
    imag_bessel_n = (imag_n*ones(1,length(kd))).*bessel_n;
    imag_bessel_nm = replicatePerOrder(imag_bessel_n);
    planewave_nm = 4*pi*(y_nm_r*ones(1,length(kd))).*imag_bessel_nm;
    numAB = zeros(size(k));
    for kk=1:length(kd)
        numAB(kk) = planewave_nm(:,kk)' * c_nm;
    end
end

% coherence
gammaAB = numAB/denAB;

end

author: polarch

repo: Spherical-Array-Processing


text:
function [one_in_K,mapping] = ml_labels2Kencoding( labels )
%ML_LABELS2KENCODING Takes class label data and transforms it to one in K
% encoding.
%
%   input -----------------------------------------------------------------
%
%       o lables : (N x 1), class labels.
%
%   output ----------------------------------------------------------------
%
%       o one_in_K : (N x K), one-in-K encoding of classes
%
%       o mapping  : (num_classes x 1)
%
%

labels      = labels(:);
N           = size(labels,1);
class_label = unique(labels);
num_classes = length(class_label);
mapping     = zeros(num_classes,1);
one_in_K    = zeros(N,num_classes);

for i=1:num_classes
    mapping(i)                           = class_label(i);
    one_in_K(class_label(i) == labels,i) = 1;
end

end

author: epfl-lasa


text:
% gauss2d() - generate a 2-dimensional gaussian matrix
%
% Usage:
%   >> [ gaussmatrix ] = gauss2d( rows, columns, ...
%                              sigmaR, sigmaC, peakR, peakC, mask)
%
% Example:
%   >> imagesc(gauss2d(50, 50)); % image a size (50,50) 2-D gaussian matrix
%
% Inputs:
%   rows    - number of rows in matrix
%   columns - number of columns in matrix
%   sigmaR  - width of standard deviation in rows (default: rows/5)
%   sigmaC  - width of standard deviation in columns (default: columns/5)
%   peakR   - location of the peak in each row (default: rows/2)
%   peakC   - location of the peak in each column (default: columns/2)
%   mask    - (0->1) portion of the matrix to mask with zeros (default: 0)
%
% Ouput:
%   gaussmatrix - 2-D gaussian matrix
%
% Author: Arnaud Delorme, CNL/Salk Institute, 2001

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

function mat = gauss2d( sizeX, sizeY, sigmaX, sigmaY, meanX, meanY, cut);

if nargin < 2
	help gauss2d
	return;
end;
if nargin < 3
	sigmaX = sizeX/5;
end;
if nargin < 4
	sigmaY = sizeY/5;
end;
if nargin < 5
	meanX = (sizeX+1)/2;
end;
if nargin < 6
	meanY = (sizeY+1)/2;
end;
if nargin < 7
	cut = 0;
end;

X = linspace(1, sizeX, sizeX)'* ones(1,sizeY);
Y = ones(1,sizeX)'   		  * linspace(1, sizeY, sizeY);
%[-sizeX/2:sizeX/2]'*ones(1,sizeX+1);
%Y = ones(1,sizeY+1)'   *[-sizeY/2:sizeY/2];

mat = exp(-0.5*(  ((X-meanX)/sigmaX).*((X-meanX)/sigmaX)...
				+((Y-meanY)/sigmaY).*((Y-meanY)/sigmaY)))...
            			/((sigmaX*sigmaY)^(0.5)*pi);

if cut > 0
	maximun = max(max(mat))*cut;
	I = find(mat < maximun);
	mat(I) = 0;
end;

return;

author: PatternRecognition

repo: OpenBMI


text:
function [cbar] = cbarEdit(cbar);
%
% [cbar] = cbarEdit([cbar]);
%
% Edit a mrVista color bar using a GUI.
%
% These color bars differ from the built-in matlab color bar tools in a few
% basic respects:
%
%   * They draw colorbars based on any colormap the user inputs, rather
%   than the color map assigned to the current figure. This is useful for
%   illustrating color codes on true color images, or images with overlaid
%   responses (like fMRI activation maps)
%
%   * In addition to being able to render vertical and horizontal color bars,
%   these tools may also render 'color wheel' colorbars, showing e.g. polar
%   angle. Other possible renderings (like rings) may be added down the
%   line.
%
% cbar is a struct with the following fields:
%
%   cmap: color map (nColors x 3) for the color bar. (Columns ar [R G B],
%   from 0-255).
%
%   nColors: # of colors to use in the cmap.
%
%   clim: color limits (aka 'clip mode'), which determines primarily
%   the labeling of the color bar. Can be set to 'auto', in which case
%   the labeling will be from 1:nColors. Otherwise, will label according to
%   the clim values (format is [min max]).
%
%   colorWheel: use a color wheel instead of a bar (e.g., to show polar
%   angle for a polar angle map).
%
%   colorWheelStart: degrees clockwise from 12-o-clock which map to beginning of
%   color map.
%
%   colorWheelDirection: direction of the color wheel. Can be:
%   1 or 'clockwise' (equivalent); or, 2 or 'counterclockwise' (equiv.)
%
%   colorWheelExtent: degrees (1 - 360) subtended by the color map, for polar
%   angle maps.
%
% ras, 08/06
if notDefined('cbar')
    cbar = cbarDefault;
end

if ischar(cbar) & isequal(lower(cbar), 'update')
    cbarEditUpdate; return;
end

% Create the GUI, attaching cmap and cbar
hfig = cbarEditGUI(cbar);

% Wait for user to respond  via GUI
uiwait;

% check if user canceled
OK = evalin('base', 'OK');
evalin('base', 'clear OK');
if OK==0
    warning('User Canceled.')
    delete(hfig);
    return
end

cbar = get(hfig, 'UserData');
delete(hfig);

cbar = rmfield(cbar, 'ui');

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function cbar = cbarEditUpdate;
% update the color map editor GUI, based on the settings in the current
% figure.
cbar = get(gcf, 'UserData');

% modify the cmap according to any callback objects
if ishandle(gcbo)
    switch get(gcbo, 'Tag')
        case 'Shift'            % shift according to slider
            delta = round( get(gcbo, 'Value') );
            cbar.cmap = circshift(cbar.cmap, delta);
            set(gcbo, 'Value', 0);   % reset slider

        case 'Flip'             % flip cmap order
            cbar.cmap = flipud(cbar.cmap);

        case 'Reset'            % reset to original cbar
            cbar = get(gcbo, 'UserData'); % stashed here

        case 'NumColors'        % set # colors in cmap
            cbar.nColors = str2num(get(gcbo, 'String'));
            while size(cbar.cmap, 1) < cbar.nColors
                cbar.cmap = [cbar.cmap; cbar.cmap];
            end
            cbar.cmap = cbar.cmap(1:cbar.nColors,:);

        case 'Preset'           % use a preset color map
            cbar.cmap = mrvColorMaps(get(gcbo, 'Label'), cbar.nColors);

        case 'ColorWheel'       % set to color wheel
            cbar.colorWheel = get(gcbo, 'Value');
            vals = {'off' 'on'};
            set(cbar.ui.colorWheelPanel, 'Visible', vals{cbar.colorWheel+1});

        case 'ColorWheelDirection' % set color wheel CW / CCW
            options = {'clockwise' 'counterclockwise'};
            cbar.colorWheelDirection = options{ get(gcbo, 'Value') };

        case 'ColorWheelStart'  % set start angle for color wheel
            cbar.colorWheelStart = str2num(get(gcbo, 'String'));

        case 'ColorWheelExtent'  % set angle subtended by cmap
			% ras 02/07: new strategy; will set things so that the wedge
			% always goes a full 360 degrees, but will pad the cmap such
			% that it only occupies [extent] degrees. This should help it
			% to work properly for all extents in [0 360].
            newExtent = str2num(get(gcbo, 'String'));

			% bounds check
			newExtent = max(1, min(360, newExtent));
			set(gcbo, 'String', num2str(newExtent));

			if newExtent ~= cbar.colorWheelExtent
				% undo any previous gray padding
				cbar.cmap = cbar.cmap(1:cbar.nColors,:);
			end

			if cbar.colorWheelExtent < 360
				prop = cbar.colorWheelExtent / 360; % proprtion of cbar filled by cmap
				totalColors = ceil(cbar.nColors / prop);
				grayPad = repmat(.5, [(totalColors - cbar.nColors) 3]);
				cbar.cmap = [cbar.cmap; grayPad];
			end

			cbar.colorWheelExtent = newExtent;

        case 'Expression'       % user inputs expression for color map
			map = cbar.cmap;
			msg = ['Enter a MATLAB expression for the color map. ' ...
				'The color map should be a matrix of 3 x ncolors, ' ...
				'with values from 0 - 255. ' ...
				'You can directly enter the expression, or assign ' ...
				'it to a variable named "map". ' ...
				'E.g.:    [gray(128); hot(128)]; ' ...
				'or:      tmp = hsv(256);  map = tmp(1:100,:); '];
			def = sprintf('gray(%i); ', cbar.nColors);
			resp = inputdlg({msg}, mfilename, 4, {def});
			if strfind(resp{1}, '=') & strfind(resp{1}, 'map')
				% presume that the expression is assigning a value to the
				% 'map' variable:
				eval(resp{1})
			else
				% the map assignment is implied
				map = eval(resp{1});
			end
			cbar.cmap = map;

    end
end

% draw the color bar
cbarDraw(cbar);
title('Edit Color Map');

% update GUI figure
set(gcf, 'UserData', cbar);

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function hfig = cbarEditGUI(cbar);
% creates the figure for cbarEdit.

%%%%% create figure, main axes
hfig = figure('Color', 'w', 'Name', 'Edit Color Map', 'MenuBar', 'none');
hax = subplot('Position', [.1 .7 .8 .16]);

%%%%% callback for all uicontrols
cb = 'cbarEdit(''update''); ';

%%%%%main slider
nC = cbar.nColors;
hs = uicontrol('Style', 'slider', 'Min', -nC/2, 'Max', nC/2, 'Value', 0, ...
   'BackgroundColor', 'w', 'ForegroundColor', [.6 .55 .6], 'Tag', 'Shift', ...
    'Units', 'normalized', 'Position', [.1 .55 .6 .08], 'Callback', cb);

% text label
uicontrol('Style', 'text', 'String', 'Amount of Rotation', ...
    'BackgroundColor', [1 1 1], 'FontSize', 12, ...
    'Units', 'normalized', 'Position', [.1 .47 .8 .08]);

% flip button
uicontrol('Style', 'pushbutton', 'String', 'Flip', 'Callback', cb, ...
    'Tag', 'Flip', ...
    'Units', 'normalized', 'Position', [.75 .55 .2 .05]);

% reset button
uicontrol('Style', 'pushbutton', 'String', 'Reset', 'Callback', cb, ...
    'Tag', 'Reset', ...
    'UserData', cbar, 'Units', 'normalized', 'Position', [.75 .5 .2 .05]);

% edit field for setting # of colors
uicontrol('Style', 'text', 'String', '# Colors in Cmap', 'Callback', cb, ...
    'BackgroundColor', 'w', 'FontSize', 10, 'FontWeight', 'bold', ...
    'UserData', cbar, 'Units', 'normalized', 'Position', [.75 .5 .2 .05]);

uicontrol('Style', 'edit', 'String', num2str(cbar.nColors), 'Callback', cb, ...
    'Tag', 'NumColors', 'Units', 'normalized', 'Position', [.75 .45 .15 .05]);

%%%%%% wedge cbar controls
% checkbox to set the colorbar as a color wheel (wedge)
uicontrol('Style', 'checkbox', 'String', 'Color Wheel (Wedge)', ...
          'Units', 'normalized', 'Position', [.1 .3 .4 .08], ...
          'BackgroundColor', 'w', 'FontSize', 10, 'FontWeight', 'bold', ...
          'Value', cbar.colorWheel, ...
          'Tag', 'ColorWheel', 'Callback', cb);

% put the subsidiary controls in a uipanel
hp = uipanel('Units', 'normalized', 'Position', [.4 .1 .5 .3], ...
             'BorderType', 'none', 'BackgroundColor', 'w');
cbar.ui.colorWheelPanel = hp;

% popup to set the colorbar direction
uicontrol('Parent', hp, 'Style', 'text', 'String', 'Direction', 'Callback', cb, ...
         'BackgroundColor', 'w', 'Units', 'normalized', 'Position', [.1 .85 .8 .1]);

uicontrol('Parent', hp, 'Style', 'popup', ...
          'String', {'clockwise' 'counterclockwise'}, ...
          'BackgroundColor', 'w', ...
          'Units', 'normalized', 'Position', [.1 .6 .8 .25], ...
          'Tag', 'ColorWheelDirection', 'Callback', cb);

% edit field to set the color wheel start angle
uicontrol('Parent', hp, 'Style', 'text', 'String', 'Start Angle', 'Callback', cb, ...
          'BackgroundColor', 'w', 'Units', 'normalized', 'Position', [.1 .35 .3 .25]);

uicontrol('Parent', hp, 'Style', 'edit', ...
          'String', num2str(cbar.colorWheelStart), ...
          'BackgroundColor', 'w', ...
          'Units', 'normalized', 'Position', [.6 .35 .3 .25], ...
          'Tag', 'ColorWheelStart', 'Callback', cb);

% edit field to set the color wheel extent
uicontrol('Parent', hp, 'Style', 'text', 'String', 'Extent', 'Callback', cb, ...
          'BackgroundColor', 'w', 'Units', 'normalized', 'Position', [.1 .05 .4 .25]);

uicontrol('Parent', hp, 'Style', 'edit', ...
          'String', num2str(cbar.colorWheelExtent), ...
          'Units', 'normalized', 'Position', [.6 .05 .4 .25], ...
          'BackgroundColor', 'w', ...
          'Tag', 'ColorWheelExtent', 'Callback', cb);

onoff = {'off' 'on'};
set(hp, 'Visible', onoff{cbar.colorWheel+1});

%%%%% also add menus to set preset color bars.
mapNames = mrvColorMaps;
mapNames = mapNames(1:end-1); % remove the 'user' option
hm = uimenu('Label','Preset Colormaps','ForegroundColor','b');
for i = 1:length(mapNames)
    uimenu(hm, 'Label', mapNames{i}, 'Tag', 'Preset', 'Callback', cb);
end

% add a callback to let the user input a MATLAB expression for the cbar
uimenu(hm, 'Label', 'Set color map using MATLAB expression...', ...
           'Tag', 'Expression', 'Callback', cb);

addFigMenuToggle(hm);

%%%%%Accept and Cancel buttons
evalin('base', 'OK=0;');
cb = 'OK = 1; close(gcf); ';
uicontrol('Style', 'pushbutton', 'String', 'Accept', 'Callback', cb, ...
    'BackgroundColor', [.1 .5 .1], 'ForegroundColor', [1 1 1], ...
    'Units', 'normalized', 'Position', [.1 .04 .3 .06]);

cb = 'OK = 0; delete(gcf); uiresume;';
uicontrol('Style', 'pushbutton', 'String', 'Cancel', 'Callback', cb, ...
    'BackgroundColor', [.7 .1 .1], 'ForegroundColor', [1 1 1], ...
    'Units', 'normalized', 'Position', [.5 .04 .3 .06]);

set(hfig, 'CloseRequestFcn', 'uiresume;', 'UserData', cbar);

cbarEditUpdate;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function eval_result=eval_holisticScencePR(predictedBbs,groundTruthBbs,thresh_iou,cls)
    % input :
    % cls : object categorys that in consideration
    % thresh_iou: iou threshold for detection with ground truth
    % ouput :
    % label_precision: measure the recall of recognition for both semantics and geometry
    % iou_recall : measure the geometric prediction recall
    % iou_precision : measure the geometric prediction precision
    if ~exist('thresh_iou','var')
        thresh_iou = [0.01,0.05:0.05:0.6,0.8:0.2:1];
    end
    if ~exist('cls','var')
        cls = {'bathtub','bed','bookshelf','box','chair','counter','desk','door','dresser','garbage_bin','lamp','monitor','night_stand','pillow','sink','sofa','table','tv','toilet'};
    end

    allOverlaps = bb3dOverlapCloseForm(predictedBbs,groundTruthBbs);
    matchpg =[];
    maxOverlaps =[];
    while sum(allOverlaps(:)>0)>0
        [os,ind] = max(allOverlaps(:));
        [i,j] = ind2sub(size(allOverlaps),ind);
        matchpg(end+1,:) =[i,j];
        maxOverlaps(end+1) = os;
        allOverlaps(i,:) =0;
        allOverlaps(:,j) =0;
    end

    iou_recall = zeros(1,length(thresh_iou));
    iou_precision = zeros(1,length(thresh_iou));
    label_precision = zeros(1,length(thresh_iou));
    label_recall = zeros(1,length(thresh_iou));
    if ~isempty(matchpg)
        for i = 1:length(thresh_iou)
            iou_match = find(maxOverlaps>thresh_iou(i));
            pclass = [predictedBbs(matchpg(iou_match,1)).classid];
            [~,gtclass] = ismember({groundTruthBbs(matchpg(iou_match,2)).classname},cls);
            if length(pclass) > 0
                pclass = double(pclass);
            end
            correctLabel = sum([pclass-gtclass]==0);
            label_precision(i) = correctLabel/length(predictedBbs);
            label_recall(i) = correctLabel/length(groundTruthBbs);
            iou_recall(i) = length(iou_match)/length(groundTruthBbs);
            iou_precision(i) = length(iou_match)/length(predictedBbs);
        end
    end
    eval_result = struct('label_precision',label_precision,'label_recall',label_recall,'iou_recall',iou_recall,'iou_precision',iou_precision);
end

author: thusiyuan


text:
function [mri, seg, hdr] = read_asa_mri(fn)

% READ_ASA_MRI reads an ASA format MRI file
%
% Use as
%   [mri, seg, hdr] = read_asa_mri(filename)
%
% The raw image data is returned, together with the position of the
% external head markers in raw image coordinates.
%
% In the ASA default PAN (pre-auricular/nasion) coordinate system
%   PointOnPositiveYAxis -> LPA
%   PointOnNegativeYAxis -> RPA
%   PointOnPositiveXAxis -> nasion

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

hdr.Nrows     = read_ini(fn, 'NumberRows=', '%d');
hdr.Ncolumns  = read_ini(fn, 'NumberColumns=', '%d');
hdr.Nslices   = read_ini(fn, 'NumberSlices=', '%d');
hdr.rows      = read_ini(fn, 'Rows=', '%s');
hdr.columns   = read_ini(fn, 'Columns=', '%s');
hdr.slices    = read_ini(fn, 'Slices=', '%s');
hdr.distance  = read_ini(fn, 'Distance=', '%f', 3);
hdr.mrifile   = read_ini(fn, 'Matrix', '%s');
hdr.segfile   = read_ini(fn, 'Segmentation', '%s');
hdr.posy      = read_ini(fn, 'PointOnPositiveYAxis', '%f', 3);
hdr.negy      = read_ini(fn, 'PointOnNegativeYAxis', '%f', 3);
hdr.posx      = read_ini(fn, 'PointOnPositiveXAxis', '%f', 3);
hdr.voxposy   = read_ini(fn, 'VoxelOnPositiveYAxis', '%f', 3);
hdr.voxnegy   = read_ini(fn, 'VoxelOnNegativeYAxis', '%f', 3);
hdr.voxposx   = read_ini(fn, 'VoxelOnPositiveXAxis', '%f', 3);
hdr.segfile   = read_ini(fn, 'Segmentation', '%s', 1);

dim = [hdr.Ncolumns hdr.Nrows hdr.Nslices];
mri = [];
seg = [];

% the data files are at the same location as the mri file, locate path
[path, name, ext] = fileparts(fn);

% this temporary needs 8x as much storage!!!
if ~isempty(hdr.mrifile)
  mrifile = fullfile(path, hdr.mrifile);
  mri = zeros(dim);
  fid = fopen_or_error(mrifile, 'rb', 'ieee-le');
  mri(:) = fread(fid, prod(dim), 'uint8');
  mri = uint8(mri);
  fclose(fid);
end

% this temporary needs 8x as much storage!!!
if ~isempty(hdr.segfile)
  segfile = fullfile(path, hdr.segfile);
  seg = zeros(dim);
  fid = fopen_or_error(segfile, 'rb', 'ieee-le');
  seg(:) = fread(fid, prod(dim), 'uint8');
  seg = uint8(seg);
  fclose(fid);
end

% flip the orientation of the MRI data, the result should be
% 'coronal(occipital-frontal)'
% 'horizontal(inferior-superior)'
% 'sagittal(right-left)'

if strcmp(hdr.columns, 'coronal(frontal-occipital)')
  hdr.columns = 'coronal(occipital-frontal)';
  mri = flipdim(mri, 1);
  seg = flipdim(seg, 1);
elseif strcmp(hdr.columns, 'horizontal(superior-inferior)')
  hdr.columns = 'horizontal(inferior-superior)';
  mri = flipdim(mri, 1);
  seg = flipdim(seg, 1);
elseif strcmp(hdr.columns, 'sagittal(left-right)')
  hdr.columns = 'sagittal(right-left)';
  mri = flipdim(mri, 1);
  seg = flipdim(seg, 1);
end

if strcmp(hdr.rows, 'coronal(frontal-occipital)')
  hdr.rows = 'coronal(occipital-frontal)';
  mri = flipdim(mri, 2);
  seg = flipdim(seg, 2);
elseif strcmp(hdr.rows, 'horizontal(superior-inferior)')
  hdr.rows = 'horizontal(inferior-superior)';
  mri = flipdim(mri, 2);
  seg = flipdim(seg, 2);
elseif strcmp(hdr.rows, 'sagittal(left-right)')
  hdr.rows = 'sagittal(right-left)';
  mri = flipdim(mri, 2);
  seg = flipdim(seg, 2);
end

if strcmp(hdr.slices, 'coronal(frontal-occipital)')
  hdr.slices = 'coronal(occipital-frontal)';
  mri = flipdim(mri, 3);
  seg = flipdim(seg, 3);
elseif strcmp(hdr.slices, 'horizontal(superior-inferior)')
  hdr.slices = 'horizontal(inferior-superior)';
  mri = flipdim(mri, 3);
  seg = flipdim(seg, 3);
elseif strcmp(hdr.slices, 'sagittal(left-right)')
  hdr.slices = 'sagittal(right-left)';
  mri = flipdim(mri, 3);
  seg = flipdim(seg, 3);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% swap the orientations of the MRI data, the result should be fixed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1st dimension corresponds to columns, which should be 'coronal(occipital-frontal)'
% 2st dimension corresponds to rows,    which should be 'sagittal(right-left)'
% 3rd dimension corresponds to slices,  which should be 'horizontal(inferior-superior)'

if     strcmp(hdr.columns, 'coronal(occipital-frontal)')
  orientation(1) = 1;
elseif strcmp(hdr.columns, 'sagittal(right-left)')
  orientation(1) = 2;
elseif strcmp(hdr.columns, 'horizontal(inferior-superior)')
  orientation(1) = 3;
end

if     strcmp(hdr.rows, 'coronal(occipital-frontal)')
  orientation(2) = 1;
elseif strcmp(hdr.rows, 'sagittal(right-left)')
  orientation(2) = 2;
elseif strcmp(hdr.rows, 'horizontal(inferior-superior)')
  orientation(2) = 3;
end

if     strcmp(hdr.slices, 'coronal(occipital-frontal)')
  orientation(3) = 1;
elseif strcmp(hdr.slices, 'sagittal(right-left)')
  orientation(3) = 2;
elseif strcmp(hdr.slices, 'horizontal(inferior-superior)')
  orientation(3) = 3;
end

mri = ipermute(mri, orientation);
seg = ipermute(seg, orientation);
hdr.rows    = 'sagittal(right-left)';
hdr.columns = 'coronal(occipital-frontal)';
hdr.slices  = 'horizontal(inferior-superior)';

% recompute the dimensions after all the swapping
hdr.Nrows    = size(mri, 1);
hdr.Ncolumns = size(mri, 2);
hdr.Nslices  = size(mri, 3);
dim = size(mri);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if possible, create the accompanying homogenous coordinate transformation matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% In case of PointOn..., ASA counts voxels from the center of the MRI
% and in case of VoxelOn..., ASA counts voxels from the corner of the MRI
% In both cases, ASA starts counting at [0 0 0], which is C convention
% whereas I want to count from the 1st voxel and number that with [1 1 1]
if ~isempty(hdr.posx) && ~isempty(hdr.negy) && ~isempty(hdr.posy)
  offset = (dim + [1 1 1])/2;
  hdr.fiducial.mri.nas = hdr.posx + offset;
  hdr.fiducial.mri.lpa = hdr.posy + offset;
  hdr.fiducial.mri.rpa = hdr.negy + offset;
else
  offset = [1 1 1];
  hdr.fiducial.mri.nas = hdr.voxposx + offset;
  hdr.fiducial.mri.lpa = hdr.voxposy + offset;
  hdr.fiducial.mri.rpa = hdr.voxnegy + offset;
end

% use the headcoordinates function (roboos/misc) to compute the transformaton matrix
hdr.transformMRI2Head = ft_headcoordinates(hdr.fiducial.mri.nas, hdr.fiducial.mri.lpa, hdr.fiducial.mri.rpa, 'asa');
hdr.transformHead2MRI = inv(hdr.transformMRI2Head);

% compute the fiducials in head coordinates
hdr.fiducial.head.nas = ft_warp_apply(hdr.transformMRI2Head, hdr.fiducial.mri.nas, 'homogenous');
hdr.fiducial.head.lpa = ft_warp_apply(hdr.transformMRI2Head, hdr.fiducial.mri.lpa, 'homogenous');
hdr.fiducial.head.rpa = ft_warp_apply(hdr.transformMRI2Head, hdr.fiducial.mri.rpa, 'homogenous');


text:
function savemphtxt(node, face, elem, filename)
%
% savemphtxt(node, face, elem, filename)
%
% save tetrahedron mesh to comsol file (.mphtxt)
%
% author: Donghyeon Kim (danielkim<at> gist.ac.kr)
% date: 2011/09/29
%
% input:
%      node: input, node list, dimension (nn,3)
%      face: input, surface face element list with label, dimension (be,4)
%      elem: input, tetrahedron element list with label, dimension (ne,5)
%      filename: input, output file name
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

n_node = size(node,1);
n_face = size(face,1);
n_elem = size(elem,1);

elem(:,1:4)=meshreorient(node(:,1:3),elem(:,1:4));

if(size(face,2)<4)
    face(:,4)=1;
elseif(min(face(:,4))==0)
    face(:,4)=face(:,4)+1;
end

if(size(elem,2)<5)
    elem(:,5)=1;
elseif(min(elem(:,5))==0)
    elem(:,5)=elem(:,5)+1;
end

fp = fopen(filename,'w');
fprintf(fp,'# Created by iso2mesh (http://iso2mesh.sf.net)\n');
fprintf(fp,'0 1\n1\n5 mesh1\n1\n3 obj\n\n');
fprintf(fp,'0 0 1\n4 Mesh\n2\n3\n%d\n1\n',n_node);

% Write Node information
fprintf(fp,'%.16f %.16f %.16f\n',node(:,1:3).');
fprintf(fp, '\n2\n\n3 tri\n');

% Write Tri information
fprintf(fp, '\n3\n');
fprintf(fp, '%d\n\n',n_face);
fprintf(fp,'%d %d %d\n',face(:,1:3).');

fprintf(fp, '\n1\n0\n');
fprintf(fp, '%d\n',n_face);
fprintf(fp,'%d\n',face(:,4).');

fprintf(fp, '\n%d\n',n_face);
fprintf(fp,'%d %d\n', zeros(0,n_face));

% Write Tet information
fprintf(fp, '\n\n3 tet\n4\n\n%d\n', n_elem);
fprintf(fp,'%d %d %d %d\n',elem(:,1:4).');

fprintf(fp, '\n4\n0\n%d\n%d\n', n_elem);
fprintf(fp,'%d\n',elem(:,5).');

fprintf(fp,'\n0\n');
fclose(fp);

author: fangq


text:
function acoustic_feat_ex( data_dir, save_dir )

rng(0);
%% Directory setting

% system(['rm -rf ', save_dir]);
%
% system(['mkdir ', save_dir]);
system(['mkdir ', save_dir, '/Normalize_Factor']);
system(['mkdir ', save_dir, '/Labels']);

%% Parameter setting

audio_sr = 16000;
split_num = 1;
name_mrcg = [save_dir, '/se_mrcg'];
name_label = [save_dir, '/se_label'];

audio_list = getAllFiles(data_dir, 'FileFilter', '\.wav$');
label_list = getAllFiles(data_dir, 'FileFilter', '\.mat$');

winlen             = ceil(audio_sr*25*0.001);	%window length (default : 25 ms)
winstep            = ceil(audio_sr*10*0.001);	%window step (default : 10 ms)

train_mean = 0;
train_std = 0;

for i = 1:1:length(audio_list)
    clc
    fprintf("MRCG extraction %d/%d ...\n", i, length(audio_list));
    %% Read audio

    noisy_speech = audioread(audio_list{i});  % noisy_speech load
    noisy_speech = noisy_speech(1:(length(noisy_speech)-mod(length(noisy_speech), split_num)));
    noisy_speech = reshape(noisy_speech, [], split_num);

    %% Caliculate MRCG
    mrcg = cell(split_num, 1);

    for j = 1:1:split_num
        mrcg{j, 1} = MRCG_features(noisy_speech(:, j), audio_sr)';
        %     imagesc(s(20000:20500,:)*1000)
    end

    mrcg_mat = cell2mat(mrcg);

    size(mrcg_mat)
    %% Save normalization factor

    temp_mean = mean(mrcg_mat,1);
    temp_std = std(mrcg_mat,1,1);
    save([save_dir, '/Normalize_Factor/normalize_factor_', sprintf('%3.3d', i)],'temp_mean', 'temp_std');
    train_mean = temp_mean + train_mean;
    train_std = temp_std + train_std;

    %% Read label
    label = cell2mat(struct2cell(load(label_list{i})));  % label load

    %% Save framed label & MRCG
    framed_label = Truelabel2Trueframe( label, winlen, winstep );
    length(framed_label)
    if (length(mrcg_mat) > length(framed_label))
        binary_saver( name_mrcg, mrcg_mat(1:length(framed_label), :), i );
        binary_saver( name_label, framed_label, i );
    else
        binary_saver( name_mrcg, mrcg_mat, i );
        binary_saver( name_label, framed_label(1:length(mrcg_mat), 1), i );
    end
end

disp('MRCG extraction done.')
%% Save global normalization factor

global_mean = train_mean / length(audio_list);
global_std = train_std / length(audio_list);
save([save_dir, '/global_normalize_factor'], 'global_mean', 'global_std');

%% Move label data

feat_list = getAllFiles(save_dir);

for i=1:1:length(feat_list)
    if ~isempty(strfind(feat_list{i}, 'label'))
        [pathstr, name, ext] = fileparts(feat_list{i});
        new_path = [pathstr, '/Labels/', name, ext];
        movefile(feat_list{i}, new_path);
    end
end

end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function s=dIOU(pred,gt)
Ibox=In(pred,gt);
I=(Ibox.b-Ibox.t)*(Ibox.r-Ibox.l);
U=(pred.b-pred.t)*(pred.r-pred.l)+(gt.b-gt.t)*(gt.r-gt.l)-I;
dI=dIbox(pred,gt);
dU=dUbox(pred,gt);
s.dt=0;
s.db=0;
s.dl=0;
s.dr=0;
if Ibox.t<Ibox.b&&Ibox.l<Ibox.r&&U>0
    s.dt=(U*dI.dt-I*dU.dt)/U^2;
    s.db=(U*dI.db-I*dU.db)/U^2;
    s.dl=(U*dI.dl-I*dU.dl)/U^2;
    s.dr=(U*dI.dr-I*dU.dr)/U^2;
else
    s.dt=0;
    s.db=0;
    s.dl=0;
    s.dr=0;
end
if pred.t<pred.b
    s.dt=s.dt;
    s.db=s.db;
else
    s.dt=s.db;
    s.db=s.dt;
end
if pred.l<pred.r
    s.dl=s.dl;
    s.dr=s.dr;
else
    s.dl=s.dr;
    s.dr=s.dl;
end

author: Zzh-tju

repo: DIoU

save_path: github-repos/MATLAB/Zzh-tju-DIoU


text:
function [ uvcoord ] = uv2coords( uv, width, height, planeID )
%UV2COORDS Summary of this function goes here
%   Detailed explanation goes here
if ~exist('planeID','var')
    planeID = 1;
end
if planeID~=1
    uv = xyz2uvN(uv2xyzN(uv, planeID), 1);
end

uvcoord = zeros(size(uv,1),2);
uvcoord(:,1) = min(round((uv(:,1)+pi)/2/pi*width+0.5), width);
uvcoord(:,2) = min(round((pi/2-uv(:,2))/pi*height+0.5), height);
end

author: zouchuhang

repo: LayoutNet


text:
function halton_test01 ( )

%*****************************************************************************80
%
%% HALTON_TEST01 tests HALTON, HALTON_STEP_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2009
%
%  Author:
%
%    John Burkardt
%
  dim_max = 3;
  test_num = 4;

  step_vec = [ 0, 5, 1000, 1000000 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HALTON_TEST01\n' );
  fprintf ( 1, '  HALTON returns the next element of a Halton sequence.\n' );
  fprintf ( 1, '  HALTON_STEP_SET sets the step.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we try several values of STEP.\n' );
  fprintf ( 1, '  We repeat the test for several dimensions.\n' );
  fprintf ( 1, '  We assume defaults for SEED, LEAP and BASE.\n' );

  for dim_num = 1 : dim_max

    for test = 1 : test_num

      halton_dim_num_set ( dim_num );
      n = 11;
      step = step_vec(test);
      halton_step_set ( step );
      seed(1:dim_num) = 0;
      halton_seed_set ( seed );
      for i = 1 : dim_num
        base(i) = prime ( i );
      end

      fprintf ( 1, '\n' );
      fprintf ( 1, '  DIM_NUM = %12d\n', dim_num );
      fprintf ( 1, '  N =    %12d\n', n );
      fprintf ( 1, '  STEP = %12d\n', step );
      i4vec_transpose_print ( dim_num, seed, '  SEED = ' );
      i4vec_transpose_print ( dim_num, base, '  BASE = ' );

      fprintf ( 1, '\n' );
      fprintf ( 1, '          STEP   Halton\n' );
      fprintf ( 1, '\n' );

      for j = 1 : n
        r = halton ( );
        fprintf ( 1, '  %12d  ', step+j-1 );
        for i = 1 : dim_num
          fprintf ( 1, '%12f  ', r(i) );
        end
        fprintf ( 1, '\n' );
      end

    end

  end

  return
end


text:
function [filesInCurDir,dirsInCurDir]=rdir(apath)
%Recursively search directories for all files and directories contained
%therein
% Example call: [filesInWorkingDir,dirsInWorkingDir]=rdir(pwd)
%
% 05/09/2013, http://www.matlabexperthelp.com/2012/08/15/recursive-directory-listing/

%% recursively get file contents of all directories
dirInfo=dir(apath);
dirInfo=dirInfo(~(strcmp({dirInfo.name},'..') | strcmp({dirInfo.name},'.' )));
for i=1:length(dirInfo)
    dirInfo(i).fullpath=fullfile(apath,dirInfo(i).name);
end
dirsInCurDir  = dirInfo([dirInfo.isdir]);
filesInCurDir = dirInfo(~[dirInfo.isdir]);

%% loop through dirs and get files and dirs inside of current dirs contents
for i=1:length(dirsInCurDir)
    % Call this function again and concatenate results
    [filesInThisDir, dirsInThisDir]=rdir(dirsInCurDir(i).fullpath);
    if ~isempty(filesInThisDir)
        filesInCurDir=[filesInCurDir; filesInThisDir ];
    end
    if ~isempty(dirsInThisDir)
        dirsInCurDir=[dirsInCurDir; dirsInThisDir];
    end
end

author: cerr


text:
% Fig. 9.9  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
% script to generate figure 9.9
% using the general simulation nonsim
clf;
N=0.4;
a=2.5;
r=0;
num=[1 1];
den=[1 0 0 ];
nf=1;
df=1;
hold on
grid on
for k=1:6
   r=r+2;
   sim('nonsim',30)
   plot(yn(:,1),yn(:,2))
end
title('Figure 9.9 Step responses of the nonlinear system')
xlabel('Time (sec)');
ylabel('y(t)');
gtext('r=2')
gtext('r=4')
gtext('r=6')
gtext('r=8')
gtext('r=10')
gtext('r=12')
hold off


text:
function spm_slice_timing(P, sliceorder, refslice, timing, prefix)
% Correct differences in slice acquisition times
% FORMAT spm_slice_timing(P, sliceorder, refslice, timing, prefix)
% P           - char array of image filenames
%               can also be a cell array of the above (multiple subjects).
% sliceorder  - slice acquisition order, a vector of integers, each
%               integer referring the slice number in the image file
%               (1=first), and the order of integers representing their
%               temporal acquisition order
%               OR vector containig the acquisition time for each slice
%               in milliseconds
% refslice    - slice for time 0
%               OR time in milliseconds for the reference slice
% timing      - additional information for sequence timing
%               timing(1) = time between slices
%                         = TA / (nslices - 1)
%               timing(2) = time between last slices and next volume
%                         = TR - TA
%               OR timing = [0 TR] when previous inputs are specified in
%               milliseconds
% prefix      - filename prefix for corrected image files, defaults to 'a'
%__________________________________________________________________________
%
%   Note: The sliceorder arg that specifies slice acquisition order is
%   a vector of N numbers, where N is the number of slices per volume.
%   Each number refers to the position of a slice within the image file.
%   The order of numbers within the vector is the temporal order in which
%   those slices were acquired.
%
%   To check the order of slices within an image file, use the SPM Display
%   option and move the crosshairs to a voxel co-ordinate of z=1.  This
%   corresponds to a point in the first slice of the volume.
%
%   The function corrects differences in slice acquisition times.
%   This routine is intended to correct for the staggered order of
%   slice acquisition that is used during echoplanar scanning. The
%   correction is necessary to make the data on each slice correspond
%   to the same point in time. Without correction, the data on one
%   slice will represent a point in time as far removed as 1/2 the TR
%   from an adjacent slice (in the case of an interleaved sequence).
%
%   This routine "shifts" a signal in time to provide an output
%   vector that represents the same (continuous) signal sampled
%   starting either later or earlier. This is accomplished by a simple
%   shift of the phase of the sines that make up the signal.
%
%   Recall that a Fourier transform allows for a representation of any
%   signal as the linear combination of sinusoids of different
%   frequencies and phases. Effectively, we will add a constant
%   to the phase of every frequency, shifting the data in time.
%
%   Shifter - This is the filter by which the signal will be convolved
%   to introduce the phase shift. It is constructed explicitly in
%   the Fourier domain. In the time domain, it may be described as
%   an impulse (delta function) that has been shifted in time the
%   amount described by TimeShift.
%
%   The correction works by lagging (shifting forward) the time-series
%   data on each slice using sinc-interpolation. This results in each
%   time series having the values that would have been obtained had
%   the slice been acquired at the same time as the reference slice.
%
%   To make this clear, consider a neural event (and ensuing hemodynamic
%   response) that occurs simultaneously on two adjacent slices. Values
%   from slice "A" are acquired starting at time zero, simultaneous to
%   the neural event, while values from slice "B" are acquired one
%   second later. Without corection, the "B" values will describe a
%   hemodynamic response that will appear to have began one second
%   EARLIER on the "B" slice than on slice "A". To correct for this,
%   the "B" values need to be shifted towards the Right, i.e., towards
%   the last value.
%
% Written by Darren Gitelman at Northwestern U., 1998
%
% Based (in large part) on ACQCORRECT.PRO from G. Aguirre and E. Zarahn
% at U. Penn.
%
% Modified by R. Henson, C. Buechel and J. Ashburner, FIL, to
% handle different reference slices and memory mapping.
%
% Modified by M. Erb, at U. Tuebingen, 1999, to ask for non-continuous
% slice timing and number of sessions.
%
% Modified by R. Henson for more general slice order and SPM2.
%
% Modified by A. Hoffmann, M. Woletz and C. Windischberger from Medical
% University of Vienna, Austria, to handle multi-band EPI sequences.
%__________________________________________________________________________
% Copyright (C) 1998-2014 Wellcome Trust Centre for Neuroimaging

% Darren Gitelman et al.
% $Id: spm_slice_timing.m 6130 2014-08-01 17:41:18Z guillaume $

SVNid = '$Rev: 6130 $';

%-Say hello
%--------------------------------------------------------------------------
SPMid = spm('FnBanner',mfilename,SVNid);

%-Parameters & Arguments
%==========================================================================
if nargin < 4, error('Not enough input arguments.'); end
if nargin < 5, prefix = 'a'; end

if ~iscell(P), P = {P}; end
nsubjects = numel(P);

% Acquisition order: 1=first slice in image
% Reference slice: 1=first slice in image, in Analyze format, slice 1 = bottom
% TR: Interscan interval (TR) {secs}
% TA: Acquisition Time (TA) {secs} [Def: TR-TR/nslices], TA <= TR
% timing(2) = TR - TA, time between last slices and next volume
% timing(1) = TA / (nslices -1), time between slices

Vin     = spm_vol(P{1}(1,:));
nslices = Vin(1).dim(3);

TR      = (nslices-1)*timing(1)+timing(2);
fprintf('%-40s: %30s\n','Number of slices is...',num2str(nslices))      %-#
fprintf('%-40s: %30s\n','Time to Repeat (TR) is...',num2str(TR))        %-#

if ~isequal(1:nslices,sort(sliceorder))
    if ~all(sliceorder >= 0 & sliceorder <= TR*1000)
        error('Input is neither slice indices nor slice times.');
    end
    unit = 'slice times (ms)';
else
    if ~ismember(refslice,sliceorder)
        error('Reference slice should contain a slice index.');
    end
    unit = 'slice indices';
end
fprintf('%-40s: %30s\n','Parameters are specified as...',unit)          %-#

if nslices ~= numel(sliceorder)
    error('Mismatch between number of slices and length of ''sliceorder'' vector.');
end

%-Slice timing correction
%==========================================================================
for subj = 1:nsubjects
    Vin   = spm_vol(P{subj});
    nimgo = numel(Vin);
    nimg  = 2^(floor(log2(nimgo))+1);
    if Vin(1).dim(3) ~= nslices
        error('Number of slices differ: %d vs %d.', nslices, Vin(1).dim(3));
    end

    % Create new header files
    Vout  = Vin;
    for k=1:nimgo
        Vout(k).fname  = spm_file(Vin(k).fname, 'prefix', prefix);
        if isfield(Vout(k),'descrip')
            desc = [Vout(k).descrip ' '];
        else
            desc = '';
        end
        Vout(k).descrip = [desc 'acq-fix ref-slice ' num2str(refslice)];
    end
    Vout = spm_create_vol(Vout);

    % Set up [time x voxels] matrix for holding image info
    slices = zeros([Vout(1).dim(1:2) nimgo]);
    stack  = zeros([nimg Vout(1).dim(1)]);

    task = sprintf('Correcting acquisition delay: session %d', subj);
    spm_progress_bar('Init',nslices,task,'planes complete');

    % Compute shifting amount from reference slice and slice order
    if isequal(unit,'slice times (ms)')
        % Compute time difference between the acquisition time of the
        % reference slice and the current slice by using slice times
        % supplied in sliceorder vector
        shiftamount = (sliceorder - refslice) / (1000 * TR);
    else
        rslice      = find(sliceorder==refslice);
        [Y, I]      = sort(sliceorder);
        shiftamount = (I - rslice) * timing(1) / TR;
    end

    % For loop to perform correction slice by slice
    for k = 1:nslices

        % Read in slice data
        B  = spm_matrix([0 0 k]);
        for m=1:nimgo
            slices(:,:,m) = spm_slice_vol(Vin(m),B,Vin(1).dim(1:2),1);
        end

        % Set up shifting variables
        len     = size(stack,1);
        phi     = zeros(1,len);

        % Check if signal is odd or even -- impacts how Phi is reflected
        %  across the Nyquist frequency. Opposite to use in pvwave.
        OffSet  = 0;
        if rem(len,2) ~= 0, OffSet = 1; end

        % Phi represents a range of phases up to the Nyquist frequency
        % Shifted phi 1 to right.
        for f = 1:len/2
            phi(f+1) = -1*shiftamount(k)*2*pi/(len/f);
        end

        % Mirror phi about the center
        % 1 is added on both sides to reflect Matlab's 1 based indices
        % Offset is opposite to program in pvwave again because indices are 1 based
        phi(len/2+1+1-OffSet:len) = -fliplr(phi(1+1:len/2+OffSet));

        % Transform phi to the frequency domain and take the complex transpose
        shifter = [cos(phi) + sin(phi)*sqrt(-1)].';
        shifter = shifter(:,ones(size(stack,2),1)); % Tony's trick

        % Loop over columns
        for i=1:Vout(1).dim(2)

            % Extract columns from slices
            stack(1:nimgo,:) = reshape(slices(:,i,:),[Vout(1).dim(1) nimgo])';

            % Fill in continous function to avoid edge effects
            for g=1:size(stack,2)
                stack(nimgo+1:end,g) = linspace(stack(nimgo,g),...
                    stack(1,g),nimg-nimgo)';
            end

            % Shift the columns
            stack = real(ifft(fft(stack,[],1).*shifter,[],1));

            % Re-insert shifted columns
            slices(:,i,:) = reshape(stack(1:nimgo,:)',[Vout(1).dim(1) 1 nimgo]);
        end

        % Write out the slice for all volumes
        for p = 1:nimgo
            Vout(p) = spm_write_plane(Vout(p),slices(:,:,p),k);
        end
        spm_progress_bar('Set',k);
    end
    spm_progress_bar('Clear');
end

fprintf('%-40s: %30s\n','Completed',spm('time'))                        %-#

author: spm


text:
function c=ref_dcti_1(f)
%REF_DCTI_1  Reference Discrete Consine Transform type I
%   Usage:  c=ref_dcti_1(f);
%
%

L=size(f,1);
W=size(f,2);

if L==1
  c=f;
  return;
end;

R=1/sqrt(2)*[eye(L);...
	     [zeros(L-2,1),flipud(eye(L-2)),zeros(L-2,1)]];

R(1,1)=1;
R(L,L)=1;

c=R'*dft(R*f);

author: ltfat

repo: ltfat


text:
function BC = createBC1D(meshvar)
% function BC = createBC1D(meshvar)
% Creates a boundary condition structure from a mesh structure
%
% SYNOPSIS:
%   BC = createBC1D(meshvar)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Define the right, and left boundary conditions
% (default = Neumann, i.e., a = 1, b = 0, c = 0)

right.a = 1;
right.b = 0;
right.c = 0;
right.periodic = 0;

left.a = 1;
left.b = 0;
left.c = 0;
left.periodic = 0;

bottom=[];
top=[];
back=[];
front=[];

BC= BoundaryCondition(meshvar, left, right, bottom, top, back, front);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool