text:
function result  = batch_omp_spr(Y, K, epsilon, sparse_output)
    if nargin < 4
        sparse_output = 1;
    end
    result = mex_batch_omp_spr(Y, K, epsilon, sparse_output);
end

author: indigits

repo: sparse-plex


text:
function f = write_mhd(filename, image, varargin)
% Write mhd
% arguments:
% write_mhd(filename,image)
% write_mhd(filename,image, param, value,...)
% params and default value:
%   'NDims'                     =   3
%   'BinaryData'                =   true
%   'BinaryDataByteOrderMSB'    =   false
%   'CompressedData'            =   false
%   'TransformMatrix'           =   1 0 0 0 1 0 0 0 1 (coherent with NDims)
%   'CenterOfRotation'          =   0 0 0 (coherent with NDims)
%   'AnatomicalOrientation'     =   RAI
%   'ElementNumberOfChannels'   =   1
%   'ElementType'               =   MET_FLOAT
%
% modifications by Thomas Kuestner, 2015
% ---------------------------------------------------------------------

% revise this
element_types=struct('uchar','MET_UCHAR','double','MET_DOUBLE','single','MET_FLOAT','logical','MET_CHAR','int8','MET_UCHAR','uint8','MET_UCHAR','int16','MET_SHORT','uint16','MET_USHORT','int32','MET_INT','uint32','MET_UINT');
% element_types=struct('double','MET_DOUBLE','int8','MET_CHAR','uint8','MET_UCHAR','int16','MET_SHORT','uint16','MET_USHORT','int32','MET_INT','uint32','MET_UINT');

%Look for arguments
NDims                       = numel(image.size);
BinaryData                  = 'True';
BinaryDataByteOrderMSB      = 'False';
CompressedData              = 'False';

TransformMatrix  = num2str(reshape(image.orientation,1,NDims*NDims)); %'1 0 0 0 1 0 0 0 1';% (coherent with NDims)
CenterOfRotation = num2str(zeros(1,NDims)) ;

AnatomicalOrientation =   'RAI';
if (isa(image,'VectorImageType') || isfield(image,'datax'))
    ElementNumberOfChannels = 3;
else
    ElementNumberOfChannels = 1;
end

%ElementType =getfield(element_types,class(image.data));
EType ='double';

for i=1:size(varargin,2)
    switch (lower(varargin{i}))
        case 'ndims'
            NDims =  varargin{i+1};
        case 'binarydata'
            BinaryData  =   varargin{i+1};
        case 'binarydatabyteordermsb'
            BinaryDataByteOrderMSB =   varargin{i+1};
        case 'compresseddata'
            CompressedData =   varargin{i+1};
        case 'transformmatrix'
            TransformMatrix =   varargin{i+1};
        case 'centerofrotation'
            CenterOfRotation =   varargin{i+1};
        case 'anatomicalorientation'
            AnatomicalOrientation =   varargin{i+1};
        case 'elementnumberofchannels'
            ElementNumberOfChannels =   varargin{i+1};
        case 'elementtype'
            EType =  varargin{i+1};
    end
end
ElementType =  element_types.(EType );
% Extract path and filenames---
path_ = regexp(filename,'/','split');
path = [];
for i=2:(size(path_,2)-1)

    path = [path '/' path_{i}];
end

raw_ = regexp(filename,'\.','split');
raw = [];
for i=1:(size(raw_,2)-1)
    raw = [raw  raw_{i} '.'];
end
rawfile = [raw 'raw'];
rawfilename_ = regexp(rawfile,'/','split');
rawfilename = rawfilename_{size(rawfilename_,2)};
% write *.mhd file
fid=fopen(filename,'w','native');
fprintf(fid,'ObjectType = Image\n');
fprintf(fid,'NDims = %d\n',NDims);
fprintf(fid,'BinaryData = %s\n',BinaryData);
fprintf(fid,'BinaryDataByteOrderMSB = %s\n',BinaryDataByteOrderMSB);
fprintf(fid,'CompressedData = %s\n',CompressedData);
fprintf(fid,'TransformMatrix = %s\n',TransformMatrix);
fprintf(fid,'CenterOfRotation = %s\n',CenterOfRotation);
fprintf(fid,'AnatomicalOrientation = %s\n',AnatomicalOrientation);
fprintf(fid,'%s\n',['Offset = ' num2str(image.origin(:)') ]);
fprintf(fid,'%s\n',['ElementSpacing = ' num2str(image.spacing(:)') ]);
fprintf(fid,'%s\n',['DimSize = ' num2str(image.size(:)') ]);

fprintf(fid,'ElementNumberOfChannels = %d\n',ElementNumberOfChannels);
fprintf(fid,'ElementType = %s\n',ElementType);
fprintf(fid,'ElementDataFile = %s\n',rawfilename);

fclose(fid);

% write *.raw file

%dataToWrite=image.data;

%if (islogical(dataToWrite))
%     dataToWrite = cast(dataToWrite,'uint8');
% else

fid=fopen(rawfile,'w','native');
if (ElementNumberOfChannels ==1)
    %fwrite(fid,image.data,class(dataToWrite) );
    dataAll = cast(image.data,EType);
elseif (ElementNumberOfChannels==2)
    %do nothing
elseif (ElementNumberOfChannels==3)
    %write point per point
    dataAll = cast(permute(cat(NDims+1,image.datax,  image.datay, image.dataz),...
        [NDims+1 1:NDims]),EType);
elseif (ElementNumberOfChannels==4)
    %TODO
end
fwrite(fid,dataAll,EType);
fclose(fid);

f=fid;

end

author: thomaskuestner


text:
function [Hhat_LMMSE,C_LMMSE,tau_p,R,H] = functionChannelEstimates_impairments(R,channelGaindB,nbrOfRealizations,M,K,L,p,f,kappatUE,kapparBS)
%Generate the channel realizations and estimates of these channels for all
%UEs in the entire network. The channels are assumed to be correlated
%Rayleigh fading. The LMMSE estimator that takes transceiver hardware
%impairments into account is used.
%
%INPUT:
%R                 = M x M x K x L x L matrix with spatial correlation
%                    matrices for all UEs in the network. R(:,:,k,j,l) is
%                    the correlation matrix for the channel between UE k
%                    in cell j and the BS in cell l. This such matrix can
%                    either include the average channel gain or can be
%                    normalized arbitrarily.
%channelGaindB     = K x L x L matrix containing the average channel gains
%                    in dB of all the channels, if these are not already
%                    included in the spatial channel correlation matrices.
%                    The product R(:,:,k,j,l)*10^(channelGaindB(k,j,l)/10)
%                    is the full spatial channel correlation matrix.
%nbrOfRealizations = Number of channel realizations
%M                 = Number of antennas per BS
%K                 = Number of UEs per cell
%L                 = Number of BSs and cells
%p                 = Uplink transmit power per UE (same for everyone)
%f                 = Pilot reuse factor
%kappatUE          = Hardware quality of the UEs' transmitters
%kapparBS          = Hardware quality of the BSs' receivers
%
%OUTPUT:
%Hhat_LMMSE   = M x nbrOfRealizations x K x L x L matrix with the LMMSE
%               channel estimates. The matrix Hhat_LMMSE(:,n,k,j,l) is the
%               n:th channel estimate of the channel between UE k in cell j
%               and the BS in cell l.
%C_LMMSE      = M x M x K x L x L matrix with estimation error correlation
%               matrices when using LMMSE estimation. The matrix is
%               organized in the same way as R.
%tau_p        = Length of pilot sequences
%R            = Scaled version of the input spatial correlation matrices R,
%               where the channel gains from channelGaindB are included
%H            = M x nbrOfRealizations x K x L x L matrix with the true
%               channel realizations. The matrix is organized in the same
%               way as Hhat_LMMSE.
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%% Generate channel realizations

%Generate uncorrelated Rayleigh fading channel realizations
H = (randn(M,nbrOfRealizations,K,L,L)+1i*randn(M,nbrOfRealizations,K,L,L));

%Prepare a matrix to save the channel gains per UE
betas = zeros(K,L,L);

%Go through all channels and apply the channel gains to the spatial
%correlation matrices
for j = 1:L

    for l = 1:L

        for k = 1:K

            %Extract channel gain in linear scale
            betas(k,j,l) = 10^(channelGaindB(k,j,l)/10);

            %Apply channel gain to correlation matrix
            R(:,:,k,j,l) = betas(k,j,l) * R(:,:,k,j,l);

            %Apply correlation to the uncorrelated channel realizations
            Rsqrt = sqrtm(R(:,:,k,j,l));
            H(:,:,k,j,l) = sqrt(0.5)*Rsqrt*H(:,:,k,j,l);

        end

    end

end

%% Perform channel estimation

%Length of pilot sequences
tau_p = f*K;

%Generate pilot pattern
if f == 1

    pilotPattern = ones(L,1);

elseif f == 2 %Only works in the running example with its 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 2; 1; 2; 1]);

elseif f == 4 %Only works in the running example with its 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 3; 4; 3; 4]);

elseif f == 16 %Only works in the running example with its 16 BSs

    pilotPattern = (1:L)';

end

%Store identity matrix of size M x M
eyeM = eye(M);

%Generate realizations of normalized noise
Np = sqrt(0.5)*(randn(M,nbrOfRealizations,K,L,f) + 1i*randn(M,nbrOfRealizations,K,L,f));

%Prepare to store LMMSE channel estimates
Hhat_LMMSE = zeros(M,nbrOfRealizations,K,L,L);

%Prepare to store estimation error correlation matrices
C_LMMSE = zeros(M,M,K,L,L);

%% Go through all f pilot groups
for g = 1:f

    %Create transmit distortion term
    etaUE = sqrt(0.5*p*tau_p*kapparBS*(1-kappatUE))*(randn(1,nbrOfRealizations,K,L,K) + 1i*randn(1,nbrOfRealizations,K,L,K));

    %Go through all cells
    for j = 1:L

        %Extract the cells that belong to pilot group g
        groupMembers = find(g==pilotPattern)';

        %Go through all UEs
        for k = 1:K

            %Create receive distortion term
            etabarBS = sqrt(0.5*p*tau_p*(1-kapparBS))*(randn(M,nbrOfRealizations,K,L) + 1i*randn(M,nbrOfRealizations,K,L));

            %Sum up transmit and receive distortion terms
            etaSum = repmat(etaUE(:,:,:,:,k),[M 1 1 1]) + etabarBS;

            %Compute processed pilot signal for all UEs that use these pilots, according to (6.22)
            ypk = sqrt(p*kappatUE*kapparBS)*tau_p*sum(H(:,:,k,g==pilotPattern,j),4) + squeeze(sum(sum(H(:,:,:,:,j).*etaSum,4),3)) + sqrt(tau_p)*Np(:,:,k,j,g);

            %Compute sum of all UEs' correlation matrices to BS j
            RjAll = sum(sum(R(:,:,:,:,j),3),4);

            %Compute the matrix that is inverted in the LMMSE estimator
            PsiInv = (p*kappatUE*kapparBS*tau_p*sum(R(:,:,k,g==pilotPattern,j),4) + p*(1-kappatUE)*kapparBS*RjAll + p*(1-kapparBS)*diag(diag(RjAll)) + eyeM);

            %Go through the cells in pilot group g
            for l = groupMembers

                %Compute LMMSE estimate of channel between BS l and UE k in
                %cell j using (6.23) in Theorem 6.1
                RPsi = R(:,:,k,l,j) / PsiInv;
                Hhat_LMMSE(:,:,k,l,j) = sqrt(p*kappatUE*kapparBS)*RPsi*ypk;

                %Compute corresponding estimation error correlation matrix, using (6.26)
                C_LMMSE(:,:,k,l,j) = R(:,:,k,l,j) - p*kappatUE*kapparBS*tau_p*RPsi*R(:,:,k,l,j);

            end

        end

    end

end

author: emilbjornson

repo: massivemimobook


text:
clear all;
fprintf('Compiling mex files ... ');

% Orthogonal matching pursuit
mex mex/mat_omp.c mex/perform_omp.c mex/invert.c mex/matrix_vector.c -o perform_omp_mex

disp('done.');

author: gpeyre

repo: matlab-toolboxes


text:
function beta = LpAnalysis(H, Ht, A, At, y, p, normfac, err, insweep, decfac, tol)

% Algorithm for solving problems of the form:
% min ||Ax||_p s.t. ||y-Hx||_2 < err

% Inputs
% A - Sparsifying Forward Transform
% At - Sparsifying Backward Transform
% H - Forward Measurement Operator
% Ht - Backward Measurement Operator
% y - collected data
% p - non-convex norm
% normfac - highest singular value of A (default 1)
% err - two norm mismatch (default 1e-6)
% insweep - number of sweeps for internal loop (default 50)
% decfac - lambda decrease factor for outside loop (default 0.5)
% tol - tolerance for convergence (default 1e-4)

% copyright (c) Angshul Majumdar 2010

c = 1.1;

if nargin < 6
    normfac=1;
end
if nargin < 7
    p = 1;
end
if nargin < 8
    err = 1e-6;
end
if nargin < 9
    insweep = 50;
end
if nargin < 10
    decfac = 0.5;
end
if nargin < 11
    tol = 1e-4;
end

alpha = 1.1*normfac;
x_initial = Ht(y);
z_initial = zeros(length(A(x_initial)),1);

x = x_initial; z = z_initial;

lambdaInit = decfac*max(abs(Ht(y))); lambda = lambdaInit;

f_current = norm(y-H(x_initial)) + lambda*norm(A(x),p);

while lambda > lambdaInit*tol
    % debug lambda
    for ins = 1:insweep
        f_previous = f_current;

        b = x + (1/alpha)*(Ht(y-H(x)));
        z = (c*z + A(b-At(z)))./((2*alpha/lambda)*(abs(A(x))).^(2-p)+c);
        x = b - At(z);

        f_current = norm(y-H(x)) + lambda*norm(A(x),p);

        if norm(f_current-f_previous)/norm(f_current + f_previous)<tol
            break;
        end
    end
    if norm(y-H(x))<err
        break;
    end
    lambda = decfac*lambda;
end
beta = x;

author: Sable

repo: mcbench-benchmarks


text:
function plaza = clear_boundary(input)
plaza = input;
[a,b] = size(plaza);
for i = 1:b
if plaza(a,i) > 0
plaza(a,i) = 0;
end
end

author: Eurus-Holmes


text:
function simulink2svg
% Show the hidden handles to get access to the simulink figures
set(0,'ShowHiddenHandles','on');
% Iterate over all children of the root
child = get(0,'children');
for i = 1:length(child)
    % Pick all scope figures
    if strcmp(get(child(i),'Tag'),'SIMULINK_SIMSCOPE_FIGURE')
        % In order to get the correct background turn inverted background
        % off. If you like a white background you should invert all labels
        % and lines.
        set(child(i),'InvertHardcopy','off');
        % Use plot2svg
        plot2svg(['Scope_' num2str(i) '.svg'], child(i))
        set(0,'ShowHiddenHandles','off')
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [removingList] = plotTracks(ax, tracks, currentFrame, params, removingList)
%% Plot bounding boxes
% Remove old lines
if nargin > 4
    if size(removingList, 1) > 0
       for i = 1:size(removingList, 1)
          delete(removingList(i));
       end
    end
end
% Important parameters for plotting the bounding boxes
removingList = [];
boundingBoxColor = params.boundingBoxColor;
% We need to draw the bounding boxes from 0 to negative because the image
% coordinate origin is in the upper left corner. Therefore the bounding box
% coordinates are mirrored.
ySign = -1;
% Now plot the bounding boxes together with annotations
for iTrack = 1:size(tracks, 2)
    track = tracks(iTrack);
    initialFrame = track.initialFrame;
    finalFrame = track.finalFrame;
    if (initialFrame <= currentFrame) && (currentFrame < finalFrame)
        % Get internal track index for the chosen frame
        currentIndex = currentFrame - initialFrame + 1;
        % Try to get the bounding box
        try
            boundingBox = track.bbox(currentIndex, :);
            velocity = track.xVelocity(currentIndex);
        catch
           continue
        end
        % Plot the bounding box
        positionBoundingBox = [boundingBox(1) ...
                               ySign*boundingBox(2)+ySign*boundingBox(4) ...
                               boundingBox(3) ...
                               boundingBox(4)];
        if params.stop == 1
            rect = rectangle(ax, 'Position', positionBoundingBox, ...
                                 'FaceColor', boundingBoxColor, ...
                                 'EdgeColor', [0 0 0], ...
                                 'PickableParts', 'visible', ...
                                 'ButtonDownFcn',{@params.onClick, track, currentFrame});
        else
            rect = rectangle(ax, 'Position', positionBoundingBox, ...
                                 'FaceColor', boundingBoxColor, ...
                                 'EdgeColor', [0 0 0], ...
                                 'PickableParts', 'visible');
        end
        removingList = [removingList; rect];

        % Plot the triangle that represents the direction of the vehicle
        if velocity < 0
            centroidSign = 1;
            front = boundingBox(1) + (boundingBox(3)*0.2);
            triangleXPosition = [front front boundingBox(1)];
        else
            centroidSign = -1;
            front = boundingBox(1) + boundingBox(3) - (boundingBox(3)*0.2);
            triangleXPosition = [front front boundingBox(1)+boundingBox(3)];
        end
        triangleYPosition = [ySign*boundingBox(2) ...
                             ySign*boundingBox(2)+ySign*boundingBox(4)...
                             ySign*boundingBox(2)+ySign*(boundingBox(4)/2)];
        triangle = fill(ax, triangleXPosition, triangleYPosition, [0 0 0]);
        removingList = [removingList; triangle];

        % Plot the text annotation
        if params.plotTextAnnotation
            velocityKmh = abs(velocity) * 3.6;
            boundingBoxAnnotationText = sprintf('%s|%.2fkm/h|ID%d', ...
                                                track.class(1), ...
                                                velocityKmh, ...
                                                track.id);
            textAnnotation = text(ax, boundingBox(1), ySign*boundingBox(2)+0.6, ...
                                  boundingBoxAnnotationText, 'FontSize',8);
            textAnnotation.BackgroundColor = [1 1 0.3];
            textAnnotation.Margin = .5;
            removingList = [removingList; textAnnotation];
        end

        % Plot the track line
        if params.plotTrackingLine
            relevantBoundingBoxes = track.bbox(1:currentIndex, :);
            centroids = [relevantBoundingBoxes(:,1) + relevantBoundingBoxes(:,3)/2, ...
                         relevantBoundingBoxes(:,2) + relevantBoundingBoxes(:,4)/2];
            plottedCentroids = line(ax, centroids(:, 1)+ centroidSign * (boundingBox(3)/2), ...
                                    ySign*centroids(:, 2));
            removingList = [removingList; plottedCentroids];
        end
    end
end
xlim(ax, [0 params.trackWidth]);
axis off;
end

author: RobertKrajewski

repo: highD-dataset


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, space, eigv, eigf, gnum, dofs_ornt] = ...
              mp_solve_maxwell_eig_2d (problem_data, method_data)

warning ('geopdes:obsolete','Function MP_SOLVE_MAXWELL_EIG_2D is obsolete. Using SOLVE_MAXWELL_EIG instead')

[geometry, msh, space, eigv, eigf, gnum, dofs_ornt] = mp_solve_maxwell_eig (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
classdef MshReader < FileReader

    properties (Access = public)

    end

    properties (Access = private)
        strcontent
        contentCell
        coord
        connec
        nElem
        nNodes
        nNodeElem
        numbers
    end

    methods (Access = public)

        function obj = MshReader(cParams)
            obj.init(cParams)
        end

        function m = read(obj)
            obj.openFile();
            obj.splitIntoLines();
            obj.getDimensions();
            obj.extractNumbers();
            obj.getCoordAndConnec();
            obj.closeFile();
            m = obj.createMesh();
        end

        function d = getDataBase(obj)
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.filePath = cParams.filePath;
        end

        function splitIntoLines(obj)
            obj.strcontent = fileread(obj.filePath);
            obj.contentCell = splitlines(obj.strcontent);
        end

        function getDimensions(obj)
            fLCell = split(obj.contentCell{1});
            fLCoor = 3;
            lLCoor = find(strcmp(obj.contentCell, 'End Coordinates'))-1;
            fLConnec = find(strcmp(obj.contentCell, 'Elements')) +1;
            lLConnec  = length(obj.contentCell) - 2;
            obj.nNodes    = lLCoor - fLCoor + 1;
            obj.nElem     = lLConnec - fLConnec + 1;
            obj.nNodeElem = str2double(fLCell{end});
        end

        function extractNumbers(obj)
            expr = '(-?0?\.?\d*)';
            match = regexp(obj.strcontent,expr, 'tokens', 'dotexceptnewline');
            cellfunres2 = cellfun(@cell2mat,match,'UniformOutput',false);
            numbersAll = str2double(cellfunres2);
            obj.numbers = numbersAll(3:end);
        end

        function getCoordAndConnec(obj)
            nN  = obj.nNodes;
            nE  = obj.nElem;
            nNE = obj.nNodeElem;
            coordRaw  = reshape(obj.numbers(1:nN*4), [4, nN])';
            connecRaw = reshape(obj.numbers(nN*4+1:end), [nNE+1, nE])';
            if isequal(coordRaw(:,4),zeros(length(coordRaw),1))
                coordRaw = coordRaw(:,1:3);
            end
            obj.coord  = coordRaw(:,2:end);
            obj.connec = connecRaw(:,2:end);
        end

        function m = createMesh(obj)
            s.coord  = obj.coord;
            s.connec = obj.connec;
            m = Mesh(s);
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% LRR | ROSL | Robust Orthonormal Subspace Learning (Shu et al. 2014)
% process_video('LRR', 'ROSL', 'dataset/demo.avi', 'output/demo_LRR-ROSL.avi');

K = 1; % The initialiation of the subspace dimension
tol = 1e-5;
maxIter = 30;
lambda = 1e-1; %2e-3;
[~,~,E_hat,A_hat] = inexact_alm_rosl(M,K,lambda,tol,maxIter);
L = A_hat;
S = E_hat;

author: andrewssobral

repo: lrslibrary


text:
function activations = virtual_electrode_activation(phis, thetas, sigmas, chanlocs, scalpmap, sph_radius)
% virtual_electrode_activation estimates the activation of electrodes in
% arbirtrary locations given a scalp map for an IC and the locations of the
% electrodes in that scalp map.
%
% Input:
% phis: Longitudes of points for which a virtual activation is desired.
%
% thetas: Latitudes of points for which a virtual activation is desired.
%
% sigmas: Standard devations of Gaussian kernels used to weight the
% original electrode activations in the calculation of the virtual
% electrode activations.
%
% chanlocs: Locations of original electrodes. Must be a struct containing
% the fields sph_phi and sph_theta. These fields should contain the
% longitudes (sph_phi) and latitudes (sph_theta) of the original
% electrodes.
%
% scalpmap: Scalp map (topography, i.e. a column of the field icawinv in
% the EEGLab data structure) for a single IC.
%
% sph_radius: Radius of the head.
%
% Output:
% activations: Topography of IC on the virtual electrode cap, i.e.
% activations of the specified virtual electrodes on this IC's scalp map.

% Copyright (C) 2013  Laura Froelich (laura.frolich@gmail.com)
%
% 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

all_electrodes_phi = cell2mat({chanlocs.sph_phi});
all_electrodes_theta = cell2mat({chanlocs.sph_theta});
nelectrodes = length(all_electrodes_theta);
nvirtual_electrodes = length(phis);

activations = NaN(nvirtual_electrodes, 1);
for ivirtual_electrode = 1:nvirtual_electrodes
    distances = spherical_distance(repmat(phis(ivirtual_electrode), 1, nelectrodes),...
        repmat(thetas(ivirtual_electrode), 1, nelectrodes),...
        all_electrodes_phi, all_electrodes_theta, sph_radius);
    distances = distances.^2;
    exp_distances = exp(-distances/(2*sigmas(ivirtual_electrode)^2));
    exp_distances=exp_distances/sum(exp_distances);
    activations(ivirtual_electrode) = sum(scalpmap.*exp_distances);

end
end

repo: REST


text:
function [legend_h,object_h,plot_h,text_strings] = columnlegend(numcolumns, str, location)
%
%   columnlegend creates a legend with a specified number of columns.
%
%   columnlegend(numcolumns, str, location)
%       numcolumns - number of columns in the legend
%       cellstr - cell array of strings for the legend
%       location - location variable for legend, default is 'NorthEast'
%                  possible values: 'NorthWest', 'NorthEast', 'SouthEast', 'SouthWest'
%
%   example:
%      legend_str = [];
%      for i=1:10,
%           x = 1:i:(10*i);
%           plot(x); hold on;
%           legend_str = [legend_str; {num2str(i)}];
%      end
%      columnlegend(3, legend_str, 'NorthWest');
%
%
%   Author: Simon Henin <shenin@gc.cuny.edu>

if nargin < 3
   location = 'NorthEast';
end

%create the legend
[legend_h,object_h,plot_h,text_strings] = legend(str);

%some variables
numlines = length(str);
numpercolumn = ceil(numlines/numcolumns);

%get old width, new width and scale factor
pos = get(legend_h, 'position');
width = numcolumns*pos(3);
rescale = pos(3)/width;

%get some old values so we can scale everything later
xdata = get(object_h(numlines+1), 'xdata');
ydata1 = get(object_h(numlines+1), 'ydata');
ydata2 = get(object_h(numlines+3), 'ydata');

%we'll use these later to align things appropriately
sheight = ydata1(1)-ydata2(1);                  % height between data lines
height = ydata1(1);                             % height of the box. Used to top margin offset
line_width = (xdata(2)-xdata(1))*rescale;   % rescaled linewidth to match original
spacer = xdata(1)*rescale;                    % rescaled spacer used for margins

%put the legend on the upper left corner to make initial adjustments easier
loci = get(gca, 'position');
set(legend_h, 'position', [loci(1) pos(2) width pos(4)]);

col = -1;
for i=1:numlines,
    if mod(i,numpercolumn)==1,
        col = col+1;
    end

    if i==1
        linenum = i+numlines;
    else
        linenum = linenum+2;
    end
    labelnum = i;

    position = mod(i,numpercolumn);
    if position == 0,
         position = numpercolumn;
    end

    %realign the labels
    set(object_h(linenum), 'ydata', [(height-(position-1)*sheight) (height-(position-1)*sheight)]);
    set(object_h(linenum), 'xdata', [col/numcolumns+spacer col/numcolumns+spacer+line_width]);

    set(object_h(linenum+1), 'ydata', [height-(position-1)*sheight height-(position-1)*sheight]);
    set(object_h(linenum+1), 'xdata', [col/numcolumns+spacer*3.5 col/numcolumns+spacer*3.5]);

    set(object_h(labelnum), 'position', [col/numcolumns+spacer*2+line_width height-(position-1)*sheight]);

end

%unfortunately, it is not possible to force the box to be smaller than the
%original height, therefore, turn it off and set background color to none
%so that it no longer appears
set(legend_h, 'Color', 'None', 'Box', 'off');

%let's put it where you want it
pos = get(legend_h, 'position');
fig_pos = get(gca, 'position');
switch lower(location),
    case {'northeast'}
        set(legend_h, 'position', [pos(1)+fig_pos(3)-pos(3) pos(2) pos(3) pos(4)]);
    case {'southeast'}
        set(legend_h, 'position', [pos(1)+fig_pos(3)-pos(3) fig_pos(2)-pos(4)/2+pos(4)/4 pos(3) pos(4)]);
    case {'southwest'}
        set(legend_h, 'position', [fig_pos(1) fig_pos(2)-pos(4)/2+pos(4)/4 pos(3) pos(4)]);
end


text:
function varargout = permute(varargin)
% Can not be permuted.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

error('file_array objects can not be permuted.');


text:
function vl_testnn(varargin)
%VL_TESTNN Run MatConvNet test suite
% VL_TESTNN('option', value, ...) takes the following options:
%  `cpu`:: true
%    Run the CPU tests.
%
%  `gpu`:: false
%    Run the GPU tests.
%
%  `single`:: true
%    Perform tests in single precision.
%
%  `double`:: false
%    Perform tests in double precision.
%
%  `command`:: `'nn'`
%    Run only tests which name starts with the specified substring.
%    E.g. `vl_testnn('command', 'nnloss') would run only the nnloss tests.
%
%  `break`:: false
%    Stop tests in case of error.
%
%  `tapFile`:: ''
%    Output the test results to a file. If a specified file does
%    exist it is overwritten.
%
%  `suiteDir`:: ''
%    Specifies the directory where the test suite files are located. If
%    left empty, the default suite is used (<MatConvNetRoot>/matlab/xtest/
%    /suite).
%
%  This function uses the Matlab unit testing framework which was
%  introduced in Matlab R2013a (v8.1).

% Copyright (C) 2015-16 Andrea Vedaldi, Karel Lenc.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.cpu = true ;
opts.gpu = false ;
opts.single = true ;
opts.double = false ;
opts.command = 'nn' ;
opts.break = false ;
opts.tapFile = '';
opts.suiteDir = '' ;
opts = vl_argparse(opts, varargin) ;

import matlab.unittest.constraints.* ;
import matlab.unittest.selectors.* ;
import matlab.unittest.plugins.TAPPlugin;
import matlab.unittest.plugins.ToFile;
addpath(fullfile(vl_rootnn, 'examples'));

% Choose which tests to run
sel = HasName(StartsWithSubstring(opts.command)) ;
if ~opts.gpu
  sel = sel & ~HasName(ContainsSubstring('device=gpu')) ;
end
if ~opts.cpu
  sel = sel & ~HasName(ContainsSubstring('device=cpu')) ;
end
if ~opts.double
  sel = sel & ~HasName(ContainsSubstring('dataType=double')) ;
end
if ~opts.single
  sel = sel & ~HasName(ContainsSubstring('dataType=single')) ;
end

root = fileparts(mfilename('fullpath')) ;
if isempty(opts.suiteDir)
  opts.suiteDir = fullfile(root, 'suite') ;
else  % any external subclasses of nntest will need it to be on the path
  addpath(fullfile(root, 'suite')) ;
end

% Run tests
suite = matlab.unittest.TestSuite.fromFolder(opts.suiteDir, sel) ;
runner = matlab.unittest.TestRunner.withTextOutput('Verbosity',3);
if opts.break
  runner.addPlugin(matlab.unittest.plugins.StopOnFailuresPlugin) ;
end
if ~isempty(opts.tapFile)
  if exist(opts.tapFile, 'file')
    delete(opts.tapFile);
  end
  runner.addPlugin(TAPPlugin.producingOriginalFormat(ToFile(opts.tapFile)));
end
result = runner.run(suite);
display(result)

repo: LapSRN


text:
function [ t, rank ] = subset_colex_successor ( n, t, rank )

%*****************************************************************************80
%
%% SUBSET_COLEX_SUCCESSOR computes the subset colexicographic successor.
%
%  Discussion:
%
%    In the original code, there is a last element with no successor.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the master set.
%    N must be positive.
%
%    Input/output, integer T(N), describes a subset.  T(I) is 0 if
%    the I-th element of the master set is not in the subset, and is
%    1 if the I-th element is part of the subset.
%    On input, T describes a subset.
%    On output, T describes the next subset in the ordering.
%    If the input T was the last in the ordering, then the output T
%    will be the first.
%
%    Input/output, integer RANK, the rank.
%    If RANK = -1 on input, then the routine understands that this is
%    the first call, and that the user wishes the routine to supply
%    the first element in the ordering, which has RANK = 0.
%    In general, the input value of RANK is increased by 1 for output,
%    unless the very last element of the ordering was input, in which
%    case the output value of RANK is 0.
%

%
%  Return the first element.
%
  if ( rank == -1 )
    t(1:n) = 0;
    rank = 0;
    return
  end
%
%  Check.
%
  subset_check ( n, t );

  for i = 1 : n

    if ( t(i) == 0 )
      t(i) = 1;
      rank = rank + 1;
      return
    else
      t(i) = 0;
    end

  end

  rank = 0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y,dzdg,dzdb] = vl_nnbnorm_old(x,g,b,varargin)
% VL_NNBNORM  CNN batch normalisation
%   Y = VL_NNBNORM(X,G,B) computes the batch normalization of the
%   input X. This is defined as:
%
%      Y(i,j,k,t) = G(k) * (X(i,j,k,t) - mu(k)) / sigma(k) + B(k)
%
%   where
%
%      mu(k) = mean_ijt X(i,j,k,t),
%      sigma(k) = sqrt(sigma2(k) + EPSILON),
%      sigma2(k) = mean_ijt (X(i,j,k,t) - mu(k))^2
%
%   are respectively the per-channel mean, standard deviation, and
%   variance of the input and G(k) and B(k) define respectively a
%   multiplicative and additive constant to scale each input
%   channel. Note that statistics are computed across all feature maps
%   in the batch packed in the 4D tensor X. Note also that the
%   constant EPSILON is used to regularize the computation of sigma(k)
%
%   [Y,DZDG,DZDB] = VL_NNBNORM(X,G,B,DZDY) computes the derviatives of
%   the output Z of the network given the derivatives with respect to
%   the output Y of this function.
%
%   VL_NNBNROM(..., 'Option', value) takes the following options:
%
%   `Epsilon`:: 1e-4
%       Specify the EPSILON constant.
%
%   See also: VL_NNNORMALIZE().

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% ISSUE - needs to store internal state, another reason for having classes?

% -------------------------------------------------------------------------
%                                                             Parse options
% -------------------------------------------------------------------------

opts.epsilon = 1e-4 ;
backMode = numel(varargin) > 0 && ~ischar(varargin{1}) ;
if backMode
  dzdy = varargin{1} ;
  opts = vl_argparse(opts, varargin(2:end)) ;
else
  opts = vl_argparse(opts, varargin) ;
end

% -------------------------------------------------------------------------
%                                                                    Do job
% -------------------------------------------------------------------------

x_size = [size(x,1), size(x,2), size(x,3), size(x,4)] ;
g_size = size(g) ;
b_size = size(b) ;
g = reshape(g, [1 x_size(3) 1]) ;
b = reshape(b, [1 x_size(3) 1]) ;
x = reshape(x, [x_size(1)*x_size(2) x_size(3) x_size(4)]) ;

mass = prod(x_size([1 2 4])) ;
mu = sum(sum(x,1),3) / mass  ;
y = bsxfun(@minus, x, mu); % y <- x_mu
sigma2 = sum(sum(y .* y,1),3) / mass + opts.epsilon ;
sigma = sqrt(sigma2) ;

if ~backMode
  y = bsxfun(@plus, bsxfun(@times, g ./ sigma, y), b) ;
else
  % remember: y contains x_mu
  dzdy = reshape(dzdy, size(x)) ;
  dzdg = sum(sum(dzdy .* y,1),3) ./ sigma ;
  dzdb = sum(sum(dzdy,1),3) ;

  muz = dzdb / mass;
  y = ...
    bsxfun(@times, g ./ sigma, bsxfun(@minus, dzdy, muz)) - ...
    bsxfun(@times, g .* dzdg ./ (sigma2 * mass), y) ;

  dzdg = reshape(dzdg, g_size) ;
  dzdb = reshape(dzdb, b_size) ;
end

y = reshape(y, x_size) ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function dstimgs = slpixlinnorm(imgs, mu, sigma)
%SLPIXLINNORM Performs linear normalization on pixel values
%
% $ Syntax $
%   - dstimgs = slpixlinnorm(imgs);
%   - dstimgs = slpixlinnorm(imgs, mu, sigma)
%
% $ Arguments $
%   - imgs:     the array of images
%   - mu:       the mean pixel value to be normalized to (default = 0)
%   - sigma:    the standard deviation relative to mean pixel (default = 1)
%
% $ Description $
%   - dstimgs = slpixlinnorm(imgs, mu, sigma) performs linear normalization
%     on the image pixels so that the average pixel value is set to mu
%     while the standard deviation is set to sigma. The normalization is
%     conducted on each page(channel) respectively.
%
%   - dstimgs = slpixlinnorm(imgs) performs linear pixel value
%     normalization using default values.
%
% $ History $
%   - Created by Dahua Lin, on Aug 8th, 2006
%

%% parse and verify input arguments

if ~isa(imgs, 'double')
    imgs = im2double(imgs);
end
[h, w, n] = size(imgs);

if nargin < 2 || isempty(mu)
    mu = 0;
end

if nargin < 3 || isempty(sigma)
    sigma = 1;
end

d = h * w;

%% perform normalization

if n == 1
    dstimgs = normalize_page(imgs, d, mu, sigma);
else
    dstimgs = zeros(size(imgs));
    for i = 1 : n
        dstimgs(:,:,i) = normalize_page(imgs(:,:,i), d, mu, sigma);
    end
end

function dstimg = normalize_page(img, d, mu, sigma)

curimg = img(:);

% compute current mean value
cur_mv = sum(curimg) / d;

% shift to zero mean
curimg = curimg - cur_mv;

% compute current standard deviation
cur_std = norm(curimg) / sqrt(d);

% normalize to specified std dev
k = sigma / cur_std;
curimg = curimg * k;

% shift to specified mean
if mu ~= 0
    curimg = curimg + mu;
end

% reshape back to origin shape
dstimg = reshape(curimg, size(img));

author: lmthang


text:
function [el, prj] = project_elec(elc, pnt, tri)

% PROJECT_ELEC projects electrodes on a triangulated surface
% and returns triangle index, la/mu parameters and distance
%
% Use as
%   [el, prj] = project_elec(elc, pnt, tri)
% which returns
%   el    = Nx4 matrix with [tri, la, mu, dist] for each electrode
%   prj   = Nx3 matrix with the projected electrode position
%
% See also TRANSFER_ELEC

% Copyright (C) 1999-2013, 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$

Nelc = size(elc,1);
el   = zeros(Nelc, 4);

% this is a work-around for http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2369
elc = double(elc);
pnt = double(pnt);
tri = double(tri);

for i=1:Nelc
  [proj,dist] = ptriprojn(pnt(tri(:,1),:), pnt(tri(:,2),:), pnt(tri(:,3),:), elc(i,:), 1);

  [mindist, minindx] = min(abs(dist));
  [la, mu] = lmoutr(pnt(tri(minindx,1),:), pnt(tri(minindx,2),:), pnt(tri(minindx,3),:), proj(minindx,:));
  smallest_dist = dist(minindx);
  smallest_tri  = minindx;
  smallest_la   = la;
  smallest_mu   = mu;

  % the following can be done faster, because the smallest_dist can be
  % directly selected
%   Ntri = size(tri,1);
%   for j=1:Ntri
%     %[proj, dist] = ptriproj(pnt(tri(j,1),:), pnt(tri(j,2),:), pnt(tri(j,3),:), elc(i,:), 1);
%     if dist(j)<smallest_dist
%       % remember the triangle index, distance and la/mu
%       [la, mu] = lmoutr(pnt(tri(j,1),:), pnt(tri(j,2),:), pnt(tri(j,3),:), proj(j,:));
%       smallest_dist = dist(j);
%       smallest_tri  = j;
%       smallest_la   = la;
%       smallest_mu   = mu;
%     end
%   end

  % store the projection for this electrode
  el(i,:) = [smallest_tri smallest_la smallest_mu smallest_dist];
end

if nargout>1
  prj = zeros(size(elc));
  for i=1:Nelc
    v1 = pnt(tri(el(i,1),1),:);
    v2 = pnt(tri(el(i,1),2),:);
    v3 = pnt(tri(el(i,1),3),:);
    la = el(i,2);
    mu = el(i,3);
    prj(i,:) = routlm(v1, v2, v3, la, mu);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION this is an alternative implementation that will also work for
% polygons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function prj = polyproj(elc,pnt)
% projects a point on a plane, e.g. an electrode on a polygon
% pnt is a Nx3 matrix with multiple vertices that span the plane
% these vertices can be slightly off the plane
center = mean(pnt,1);
% shift the vertices to have zero mean
pnt(:,1) = pnt(:,1) - center(1);
pnt(:,2) = pnt(:,2) - center(2);
pnt(:,3) = pnt(:,3) - center(3);
elc(:,1) = elc(:,1) - center(1);
elc(:,2) = elc(:,2) - center(2);
elc(:,3) = elc(:,3) - center(3);
pnt = pnt';
elc = elc';
[u, s, v] = svd(pnt);
% The vertices are assumed to ly in plane, at least reasonably. That means
% that from the three eigenvectors there is one which is very small, i.e.
% the one orthogonal to the plane. Project the electrodes along that
% direction.
u(:,3) = 0;
prj = u * u' * elc;
prj = prj';
prj(:,1) = prj(:,1) + center(1);
prj(:,2) = prj(:,2) + center(2);
prj(:,3) = prj(:,3) + center(3);


text:
function r = Tsqrt(x)

%function r = Tsqrt(x)
%
%The same as CHOL, only
%R*R' = X;
%instead of
%R'*R = X.
%
%See also CHOL.

r = chol(x)';


text:
% drawobject ... draws a (calibration) object of a specific type
%
% [fig] = drawobject(X, ctype, fig, {color})
%
% X ...... 3xn matrix containg the object's corner points
%          if X is 4xn only the first 3 rows are used
% ctype .. a string specifying the type of the object
%          according to the type the appropriate drawing functions is called
%          supported types are:
%          'cube'    - a cube
%          'octagon' - a planar octagon
%          'cloud'   - a point cloud
% fig .... figure handle
% color .. color of plotting; defaults to blue
%
% fig .... returns the figure handle
%
% $Id: drawobject.m,v 2.0 2003/06/19 12:07:12 svoboda Exp $

function [fig] = drawobject(X,ctype,fig,color)

if nargin < 4
  color = 'b';	% default color
end

switch ctype
case 'cube',
    drawcube(X,fig,color);
case 'octagon',
    drawoctagon(X,fig,color);
case 'cloud',
    drawcloud(X,fig,color);
otherwise,
    error('unknown object type: ', ctype)
end

return

author: strawlab

repo: MultiCamSelfCal


text:
% Publicly accessible dependent properties of an object.
%
% See also: meta.property

% Copyright 2010 Levente Hunyadi
function dependent = getdependentproperties(obj)

dependent = {};
if isstruct(obj)  % structures have no dependent properties
    return;
end

try
    clazz = metaclass(obj);
catch %#ok<CTCH>
    return;  % old-style class (i.e. not defined with the classdef keyword) have no dependent properties
end

k = 0;  % number of dependent properties found
n = numel(clazz.Properties);  % maximum number of properties
dependent = cell(n, 1);
for i = 1 : n
    property = clazz.Properties{i};
    if property.Abstract || property.Hidden || ~strcmp(property.GetAccess, 'public') || ~property.Dependent
        continue;  % skip abstract, hidden, inaccessible and independent properties
    end

    k = k + 1;
    dependent{k} = property.Name;
end
dependent(k+1:end) = [];  % drop unused cells

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function hole_num = p07_hole_num ( )

%*****************************************************************************80
%
%% P07_HOLE_NUM counts the holes in problem 07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer HOLE_NUM, the number of holes.
%
  hole_num = 0;

  return
end


text:
function [XS,XF,O,H,L,C,V] = filterdata_lookback(XS,XF,D,O,H,L,C,V)

ndays = length(D);

% Take out the days that are NaN
index_nonNaN_linear = find(~isnan(O));
index_nonNaN_sub = zeros(size(O));
index_nonNaN_sub(index_nonNaN_linear(:)) = 1; % matrix of values with non-NaN
ndays_valid = sum(index_nonNaN_sub,1);

% Now we filter by stocks, then by days
ndays_filter = mode(ndays_valid); % use the charts with the modal number of days
%index_stocks_valid = find(ndays_valid == ndays_filter); % filter out stocks with not many days
index_stocks_invalid = find(ndays_valid ~= ndays_filter); % filter out stocks with not many days

% Filter out stocks that don't have modal number of available days
%ntickers = size(XS(:,1),1);
O(:,index_stocks_invalid)=[]; H(:,index_stocks_invalid)=[]; L(:,index_stocks_invalid)=[];
C(:,index_stocks_invalid)=[]; V(:,index_stocks_invalid)=[];
XS(index_stocks_invalid,:)=[]; XF(index_stocks_invalid,:)=[];

fprintf('Filtering modal days... [%d days, %d stocks]\n',size(O,1),size(O,2))

% Now take out all rows with nans (corresponds to holidays and weekends)
[rows, columns] = find(isnan(O)); % finds the rows and columns of elements that have nan
index_weekendsholidays = unique(rows); % find the unique rows
O(index_weekendsholidays,:)=[]; H(index_weekendsholidays,:)=[]; L(index_weekendsholidays,:)=[];
C(index_weekendsholidays,:)=[]; V(index_weekendsholidays,:)=[];
fprintf('Filtering NaN days... [%d days, %d stocks]\n',size(O,1),size(O,2));

% Make sure there are no more NaNs
if( length(find(isnan(O)))>0 )
  error('NaNs are still in the data. Remove them!')
end
end

author: gudbrandtandberg


text:
function phaseChangeEnthalpyMethodExample( varargin )
%phaseChangeEnthalpyMethodExample() demonstrates the 2D phase change using
% the enthalpy method
%   In this example, the heat equation is solved on a 2D initially solid
%   rectangular domain. This file contains three test cases in total.
%   By setting COMPUTE_FLAG to either 1, 2 or 3, different kinds of phase
%   change problems for water-ice are solved. Depending on COMPUTE_FLAG,
%   the visualization is either a plot comparing the numerical solution
%   with the analytical solution, a full 2D visualization of the
%   temperature distribution or a full 2D visualization of the phase
%   distribution within the domain.
%
%   See also StefanProblemAnalyticalSolution.
%
%   Date:   27.09.2015
%   Author: Kai Schueller

close all; % closes all figures

addpath('Functions'); % add Functions folder to the search path

%% Set the flag for the computation:
% COMPUTE_FLAG=1; % Solves the two-phase Stefan problem and plots the
                  % solution together with the analytical solution.
                  % Please note: the analytical solution is for a
                  % semi-infinite domain, while the numerical solution is
                  % for a finite domain. Therefore, the difference between
                  % the solutions increases with increasing time at the
                  % right side of the domain
% COMPUTE_FLAG=2; % Same as COMPUTE_FLAG=1, but 2D visualization without
                  % comparison with analytical solution and Neumann BC
% COMPUTE_FLAG=3; % Same as COMPUTE_FLAG=2, but with convection in the
                  % direction of the heat source (i.e. in negative x) so
                  % that there will be a steady state after some time
COMPUTE_FLAG=1;

%% User Input
% Create the mesh
m=createMesh2D(50,3,2,2);
 visualizeMesh2D(m); % uncomment if mesh should be displayed

% Thermophysical properties of the phase change material (PCM)
rho_L=1000;         % density of the liquid phase [kg/m^3]
rho_S=rho_L;%920;   % density of the solid phase [kg/m^3]
h_melt=333400;      % latent heat of melting [J/kg]
T_m=0;              % melting temperature [degC]
T_L=0.01;           % liquidus temperature [degC]
T_S=-0.01;          % solidus temperature [degC]
k_L=0.6;            % thermal conductivity at liquid phase [W/(mK)]
k_S=2.3;            % thermal conductivity at solid phase [W/(mK)]
c_L=4200;           % Heat capacity liquid Phase [J/(kgK)]
c_S=2000;           % Heat capacity of the solid phase [J/(kgK)]
T_init = -5;

% Create the boundary condition structure
BC = createBC(m);

if COMPUTE_FLAG==1
    % We will later extract a 1D solution, so there is no need for many
    % y-cells
    m=createMesh2D(m.dims(1),2,m.facecenters.x(end),m.facecenters.y(end));
    BC = createBC(m);

    u = [0,0];
    T_Solid=-5;     % initial temperature of the solid phase [degC]
    T_Liquid=5;     % initial temperature of the liquid phase [degC]
    analSol.x=0:0.01:2; % x positions for analytical solution

    % Assign boundary conditions
    BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=T_Liquid;   % Dirichlet for the left boundary
    BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=T_Solid; % Dirichlet for the right boundary
end
if COMPUTE_FLAG==2
    u = [0,0];
    q=-1000; % heat flux [W/m^2]
    % Assign boundary conditions
    BC.left.a(:) = 1; BC.left.b(:)=0; BC.left.c(:)=q;     % Dirichlet for the left boundary
    BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=-5; % Dirichlet for the right boundary
end
if COMPUTE_FLAG==3
    u = [-0.005,0];      % velocity of the phase change material in x and y
    q=-1000; % heat flux [W/m^2]
    % Assign boundary conditions
    BC.left.a(:) = 1; BC.left.b(:)=0; BC.left.c(:)=q;     % Dirichlet for the left boundary
    BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=-5; % Dirichlet for the right boundary
end

% Solver settings
dt = 36000; % time step size
final_t = 2000000; % final time
convergence_tolerance=1e-8;
relaxation_coeff=0.7;

%% Set initial values
liquidFraction=createCellVariable(m,0);   % initialize liquid fraction,
                                          % 0 = initially solid
Temp = createCellVariable(m, T_init, BC); % initial temperature

RHS_liquidFraction=constantSourceTerm(liquidFraction);

%% main loop

% we use constant boundary conditions, so there is no need to change them
% within the main loop
[Mbc, RHSbc] = boundaryCondition2D(BC);

% the convection velocity does also not change
u_face = createFaceVariable(m, u);

for t=dt:dt:final_t
    % The transient term is calculated outside of the liquid fraction
    % update loop for each time step. So the liquid fraction of the last
    % time step is used.
    rho_c_mix=createMixCellVar(m,liquidFraction,rho_L.*c_L,rho_S.*c_S);
    [M_trans, RHS_trans] = transientTerm(Temp, 1, rho_c_mix);

    RHS_liquidFraction_old=RHS_liquidFraction;

    iterations=0; % set the iteration counter to zero

    % set the error higher than tolerance to enter the while-loop for the
    % first time
    error=convergence_tolerance+1;

    while error>convergence_tolerance || iterations>10000
        iterations=iterations+1;

        % create mixture values
        k_mix=createMixCellVar(m,liquidFraction,k_L,k_S);
        k_mix_face = harmonicMean(k_mix);
        rho_mix=createMixCellVar(m,liquidFraction,rho_L,rho_S);
        rho_mix_face = harmonicMean(rho_mix);
        c_mix=createMixCellVar(m,liquidFraction,c_L,c_S);
        u_rho_mix_face=u_face;
        u_rho_mix_face.xvalue=rho_mix_face.xvalue.*u_face.xvalue;
        u_rho_mix_face.yvalue=rho_mix_face.yvalue.*u_face.yvalue;

        % calculate the matrix M
        Mdiff = dt*diffusionTerm(k_mix_face);
        Mconv =  dt*convectionTerm(u_rho_mix_face);
        M = M_trans-Mdiff+Mbc+Mconv;

        % calculate the right hand side vector
        RHS_gammaTerm=constantSourceTerm(rho_mix).*h_melt.*(RHS_liquidFraction_old-RHS_liquidFraction);
        RHS = RHS_trans+RHSbc+RHS_gammaTerm;

        Temp_old=Temp; % store the old temperature distribution

        % solve the heat equation to obtain the new temperature
        Temp = solvePDE(m,M, RHS);

        % calculate maximum error between current and previous temperatures
        error=max(max(abs(Temp.value-Temp_old.value)));

        % update the liquid fraction based on the result of heat eq.
        liquidFraction=updateLiquidFraction(liquidFraction,...
                                                Temp,...
                                                c_mix,...
                                                h_melt,...
                                                T_L,...
                                                T_S,...
                                                relaxation_coeff);

        RHS_liquidFraction=constantSourceTerm(liquidFraction);
    end

    if COMPUTE_FLAG==1

        % extract 1D data from solution for the plot
        numSol.T=liquidFraction.value(:,3);

        % now interpolate over the ghost cell and first internal cell
        values_temp(1)=(numSol.T(1)+numSol.T(2))/2;
        values_temp(2)=(numSol.T(end)+numSol.T(end-1))/2;
        numSol.T=[values_temp(1); numSol.T(2:end-1); values_temp(2)];

        % extracs the x-positions
        numSol.x = [m.facecenters.x(1); m.cellcenters.x; m.facecenters.x(end)];

        % calculate analytical solution
        [analSol.T,analSol.InterfacePos] = StefanProblemAnalyticalSolution(...
                                             T_Liquid,...
                                             T_Solid,...
                                             T_m,...
                                             h_melt,...
                                             rho_L,...
                                             k_L,...
                                             k_S,...
                                             c_L,...
                                             c_S,...
                                             analSol.x,...
                                             t);

        % plot the analytical and numerical solution
        if t==dt
            h_fig=figure('NumberTitle','Off');
            set(h_fig,'Name',sprintf('Stefan Problem: %d s (%.2f h)',t,t/3600));
            h1=plot(numSol.x,numSol.T,'-x','color',[0 0 1]);
            hold on;
            %h2=plot(analSol.x,analSol.T,'-','color',[1 0 0]);
            h3=plot([analSol.InterfacePos analSol.InterfacePos],[0 1],'--','color',[0 0 0]);
            legend([h1 h3],'Numerical Solution: Liquid fraction','Analytical Solution: Phase Interface');
            xlabel('x [m]');
            ylabel('Liquid fraction');
            hold off;
        else
            set(h1,'XData',numSol.x,'YData',numSol.T);
            %set(h2,'XData',analSol.x,'YData',analSol.T);
            set(h3,'XData',[analSol.InterfacePos analSol.InterfacePos]);
            set(h_fig,'Name',sprintf('Stefan Problem: %d s (%.2f h)',t,t/3600));
        end
        xlim([0 0.5]);
        saveas(h_fig,['figureSaveGamma/fig_',num2str(t),'.eps'],'epsc');
        drawnow;
    elseif COMPUTE_FLAG==2 || COMPUTE_FLAG==3
        if COMPUTE_FLAG==2
            visualizeCells(Temp);
            shading interp;
        else
            visualizeCells(liquidFraction);
        end

        colormap(jet);
        drawnow;
    end
end

end

%% Additional Functions

% createMixCellVar() calculates the mixture value of a given property by
% using the liquid fraction of the mixture
function mixCellVar=createMixCellVar(m,liquidFraction,var_L,var_S)
    mixCellVar=liquidFraction.value*var_L+(1-liquidFraction.value)*var_S;
    mixCellVar=createCellVariable(m,mixCellVar(2:end-1,2:end-1));
end

% updateGamma() calculates a new liquid fraction that better fits the
% mixture enthalpy
function helpVar = updateLiquidFraction( liquidFraction_old,T,c,h_melt,T_L,T_S,relaxation_coefficient )

    T_corr=liquidFraction_old.value*(T_L-T_S)+T_S;

    liquidFraction_new=min(max(liquidFraction_old.value+...
        relaxation_coefficient*c.value./h_melt.*(T.value-T_corr),0),1);

    % set liquid fraction of corner cells to zero
    liquidFraction_new([1 end],[1 end])=0;
    helpVar.value=liquidFraction_new;
    helpVar.domain=liquidFraction_old.domain;

end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function data = read_asa_msr(fn);

% READ_ASA_MSR reads EEG or MEG data from an ASA data file
% converting the units to uV or fT

% Copyright (C) 2002, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_asa_msr.m 945 2010-04-21 17:41:20Z roboos $

Npnt      = read_asa(fn, 'NumberPositions=', '%d');
Ntime     = read_asa(fn, 'NumberTimesteps=', '%d');
UnitT     = read_asa(fn, 'UnitTime', '%s');
UnitM     = read_asa(fn, 'UnitMeas', '%s');
Timesteps = read_asa(fn, 'Timesteps', '%s');
lab       = read_asa(fn, 'Labels', '%s', Npnt);

val = read_asa(fn, 'Values', '%f');
if any(size(val)~=[Npnt,Ntime])
  msm_file = read_asa(fn, 'Values', '%s');
  [path, name, ext] = fileparts(fn);
  fid = fopen(fullfile(path, msm_file), 'rb', 'ieee-le');
  val = fread(fid, [Ntime, Npnt], 'float32')';
  fclose(fid);
end

tmp = sscanf(Timesteps, '%f(%f)%f');
time = linspace(tmp(1), tmp(3), Ntime);

if strcmpi(UnitT,'ms')
  time = 1*time;
elseif strcmpi(UnitT,'s')
  time = 1000*time;
elseif ~isempty(UnitT)
  error(sprintf('Unknown unit of time (%s)', UnitT));
end

if strcmpi(UnitM,'uv')
  val = 1*val;
elseif strcmpi(UnitM,'?v')
  val = 1*val;
elseif strcmpi(UnitM,'mv')
  val = 1000*val;
elseif strcmpi(UnitM,'v')
  val = 1000000*val;
elseif strcmpi(UnitM,'ft')
  val = 1*val;
elseif strcmpi(UnitM,'pt')
  val = 1000*val;
elseif ~isempty(UnitM)
  error(sprintf('Unknown unit of measurement (%s)', UnitM));
end

if length(size(lab))==2
  lab = tokenize(lab{1});
end

data.time  = time;
data.data  = val;
data.label = lab;

author: PatternRecognition

repo: OpenBMI


text:
function RDr=scd_model_GPD_RDr(d,Delta,delta,Dr)
% RDr=scd_model_GPD_RDr(d,Delta,delta,Dr)

alpha=4*delta*Dr/d^2; beta=4*Delta*Dr/d^2;
a=[1.84118378134065	5.33144277352503	8.53631636634628	11.7060049025920	14.8635886339090].^2;% roots of bessel J1'
td=Delta-delta/3;

k=zeros(length(delta),5);
for m=1:5;
    numerator=2*alpha*a(m) - 2 + 2*exp(-alpha*a(m)) + (2-exp(alpha*a(m))-exp(-alpha*a(m))).*exp(-beta*a(m));
    denominator=alpha.^2*a(m).^3.*(a(m)-1);
    k(:,m)=numerator(:)./denominator(:);
end
k(isnan(k))=0;
RDr=sum(k,2).*d.^2./(2*td);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function title = p03_title ( )

%*****************************************************************************80
%
%% P03_TITLE returns the title of problem p03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Corner Peak';

  return
end


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

figure; % open a new figure window

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

plot(x, y, 'rx', 'MarkerSize', 10); % Plot the data
ylabel('Profit in $10,000s'); % Set the y􀀀axis label
xlabel('Population of City in 10,000s'); % Set the x􀀀axis label

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

end

author: vkosuri

repo: CourseraMachineLearning


text:
function value = bmi_english ( w_lb, h_ft, h_in )

%*****************************************************************************80
%
%% BMI_ENGLISH computes the body mass index given English measurements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real W_LB, the body weight in pounds.
%
%    Input, real H_FT, H_IN, the body height in feet and inches
%
%    Output, real VALUE, the body mass index.
%
  w_kg = pounds_to_kilograms ( w_lb );

  h_m = feet_to_meters ( h_ft + ( h_in / 12.0 ) );

  value = bmi_metric ( w_kg, h_m );

  return
end


text:
function [n maxLflag] = loadfactor(tState)

if tState.omega >= 0 % if level flight or sustained turn
    n = sqrt((tState.omega*tState.V/g0)^2+1);
    maxLflag = false;
elseif tState.omega < 0 % if instantaneous turn indicated by negative omega
    % NB: rotor itertia doesn't help an inst. turn
    n = abs(tState.omega).*tState.V/g0;
    maxLflag = true;
else
    disp('Ill-defined turning condition')
end

end

author: Sable


text:
function display(m,varargin)
% standard output

if check_option(varargin,'onlyShowMiller')
  eps = 1e4;
  % extract coordinates in the correct form
  d = round(m.coordinates * eps)./eps;
  % set up coordinate names
  columnNames = vec2cell(char(m.dispStyle));
  cprintf(d,'-L','  ','-Lc',columnNames);
  return
end

displayClass(m,inputname(1),'moreInfo',char(m.CS,'compact'),varargin{:});

display@vector3d(m,'skipHeader', 'skipCoordinates');

% display coordinates
if length(m) < 20 && ~isempty(m)

  display(m,'onlyShowMiller')

elseif ~getMTEXpref('generatingHelpMode') && ~isempty(m)

  disp(' ')
  s = setappdata(0,'data2beDisplayed',m);
  disp(['  show Miller'])
  disp(' ')

end

author: mtex-toolbox

repo: mtex

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


text:
function output = callcbc(interfacedata)

% Standard input interface
model = yalmip2cbc(interfacedata);

options = interfacedata.options;
model.opts = options.cbc;
model.opts.display = options.verbose;

if options.savedebug
    save cbcdebug model
end

showprogress('Calling CBC',options.showprogress);
solvertime = tic;
if nnz(model.H)==0
    [x,fval,exitflag,nodes,xc] = cbc([],model.f, model.A, model.rl, model.ru, model.lb, model.ub, model.xtype,model.sos,model.x0,model.opts);
else
    [x,fval,exitflag,nodes,xc] = cbc(model.H,model.f, model.A, model.rl, model.ru, model.lb, model.ub, model.xtype,model.sos,model.x0,model.opts);
end
solvertime = toc(solvertime);

% No duals
D_struc = [];
if isempty(x)
    x = zeros(length(c),1);
end

% Check, currently not exhaustive...
problem = 0;
switch exitflag
    case {0,2,6}
        problem = 0;
    case {1,8}
        problem = 1;
    case {3,4}
        problem = 3;
    case 5
        problem = 16;
    case 7
        problem = 2;
    otherwise
        problem = -1;
end

infostr = yalmiperror(problem,'CBC');

% Save all data sent to solver?
if options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fval = fval;
    solveroutput.exitflag = exitflag;
    solveroutput.nodes = nodes;
    solveroutput.xc = xc;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x(:),[],[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% std_readitc()  - load ITC measures for data channels or
%                  for all components of a specified cluster.
% Usage:
%         >> [STUDY, itcdata, times, freqs] = ...
%                   std_readitc(STUDY, ALLEEG, varargin);
%
% Note: this function is a helper function that contains a call to the
% std_readersp function that reads all 2-D data matrices for EEGLAB STUDY.
% See the std_readersp help message for more information.
%
% Author: Arnaud Delorme, CERCO, 2006-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2004, arno@sccn.ucsd.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 [STUDY, erspdata, alltimes, allfreqs, erspbase] = std_readersp(STUDY, ALLEEG, varargin);

[STUDY, erspdata, alltimes, allfreqs] = std_readersp(STUDY, ALLEEG, 'infotype','itc', varargin{:});

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function G_vctr_sym = G_vctr_fcn(in1,in2)
%G_VCTR_FCN
%    G_VCTR_SYM = G_VCTR_FCN(IN1,IN2)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    19-Oct-2021 16:04:40

pi171 = in2(17,:);
pi191 = in2(19,:);
pi271 = in2(27,:);
pi291 = in2(29,:);
pi301 = in2(30,:);
pi371 = in2(37,:);
pi391 = in2(39,:);
pi401 = in2(40,:);
pi471 = in2(47,:);
pi481 = in2(48,:);
pi491 = in2(49,:);
pi501 = in2(50,:);
pi571 = in2(57,:);
pi581 = in2(58,:);
pi591 = in2(59,:);
pi601 = in2(60,:);
q2 = in1(2,:);
q3 = in1(3,:);
q4 = in1(4,:);
q5 = in1(5,:);
q6 = in1(6,:);
t2 = sin(1.57079632679);
t3 = cos(q3);
t4 = sin(q2);
t5 = cos(q2);
t6 = sin(q3);
t7 = t2.*t4.*(1.81e2./1.0e3);
t8 = t7-6.13e2./1.0e3;
t9 = t2.*t3.*t5.*(9.81e2./1.0e2);
t15 = t2.*t4.*t6.*(9.81e2./1.0e2);
t10 = t9-t15;
t11 = t2.*t5.*t6.*(1.81e2./1.0e3);
t12 = t2.*t5.*t6.*(9.81e2./1.0e2);
t13 = t2.*t3.*t4.*(9.81e2./1.0e2);
t14 = t12+t13;
t16 = cos(q4);
t17 = sin(q4);
t18 = t2.*t14.*t16;
t19 = t2.*t10.*t17;
t20 = t18+t19;
t21 = t3.*t8;
t22 = t6.*t8;
t23 = t2.*t3.*t5.*(1.81e2./1.0e3);
t24 = t22-t23;
t25 = t11+t21-5.71e2./1.0e3;
t26 = cos(q6);
t27 = cos(q5);
t28 = sin(q6);
t29 = t2.*t10.*t16;
t40 = t2.*t14.*t17;
t30 = t29-t40;
t31 = sin(q5);
t32 = t31.*(8.7e1./5.0e2);
t33 = t2.*t16.*t25;
t50 = t2.*t17.*t24;
t34 = t33-t50;
t35 = t27.*t34;
t36 = t32+t35;
t37 = t2.*t16.*t24;
t38 = t2.*t17.*t25;
t39 = t37+t38+3.0./2.5e1;
t41 = t2.*t3.*t4.*(1.81e2./1.0e3);
t42 = t11+t41;
t44 = t2.*t4.*t6.*(1.81e2./1.0e3);
t43 = t23-t44;
t45 = t2.*t16.*t42;
t46 = t2.*t17.*t43;
t47 = t45+t46;
t48 = t2.*t16.*t43;
t51 = t2.*t17.*t42;
t49 = t48-t51;
t52 = t27.*(8.7e1./5.0e2);
t68 = t31.*t34;
t53 = t52-t68;
t54 = t30.*t31.*t53;
t55 = t31.^2;
t56 = t20.*t28;
t72 = t26.*t27.*t30;
t57 = t56-t72;
t58 = t20.*t26;
t59 = t27.*t28.*t30;
t60 = t58+t59;
t61 = t37+t38;
t62 = t20.*t61;
t63 = t11+t21;
t64 = t27.^2;
t65 = t2.*t17.*t63;
t66 = t37+t65;
t67 = t27.*t30.*t36;
t69 = pi371.*t30;
t70 = pi391.*t20;
t71 = pi491.*t20;
t73 = pi591.*t60;
t74 = pi271.*t14;
t77 = t2.*t16.*t63;
t75 = t50-t77;
t76 = t30.*t34;
t78 = t28.*t30;
t79 = t20.*t26.*t27;
t80 = t78+t79;
t81 = t26.*t30;
t91 = t20.*t27.*t28;
t82 = t81-t91;
t83 = t28.*t36;
t84 = t26.*t39;
t85 = t83+t84;
t86 = t26.*t36;
t92 = t28.*t39;
t87 = t86-t92;
t88 = t57.*t87;
t89 = t20.*t55.*t66;
t90 = pi471.*t27.*t30;
t93 = t20.*t55.*t61;
G_vctr_sym = [0.0;t69+t70+t71+t73+t74+t90-pi291.*t10-pi571.*t57+pi301.*(-t10.*t42+t14.*t43+t10.*t63+t14.*(t22-t2.*t3.*t5.*(1.81e2./1.0e3)))-pi601.*(t54+t88+t80.*(t28.*t47-t26.*t27.*t49)+t82.*(t26.*t47+t27.*t28.*t49)-t60.*t85-t20.*t49.*t55)+pi401.*(t62+t76+t20.*t49-t30.*t47)+pi501.*(-t54+t62+t67-t30.*t47+t20.*t49.*t55+t20.*t49.*t64)+pi171.*t2.*t4.*(9.81e2./1.0e2)-pi191.*t2.*t5.*(9.81e2./1.0e2)-pi481.*t30.*t31-pi581.*t30.*t31;t69+t70+t71+t73+t74+t90-pi601.*(t54+t88+t89-t80.*(t28.*t75-t26.*t27.*t66)-t82.*(t26.*t75+t27.*t28.*t66)-t60.*t85)-pi291.*t10-pi571.*t57+pi401.*(t62+t76-t20.*t66+t30.*t75)-pi501.*(t54-t62-t67+t89-t30.*t75+t20.*t64.*t66)-pi481.*t30.*t31-pi581.*t30.*t31;t69+t70+t71+t73+t90-pi601.*(t54+t88+t93+t80.*(t28.*t34+t26.*t27.*t61)+t82.*(t26.*t34-t27.*t28.*t61)-t60.*t85)-pi571.*t57-pi501.*(t54-t62-t67+t76+t93+t20.*t61.*t64)-pi481.*t30.*t31-pi581.*t30.*t31;-pi601.*(-t20.*t31.*t36+t20.*t27.*t53-t26.*t53.*t80+t28.*t53.*t82+t20.*t26.*t31.*t87+t20.*t28.*t31.*t85)-pi471.*t20.*t31-pi481.*t20.*t27-pi581.*t20.*t27-pi571.*t20.*t26.*t31-pi591.*t20.*t28.*t31;pi571.*t82+pi591.*t80];

author: shamilmamedov


text:
function [ error_per_image ] = compute_eye_corner_error( ground_truth_all, detected_points_all )
num_of_images = size(ground_truth_all,3);
num_of_points = size(ground_truth_all,1);

error_per_image = zeros(num_of_images,1);

for i =1:num_of_images
    detected_points      = detected_points_all(:,:,i);
    ground_truth_points  = ground_truth_all(:,:,i);
    interocular_distance = norm(ground_truth_points(37,:)-ground_truth_points(46,:));
    %interocular_distance = norm(mean(ground_truth_points(37:42,:))-mean(ground_truth_points(43:48,:)));
    sum=0;
    for j=1:num_of_points
        sum = sum+norm(detected_points(j,:)-ground_truth_points(j,:));
    end
    error_per_image(i) = sum/(num_of_points*interocular_distance);
end
end

author: jiankangdeng

repo: MenpoBenchmark


text:
% Helper function for decom.m, run DECOM method at regular time intervals.
%
% Octave compatible
%
% (please see decom.m for documentation and HOWTO_egg.m for a usage example)
%
%
% Copyright (c) 2013 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [Oqeggmess, f0mess, npicferms, npicouvers, atimes] = decom_sig(s, fs, f0s, method)

    if nargin<3; f0s=[]; f0=[]; end
    if nargin<4; method = []; end

    winLen=2*round(20/1000*fs /2)+1;
    winShift=round(5/1000*fs);

    % Do processing
    start=1;
    stop=start+winLen-1;
    cnt=1;

    while stop <= length(s)

        % Get windowed frame
        x_frame = s(start:stop);
        atimes(cnt) = (0.5*(stop+start))/fs;

        if ~isempty(f0s) && size(f0s,2)>1
            f0 = interp1(f0s(:,1), f0s(:,2), atimes(cnt), 'nearest', 'extrap');
        end

        [Oqeggmess(cnt), f0mess(cnt), npicferms(cnt), npicouvers(cnt)] = decom(x_frame, fs, f0, method);

        % Increment
        start=start+winShift-1;
        stop=start+winLen-1;
        cnt=cnt+1;
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function s=dUbox(pred,gt)
dX=dXbox(pred);
dI=dIbox(pred,gt);
s.dt=dX.dt-dI.dt;
s.db=dX.db-dI.db;
s.dl=dX.dl-dI.dl;
s.dr=dX.dr-dI.dr;

author: Zzh-tju

repo: DIoU

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


text:
function [ff,f]=lpcar2ff(ar,np)
%LPCAR2FF LPC: Convert AR coefs to complex spectrum FF=(AR,NP)
%
%  Inputs: ar(nf,n)     AR coefficients, one frame per row
%          np           Size of output spectrum is np+1 [n]
%
% Outputs: ff(nf,np+1)  Complex spectrum from DC to Nyquist
%          f(1,np+1)    Normalized frequencies (0 to 0.5)
%
% For high speed make np equal to a power of 2

%      Copyright (C) Mike Brookes 1998-2014
%      Version: $Id: lpcar2ff.m 5026 2014-08-22 17:47:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,p1]=size(ar);
if nargin<2
    if nargout
        np=p1-1;
    else
        np=128;
    end
end
ff=(rfft(ar.',2*np).').^(-1);
f=(0:np)/(2*np);
if ~nargout
    subplot(2,1,2);
    plot(f,unwrap(angle(ff)));
    xlabel('Normalized frequency f/f_s');
    ylabel('Phase (rad)');
    subplot(2,1,1);
    plot(f,db(abs(ff)));
    xlabel('Normalized frequency f/f_s');
    ylabel('Gain (dB)');
    title('LPC Spectrum');
end

author: jtkim-kaist

repo: VAD

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


text:
%% Copyright (C) 2014-2018 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

function convert_comments (basedir, subdir, dirout)
% this slightly strange way of doing things (basedir, subdir) is
% b/c I must "chdir" into base, but get_first_help_sentence() must
% not be in the class dir...

  %basedir, subdir, dirout
  files = dir([basedir subdir]);
  chdir(basedir)

  for i=1:length(files)
    if (~files(i).isdir)
      [dir, name, ext] = fileparts(files(i).name);
      if (strcmp(ext, '.m'))
        if isempty(subdir)
          octname = [name ext];
        else
          octname = [subdir '/' name ext];
        end
        fprintf('Converting texinfo to Matlab-style documentation: %s\n', octname)
        r = convert_oct_2_ml (octname, [dirout octname]);
        if ~r
          [status, msg, msgid] = copyfile (octname, [dirout octname], 'f');
          if (status ~= 1)
            error(msg)
          end
          fprintf('**** COPYING %s UNMODIFIED ****\n', octname)
        end
      end
    end
  end
end

function success = convert_oct_2_ml (fname, foutname)

  [dir, fcn, ext] = fileparts(fname);

  newl = sprintf('\n');

  [fi,msg] = fopen(fname, 'r');
  if (fi < 0)
    error(msg)
  end

  ins = {}; i = 0;
  while (1)
    temp = fgets(fi);
    if ~ischar(temp) && temp == -1
      break
    end
    i = i + 1;
    ins{i} = temp;
    % todo, possible strip newl
  end

  fclose(fi);

  % trim newlines
  ins = deblank(ins);

  %% find the actual function [] = ... line
  Nfcn = [];
  for i = 1:length(ins)
    I = strfind (ins{i}, 'function');
    if ~isempty(I) && I(1) == 1
      %disp ('found function header')
      Nfcn = i;
      break
    end
  end
  if isempty(Nfcn)
    disp('AFAICT, this is a script, not a function')
    success = false;
    return
  end

  %% copyright block
  [cr,N] = findblock(ins, 1);
  if (Nfcn < N)
    warning('function header in first block (where copyright block should be), not converting')
    success = false;
    return
  end
  cr = ltrim(cr, 3);

  % cut 2nd line if empty
  if isempty(cr{2})
    cr2 = cell(1,length(cr)-1);
    cr2(1) = cr(1);
    cr2(2:end) = cr(3:end);
    cr = cr2;
  end

  cr = prepend_each_line(cr, '%', ' ');
  cr{1} = ['%' cr{1}];
  copyright_summary = 'This is free software under the GPL, see .m file for full details.';
  % alternative for BSD:
  %copyright_summary = 'This is free software, see .m file for license.';

  %% use block
  % we don't parse this, just call get_help_text
  temp = ins{N};
  if ~strcmp(temp, '%% -*- texinfo -*-')
    error('can''t find the texinfo line, aborting')
    %success = false;
    %return
  end

  %% the "lookfor" line
  lookforstr = get_first_help_sentence (fname);
  if (~isempty(strfind(lookforstr, newl)))
    lookforstr
    error('lookfor string contains newline: missing period? too long? some other issue?')
    %success = false;
    %return
  end
  if (length(lookforstr) > 76)
    error(sprintf('lookfor string of length %d deemed too long', length(lookforstr)))
  end

  %% get the texinfo source, and format it
  [text, form] = get_help_text(fname);
  if ~strcmp(form, 'texinfo')
    text
    form
    error('formatted incorrectly, help text not texinfo')
  end
  usestr = __makeinfo__(text, 'plain text');

  %% remove the lookforstr from the text
  I = strfind(usestr, lookforstr);
  if length(I) ~= 1
    I
    lookforstr
    usestr
    error('too many lookfor lines?')
  end
  len = length(lookforstr);
  J = I + len;

  % if usestr has only a lookfor line then no need to see what's next
  if (J < length(usestr))
    % find next non-empty char
    %while isspace(usestr(J))
    %  J = J + 1;
    %end

    % let's be more conservative trim newline in usual case:
    if ~isspace(usestr(J))
      error('no space or newline after lookfor line?');
    end
    J = J + 1;
  end

  usestr = usestr([1:(I-1) J:end]);

  use = strsplit(usestr, newl, 'CollapseDelimiters', false);

  %% 2017: we have more choices than "-- Function File"
  % just trim "--" so it doesn't break indents
  def_line_has_function_name = false;
  for i=1:length(use)
    if regexp(use{i}, '^ -- ')
      if (~ isempty(strfind(use{i}, [' ' fcn])))
        def_line_has_function_name = true;
      end
    end
    use{i} = regexprep(use{i}, '^ -- ', '     ');
  end
  %usestr = strrep(usestr, lookforstr, '');
  if (~ def_line_has_function_name)
    warning('function @deftypefn line may not include function name:')
    use{i}
  end

  use = ltrim(use, 2);
  while isempty(use{end})
    use = use(1:end-1);
  end

  %% the rest
  N = Nfcn;
  fcn_line = ins{N};

  % sanity checks
  I = strfind(ins{N+1}, '%');
  if ~isempty(I) && I(1) == 1
    ins{N}
    ins{N+1}
    error('possible duplicate comment header following function')
  end

  therest = ins(N+1:end);

  %% Output
  [f, msg] = fopen(foutname, 'w');
  if (f < 0)
    msg
    error('failed to open file for writing')
  end

  fdisp(f, fcn_line)

  fprintf(f, '%%%s  %s\n', upper(fcn), lookforstr)

  for i=1:length(use)
    fprintf(f, '%%%s\n', use{i});
  end

  fdisp(f, '%');
  fprintf(f, '%%   %s\n', copyright_summary);

  %fdisp(f, '%');
  %fdisp(f, '%   [Genereated from a GNU Octave .m file, edit that instead.]');

  %fprintf(f,(s)

  fdisp(f, '');
  fdisp(f, '%% Note for developers');
  fdisp(f, '% This file is autogenerated from a GNU Octave .m file.');
  fdisp(f, '% If you want to edit, please make changes to the original instead');

  fdisp(f, '');
  for i=1:length(cr)
    fprintf(f, '%s\n', cr{i});
  end

  fdisp(f, '');

  for i=1:length(therest)
    fprintf(f, '%s\n', therest{i});
  end

  fclose(f);

  success = true;

end

function [block,endl] = findblock(f, j)
  block = {}; c = 0;
  %newl = sprintf('\n');
  for i = j:length(f)
    temp = f{i};
    %if (strcmp(temp, newl))
    if (isempty(temp))
      endl = i + 1;
      break
    end
    c = c + 1;
    block{c} = temp;
  end
end

function g = ltrim(f, n)
  g = {};
  for i = 1:length(f)
    temp = f{i};
    if length(temp) < n
      g{i} = '';
    else
      g{i} = substr(temp, n+1);
    end
  end
end

function g = prepend_each_line(f, pre, pad)
  g = {};
  for i = 1:length(f)
    temp = f{i};
    if isempty(temp)
      g{i} = pre;
    else
      g{i} = [pre pad temp];
    end
  end
end

repo: octsympy


text:
classdef Slice < dagnn.ElementWise
  properties
    dim = 3
    slicePoint = 3
  end

  properties (Transient)
    inputSizes = {}
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs = vl_nnslice(inputs{1}, obj.dim, obj.slicePoint, []) ;
      obj.inputSizes = cellfun(@size, inputs, 'UniformOutput', false) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs = vl_nnslice(inputs, obj.dim, obj.slicePoint, ...
                            derOutputs, 'inputSizes', obj.inputSizes) ;
      derParams = {} ;
    end

    function reset(obj)
      obj.inputSizes = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      sz = inputSizes{1} ;
      for k = 2:numel(inputSizes)
        sz(obj.dim) = sz(obj.dim) + inputSizes{k}(obj.dim) ;
      end
      outputSizes{1} = sz ;
    end

    function rfs = getReceptiveFields(obj)
      numInputs = numel(obj.net.layers(obj.layerIndex).inputs) ;
      if obj.dim == 3 || obj.dim == 4
        rfs = getReceptiveFields@dagnn.ElementWise(obj) ;
        rfs = repmat(rfs, numInputs, 1) ;
      else
        for i = 1:numInputs
          rfs(i,1).size = [NaN NaN] ;
          rfs(i,1).stride = [NaN NaN] ;
          rfs(i,1).offset = [NaN NaN] ;
        end
      end
    end

    function load(obj, varargin)
      s = dagnn.Layer.argsToStruct(varargin{:}) ;
      % backward file compatibility
      if isfield(s, 'numInputs'), s = rmfield(s, 'numInputs') ; end
      load@dagnn.Layer(obj, s) ;
    end

    function obj = Slice(varargin)
      obj.load(varargin{:}) ;
    end
  end
end

repo: MFT


text:
function handles = dtiDeleteROI(roiNum,handles)
%
%   handles = dtiDeleteROI(roiNum,handles)
%
%Author: Wandell
%Purpose:
%
% Example:
%    handles = dtiDeleteROI(handles.curRoi,handles)

if ieNotDefined('roiNum'), error('roi number required.'); end
if isempty(handles.rois), warning('No ROIS to delete'); return;  end

nRois = length(handles.rois);
if (roiNum < 1) | (roiNum > nRois), error('Bad roi number.'); end

handles.rois(roiNum) = [];
if (roiNum <= handles.curRoi), handles.curRoi = handles.curRoi - 1; end

if isempty(handles.rois), handles.curRoi = 0;
elseif handles.curRoi < 1, handles.curRoi = 1; end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
 function [coef, arg] = nufft_coef(om, J, K, kernel)
%function [coef, arg] = nufft_coef(om, J, K, kernel)
%|
%| NUFFT interpolation coefficient vectors for a given kernel function.
%|
%| in
%|	om	[M 1]	digital frequency omega in radians
%|	J		# of neighbors used per frequency location
%|	K		FFT size (should be >= N, the signal_length)
%|	kernel		kernel function handle
%|
%| out
%|	coef	[J M]	coef vector for each frequency
%|	arg	[J M]	kernel argument
%|
%| Copyright 2002-4-11, Jeff Fessler, University of Michigan

if nargin < 4, ir_usage, end

M = length(om);
gam = 2*pi/K;
dk = om / gam - nufft_offset(om, J, K);		% [M 1]
arg = outer_sum(-[1:J]', dk');			% [J M] kernel arg

coef = feval(kernel, arg, J);			% [J M]

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function result = compHolisticRanking( aid )
%COMPHOLISTICRANKING Summary of this function goes here
%   Detailed explanation goes here
global config
bufname = config.bufname;

did = config.valid_data_id(config.valid_anno_id==aid);

load(config.globalSceneSVMFile);
scenesvm = MODEL{config.globalSceneSVMiterID}.initModel;
clear MODEL

load(sprintf([bufname config.sceneFeatFile], 1));
SEED_HYPS = ALL_HYPS(1);
names = fieldnames(SEED_HYPS);
for i = 1:length(names)
    SEED_HYPS.(names{i}) = [];
end

HYPS = repmat(SEED_HYPS, 5000, 1);
ALL_COST = -10*ones(5000,1);
for rid = 1:config.sampleroomnum
    fprintf('>>%d:%d ', aid, rid);
    try
        load(sprintf([bufname '/' config.sceneFeatFile], rid));
    catch
        fprintf('no data\n');
        continue;
    end
    if isempty(ALL_HYPS)
        continue;
    end
    HYPS((rid-1)*100+1:rid*100) = ALL_HYPS;
    COST((rid-1)*100+1:rid*100) = vertcat(ALL_HYPS.COST);

    [ ~, hypScore] = sceneHypsEvaluation( vertcat(ALL_HYPS.sceneImgFea), ...
            vertcat(ALL_HYPS.MINSCORE), scenesvm, did );
    ALL_COST((rid-1)*100+1:rid*100) = hypScore;
    ROOMID((rid-1)*100+1:rid*100) = rid;
end
fprintf('\n');

load([bufname config.objHypoFile3D]);

%% get our detection performance
[~,I] = sort(ALL_COST,'descend');
%keyboard
scene = packupScene(HYPS(I(1)), CAN_POOL{ROOMID(I(1))});
DET = scene{1};

result.Hyps = HYPS(I(1));
result.GLOBAL = ALL_COST(I(1));
result.ROOMID = ROOMID(I(1));
result.DET = DET;

rot_file_name = [bufname config.IMGLIST(aid).name config.rotateImageSuffix];
rotImg = imread(rot_file_name);
figure; imshow(visualObject2D( rotImg, DET, 512, 1024 ));
figure; visualObject3D( DET, 4, 1 );

end

author: zouchuhang

repo: LayoutNet


text:
function frf=dfracft(f,a,varargin)
%DFRACFT  Discrete Fractional Fourier transform
%   Usage:  V=dfracft(f,a,p);
%           V=dfracft(f,a);
%
%   `dfracft(f,a)` computes the discrete fractional Fourier Transform of the
%   signal *f* to the power *a*. For *a=1* it corresponds to the ordinary
%   discrete Fourier Transform. If *f* is multi-dimensional, the
%   transformation is applied along the first non-singleton dimension.
%
%   `dfracft(f,a,dim)` does the same along dimension *dim*.
%
%   `dfracft(f,a,[],p)` or `dfracft(f,a,dim,p)` allows to choose the order
%   of approximation of the second difference operator (default: *p=2*).
%
%   See also:  ffracft, dft, hermbasis, pherm
%
%   References: ozzaku01,buma04

%   AUTHOR : Christoph Wiesmeyr
%   TESTING: TEST_HERMBASIS
%   REFERENCE: OK

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.p   = 2;
definput.keyvals.dim = [];
[flags,keyvals,dim,p]=ltfatarghelper({'dim','p'},definput,varargin);

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],dim,upper(mfilename));

H = hermbasis(L,p);

% set up the eigenvalues
k=0:L-1;
lam = exp(-1i*k*a*pi/2);
lam=lam(:);

% correction for even signal lengths
if ~rem(L,2)
    lam(end)=exp(-1i*L*a*pi/2);
end

% shuffle the eigenvalues in the right order
even=~mod(L,2);
cor=2*floor(L/4)+1;
for k=(cor+1):2:(L-even)
    lam([k,k+1])=lam([k+1,k]);
end

frf =H*(bsxfun(@times,lam,H'*f));

frf=assert_sigreshape_post(frf,dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [peakVals,isMax,peakIdx]=findMaxMinPoints(data)
%%FINDMAXMINPOINTS Given a grid of points in one or more dimensions, this
%                  function finds all maxima and minima. A maximum/minimum
%                  is declared if the point in question is higher/lower
%                  than all neighboring points (including those diagonally
%                  offset). Points at the edge of the matrix are not
%                  considered to be maxima or minima.
%
%INPUTS: data An n1Xn2X...Xnk dimensional matrix.  Singleton dimensions
%             (ni=1) are ignored. Non-singleton dimensions must be > 2. to
%             consider points that are not on the edge. data must be a real
%             matrix.
%
%OUTPUTS: peakVals A numValsX1 vector of the maximum and minimum values
%                  found. If no values are found, this is an empty matrix.
%            isMax A numValsX1 vector indicating whether each value in
%                  peakVals is a maximum or a minimum. 1 indicates maximum,
%                  0 indicated minimum.
%          peakIdx A numValsX1 vector of indices of the peaks such that
%                  data(peakIdx(i)) provides the ith peak value.
%
%EXAMPLE:
%This function is useful for finding maxima and minima of cost functions.
%Here, we find all maxima and minima of a surface.
% numPoints=100;
% pts=linspace(0,2*pi,numPoints);
% [X,Y]=meshgrid(pts,pts);
% Z=sin(2*X+3*Y)+cos(X-2*Y);
% %Display the surface whose peaks are desired.
% figure(1)
% clf
% surface(X,Y,Z,'EdgeColor','None')
% [peakVals,isMax,peakIdx]=findMaxMinPoints(Z)
%One sees that all 13 of the maxma and minima (-1 and 1)  that are found.
%Values on the edge are not counted for maxma or minima.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

dims=size(data);

%Get rid of singleton dimensions.
dims=dims(dims~=1);
if(isempty(dims)||isscalar(dims)&&dims==1||any(dims==2))
   peakVals=[];
   peakIdx=[];
   isMax=[];
   return;
end

numDims=length(dims);
switch(numDims)
    case 1%Find peaks for linear values.
        peakVals=[];
        peakIdx=[];
        isMax=[];
        for idx1=2:(dims(1)-1)
            curVal=data(idx1);

            vals=[curVal-data(idx1-1);
                  curVal-data(idx1+1)];
            if(all(vals>0))
                peakVals=[peakVals;curVal];
                peakIdx=[peakIdx;idx1];
                isMax=[isMax;1];
            elseif(all(vals<0))
                peakVals=[peakVals;curVal];
                peakIdx=[peakIdx;idx1];
                isMax=[isMax;0];
            end
        end
    case 2
        data=reshape(data,dims(1),dims(2));

        peakVals=[];
        peakIdx=[];
        isMax=[];

        for idx1=2:(dims(1)-1)
            for idx2=2:(dims(2)-1)
                curVal=data(idx1,idx2);

                vals=[curVal-data(idx1,idx2-1);
                      curVal-data(idx1-1,idx2-1);
                      curVal-data(idx1+1,idx2-1);
                      curVal-data(idx1,idx2+1);
                      curVal-data(idx1-1,idx2+1);
                      curVal-data(idx1+1,idx2+1);
                      curVal-data(idx1-1,idx2);
                      curVal-data(idx1+1,idx2)];

                %If we have found a maximum or a minimum
                if(all(vals>0))
                    peakVals=[peakVals;curVal];
                    idx=sub2ind(dims,idx1,idx2);
                    peakIdx=[peakIdx;idx];
                    isMax=[isMax;1];
                elseif(all(vals<0))
                    peakVals=[peakVals;curVal];
                    idx=sub2ind(dims,idx1,idx2);
                    peakIdx=[peakIdx;idx];
                    isMax=[isMax;0];
                end
            end
        end
    otherwise%3 or more dimensions
        maxTupleVals=dims-3;

        peakVals=[];
        peakIdx=[];
        isMax=[];

        idxList=getNextTuple(numDims);
        while(~isempty(idxList))
            %Points that will be compared to their neighbors are never the
            %edge points. Thus, the minimum value possible for an element
            %is 2 and the maximum value is one less than the number of
            %things in that dimension.
            shiftIdxList=idxList+2;

            curIdx=nDim2Index(dims,shiftIdxList);
            curVal=data(curIdx);

            %Now, we look at all neighbors of the current point. This means
            %that all points go through all combinations of +1,-1, and 0,
            %except for the all zero case. We can get the values using
            %tuples.
            maxSignVals=2*ones(numDims,1);

            signList=getNextTuple(zeros(numDims,1),maxSignVals);

            isMax=true;
            isMin=true;
            while(~isempty(signList)&&(isMax==true||isMin==true))
                temp=signList;
                temp(temp==2)=-1;

                shiftedIdx=nDim2Index(dims,shiftIdxList+temp);
                compVal=data(shiftedIdx);

                if(compVal>=curVal)
                   isMax=false;
                end

                if(compVal<=curVal)
                   isMin=false;
                end

                signList=getNextTuple(signList,maxSignVals);
            end

            if(~(isMax&&isMin))
                if(isMax)
                    peakVals=[peakVals;curVal];
                    peakIdx=[peakIdx;curIdx];
                    isMax=[isMax;1];
                elseif(isMin)
                    peakVals=[peakVals;curVal];
                    peakIdx=[peakIdx;curIdx];
                    isMax=[isMax;0];
                end
            end

            idxList=getNextTuple(idxList,maxTupleVals);
        end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Show basic use of spm_mlm_bayes

clear all
close all

N=100;
d=2;
p=3;

y=randn(N,d);
x=randn(N,p);
verbose=1;

mlm = spm_mlm_bayes (y,x,'input',verbose);

figure
imagesc(mlm.wmean);
colorbar
ylabel('Inputs');
xlabel('Outputs');
colormap(gray);

disp(sprintf('Model evidence = %1.2f', mlm.fm));

author: spm


text:
function [ IDX, C, d ] = kmeans2( X, k, varargin )
% Fast version of kmeans clustering.
%
% Cluster the N x p matrix X into k clusters using the kmeans algorithm. It
% returns the cluster memberships for each data point in the N x 1 vector
% IDX and the K x p matrix of cluster means in C.
%
% This function is in some ways less general than Matlab's kmeans.m (for
% example it only uses euclidian distance), but it has some options that
% the Matlab version does not (for example, it has a notion of outliers and
% min-cluster size).  It is also many times faster than matlab's kmeans.
% General kmeans help can be found in help for the matlab implementation of
% kmeans. Note that the although the names and conventions for this
% algorithm are taken from Matlab's implementation, there are slight
% alterations (for example, IDX==-1 is used to indicate outliers).
%
% IDX is a n-by-1 vector used to indicated cluster membership.  Let X be a
% set of n points.  Then the ID of X - or IDX is a column vector of length
% n, where each element is an integer indicating the cluster membership of
% the corresponding element in X.  IDX(i)=c indicates that the ith point in
% X belongs to cluster c. Cluster labels range from 1 to k, and thus
% k=max(IDX) is typically the number of clusters IDX divides X into. The
% cluster label "-1" is reserved for outliers. IDX(i)==-1 indicates that
% the given point does not belong to any of the discovered clusters. Note
% that matlab's version of kmeans does not have outliers.
%
% USAGE
%  [ IDX, C, d ] = kmeans2( X, k, [varargin] )
%
% INPUTS
%  X       - [n x p] matrix of n p-dim vectors.
%  k       - maximum nuber of clusters (actual number may be smaller)
%  prm     - additional params (struct or name/value pairs)
%   .k         - [] alternate way of specifying k (if not given above)
%   .nTrial    - [1] number random restarts
%   .maxIter   - [100] max number of iterations
%   .display   - [0] Whether or not to display algorithm status
%   .rndSeed   - [] random seed for kmeans; useful for replicability
%   .outFrac   - [0] max frac points that can be treated as outliers
%   .minCl     - [1] min cluster size (smaller clusters get eliminated)
%   .metric    - [] metric for pdist2
%   .C0        - [] initial cluster centers for first trial
%
% OUTPUTS
%  IDX    - [n x 1] cluster membership (see above)
%  C      - [k x p] matrix of centroid locations C(j,:) = mean(X(IDX==j,:))
%  d      - [1 x k] d(j) is sum of distances from X(IDX==j,:) to C(j,:)
%           sum(d) is a typical measure of the quality of a clustering
%
% EXAMPLE
%
% See also DEMOCLUSTER
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.24
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get input args
dfs = {'nTrial',1, 'maxIter',100, 'display',0, 'rndSeed',[],...
  'outFrac',0, 'minCl',1, 'metric',[], 'C0',[],'k',k };
[nTrial,maxt,dsp,rndSeed,outFrac,minCl,metric,C0,k] = ...
  getPrmDflt(varargin,dfs); assert(~isempty(k) && k>0);

% error checking
if(k<1); error('k must be greater than 1'); end
if(~ismatrix(X) || any(size(X)==0)); error('Illegal X'); end
if(outFrac<0 || outFrac>=1), error('outFrac must be in [0,1)'); end
nOut = floor( size(X,1)*outFrac );

% initialize random seed if specified
if(~isempty(rndSeed)); rand('state',rndSeed); end; %#ok<RAND>

% run kmeans2main nTrial times
bd=inf; t0=clock;
for i=1:nTrial, t1=clock; if(i>1), C0=[]; end
  if(dsp), fprintf('kmeans2 iter %i/%i step: ',i,nTrial); end
  [IDX,C,d]=kmeans2main(X,k,nOut,minCl,maxt,dsp,metric,C0);
  if(sum(d)<sum(bd)), bIDX=IDX; bC=C; bd=d; end
  if(dsp), fprintf('  d=%f  t=%fs\n',sum(d),etime(clock,t1)); end
end
IDX=bIDX; C=bC; d=bd; k=max(IDX);
if(dsp), fprintf('k=%i  d=%f  t=%fs\n',k,sum(d),etime(clock,t0)); end

% sort IDX to have biggest clusters have lower indicies
cnts = zeros(1,k); for i=1:k; cnts(i) = sum( IDX==i ); end
[~,order] = sort( -cnts ); C = C(order,:); d = d(order);
IDX2=IDX; for i=1:k; IDX2(IDX==order(i))=i; end; IDX = IDX2;

end

function [IDX,C,d] = kmeans2main( X, k, nOut, minCl, maxt, dsp, metric, C )

% initialize cluster centers to be k random X points
[N,p] = size(X); k = min(k,N); t=0;
IDX = ones(N,1); oldIDX = zeros(N,1);
if(isempty(C)), C = X(randperm(N,k),:)+randn(k,p)/1e5; end

% MAIN LOOP: loop until the cluster assigments do not change
if(dsp), nDg=ceil(log10(maxt-1)); fprintf(int2str2(0,nDg)); end
while( any(oldIDX~=IDX) && t<maxt )
  % assign each point to closest cluster center
  oldIDX=IDX; D=pdist2(X,C,metric); [mind,IDX]=min(D,[],2);

  % do not use most distant nOut elements in computation of centers
  mind1=sort(mind); thr=mind1(end-nOut); IDX(mind>thr)=-1;

  % Recalculate means based on new assignment, discard small clusters
  k0=0; C=zeros(k,p);
  for IDx=1:k
    ids=find(IDX==IDx); nCl=size(ids,1);
    if( nCl<minCl ), IDX(ids)=-1; continue; end
    k0=k0+1; IDX(ids)=k0; C(k0,:)=sum(X(ids,:),1)/nCl;
  end
  if(k0>0), k=k0; C=C(1:k,:); else k=1; C=X(randint2(1,1,[1 N]),:); end
  t=t+1; if(dsp), fprintf([repmat('\b',[1 nDg]) int2str2(t,nDg)]); end
end

% record within-cluster sums of point-to-centroid distances
d=zeros(1,k); for i=1:k, d(i)=sum(mind(IDX==i)); end

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox