text:
function [vertices, keypoints, smpl_parameter] = read_body_recon(reconstruction_dir)
    vertices = import_vertices_keypoints([reconstruction_dir, '/smpl_vertex.txt']);
    keypoints = import_vertices_keypoints([reconstruction_dir, '/keypoints.txt']);
    smpl_parameter = import_params([reconstruction_dir, '/smpl_parameter.txt']);
end

function vertices = import_vertices_keypoints(filename)
dataLines = [1, Inf];
opts = delimitedTextImportOptions("NumVariables", 3);
opts.DataLines = dataLines;
opts.Delimiter = " ";
opts.VariableNames = ["VarName1", "VarName2", "VarName3"];
opts.VariableTypes = ["double", "double", "double"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
opts.ConsecutiveDelimitersRule = "join";
opts.LeadingDelimitersRule = "ignore";
vertices = readtable(filename, opts);
vertices = table2array(vertices);
end

function params = import_params(filename, dataLines)
if nargin < 2
    dataLines = [1, Inf];
end
opts = delimitedTextImportOptions("NumVariables", 1);
opts.DataLines = dataLines;
opts.Delimiter = ",";
opts.VariableNames = "VarName1";
opts.VariableTypes = "double";
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
params = readtable(filename, opts);
params = table2array(params);
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
classdef CumMomentAnalysisWindow < AnalysisWindow
    % CUMMOMENTANALYSISWINDOW shows cumulative moment release
    properties
    end
    methods
        function obj=CumMomentAnalysisWindow(ax)
            obj@AnalysisWindow(ax);
        end

        function prepare_axes(obj)
            % prepare the moment release axes
            if isempty(obj.ax.Tag)
                obj.ax.Tag = 'dvMoment';
            end
            obj.ax.Title.String='Cum Moment Release';
            obj.ax.XLabel.String='time';
            obj.ax.YLabel.String='Cumulative Moment [N m]'; %units as per calc_moment
        end

        function [x,y]=calculate(~,catalog)
            if ~isa(catalog,'ZmapCatalog')
                catalog=catalog.Catalog;
            end
            % return the datetime and cumulative-moment for each event in the catalog
            x=catalog.Date;
            [~, y, ~] = calc_moment(catalog);
        end

    end
end

author: CelsoReyes


text:
function x0 = addfixparameters(x0,x,fixedparam)
x0(~fixedparam) = x;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% SYNTAX:
% [data_d, svs, nSV] = hmrR_MotionCorrectPCA(data_d, mlActMan, mlActAuto, tIncMan, tIncAuto, nSV)
%
% UI NAME:
% Motion_Correct_PCA
%
% This function uses PCA to filter only the segments of data identified as
% a motion artifact. The motion artifacts are indicated in the tInc vector
% by the value of 0.
%
%
% INPUTS
% data_d:  SNIRF object containing data matrix, timepoints x sd pairs
% mlActMan: Cell array of vectors, one for each time base in data, specifying
%        active/inactive channels with 1 meaning active, 0 meaning inactive
% mlActAuto: Cell array of vectors, one for each time base in data, specifying
%        active/inactive channels with 1 meaning active, 0 meaning inactive
% tIncMan: Cell array for eqach data block with vectors of length time points
%          where 1's indicating data included and 0's indicating motion artifact
% tIncAuto: Cell array for eqach data block with vectors of length time points
%          where 1's indicating data included and 0's indicating motion artifact
% nSV: Cell array for each data block with the number of principal components to remove
%      from each data block. If this number is less than 1, then the filter removes the first n
%      components of the data that removes a fraction of the variance up to nSV. Percent variance
%      you want to remove, or give an integer for number ofcomponents to remove
%
%
% OUTPUTS
% data_dN: SNIRF object containing the motion corrected data matrix.
% svs: Cell array for each data block with the singular values of the PCA
% nSV: Cell array for each data block with the number of singular values removed from the data.
%
%
% USAGE OPTIONS:
% Motion_Correct_PCA:  [dod, svs, nSV] = hmrR_MotionCorrectPCA(dod, mlActMan, mlActAuto, tIncMan, tIncAuto, nSV)
%
% PARAMETERS:
% nSV: 0.0
%
% PREREQUISITES:
% Intensity_to_Delta_OD: dod = hmrR_Intensity2OD( intensity )
%
function [data_dN, svs, nSV] = hmrR_MotionCorrectPCA(data_d,  mlActMan, mlActAuto, tIncMan, tIncAuto, nSV)

% Init output
data_dN = DataClass.empty();
svs = cell(length(data_d),1);
if ~iscell(nSV)
    nSV = repmat({nSV}, length(data_d),1);
end

% Check input args
if isempty(tIncMan)
    tIncMan = cell(length(data_d),1);
end
if isempty(tIncAuto)
    tIncAuto = cell(length(data_d),1);
end
if isempty(mlActMan)
    mlActMan = cell(length(data_d),1);
end
if isempty(mlActAuto)
    mlActAuto = cell(length(data_d),1);
end

for iBlk=1:length(data_d)
    data_dN(iBlk) = DataClass(data_d(iBlk));

    d           = data_d(iBlk).GetDataTimeSeries();
    MeasList    = data_d(iBlk).GetMeasList();
    mlActMan{iBlk} = mlAct_Initialize(mlActMan{iBlk}, MeasList);
    mlActAuto{iBlk} = mlAct_Initialize(mlActAuto{iBlk}, MeasList);
    mlAct = mlActMan{iBlk}(:,3) & mlActAuto{iBlk}(:,3);

    if isempty(tIncMan{iBlk})
        tIncMan{iBlk} = ones(size(d,1),1);
    end
    if isempty(tIncAuto{iBlk})
        tIncAuto{iBlk} = ones(size(d,1),1);
    end
    tInc = tIncMan{iBlk} & tIncAuto{iBlk};

    lstNoInc = find(tInc==0);
    lstAct = find(mlAct==1);

    if isempty(lstNoInc)
        nSV{iBlk} = 0;
        continue;
    end

    %
    % Do PCA
    %
    y = d(lstNoInc,lstAct);
    yo = y;

    c = y.' * y;
    [V,St,foo] = svd(c);
    svs{iBlk} = diag(St) / sum(diag(St));

    svsc = svs{iBlk};
    for idx = 2:size(svs{iBlk},1)
        svsc(idx) = svsc(idx-1) + svs{iBlk}(idx);
    end
    if nSV{iBlk}<1 & nSV{iBlk}>0 % find number of SV to get variance up to nSV{iBlk}
        ev = diag(svsc<nSV{iBlk});
        nSV{iBlk} = find(diag(ev)==0,1)-1;
    end

    ev = zeros(size(svs{iBlk},1),1);
    ev(1:nSV{iBlk}) = 1;
    ev = diag(ev);

    yc = yo - y*V*ev*V';

    %
    % splice the segments of data together
    %
    lstMs = find(diff(tInc)==-1);
    lstMf = find(diff(tInc)==1);
    if isempty(lstMf)
        lstMf = length(tInc);
    end
    if isempty(lstMs)
        lstMs = 1;
    end
    if lstMs(1)>lstMf(1)
        lstMs = [1;lstMs];
    end
    if lstMs(end)>lstMf(end)
        lstMf(end+1,1) = length(tInc);
    end
    lstMb = lstMf-lstMs;
    for ii=2:length(lstMb)
        lstMb(ii) = lstMb(ii-1) + lstMb(ii);
    end

    dN = d;

    for ii=1:length(lstAct)
        jj = lstAct(ii);

        lst = (lstMs(1)):(lstMf(1)-1);
        if lstMs(1)>1
            dN(lst,jj) = yc(1:lstMb(1),ii) - yc(1,ii) + dN(lst(1),jj);
        else
            dN(lst,jj) = yc(1:lstMb(1),ii) - yc(lstMb(1),ii) + dN(lst(end),jj);
        end

        for kk=1:(length(lstMf)-1)
            lst = (lstMf(kk)-1):lstMs(kk+1);
            dN(lst,jj) = d(lst,jj) - d(lst(1),jj) + dN(lst(1),jj);

            lst = (lstMs(kk+1)):(lstMf(kk+1)-1);
            dN(lst,jj) = yc((lstMb(kk)+1):lstMb(kk+1),ii) - yc(lstMb(kk)+1,ii) + dN(lst(1),jj);
        end

        if lstMf(end)<length(d)
            lst = (lstMf(end)-1):length(d);
            dN(lst,jj) = d(lst,jj) - d(lst(1),jj) + dN(lst(1),jj);
        end
    end

    data_dN(iBlk).SetDataTimeSeries(dN);
end

author: BUNPC


text:
% Group (split) position data by time or in specified number of partitions
%
% This function divides position data into groups based on pair of values (groupType, groupingValue).
% You can use it to split position data into predefined number of groups (2, 4, e.t.c.), or into
% groups of equal time.
%
%  USAGE
%   [groupIndices, edges] = general.groupPositions(pos, groupType, groupingValue)
%   pos             Position samples. Matrix of size at least Nx1.
%   groupType       String that defines type of grouping. Affects how groupingValue
%                   argument is used. 'time' defines group by time interval.
%                   'num' groups positions in specified number of groups.
%   groupingValue   Integer that specify value for grouping. If groupType == 'time'
%                   this argument should be a value of seconds.
%                   If groupType == 'num' this argument is the number of groups.
%
%   groupIndices    Vector of length equal to number of position samples. Each value
%                   indicates group belongingness. If groupIndices(3) == 1, means that
%                   position sample number 3 belongs to group number 1.
%   edges           Vector of edges that define groups. For example, if data is split
%                   in halves, then edges = [0 <middle_point> Inf].
%
%  EXAMPLE
%
%   pos = data.getPositions();
%   groupIndices = general.groupPositions(pos, 'num', 2); % divide data into halves
%   positionsOfFirstGroup = pos(groupIndices == 1, :);
%
%   general.groupPositions(pos, 'time', 1); % divide data into groups each duration of 1 second
%
function [groupIndices, edges] = groupPositions(pos, groupType, groupingValue)
    if nargin < 3
        error('BNT:numArgs', 'Incorrect number of parameters (type ''help general.groupPositions'' for details).');
    end

    if size(pos, 1) == 1 && size(pos, 2) > 2
        % this is probably a row vector, transpose it
        pos = pos';
    end

    if size(pos, 2) < 1
        error('Incorrect value for argument ''pos'' (type ''help general.groupPositions'' for details).');
    end

    if ~helpers.isstring(groupType, 'time', 'num')
        error('Incorrect value for argument ''groupType'' (type ''help general.groupPositions'' for details).');
    end

    if ~helpers.isdscalar(groupingValue, '>0')
        error('Incorrect value for argument ''groupingValue'', it should be >= 1 (type ''help general.groupPositions'' for details).');
    end

    post = pos(:, 1);
    duration = max(post) - min(post);

    if strcmpi(groupType, 'num')
        interval = duration / groupingValue;
        edges = 0:interval:duration;
        edges = min(post) + edges; % shift if for case min(post) ~= 0
        edges(end) = Inf;
    else
        interval = groupingValue;
        edges = 0:interval:duration;
        edges = min(post) + edges; % shift if for case min(post) ~= 0
        edges(end+1) = Inf;
    end

    [~, groupIndices] = histc(post, edges);
end

author: kavli-ntnu


text:
% This function dynamically generates distorted speech signals. It is
% usually used for speech enhancement or robust ASR tasks.
%
% Authors: Xiong Xiao, NTU, Singapore
% Last Modified: 28 Feb 2017
%
function [data, paraDD, data_log] = GenDynamicDistortion(base_data, paraDD)

clean_data = base_data(1).data;
rir_data = base_data(2).data;
noise_data = base_data(3).data;

% get settings
useSoftMask = ReturnFieldWithDefaultValue(paraDD, 'useSoftMask', 0);
randomizeFileOrder = ReturnFieldWithDefaultValue(paraDD, 'randomizeFileOrder', 1);
gainNorm = ReturnFieldWithDefaultValue(paraDD, 'gainNorm', 1);
vadStream = ReturnFieldWithDefaultValue(paraDD, 'vadStream', []);
if ~isempty(vadStream)
    vadData = base_data(vadStream);
end

[clean_data_loaded, clean_idx] = LoadOriginalWave(clean_data, paraDD, randomizeFileOrder,1);
if ~isempty(rir_data)
    [rir_data_loaded, rir_idx] = LoadOriginalWave(rir_data, paraDD, randomizeFileOrder,2);
end
if ~isempty(noise_data)
    [noise_data_loaded, noise_idx] = LoadOriginalWave(noise_data, paraDD, randomizeFileOrder,3);
end

if length(base_data)>3  % sometimes, we have other data, such as frame phone label. We just copy that to the output data.
    hasExtraData = 0; extraData = [];
    for i=4:length(base_data)
        if ~isempty(vadStream) && i == vadStream; continue; end
        extraData(end+1).data = base_data(i).data(clean_idx);  % note that the extra data streams are binded with the clean data.
        hasExtraData = 1;
    end
end

switch lower(paraDD.SNR_PDF)
    case 'uniform'
        SNR = rand(1,paraDD.nUtt4Iteration) * (paraDD.SNR_para(2) - paraDD.SNR_para(1)) + paraDD.SNR_para(1);
    case 'normal'
        SNR = randn(1,paraDD.nUtt4Iteration) * paraDD.SNR_para(2) + paraDD.SNR_para(1);
        SNR = max(-50, min(50, SNR));   % do not allow too low or high SNRs
end

if isfield(paraDD, 'seglen') && paraDD.seglen>0
    doSegmentation = 1;
    seglen = paraDD.seglen;
    segshift = paraDD.segshift;
    frame_len = paraDD.frame_len;
    frame_shift = paraDD.frame_shift;
else
    doSegmentation = 0;
end
if isfield(paraDD, 'outputDataIdxMask') && ~isempty(paraDD.outputDataIdxMask) && paraDD.outputDataIdxMask>0
    genMask = 1;
else
    genMask = 0;
end

if ~isfield(paraDD, 'outputDataIdxClean')   % by default, the clean speech should be in stream 2
    cleanStreamIdx = 2;
elseif ~isempty(paraDD.outputDataIdxClean) && paraDD.outputDataIdxClean>0   % otherwise, put it in the defined index
    cleanStreamIdx = paraDD.outputDataIdxClean;
else
    cleanStreamIdx = -1;
end

distorted = {}; clean_data_aligned = {}; mask = {};
for j=1:length(extraData)
    extraDataSeg(j).data = {};
end
for i=1:length(clean_idx)
    data_log{i} = [clean_data{clean_idx(i)}];
    if paraDD.singlePrecision
        curr_clean_wav = single(clean_data_loaded{clean_idx(i)})';
        if ~isempty(rir_data); curr_rir = single(rir_data_loaded{rir_idx(i)})'; else; curr_rir = []; end
        if ~isempty(noise_data); curr_noise = single(noise_data_loaded{noise_idx(i)})'; else; curr_noise = [];  end
    else
        curr_clean_wav = double(clean_data_loaded{clean_idx(i)})';
        if ~isempty(rir_data); curr_rir = double(rir_data_loaded{rir_idx(i)})'; else; curr_rir = []; end
        if ~isempty(noise_data); curr_noise = double(noise_data_loaded{noise_idx(i)})'; else; curr_noise = [];  end
    end
    if ~isempty(curr_rir)   % we want the direct sound to have the gain of 1
        curr_rir = curr_rir/max(curr_rir(:));
        data_log{i} = sprintf('%s\t%s', data_log{i}, rir_data{rir_idx(i)});
    else
        data_log{i} = sprintf('%s\tNULL', data_log{i});
    end
    if ~isempty(noise_data)
        data_log{i} = sprintf('%s\t%s', data_log{i}, noise_data{noise_idx(i)});
    else
        data_log{i} = sprintf('%s\tNULL', data_log{i});
    end
    data_log{i} = sprintf('%s\t%fdB', data_log{i}, SNR(i));

    [curr_distorted, curr_reverb, curr_direct] = ApplyConstRirNoise(curr_clean_wav, paraDD.fs, curr_rir, curr_noise, SNR(i), paraDD.useGPU);
    curr_distorted = gather(curr_distorted)';
    curr_distorted(:, length(curr_clean_wav)+1:end) = [];
    curr_direct = gather(curr_direct)';
    curr_direct(:, length(curr_clean_wav)+1:end) = [];
    %curr_reverb = gather(curr_reverb)';
    %curr_reverb(:, length(curr_clean_wav)+1:end) = [];

    if genMask
        if 0
            [curr_mask, curr_SNR] = genMaskFromParallelData(curr_clean_wav, curr_distorted', paraDD.fs, 0);
        else    % it's better to use curr_direct, which contains the early reflection up to 50ms after the direct sound, as the clean reference.
            if isempty(vadStream)
                [curr_mask] = genMaskFromParallelData(curr_clean_wav, curr_direct', curr_distorted', [], paraDD.fs, useSoftMask, 0);
            else
                [curr_mask] = genMaskFromParallelData(curr_clean_wav, curr_direct', curr_distorted', vadData.data{clean_idx(i)}, paraDD.fs, useSoftMask, 0);
            end
        end
    end
    if gainNorm
        curr_distorted = curr_distorted / max(abs(curr_distorted(:)));
    end
    switch class(gather(clean_data_loaded{clean_idx(i)}(1)))
        case 'int16'
            curr_distorted = StoreWavInt16(curr_distorted);
        case {'single', 'float'}
            curr_distorted = single(curr_distorted);
        otherwise
            curr_distorted = double(curr_distorted);
    end

    if doSegmentation
        curr_distorted_seg = DivideSent2Segments(curr_distorted, (seglen-1)*frame_shift+frame_len, segshift*frame_shift, 1);
        distorted = [distorted curr_distorted_seg'];
        curr_clean_wav_seg = DivideSent2Segments(clean_data_loaded{clean_idx(i)}, (seglen-1)*frame_shift+frame_len, segshift*frame_shift, 1);
        if cleanStreamIdx>0
            clean_data_aligned = [clean_data_aligned curr_clean_wav_seg'];
        end
        if genMask
            curr_mask_seg = DivideSent2Segments(curr_mask, seglen, segshift, 1);
            mask = [mask curr_mask_seg'];
        end
        for j=1:length(extraData)
            tmpDataSeg = DivideSent2Segments(extraData(j).data{i}, seglen, segshift, 1);    % we assume that the extra data are having 100Hz frame rate
            extraDataSeg(j).data = [extraDataSeg(j).data tmpDataSeg'];
        end
    else
        distorted{i} = curr_distorted;
        if cleanStreamIdx>0; clean_data_aligned{i} = clean_data_loaded{clean_idx(i)}; end
        if genMask; mask{i} = curr_mask; end
    end

    PrintProgress(i, length(clean_idx), max(10, round(length(clean_idx)/10)), ...
            'GenDynamicDistortion->Mixing clean speech with RIR and noise');
end

nOutputStream = 3;
if ~isfield(paraDD, 'outputDataIdxDistorted')   % by default, the distorted speech should be in stream 1
    data(1).data = distorted;
elseif ~isempty(paraDD.outputDataIdxDistorted) && paraDD.outputDataIdxDistorted>0   % otherwise, put it in the defined index
    data(paraDD.outputDataIdxDistorted).data = distorted;
else
    nOutputStream = nOutputStream -1;
    % do not output distorted speech
end
if cleanStreamIdx>0
    data(cleanStreamIdx).data = clean_data_aligned;
else
    nOutputStream = nOutputStream -1;
    % do not output clean speech
end
if genMask
    data(paraDD.outputDataIdxMask).data = mask;
else
    nOutputStream = nOutputStream -1;
end

if hasExtraData
    if doSegmentation
        data = [data extraDataSeg];
    else
        data = [data extraData];
    end
end

end

function [data, idx] = LoadOriginalWave(wavlist, paraDD, randomizeFileOrder, streamIdx)
if isempty(wavlist)
    return;
end
if randomizeFileOrder
    idx = randperm(length(wavlist));
else
    idx = 1:length(wavlist);
end
nRepeat = ceil(paraDD.nUtt4Iteration/length(idx));
idx = repmat(idx, 1, nRepeat);
idx(paraDD.nUtt4Iteration+1:end) = [];

if paraDD.inputFeature(streamIdx)
    data = wavlist;
else
    need2load = unique(idx);
    for i=1:length(need2load)
        PrintProgress(i, length(need2load), max(100, round(length(need2load)/10)), ...
            sprintf('GenDynamicDistortion->LoadOriginalWav->Stream %d', streamIdx));
        % data{need2load(i)} = InputReader(wavlist{need2load(i)}, paraDD.fileReader(streamIdx));
        data{need2load(i)} = paraDD.fileReader{streamIdx}.read(wavlist{need2load(i)});
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function varargout = waterfall(varargin)
%WATERFALL   Waterfall plot of a DISKFUN.
%   WATERFALL(F) displays the waterfall plot of F.
%
%   WATERFALL(F, S) displays the column and row chebfuns of F that are used
%   for its approximation. This is a 3D version of plot(f, S), where S is a
%   string (see PLOT).
%
%   WATERFALL(F, S, 'nslices', N) displays the first min(N, length(f))
%   column and rows.
%
%   WATERFALL supports passing options to the plot, similar to standard
%   Matlab plot commands. The options supported are:
%       'color':      Color of lines and markers plotted.
%       'marker':     Marker for pivot points.
%       'markersize': Size of markers plotted at pivot points.
%
%   H = WATERFALL(...) returns a handle to a waterfall plot object.
%
% See also DISKFUN/PLOT.

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

%switch to polar
f = varargin{1};
f = cart2pol(f, 'cdr');

[varargout{1:nargout}] = waterfall@separableApprox(f,varargin{2:end});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Matdetec2 = findMatdetecFootprint(DetectFootPrint,XPsizeOut,YPsizeOut)
% Fixed the calculation of footprint of Sentinel-2 detector based on old
% version of Sentinel-2 data (17/12/2021 Shi)

clear Matdetec k
% sorted foot print from 1 to 4
[DetectFootPrint.Ndect, idnews] = sort(DetectFootPrint.Ndect);
DetectFootPrint.Nband = DetectFootPrint.Nband(idnews);
DetectFootPrint.Data = DetectFootPrint.Data(idnews);

for i = 1:length(DetectFootPrint.Ndect)

  IDX = knnsearch(XPsizeOut',DetectFootPrint.Data{i,1}(:,1));
  IDY = knnsearch(YPsizeOut,DetectFootPrint.Data{i,1}(:,2));

  dum2 = single(poly2mask(double(IDX), double(IDY),length(XPsizeOut),length(XPsizeOut))) ;
  clear IDX IDY;
  dum2 = conv2(dum2,ones(3),'same')>0; % to fill boundary
  Matdetec(:,:,i) = dum2;
  clear dum*

  %   find orientation of detect + slope for computing perpendicular kernel
  I=nan(size(Matdetec,2),1);
  for ii=1:size(Matdetec,1)
    dum=find(Matdetec(ii,:,i)==1,1,'first');
    if ~isempty(dum)
      I(ii,1)=dum;
    end
  end
  clear dum;
  J = [1:size(Matdetec,1)]' ;
  test = ~isnan(I) & I > 1 & I < size(Matdetec,2);
  warning off all % if warning => not enough point => slope=0 => good because tile boundary
  k{i,1} = polyfit(J(test),I(test),1);
  clear test;

  I=nan(size(Matdetec,2),1);
  for ii=1:size(Matdetec,1)
    dum=find(Matdetec(ii,:,i)==1,1,'last');
    if ~isempty(dum)
      I(ii,1)=dum;
    end
  end
  J = [1:size(Matdetec,1)]' ;
  test = ~isnan(I) & I > 1 & I < size(Matdetec,2);
  k{i,2} = polyfit(J(test),I(test),1);
  clear test;
  warning on all

end

% mediane
for i = 1:length(DetectFootPrint.Ndect)-1
  mediane = mean( [k{i,2} ; k{i+1,1}] ) ;

  k{i,2} = mediane ;
  k{i+1,1} = mediane ;
  clear mediane;
end
J = [1:size(Matdetec,1)]' ;
I = [1:size(Matdetec,2)] ;

[Jmat Imat] = meshgrid(I,J);
clear I J;

Matdetec2 = nan(size(Matdetec,1),size(Matdetec,2));
clear Matdetec;
for i = 1:length(DetectFootPrint.Ndect)

  liminf = polyval(k{i,1},Jmat);
  limsup = polyval(k{i,2},Jmat);

  if sum(k{i,2}) == 0 % footprint at low-right corner
      Matdetec2(Imat>=liminf) = DetectFootPrint.Ndect(i) ;
  elseif sum(k{i,1}) == 0 % footprint at up-left corner
     Matdetec2(Imat<=limsup) = DetectFootPrint.Ndect(i) ;
  else
     Matdetec2(Imat>=liminf & Imat<=limsup) = DetectFootPrint.Ndect(i) ;
  end
  clear liminf limsup;
end
clear Imat ImatJ k;

Matdetec2 = Matdetec2';

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
% StackExchange Mathematics Q3892375
% https://math.stackexchange.com/questions/3892375
% Solve Linear Least Squares with L1 Norm Regularization with Linear
% Equality and Non Negativity Constraints.
% References:
%   1.
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     27/12/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 0;2165;42; %<! Set to 0 for Random

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Parameters

numRows         = 100;
numCols         = 20;
numRowsC        = 80;
numRowsD        = 20;
paramLambda     = 0;

% Solvers parameters
numIterations       = 5000;
stepSizeGd          = 5e-5;
stepSizeMomentum    = 0.8;
stepSizeAccel       = 0.4;

%% Load / Generate Data

mA = randn(numRows, numCols);
mC = randn(numRows, numCols);
mD = randn(numRows, numCols);
vB = randn(numRows, 1);

hObjFun = @(vX) 0.5 * sum((mA * vX - vB) .^ 2) + (paramLambda * sum(abs(mC * vX)));

mAA = mA.' * mA;
vAb = mA.' * vB;
mDD = mD * mD.';
mInvDD = pinv(mDD);

hG = @(vX) (mAA * vX) - vAb + (paramLambda * mD.' * sign(mD * vX)); %<! Gradient of the Objection Function

hP1 = @(vX) vX - (mD.' * mInvDD * (mD * vX)); %<! Projection onto the constraint set
hP2 = @(vX) max(vX, 0); %<! Projection onto the constraint set

hP = @(vX) OrthogonalProjectionOntoConvexSets({hP1; hP2}, vX, 100, 1e-6);

solverIdx       = 0;
cMethodString   = {};

mObjFunValMse   = zeros([numIterations, 1]);
mSolMse         = zeros([numIterations, 1]);

%% Solution by CVX

solverString = 'CVX';

% cvx_solver('SDPT3'); %<! Default, Keep numRows low
% cvx_solver('SeDuMi');
% cvx_solver('Mosek'); %<! Can handle numRows > 500, Very Good!
% cvx_solver('Gurobi');

hRunTime = tic();

cvx_begin('quiet')
% cvx_begin()
    % cvx_precision('best');
    variable vX(numCols, 1);
    minimize( 0.5 * sum_square(mA * vX - vB) + (paramLambda * norm(mC * vX, 1)) );
    subject to
        mD * vX == 0;
        vX >= 0;
cvx_end

runTime = toc(hRunTime);

% vX = mX(:);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The ', solverString, ' Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

sCvxSol.vXCvx     = vX;
sCvxSol.cvxOptVal = hObjFun(vX);

%% Solution by Projected Gradient Descent

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Method'];

hRunTime = tic();

[vX, mX] = ProjectedGd(zeros(numCols, 1), hG, hP, numIterations, stepSizeGd);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Solution by Projected Gradient Descent with Momentum

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Descent with Momentum'];

hRunTime = tic();

[vX, ~, mX] = ProjectedGdMomentum(zeros(numCols, 1), zeros(numCols, 1), hG, hP, numIterations, stepSizeGd, stepSizeMomentum);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Solution by Projected Gradient Descent with Nesterov Acceleration

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Descent with Nesterov Acceleration'];

hRunTime = tic();

[vX, ~, mX] = ProjectedGdAccel(zeros(numCols, 1), zeros(numCols, 1), hG, hP, numIterations, stepSizeGd, stepSizeAccel);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Solution by Projected Gradient Descent with FISTA

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Descent with FISTA Acceleration'];

hRunTime = tic();

[vX, ~, mX] = ProjectedGdFista(zeros(numCols, 1), zeros(numCols, 1), hG, hP, numIterations, stepSizeGd);
% [vX, mX] = SolveLsFista(zeros(numCols, 1), mA, vB, paramLambda, numIterations, stepSizeGd);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Display Results

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);

hAxes       = subplot(2, 1, 1);
hLineSeries = plot(1:numIterations, 10 * log10(mObjFunValMse));
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', ['Objective Function Value vs. Optimal Value (CVX)'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', '$ 10 \log_{10} {\left( \left| f \left( x \right) - f \left( {x}_{CVX} \right) \right| \right)}^{2} $', ...
    'FontSize', fontSizeAxis, 'Interpreter', 'latex');
set(hAxes, 'XLim', [1, numIterations]);
hLegend = ClickableLegend(cLegendString);

hAxes       = subplot(2, 1, 2);
hLineSeries = plot(1:numIterations, 10 * log10(mSolMse));
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', ['Solution Error Norm'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', '$ 10 \log_{10} \left( {\left\| x - {x}_{CVX} \right\|}_{2}^{2} \right) $', ...
    'FontSize', fontSizeAxis, 'Interpreter', 'latex');
set(hAxes, 'XLim', [1, numIterations]);
hLegend = ClickableLegend(cLegendString);

if(generateFigures == ON)
    % saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
% adjustlocs() - read neuroscan polar location file (.asc)
%
% Usage:
%   >> chanlocs = adjustlocs( chanlocs );
%   >> chanlocs = adjustlocs( chanlocs, 'key1', val1, 'key2', val2, ...);
%
% Inputs:
%   chanlocs       - EEGLAB channel location data structure. See
%                    help readlocs()
%
% Optional inputs:
%   'center'       - [cell array] one or several electrode names to compute
%                    center location (if several electrodes are provided, the
%                    function use their iso-barycenter). Ex: { 'cz' } or
%                    { 'fz' 'pz' }.
%   'autocenter'   - ['on'|'off']  attempt to automatically detect all symetrical
%                    electrode in the 10-20 system to compute the center.
%                    Default: 'on'.
%   'rotate'       - [cell array] name and planar angle of landmark electrodes
%                    to use as template planar rotation. Ex: { 'c3', 90 }.
%   'autorotate'   - ['on'|'off'] attempt to automatically detect
%                    electrode in the 10-20 system to compute the average
%                    planar rotation. Default 'on'.
%   'scale'        - [cell array] name and phi angle of electrodes along
%                    horizontal central line. Ex: { 'c3' 44 }. Scale uniformly
%                    all direction. Use 'hscale' and 'vscale' for scaling x and
%                    y axis independently.
%   'hscale'       - [cell array] name and phi angle of electrodes along
%                    honrizontal central line. Ex: { 'c3' 44 }.
%   'vscale'       - [cell array] name and phi angle of one electrodes along
%                    central vertical axis. Ex: { 'fz' 44 }.
%   'autoscale'    - ['on'|'off'] automatic scaling with 10-20 system used as
%                    reference. Default is 'on'.
%   'uniform'      - ['on'|'off'] force the scaling to be uniform along the X
%                    and the Y axis. Default is 'on'.
%   'coordinates'  - ['pol'|'sph'|'cart'] use polar coordinates ('pol'), sperical
%                    coordinates ('sph') or cartesian ('cart'). Default is 'sph'.
%                    (Note that using polar and spherical coordinates have the
%                    same effect, except that units are different).
%                    Default is 'sph'.
%
% Outputs:
%   chanlocs       - EEGLAB channel location data structure. See
%                    help readlocs()
%
% Note: operations are performed in the following order, first re-centering
%       then planar rotation and finally sperical re-scaling.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 1 Dec 2003
%
% See also: readlocs()

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

function chanlocs = adjustlocs( chanlocs, varargin)

    if nargin < 1
        help adjustlocs;
        return;
    end;

    % check input parameters
    % ----------------------
    g = finputcheck( varargin, { 'hscale'     'cell'  []   {};
                                 'vscale'     'cell'  []   {};
                                 'scale'      'cell'  []   {};
                                 'center'     'cell'  []   {};
                                 'rotate'     'cell'  []   {};
                                 'autoscale'  'string'  { 'on' 'off' }    'off';
                                 'autocenter' 'string'  { 'on' 'off' }    'on';
                                 'autorotate' 'string'  { 'on' 'off' }    'on';
                                 'uniform'    'string'  { 'on' 'off' }    'on';
                                 'coordinates' 'string'  { 'pol' 'sph' 'cart' } 'sph' });
    if ischar(g), error(g); end;

    names = { chanlocs.labels };

    % auto center
    % -----------
    if strcmpi(g.autocenter, 'on') & isempty(g.center)
        disp('Reading template 10-20 file');
        locs1020 = readlocs('eeglab1020.ced');

        % scan electrodes for horiz pos
        % -----------------------------
        tmpnames      = lower(names);
        tmpnames1020  = lower({ locs1020.labels });
        [tmp indelec] = intersect(tmpnames1020, tmpnames);

        % remove non-symetrical electrodes
        % --------------------------------
        if ~isempty(indelec)
            if find(indelec == 79), indelec(end+1) = 80; end; % for Cz
            ind2remove = [];
            for index = 1:length(indelec)
                if mod(indelec(index),2)
                    if ~ismember(indelec(index)+1, indelec)
                        ind2remove = [ ind2remove index ];
                    end;
                else
                    if ~ismember(indelec(index)-1, indelec)
                        ind2remove = [ ind2remove index ];
                    end;
                end;
            end;
            indelec(ind2remove) = [];
            if find(indelec == 80), indelec(end) = []; end; % for Cz

            g.center = tmpnames1020(indelec);
        end;
        if isempty(g.center)
            disp('No electrodes found for auto-centering')
        else
            disp([ num2str(length(g.center)) ' landmark electrodes found for position auto-centering' ])
        end;
    end;

    % auto rotate
    % -----------
    if strcmpi(g.autorotate, 'on') & isempty(g.rotate)
        if exist('locs1020') ~= 1
            disp('Reading template 10-20 file');
            locs1020 = readlocs('eeglab1020.ced');
        end;

        % scan electrodes for horiz pos
        % -----------------------------
        tmpnames      = lower(names);
        tmpnames1020  = lower({ locs1020.labels });
        tmptheta1020 = { locs1020.theta };
        [tmp indelec] = intersect(tmpnames1020(1:end-1), tmpnames); % do not use cz

        g.rotate(1:2:2*length(indelec))   = tmpnames1020 (indelec);
        g.rotate(2:2:2*length(indelec)+1) = tmptheta1020(indelec);

        if isempty(g.rotate)
            disp('No electrodes found for auto planar rotation')
        else
            disp([ num2str(length(g.rotate)/2) ' landmark electrodes found for auto planar rotation' ])
        end;
    end;

    % auto scale
    % ----------
    if strcmpi(g.autoscale, 'on') & isempty(g.hscale) & isempty(g.vscale)
        if exist('locs1020') ~= 1
            disp('Reading template 10-20 file');
            locs1020 = readlocs('eeglab1020.ced');
        end;

        if strcmpi(g.uniform, 'off')
            % remove all vertical electrodes for horizontal scaling
            % -----------------------------------------------------
            theta  = [ locs1020.theta ];
            indxh  = find(abs(theta) == 90);
            indxv  = union(find(theta == 0) , find(theta == 180));
            locs1020horiz = locs1020(indxh);
            locs1020vert  = locs1020(indxv);

            % scan electrodes for horiz pos
            % -----------------------------
            tmpnames      = lower(names);
            tmpnames1020  = lower({ locs1020horiz.labels });
            tmpradius1020 = { locs1020horiz.radius };
            [tmp indelec] = intersect(tmpnames1020, tmpnames);

            if isempty(indelec)
                disp('No electrodes found for horiz. position spherical re-scaling')
            else
                disp([ num2str(length(indelec)) ' landmark electrodes found for horiz. position spherical re-scaling' ]);
                if strcmpi(g.coordinates, 'cart')
                    g.hscale(2:2:2*length(indelec)+1) = [ locs1020horiz.Y ];
                else
                    g.hscale(2:2:2*length(indelec)+1) = tmpradius1020(indelec);
                end;
            end;

            % scan electrodes for vert. pos
            % -----------------------------
            tmpnames1020  = lower({ locs1020vert.labels });
            tmpradius1020 = { locs1020vert.radius };
            [tmp indelec] = intersect(tmpnames1020, tmpnames);

            if isempty(indelec)
                disp('No electrodes found for vertical position spherical re-scaling')
            else
                disp([ num2str(length(indelec)) ' landmark electrodes found for vertical spherical re-scaling' ]);
                g.vscale(1:2:2*length(indelec))   = tmpnames1020 (indelec);
                if strcmpi(g.coordinates, 'cart')
                    g.vscale(2:2:2*length(indelec)+1) = [ locs1020vert.X ];
                else
                    g.vscale(2:2:2*length(indelec)+1) = tmpradius1020(indelec);
                end;
            end;
        else
            % uniform scaling
            % ---------------
            tmpnames      = lower(names);
            tmpnames1020  = lower({ locs1020.labels });
            tmpradius1020 = { locs1020.radius };
            [tmp indelec] = intersect(tmpnames1020, tmpnames);

            if isempty(indelec)
                disp('No electrodes found for uniform spherical re-scaling')
            else
                disp([ num2str(length(indelec)) ' landmark electrodes found for uniform spherical re-scaling' ]);
                g.scale(1:2:2*length(indelec))   = tmpnames1020 (indelec);
                if strcmpi(g.coordinates, 'cart')
                    tmpabsxyz = mattocell(abs([ locs1020.X ]+j*[ locs1020.Y ]));
                    g.scale(2:2:2*length(indelec)+1) = tmpabsxyz(indelec);
                else
                    g.scale(2:2:2*length(indelec)+1) = tmpradius1020(indelec);
                end;
            end;
        end;
        if strcmpi(g.coordinates, 'sph')
            g.coordinates = 'pol'; % use polar coordinates for scaling
        end;
    end;

    % get X and Y coordinates
    % -----------------------
    if strcmpi(g.coordinates, 'sph') | strcmpi(g.coordinates, 'pol')
        [X Y] = pol2cart( [ chanlocs.theta ]/180*pi, [ chanlocs.radius ]); Z = 1;
        if strcmpi(g.coordinates, 'sph')
            X = X/0.25*46;
            Y = Y/0.25*46;
        end;
    else
        X = [ chanlocs.X ];
        Y = [ chanlocs.Y ];
        Z = [ chanlocs.Z ];
    end;

    % recenter
    % --------
    if ~isempty(g.center)
        for index = 1:length(g.center)
            tmpindex = strmatch( lower(g.center{index}), lower(names), 'exact' );
            if isempty(tmpindex)
                error(['Electrode ''' g.center{index} ''' not found for re-centering']);
            end;
            indexelec(index) = tmpindex;
        end;
        showmsg('Using electrode', 'for re-centering', g.center);
        centerx = mean(X(indexelec));
        centery = mean(Y(indexelec));
        X = X - centerx;
        Y = Y - centery;
    end;

    % planar rotation
    % ---------------
    if ~isempty(g.rotate)
        % find electrodes
        % ---------------
        clear elec;
        for index = 1:2:length(g.rotate)
            tmpindex = strmatch( lower(g.rotate{index}), lower(names), 'exact' );
            if isempty(tmpindex)
                error(['Electrode ''' g.rotate{index} ''' not found for left-right scaling']);
            end;
            elec((index+1)/2) = tmpindex;
        end;
        vals = [ g.rotate{2:2:end} ];

        % compute average scaling factor
        % ------------------------------
        [ allangles tmp ] = cart2pol(X(elec), Y(elec));
        allangles = allangles/pi*180;
        diffangle = allangles - vals;
        %diffangle2 = allangles + vals;
        %if abs(diffangle1) > abs(diffangle2), diffangle = diffangle2;
        %else                                  diffangle = diffangle1;
        %end;
        tmpind    = find(diffangle >  180); diffangle(tmpind) = diffangle(tmpind)-360;
        tmpind    = find(diffangle < -180); diffangle(tmpind) = diffangle(tmpind)+360;
        anglerot  = mean(diffangle);
        tmpcplx = (X+j*Y)*exp(-j*anglerot/180*pi);
        X = real(tmpcplx);
        Y = imag(tmpcplx);
        showmsg('Using electrode', ['for planar rotation (' num2str(anglerot,2) ' degrees)'], g.rotate(1:2:end));
    end;

    % computing scaling factors
    % -------------------------
    if ~isempty(g.scale)
        % find electrodes
        % ---------------
        clear elec;
        for index = 1:2:length(g.scale)
            tmpindex = strmatch( lower(g.scale{index}), lower(names), 'exact' );
            if isempty(tmpindex)
                error(['Electrode ''' g.scale{index} ''' not found for left-right scaling']);
            end;
            elec((index+1)/2) = tmpindex;
        end;
        vals = [ g.scale{2:2:end} ];

        % compute average scaling factor
        % ------------------------------
        nonzero = find(vals > 0);
        hscalefact = mean(abs(Y(elec(nonzero))+j*X(elec(nonzero)))./vals(nonzero)); % *46/0.25; %/44/0.25;
        vscalefact = hscalefact;
        showmsg('Using electrode', ['for uniform spherical re-scaling (x' num2str(1/hscalefact,4) ')'], g.scale(1:2:end));
    else
        if ~isempty(g.hscale)
            % find electrodes
            % ---------------
            clear elec;
            for index = 1:2:length(g.hscale)
                tmpindex = strmatch( lower(g.hscale{index}), lower(names), 'exact' );
                if isempty(tmpindex)
                    error(['Electrode ''' g.hscale{index} ''' not found for left-right scaling']);
                end;
                elec((index+1)/2) = tmpindex;
            end;
            vals = [ g.hscale{2:2:end} ];
            showmsg('Using electrode', [ 'for left-right spherical re-scaling (x' num2str(1/hscalefact,4) ')'], g.hscale(1:2:end));

            % compute average scaling factor
            % ------------------------------
            hscalefact = mean(abs(Y(elec))./vals); % *46/0.25; %/44/0.25;
            if isempty(g.vscale)
                vscalefact =  hscalefact;
            end;
        end;
        if ~isempty(g.vscale)
            % find electrodes
            % ---------------
            clear elec;
            for index = 1:2:length(g.vscale)
                tmpindex = strmatch( lower(g.vscale{index}), lower(names), 'exact' );
                if isempty(tmpindex)
                    error(['Electrode ''' g.vscale{index} ''' not found for rear-front scaling']);
                end;
                elec((index+1)/2) = tmpindex;
            end;
            vals = [ g.vscale{2:2:end} ];
            showmsg('Using electrode', ['for rear-front spherical re-scaling (x' num2str(1/vscalefact,4) ')'], g.vscale(1:2:end));

            % compute average scaling factor
            % ------------------------------
            vscalefact = mean(abs(X(elec))./vals); % *46/0.25; %/44/0.25;
            if isempty(g.vscale)
                hscalefact =  vscalefact;
            end;
        end;
    end;

    % uniform?
    % --------
    if strcmpi(g.uniform, 'on') & ( ~isempty(g.vscale) | ~isempty(g.hscale))
        disp('uniform scaling: averaging left-right and rear-front scaling factor');
        hscalefact = mean([hscalefact vscalefact]);
        vscalefact = hscalefact;
    end;

    % scaling data
    % ------------
    if ~isempty(g.vscale) | ~isempty(g.hscale) | ~isempty(g.scale)
        Y = Y/hscalefact;
        X = X/vscalefact;
        Z = Z/((hscalefact+vscalefact)/2);
    end;

    % updating structure
    % ------------------
    if strcmpi(g.coordinates, 'sph') |  strcmpi(g.coordinates, 'pol')
        [phi,theta] = cart2pol(Y, X);
        phi = phi/pi*180;
        if strcmpi(g.coordinates, 'pol')
            theta = theta/0.25*46;
        end;

        % convert to other types of coordinates
        % -------------------------------------
        labels = names';
        chanlocs = struct('labels', names, 'sph_theta_besa', mattocell(theta), ...
                'sph_phi_besa', mattocell(phi), 'sph_radius', { chanlocs.sph_radius });
        chanlocs = convertlocs( chanlocs, 'sphbesa2all');
    else
        for index = 1:length(chanlocs)
            chanlocs(index).X = X(index);
            chanlocs(index).Y = Y(index);
            chanlocs(index).Z = Z(index);
        end;
        chanlocs = convertlocs(chanlocs, 'cart2all');
    end;

function showmsg(begmsg, endmsg, struct);
    if length(struct) <= 1
        disp([ begmsg ' ''' struct{1} ''' ' endmsg]);
    elseif length(struct) <= 2
        disp([ begmsg ' ''' struct{1} ''' and ''' struct{2}  ''' ' endmsg]);
    elseif length(struct) <= 3
        disp([ begmsg ' ''' struct{1} ''', ''' struct{2} ''' and ''' struct{3}  ''' ' endmsg]);
    else
        disp([ begmsg ' ''' struct{1} ''', ''' struct{2} ''', ''' struct{3}  ''' ... ' endmsg]);
    end;

author: PatternRecognition

repo: OpenBMI


text:
% ======================================================== %
% Files of the Matlab programs included in the book:       %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,  %
% Second Edition, Luniver Press, (2010).   www.luniver.com %
% ======================================================== %

% -------------------------------------------------------- %
% Firefly Algorithm for constrained optimization using     %
% for the design of a spring (benchmark)                   %
% by Xin-She Yang (Cambridge University) Copyright @2009   %
% -------------------------------------------------------- %

function fa_ndim
% parameters [n N_iteration alpha betamin gamma]
para=[20 500 0.5 0.2 1];

help fa_ndim.m

% Simple bounds/limits for d-dimensional problems
d=15;
Lb=zeros(1,d);
Ub=2*ones(1,d);

% Initial random guess
u0=Lb+(Ub-Lb).*rand(1,d);

[u,fval,NumEval]=ffa_mincon(@cost,u0,Lb,Ub,para);

% Display results
bestsolution=u
bestojb=fval
total_number_of_function_evaluations=NumEval

%%% Put your own cost/objective function here --------%%%
%% Cost or Objective function
 function z=cost(x)
% Exact solutions should be (1,1,...,1)
z=sum((x-1).^2);

%%% End of the part to be modified -------------------%%%

%%% --------------------------------------------------%%%
%%% Do not modify the following codes unless you want %%%
%%% to improve its performance etc                    %%%
% -------------------------------------------------------
% ===Start of the Firefly Algorithm Implementation ======
%         Lb = lower bounds/limits
%         Ub = upper bounds/limits
%   para == optional (to control the Firefly algorithm)
% Outputs: nbest   = the best solution found so far
%          fbest   = the best objective value
%      NumEval = number of evaluations: n*MaxGeneration
% Optional:
% The alpha can be reduced (as to reduce the randomness)
% ---------------------------------------------------------

% Start FA
function [nbest,fbest,NumEval]...
           =ffa_mincon(fhandle,u0, Lb, Ub, para)
% Check input parameters (otherwise set as default values)
if nargin<5, para=[20 500 0.25 0.20 1]; end
if nargin<4, Ub=[]; end
if nargin<3, Lb=[]; end
if nargin<2,
disp('Usuage: FA_mincon(@cost,u0,Lb,Ub,para)');
end

% n=number of fireflies
% MaxGeneration=number of pseudo time steps
% ------------------------------------------------
% alpha=0.25;      % Randomness 0--1 (highly random)
% betamn=0.20;     % minimum value of beta
% gamma=1;         % Absorption coefficient
% ------------------------------------------------
n=para(1);  MaxGeneration=para(2);
alpha=para(3); betamin=para(4); gamma=para(5);

% Total number of function evaluations
NumEval=n*MaxGeneration;

% Check if the upper bound & lower bound are the same size
if length(Lb) ~=length(Ub),
    disp('Simple bounds/limits are improper!');
    return
end

% Calcualte dimension
d=length(u0);

% Initial values of an array
zn=ones(n,1)*10^100;
% ------------------------------------------------
% generating the initial locations of n fireflies
[ns,Lightn]=init_ffa(n,d,Lb,Ub,u0);

% Iterations or pseudo time marching
for k=1:MaxGeneration,     %%%%% start iterations

% This line of reducing alpha is optional
 alpha=alpha_new(alpha,MaxGeneration);

% Evaluate new solutions (for all n fireflies)
for i=1:n,
   zn(i)=fhandle(ns(i,:));
   Lightn(i)=zn(i);
end

% Ranking fireflies by their light intensity/objectives
[Lightn,Index]=sort(zn);
ns_tmp=ns;
for i=1:n,
 ns(i,:)=ns_tmp(Index(i),:);
end

%% Find the current best
nso=ns; Lighto=Lightn;
nbest=ns(1,:); Lightbest=Lightn(1);

% For output only
fbest=Lightbest;

% Move all fireflies to the better locations
[ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,nbest,...
      Lightbest,alpha,betamin,gamma,Lb,Ub);

end   %%%%% end of iterations

% -------------------------------------------------------
% ----- All the subfunctions are listed here ------------
% The initial locations of n fireflies
function [ns,Lightn]=init_ffa(n,d,Lb,Ub,u0)
  % if there are bounds/limits,
if length(Lb)>0,
   for i=1:n,
   ns(i,:)=Lb+(Ub-Lb).*rand(1,d);
   end
else
   % generate solutions around the random guess
   for i=1:n,
   ns(i,:)=u0+randn(1,d);
   end
end

% initial value before function evaluations
Lightn=ones(n,1)*10^100;

% Move all fireflies toward brighter ones
function [ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,...
             nbest,Lightbest,alpha,betamin,gamma,Lb,Ub)
% Scaling of the system
scale=abs(Ub-Lb);

% Updating fireflies
for i=1:n,
% The attractiveness parameter beta=exp(-gamma*r)
   for j=1:n,
      r=sqrt(sum((ns(i,:)-ns(j,:)).^2));
      % Update moves
if Lightn(i)>Lighto(j), % Brighter and more attractive
   beta0=1; beta=(beta0-betamin)*exp(-gamma*r.^2)+betamin;
   tmpf=alpha.*(rand(1,d)-0.5).*scale;
   ns(i,:)=ns(i,:).*(1-beta)+nso(j,:).*beta+tmpf;
      end
   end % end for j

end % end for i

% Check if the updated solutions/locations are within limits
[ns]=findlimits(n,ns,Lb,Ub);

% This function is optional, as it is not in the original FA
% The idea to reduce randomness is to increase the convergence,
% however, if you reduce randomness too quickly, then premature
% convergence can occur. So use with care.
function alpha=alpha_new(alpha,NGen)
% alpha_n=alpha_0(1-delta)^NGen=10^(-4);
% alpha_0=0.9
delta=1-(10^(-4)/0.9)^(1/NGen);
alpha=(1-delta)*alpha;

% Make sure the fireflies are within the bounds/limits
function [ns]=findlimits(n,ns,Lb,Ub)
for i=1:n,
     % Apply the lower bound
  ns_tmp=ns(i,:);
  I=ns_tmp<Lb;
  ns_tmp(I)=Lb(I);

  % Apply the upper bounds
  J=ns_tmp>Ub;
  ns_tmp(J)=Ub(J);
  % Update this new move
  ns(i,:)=ns_tmp;
end

%% ==== End of Firefly Algorithm implementation ======

author: Sable

repo: mcbench-benchmarks


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the column row of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

epsilon_init = (6/(L_in))^0.5;

W = rand(L_out,1+L_in)*2*epsilon_init-epsilon_init;

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

end

repo: MachineLearningEx


text:
%NULIBSVC Trainable classifier: LIBSVM, nu-algorithme
%
% 	[W,J,NU] = NULIBSVC(A,KERNEL,NU)
% 	[W,J,NU] = A*NULIBSVC([],KERNEL,NU)
% 	[W,J,NU] = A*NULIBSVC(KERNEL,NU)
%
% INPUT
%   A	      Dataset
%   KERNEL  Mapping to compute kernel by A*MAP(A,KERNEL)
%           or string to compute kernel by FEVAL(KERNEL,A,A)
%           or cell array with strings and parameters to compute kernel by
%           FEVAL(KERNEL{1},A,A,KERNEL{2:END})
%           Default: linear kernel.
%   NU      nu value, upperbound error.
%           Default NU is derived from 1-NN error.
%
% OUTPUT
%   W       Mapping: Support Vector Classifier
%   J       Object idences of support objects. Can be also obtained as W{4}
%   NU      Actual nu_value used
%
% DESCRIPTION
% Optimizes a support vector classifier for the dataset A by the libsvm
% package, see http://www.csie.ntu.edu.tw/~cjlin/libsvm/. LIBSVC calls the
% svmtrain routine of libsvm for training. Classifier execution for a
% test dataset B may be done by D = B*W; In D posterior probabilities are
% given as computed by svmpredict using the '-b 1' option.
%
% The kernel may be supplied in KERNEL by
% - an untrained mapping, e.g. a call to PROXM like W = LIBSVC(A,PROXM([],'R',1))
% - a string with the name of the routine to compute the kernel from A
% - a cell-array with this name and additional parameters.
% This will be used for the evaluation of a dataset B by B*W or PRMAP(B,W) as
% well.
%
% If KERNEL = 0 (or not given) it is assumed that A is already the
% kernelmatrix (square). In this also a kernel matrix should be supplied at
% evaluation by B*W or PRMAP(B,W). However, the kernel has to be computed with
% respect to support objects listed in J (the order of objects in J does matter).
%
% REFERENCES
% R.-E. Fan, P.-H. Chen, and C.-J. Lin. Working set selection using the second order
% information for training SVM. Journal of Machine Learning Research 6, 1889-1918, 2005
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, LIBSVC, SVC, PROXM

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

function [W,J,NU] = nulibsvc(varargin)

	checktoolbox('libsvm');

  mapname = 'nuLIBSVM';
  argin = shiftargin(varargin,{'prmapping','char','cell'});
  argin = setdefaults(argin,[],proxm([],'p',1),[],[]);

  if mapping_task(argin,'definition')

    W = define_mapping(argin,'untrained',mapname);

	elseif mapping_task(argin,'training')			% Train a mapping.

    [a,kernel,NU] = check_for_old_call(argin);
		if isempty(NU), NU = 2*min(max(testk(a,1),0.01),(0.8*min(classsizes(a))/size(a,1))); end
		%opt = ['-s 0 -t 4 -b 1 -c ',num2str(NU)];
    opt = ['-s 1 -t 4 -b 1 -n ',num2str(NU), ' -q'];
		islabtype(a,'crisp');
		isvaldset(a,1,2); % at least 1 object per class, 2 classes
		[m,k,c] = getsize(a);
		nlab = getnlab(a);

		K = compute_kernel(a,a,kernel);
		K = min(K,K');   % make sure kernel is symmetric
		K = [[1:m]' K];  % as libsvm wants it
	                   % call libsvm
		u = svmtrain(nlab,K,opt);
		if isempty(u)
			prwarning(1,'nulibsvc: no solution for SVM, pseudo-inverse will be used')
			W = lkc(prdataset(K(:,2:end),getlabels(a)),0);
			J = [1:m]';
			return
		end
		                 % Store the results:
   	J = full(u.SVs);
    if isequal(kernel,0)
      s = [];
      in_size = length(J);
      % in_size = 0; % to allow old and new style calls
    else
      s = a(J,:);
      in_size = k;
    end

    lablist = getlablist(a);
    W = prmapping(mfilename,'trained',{u,s,kernel,J,opt},lablist(u.Label,:),in_size,c);

		W = setname(W,mapname);
		W = setcost(W,a);

  else % Evaluation

    [a,W] = deal(argin{1:2});
    [u,s,kernel,J,opt] = getdata(W);
		m = size(a,1);

		K = compute_kernel(a,s,kernel);
    k = size(K,2);
    if k ~= length(J)
      if isequal(kernel,0)
        if (k > length(J)) &  (k >= max(J))
          % precomputed kernel; old style call
          prwarning(1,'Old style execution call: The precomputed kernel was calculated on a test set and the whole training set!')
        else
          error('Inappropriate precomputed kernel!\nFor the execution the kernel matrix should be computed on a test set and the set of support objects');
        end
      else
        error('Kernel matrix has the wrong number of columns');
      end
    else
      % kernel was computed with respect to the support objects
      % we make an approprite correction in the libsvm structure
      u.SVs = sparse((1:length(J))');
    end
    K = [[1:m]' K];  % as libsvm wants it
    %[lab,acc,d] = svmpredict(getnlab(a),K,u,'-b 1');
    [lab,acc,d] = svmpredict(ones(m,1),K,u,'-b 1');
		W = setdat(a,d,W);

	end

 return;

function K = compute_kernel(a,s,kernel)

	% compute a kernel matrix for the objects a w.r.t. the support objects s
	% given a kernel description

	if  isstr(kernel) % routine supplied to compute kernel
		K = feval(kernel,a,s);
	elseif iscell(kernel)
		K = feval(kernel{1},a,s,kernel{2:end});
	elseif ismapping(kernel)
		K = a*prmap(s,kernel);
	elseif kernel == 0 % we have already a kernel
		K = a;
	else
		error('Do not know how to compute kernel matrix')
	end

	K = +K;

return

function [a,kernel,NU] = check_for_old_call(argin)

[a,kernel,NU,par] = deal(argin{:});
if ischar(kernel) && exist(kernel,'file') ~= 2
  kernel = proxm(kernel,NU);
  NU = par;
end

author: marianux

repo: ecg-kit

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


text:
function h=joint_histogram(x,y)
%
% Takes a pair of images of equal size and returns the 2d joint histogram.
% used for MI calculation
%
% written by Amir Pasha Mahmoudzadeh
% Wright State University
% Biomedical Imaging Lab

rows=size(x,1);
cols=size(y,2);
N=256;

h=zeros(N,N);

for i=1:rows;
  for j=1:cols;
    h(x(i,j)+1,y(i,j)+1)= h(x(i,j)+1,y(i,j)+1)+1;
  end
end

imshow(h)

end


text:
function varargout = contourf( v, data, varargin )
% spherical filled contour plot
%
% Syntax
%   contourf(v,data)
%
% Input
%  v - @vector3d
%  data - double
%
% Options
%  contours - number of contours
%
% See also
% vector3d/plot vector3d/contour

if nargin == 1, data = []; end
if ischar(data)
  varargin = [data,varargin];
  data = [];
end

% in older matlab version we have to plot contour and countour lines
% seperately to avoid artefacts
if verLessThan('matlab','8.5')

  [varargout{1:nargout}] = v.smooth(data,'contours',10,'LineStyle','none',varargin{:});

  v.smooth(data,'contours',10,'fill','off','LineStyle','-','LineColor','k','hold',varargin{:});
else
  [varargout{1:nargout}] = v.smooth(data,'contours',10,'LineStyle','-','LineColor','k',varargin{:});
end

if nargout == 0, clear h; end

% TODO: data may not set

author: mtex-toolbox

repo: mtex

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


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverse dynamics for the 2dof planar robot
%
%   tau = inversedynamics_2dofplanar_fast(robot, q, qd, qdd, fe)
%
%   Where robot stores the kinematic and dynamic parameters for this robot.
%   q: joint positions.
%   qd: joint velocities.
%   qdd: joint accelerations.
%   fe: vector of external forces. (unused)
%
%   This function just executes the inverse dynamic model for this robot.
%   The equations to compute this dynamic model can be found in:
%   "ROBOT ANALYSIS. The mechanics of Serial and Parallel
%        manipulators". Lung Weng Tsai. John Wiley and Sons, inc. ISBN:
%        0-471-32593-7. page 405.
%
%   Author: Arturo Gil Aparicio arturo.gil@umh.es
%   Date: 08/03/2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

fprintf('\nComputing inverse dynamics for the %s robot', robot.name);

a = eval(robot.DH.a);
a1=a(1);
a2=a(2);

gc=9.81;

m1=robot.dynamics.masses(1);
m2=robot.dynamics.masses(2);

%express it as a general dynamic function
%M*qdd + V + G = Q, where Q is a vector of generalized forces or moments
% in this case, it is done just to represent the equations more clearly.
% have a look at the directdynamics_2dofplanar_fast function to observe
% this expression used to compute the direct dynamic model

%M is a 2x2 manipulator inertia matrix
M = [(m1*(1/3) + m2)*a1^2 + m2*a1*a2*cos(q(2)) + (1/3)*m2*a2^2 ...
      (1/2)*m2*a1*a2*cos(q(2))+(1/3)*m2*a2^2 ;
      (1/2)*m2*a1*a2*cos(q(2))+(1/3)*m2*a2^2  (1/3)*m2*a2^2];

V=[-m2*a1*a2*sin(q(2))*(qd(1)*qd(2)+(1/2)*qd(2)^2);
    (1/2)*m2*a1*a2*sin(q(2))*qd(1)^2];

G=[gc*(((1/2)*m1+m2)*a1*cos(q(1)) + (1/2)*m2*a2*cos(q(1)+q(2)));
    (1/2)*m2*gc*a2*cos(q(1)+q(2))];

tau = M*qdd + V + G;

repo: ARTE


text:
function [fnames, defaults] = mysubs_fields

% function [fnames, defaults] = mysubs_fields
% Additional fields for class cfg_choice. See help of
% @cfg_item/subs_fields for general help about this function.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: mysubs_fields.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

fnames = {'values'};
defaults = {{}};

author: spm


text:
function i = mytsearch ( x, y, t, xi, yi, i )

%*****************************************************************************80
%
%% MYTSEARCH: Find the enclosing triangle for points in a 2D plane.
%
%  Discussion:
%
%    The indices of the triangles enclosing the points in [XI,YI] are
%    returned.  The triangulation T of [X,Y] must be convex.  Points lying
%    outside the triangulation are assigned a NaN index.
%
%    I is an optional initial guess for the indices. A full search is
%    done using the standard TSEARCH function for points with an invalid
%    initial guess.
%
%  Author:
%
%    Darren Engwirda
%.

% I/O and error checks
if nargin<6
   i = [];
   if nargin<5
      error('Wrong number of inputs');
   end
elseif nargin>6
   error('Wrong number of inputs');
end
if nargout>1
   error('Wrong number of outputs');
end
ni = size(xi,1);
if (numel(xi)~=ni) || (numel(yi)~=ni) || ...
      (numel(x)~=numel(y)) || (~isempty(i) && (numel(i)~=ni))
   error('Wrong input dimensions');
end

% Translate to the origin and scale the min xy range onto [-1,1]
% This is absolutely critical to avoid precision issues for large problems!
maxxy = max([x,y]);
minxy = min([x,y]);
den = 0.5*min(maxxy-minxy);

x  = ( x-0.5*(minxy(1)+maxxy(1))) / den;
y  = ( y-0.5*(minxy(2)+maxxy(2))) / den;
xi = (xi-0.5*(minxy(1)+maxxy(1))) / den;
yi = (yi-0.5*(minxy(2)+maxxy(2))) / den;

% Check initial guess
if ( ~isempty(i) )
   k = find(i>0 & ~isnan(i));

   tri = i(k);

   n1 = t(tri,1);
   n2 = t(tri,2);
   n3 = t(tri,3);

   ok = sameside(x(n1),y(n1),x(n2),y(n2),xi(k),yi(k),x(n3),y(n3)) & ...
        sameside(x(n2),y(n2),x(n3),y(n3),xi(k),yi(k),x(n1),y(n1)) & ...
        sameside(x(n3),y(n3),x(n1),y(n1),xi(k),yi(k),x(n2),y(n2));

   j = true(ni,1);
   j(k(ok)) = false;
else
   j = true(ni,1);
end
%
% Do a full search for points that failed
%
if any(j)
   i(j) = tsearch(x,y,t,xi(j),yi(j));
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function i = sameside(xa,ya,xb,yb,x1,y1,x2,y2)

% Test if [x1(i),y1(i)] and [x2(i),y2(i)] lie on the same side of the line
% AB(i).

  dx = xb-xa;
  dy = yb-ya;
  a1 = (x1-xa).*dy-(y1-ya).*dx;
  a2 = (x2-xa).*dy-(y2-ya).*dx;

% If sign(a1)=sign(a2) the points lie on the same side
  i = false(length(xa),1);
  i(a1.*a2>=0.0) = true;

end

author: johannesgerer

repo: jburkardt-m


text:
function res = fft2c(x)

S = size(x);
fctr = S(1)*S(2);

x = reshape(x,S(1),S(2),prod(S(3:end)));

res = zeros(size(x));
for n=1:size(x,3)
	res(:,:,n) = 1/sqrt(fctr)*fftshift(fft2(ifftshift(x(:,:,n))));
end

res = reshape(res,S);

author: thomaskuestner


text:
function [ H, f, c ] = trifbank( M, K, R, fs, h2w, w2h )
% TRIFBANK Triangular filterbank.
%
%   [H,F,C]=TRIFBANK(M,K,R,FS,H2W,W2H) returns matrix of M triangular filters
%   (one per row), each K coefficients long along with a K coefficient long
%   frequency vector F and M+2 coefficient long cutoff frequency vector C.
%   The triangular filters are between limits given in R (Hz) and are
%   uniformly spaced on a warped scale defined by forward (H2W) and backward
%   (W2H) warping functions.
%
%   Inputs
%           M is the number of filters, i.e., number of rows of H
%
%           K is the length of frequency response of each filter
%             i.e., number of columns of H
%
%           R is a two element vector that specifies frequency limits (Hz),
%             i.e., R = [ low_frequency high_frequency ];
%
%           FS is the sampling frequency (Hz)
%
%           H2W is a Hertz scale to warped scale function handle
%
%           W2H is a wared scale to Hertz scale function handle
%
%   Outputs
%           H is a M by K triangular filterbank matrix (one filter per row)
%
%           F is a frequency vector (Hz) of 1xK dimension
%
%           C is a vector of filter cutoff frequencies (Hz),
%             note that C(2:end) also represents filter center frequencies,
%             and the dimension of C is 1x(M+2)
%
%   Example
%           fs = 16000;               % sampling frequency (Hz)
%           nfft = 2^12;              % fft size (number of frequency bins)
%           K = nfft/2+1;             % length of each filter
%           M = 23;                   % number of filters
%
%           hz2mel = @(hz)(1127*log(1+hz/700)); % Hertz to mel warping function
%           mel2hz = @(mel)(700*exp(mel/1127)-700); % mel to Hertz warping function
%
%           % Design mel filterbank of M filters each K coefficients long,
%           % filters are uniformly spaced on the mel scale between 0 and Fs/2 Hz
%           [ H1, freq ] = trifbank( M, K, [0 fs/2], fs, hz2mel, mel2hz );
%
%           % Design mel filterbank of M filters each K coefficients long,
%           % filters are uniformly spaced on the mel scale between 300 and 3750 Hz
%           [ H2, freq ] = trifbank( M, K, [300 3750], fs, hz2mel, mel2hz );
%
%           % Design mel filterbank of 18 filters each K coefficients long,
%           % filters are uniformly spaced on the Hertz scale between 4 and 6 kHz
%           [ H3, freq ] = trifbank( 18, K, [4 6]*1E3, fs, @(h)(h), @(h)(h) );
%
%            hfig = figure('Position', [25 100 800 600], 'PaperPositionMode', ...
%                              'auto', 'Visible', 'on', 'color', 'w'); hold on;
%           subplot( 3,1,1 );
%           plot( freq, H1 );
%           xlabel( 'Frequency (Hz)' ); ylabel( 'Weight' ); set( gca, 'box', 'off' );
%
%           subplot( 3,1,2 );
%           plot( freq, H2 );
%           xlabel( 'Frequency (Hz)' ); ylabel( 'Weight' ); set( gca, 'box', 'off' );
%
%           subplot( 3,1,3 );
%           plot( freq, H3 );
%           xlabel( 'Frequency (Hz)' ); ylabel( 'Weight' ); set( gca, 'box', 'off' );
%
%   Reference
%           [1] Huang, X., Acero, A., Hon, H., 2001. Spoken Language Processing:
%               A guide to theory, algorithm, and system development.
%               Prentice Hall, Upper Saddle River, NJ, USA (pp. 314-315).

%   Author  Kamil Wojcicki, UTD, June 2011

    if( nargin~= 6 ), help trifbank; return; end; % very lite input validation

    f_min = 0;          % filter coefficients start at this frequency (Hz)
    f_low = R(1);       % lower cutoff frequency (Hz) for the filterbank
    f_high = R(2);      % upper cutoff frequency (Hz) for the filterbank
    f_max = 0.5*fs;     % filter coefficients end at this frequency (Hz)
    f = linspace( f_min, f_max, K ); % frequency range (Hz), size 1xK
    fw = h2w( f );

    % filter cutoff frequencies (Hz) for all filters, size 1x(M+2)
    c = w2h( h2w(f_low)+[0:M+1]*((h2w(f_high)-h2w(f_low))/(M+1)) );
    cw = h2w( c );

    H = zeros( M, K );                  % zero otherwise
    for m = 1:M

        % implements Eq. (6.140) on page 314 of [1]
        % k = f>=c(m)&f<=c(m+1); % up-slope
        % H(m,k) = 2*(f(k)-c(m)) / ((c(m+2)-c(m))*(c(m+1)-c(m)));
        % k = f>=c(m+1)&f<=c(m+2); % down-slope
        % H(m,k) = 2*(c(m+2)-f(k)) / ((c(m+2)-c(m))*(c(m+2)-c(m+1)));

        % implements Eq. (6.141) on page 315 of [1]
        k = f>=c(m)&f<=c(m+1); % up-slope
        H(m,k) = (f(k)-c(m))/(c(m+1)-c(m));
        k = f>=c(m+1)&f<=c(m+2); % down-slope
        H(m,k) = (c(m+2)-f(k))/(c(m+2)-c(m+1));

   end

   % H = H./repmat(max(H,[],2),1,K);  % normalize to unit height (inherently done)
   % H = H./repmat(trapz(f,H,2),1,K); % normalize to unit area

% EOF

author: a-nagrani

repo: VGGVox

save_path: github-repos/MATLAB/a-nagrani-VGGVox


text:
clear

%%

load('results/results_clnf_cross-data_general.mat');
load('results/menpo_labels.mat');
[clnf_error, frontal_ids] = compute_error_menpo_1( labels,  experiment.shapes);
clnf_error_frontal = clnf_error(frontal_ids);
clnf_error_profile = clnf_error(~frontal_ids);

load('results/results_ceclm_cross-data.mat');

[ceclm_error, frontal_ids] = compute_error_menpo_1( labels,  experiment.shapes);
ceclm_error_frontal = ceclm_error(frontal_ids);
ceclm_error_profile = ceclm_error(~frontal_ids);

load('results/CFAN_menpo_train.mat');
for i = 1:numel(shapes)
    shapes{i} = shapes{i}-0.5;
end

[cfan_error, frontal_ids] = compute_error_menpo_1(labels, shapes);
cfan_error_frontal = cfan_error(frontal_ids);
cfan_error_profile = cfan_error(~frontal_ids);

load('results/tcdcn_menpo.mat');
for i = 1:numel(shapes)
    shapes{i} = shapes{i};
end

[tcdcn_error, frontal_ids] = compute_error_menpo_1(labels, shapes);
tcdcn_error_frontal = tcdcn_error(frontal_ids);
tcdcn_error_profile = tcdcn_error(~frontal_ids);

load('results/menpo_train_3DDFA.mat');
for i = 1:numel(shapes)
    shapes{i} = shapes{i}-0.5;
end

[error_3ddfa, frontal_ids] = compute_error_menpo_1(labels, shapes);
error_3ddfa_frontal = error_3ddfa(frontal_ids);
error_3ddfa_profile = error_3ddfa(~frontal_ids);

load('results/Menpo-CFSS_train.mat');
shapes = cell(size(estimatedPoseFull,1),1);

for i = 1:numel(shapes)
    shape = cat(2, estimatedPoseFull(i,1:68)', estimatedPoseFull(i,69:end)');
    shapes{i} = shape-0.5;
end

[cfss_error, frontal_ids] = compute_error_menpo_1(labels, shapes);
cfss_error_frontal = cfss_error(frontal_ids);
cfss_error_profile = cfss_error(~frontal_ids);

author: TadasBaltrusaitis

repo: OpenFace


text:
function alpha = anglePoints3d(varargin)
%ANGLEPOINTS3D Compute angle between three 3D points.
%
%   ALPHA = anglePoints3d(P1, P2)
%   Computes angle (P1, O, P2), in radians, between 0 and PI.
%
%   ALPHA = anglePoints3d(P1, P2, P3)
%   Computes angle (P1, P2, P3), in radians, between 0 and PI.
%
%   ALPHA = anglePoints3d(PTS)
%   PTS is a 3x3 or 2x3 array containing coordinate of points.
%
%   Example
%       rad2deg(anglePoints3d([0 0 1],[1 1 0]))
%
%   See also
%   points3d, angles3d

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-21
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

p2 = [0 0 0];
if length(varargin) == 1
    pts = varargin{1};
    if size(pts, 1)==2
        p1 = pts(1,:);
        p0 = [0 0 0];
        p2 = pts(2,:);
    else
        p1 = pts(1,:);
        p0 = pts(2,:);
        p2 = pts(3,:);
    end

elseif length(varargin) == 2
    p1 = varargin{1};
    p0 = [0 0 0];
    p2 = varargin{2};

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

% ensure all data have same size
n1 = size(p1, 1);
n2 = size(p2, 1);
n0 = size(p0, 1);

if n1 ~= n0
    if n1 == 1
        p1 = repmat(p1, [n0 1]);
        n1 = n0;
    elseif n0==1
        p0 = repmat(p0, [n1 1]);
    else
        error('Arguments P1 and P0 must have the same size');
    end
end

if n1 ~= n2
    if n1 == 1
        p1 = repmat(p1, [n2 1]);
    elseif n2 == 1
        p2 = repmat(p2, [n1 1]);
    else
        error('Arguments P1 and P2 must have the same size');
    end
end

% normalized vectors
p1 = normalizeVector3d(p1 - p0);
p2 = normalizeVector3d(p2 - p0);

% compute angle
alpha = acos(dot(p1, p2, 2));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This is a testing environment for the files in the folder kernel/viewer
% 1. Based on data/contents, a list of required files is generated and it
%    is verified, that all files are present; additional files are listed.
% 2. All c-files are compiled.
% 3. All files are executed.
%==============================================================================

FAIRcheckFiles(mfilename);

viewImage('reset','viewImage','viewImage2D','colormap','gray(256)');
viewImage('disp');
[scheme,parameter] = viewImage
viewImage('clear');
viewImage('disp');
testEnd

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


text:
function v = head(s)
% HEAD  As streamCar
v = streamCar(s);


text:
clear all, close all, clc

xC = [2; 1;];                           % Center of data (mean)
sig = [2; .5;];                         % Principal axes

theta = pi/3;                           % Rotate cloud by pi/3
R = [cos(theta) -sin(theta);            % Rotation matrix
    sin(theta) cos(theta)];

nPoints = 10000;                        % Create 10,000 points
X = R*diag(sig)*randn(2,nPoints) + diag(xC)*ones(2,nPoints);

subplot(1,2,1)                          % Plot cloud of noisy data
scatter(X(1,:),X(2,:),'k.','LineWidth',2)
hold on, box on, grid on
axis([-6 8 -6 8])

%% f_ch01_ex03_1b

Xavg = mean(X,2);                       % Compute mean
B = X - Xavg*ones(1,nPoints);           % Mean-subtracted Data
[U,S,V] = svd(B/sqrt(nPoints),'econ');  % Find principal components (SVD)

subplot(1,2,2)
scatter(X(1,:),X(2,:),'k.','LineWidth',2)  % Plot data to overlay PCA
hold on, box on, grid on
axis([-6 8 -6 8])

theta = (0:.01:1)*2*pi;
[Xstd] = U*S*[cos(theta); sin(theta)];  % 1-std confidence interval
plot(Xavg(1)+Xstd(1,:),Xavg(2) + Xstd(2,:),'r-','LineWidth',1.5)
plot(Xavg(1)+2*Xstd(1,:),Xavg(2) + 2*Xstd(2,:),'r-','LineWidth',1.5)
plot(Xavg(1)+3*Xstd(1,:),Xavg(2) + 3*Xstd(2,:),'r-','LineWidth',1.5)

% Plot principal components U(:,1)S(1,1) and U(:,2)S(2,2)
plot([Xavg(1) Xavg(1)+U(1,1)*S(1,1)],[Xavg(2) Xavg(2)+U(2,1)*S(1,1)],'c-','LineWidth',2)
plot([Xavg(1) Xavg(1)+U(1,2)*S(2,2)],[Xavg(2) Xavg(2)+U(2,2)*S(2,2)],'c-','LineWidth',2)

set(gcf,'Position',[100 100 600 300])

author: dynamicslab


text:
function output = fetch_output(handle)
    if isa(handle, 'caffe.Net')
        output = fetch_net_output(handle);
    elseif isa(handle, 'caffe.Solver')
        cnt = zeros(20,1);
        for i = 1 : length(handle.nets)
            output_t = fetch_net_output(handle.nets{i});
            if i == 1
                output = output_t;
                for j = 1 : length(output_t)
                    if output_t(j).value > 0
                        cnt(j)=cnt(j)+1;
                    end
                end
            else
                for j = 1 : length(output_t)
                    output(j).value = output(j).value + output_t(j).value;
                    if output_t(j).value > 0
                        cnt(j)=cnt(j)+1;
                    end
                end
            end
        end
        for j = 1 : length(output)
            output(j).value = output(j).value / cnt(j);
        end
    else
        error('only accept net or solver handle');
    end
end

function output = fetch_net_output(handle)
    outputs = handle.outputs;
    for i = 1 : length(outputs)
        output(i).name = outputs{i};
        output(i).value = handle.blob_vec(handle.name2blob_index(output(i).name)).get_data();
    end
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function [pkgNameStr] = pkgName(varargin)
	% Get name for the package (CIAtah by default).
	% Biafra Ahanonu
	% started: updated: 2021.03.16 [14:39:41]
	% inputs
		%
	% outputs
		%

	% changelog
		%
	% TODO
		%

	%========================
	% DESCRIPTION
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		pkgNameStr = 'CIAtah';
	catch err
		pkgNameStr = 'CIAtah';
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%% Optimise parameters of GMM through stochastic gradient descent
clear all;
%% Create ground truth GMM

Priors = ones(1,3)./3;
Mu     = [-8 2;0 0; 3 -3]';
Sigma(:,:,1) = [2 0.5; 0.5 1];
Sigma(:,:,2) = [1 0;0  5];
Sigma(:,:,3) = [4 0;0 1];

X = gmm_sample(500,Priors,Mu,Sigma);

%% Plot groud truth GMM & Data

close all;box on; grid on;
plot_gmm_contour(gca,Priors,Mu,Sigma,[0 0 1],1);
plot(X(1,:),X(2,:),'r+');

%% Initial parameters of GMM to train & Data

Priors_I = ones(1,3)./3;
Mu_I     = [ [-8;-2],[8;0],[0;8] ];
for i=1:3
    Sigma_I(:,:,i) = eye(2,2);
end

%% Plot Initial parameters & Data

close all;
figure; hold on;grid on; box on;
plot_gmm_contour(gca,Priors_I,Mu_I,Sigma_I,[0 0 1],1);
plot(X(1,:),X(2,:),'r+');
axis equal;

%% Gradient descent on Likelihood

W           = log_gmm_div_w(X,Priors_I,Mu_I,Sigma_I);
[dMu,dMu_n] = log_gmm_div_mu(X,W,Mu_I,Sigma_I);
dPriors     = log_gmm_div_prior(W,Priors);
dSigma      = log_gmm_div_sigma(X,W,Mu_I + dMu,Sigma_I);

dTSigma = dSigma;

%% Plot gradient direction of Mu

colors = [1 0 0;
          0 1 0;
          0 0 1];
close all;
figure; hold on;grid on; box on;
plot_gmm_contour(gca,Priors_I,Mu_I,Sigma_I,colors,3);

plot_gmm_contour(gca,Priors_I,Mu_I,dTSigma,[0 0 0],3);

%      (3 x 3)
color = (colors * W')';
scatter(X(1,:),X(2,:),10,color,'filled');

for k=1:3
    %(D x N)
    x = repmat(Mu_I(:,k),1,size(dMu_n,2));
    quiver(x(1,:),x(2,:),dMu_n(1,:,k),dMu_n(2,:,k));
    quiver(x(1,1),x(2,1),dMu(1,k),dMu(2,k),'r');

end

axis equal;

%% Optimise over mean

Priors_I = ones(1,3)./3;
Mu_I     = [ [-8;-2],[8;0],[0;8] ];
for i=1:3
    Sigma_I(:,:,i) = eye(2,2);
end

Priors_k = Priors_I;
Mu_k     = Mu_I;
Sigma_k  = Sigma_I;

W        = log_gmm_div_w(X,Priors_I,Mu_I,Sigma_I);
color = (colors * W')';
scatter(X(1,:),X(2,:),10,color,'filled');

K       = 1;
logliks = [];
bVision = true;

if bVision
    close all;
    hf = figure; hold on;grid on; box on;
    ghandle = plot_gmm_contour(gca,Priors_k,Mu_k,Sigma_k,colors,3);
    hs      = scatter(X(1,:),X(2,:),10,color,'filled');
    axis equal;
    axis([get(gca,'XLim'),get(gca,'YLim')]);
end

%% Optimisation

K = 1000;

for k=1:K

    % compute log likelihood
    logliks(k) =  LogLikelihood_gmm(X,Priors_k,Mu_k,Sigma_k);

    %(D x N x K)
    W       = log_gmm_div_w(X,Priors_k,Mu_k,Sigma_k);
    dPriors = log_gmm_div_prior(W,Priors);
    dMu     = log_gmm_div_mu(X,W,Mu_k,Sigma_k);
%    dSigma  = log_gmm_div_sigma(X,W, Mu_k + dMu,Sigma_k);

    Priors_k = Priors_k + dPriors;
    Priors_k = Priors_k./sum(Priors_k);
    Mu_k     = Mu_k + dMu;
%    Sigma_k  = Sigma_k + 0.001 .* dSigma;

    disp(['it(' num2str(k) ') loglik: ' num2str(logliks(k))]);

    if bVision
        plot_gmm_contour(gca,Priors_k,Mu_k,Sigma_k,[0 0 1],1,ghandle);
        color = (colors * W')';
        set(hs,'CData',color);
        pause(0.01);
    end

    if k > 1
       if abs(logliks(k) - logliks(k-1)) < 1e-8
            k = K +1;
            disp(['Optimisation finished diff: ' num2str(abs(logliks(k) - logliks(k-1)))]);
       end
    end

end

%% Plot Likelihood

figure; hold on;
plot(-logliks);
title('Log lik vs epoch');
xlabel('epoch');

%% Plot result Priors_k, Mu_k, Sigma_k for k=END

close all;
figure; hold on;grid on; box on;
plot_gmm_contour(gca,Priors_k,Mu_k,Sigma_k,[0 0 1],1);
plot(X(1,:),X(2,:),'r+');

for k=1:3
        %(D x N)
    x = repmat(Mu_k(:,k),1,size(dMu,2));
end

axis equal;

%% EM solution

EM(X, Priors_I, Mu_I, Sigma_I,'full');

%%

Priors_I = ones(1,3)./3;
Mu_I     = [ [-8;-2],[8;0],[0;8] ];
for i=1:3
    Sigma_I(:,:,i) = eye(2,2);
end

K       = 50;
options = statset('Display','final','MaxIter',K);
S       = struct('mu',Mu_I','Sigma',Sigma_I,'ComponentProportion',Priors_I);
GMModel = fitgmdist(X',3,'Options',options,'CovarianceType','full','Start',S);

Priors_k        = GMModel.ComponentProportion;
Mu_k            = GMModel.mu';
Sigma_k         = GMModel.Sigma;

author: epfl-lasa


text:
function out = distfcm(center, data)
%DISTFCM Distance measure in fuzzy c-mean clustering.
%	OUT = DISTFCM(CENTER, DATA) calculates the Euclidean distance
%	between each row in CENTER and each row in DATA, and returns a
%	distance matrix OUT of size M by N, where M and N are row
%	dimensions of CENTER and DATA, respectively, and OUT(I, J) is
%	the distance between CENTER(I,:) and DATA(J,:).
%
%       See also FCMDEMO, INITFCM, IRISFCM, STEPFCM, and FCM.

%	Roger Jang, 11-22-94, 6-27-95.
%       Copyright 1994-2016 The MathWorks, Inc.

out = zeros(size(center, 1), size(data, 1));

% fill the output matrix

if size(center, 2) > 1
    for k = 1:size(center, 1)
	out(k, :) = sqrt(sum(((data-ones(size(data, 1), 1)*center(k, :)).^2), 2));
    end
else	% 1-D data
    for k = 1:size(center, 1)
	out(k, :) = abs(center(k)-data)';
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [lmval,indd]=lmax(xx,filt)
%LMAX 	[lmval, indd]=lmax(xx,filt). Find local maxima in vector XX,where
%	LMVAL is the output vector with maxima values, INDD  is the
%	corresponding indexes, FILT is the number of passes of the small
%	running average filter in order to get rid of small peaks.  Default
%	value FILT =0 (no filtering). FILT in the range from 1 to 3 is
%	usially sufficient to remove most of a small peaks
%	For example:
%	xx=0:0.01:35; y=sin(xx) + cos(xx ./3);
%	plot(xx,y); grid; hold on;
%	[b,a]=lmax(y,2)
%	 plot(xx(a),y(a),'r+')
%	see also LMIN, MAX, MIN

%**************************************************|
% 	Serge Koptenko, Guigne International Ltd., |
%	phone (709)895-3819, fax (709)895-3822     |
%--------------06/03/97----------------------------|

x=xx;
len_x = length(x);
	fltr=[1 1 1]/3;
  if nargin <2, filt=0;
	else
x1=x(1); x2=x(len_x);
	for jj=1:filt,
	c=conv(fltr,x);
	x=c(2:len_x+1);
	x(1)=x1;
        x(len_x)=x2;
	end
  end
lmval=[]; indd=[];
i=2;		% start at second data point in time series
    while i < len_x-1,
	if x(i) > x(i-1)
	   if x(i) > x(i+1)	% definite max
lmval =[lmval x(i)];
indd = [ indd i];
	   elseif x(i)==x(i+1)&x(i)==x(i+2)	% 'long' flat spot
%lmval =[lmval x(i)];  	%1   comment these two lines for strict case
%indd = [ indd i];	%2 when only  definite max included
i = i + 2;  		% skip 2 points
	   elseif x(i)==x(i+1)	% 'short' flat spot
%lmval =[lmval x(i)];	%1   comment these two lines for strict case
%indd = [ indd i];	%2 when only  definite max included
i = i + 1;		% skip one point
	   end
	end
	i = i + 1;
    end
if filt>0 & ~isempty(indd),
	if (indd(1)<= 3)|(indd(length(indd))+2>length(xx)),
	   rng=1;	%check if index too close to the edge
	else rng=2;
	end
	  for ii=1:length(indd), 	% Find the real maximum value
	    [val(ii) iind(ii)] = max(xx(indd(ii) -rng:indd(ii) +rng));
	    iind(ii)=indd(ii) + iind(ii)  -rng-1;
	  end
  indd=iind; lmval=val;
else
end


text:
function feature = import_idt(file,tra_len)
if nargin < 2
    tra_len = 15;
end
    fid = fopen(file,'rb');
    feat = fread(fid,[10+4*tra_len+96*3+108,inf],'float');
	feature = struct('info',[],'tra',[],'tra_shape',[],'hog',[],'hof',[],'mbhx',[],'mbhy',[]);
	if ~isempty(feat)
		feature.info = feat(1:10,:);
		feature.tra = feat(11:10+tra_len*2,:);
        feature.tra_shape = feat(11+tra_len*2:10+tra_len*4,:);
        ind = 10+tra_len*4;
		feature.hog = feat(ind+1:ind+96,:);
		feature.hof = feat(ind+97:ind+204,:);
		feature.mbhx = feat(ind+205:ind+300,:);
		feature.mbhy = feat(ind+301:end,:);
	end
    fclose(fid);
end

author: wanglimin

repo: TDD

save_path: github-repos/MATLAB/wanglimin-TDD


text:
classdef Create < handle
    properties
        serPort
        delay
        cur_mode           % current mode: passive, safe, full
    end

    methods
        function robot = Create(varargin)
        %[serPort] = RoombaInit(my_COM)
        % initializes serial port for use with Roomba
        % COMMport is the number of the comm port
        % ex. RoombaInit(1) sets port = 'COM1'
        % note that it sets baudrate to a default of 57600
        % can be changed (see SCI).
        % An optional time delay can be added after all commands
        % if your code crashes frequently.  15 ms is recommended by irobot
        % By; Joel Esposito, US Naval Academy, 2011

        % This code puts the robot in CONTROL(132) mode, which means does NOT stop
        % when cliff sensors or wheel drops are true; can also run while plugged into charger

            Contrl = 132;

            % Esposito 9/2008

            warning off

            opt.port = '';
            opt.delay = 0.015;

            opt = tb_optparse(opt, varargin);

            robot.delay = opt.delay;

            % clean up any previous instances of the port, can happen...
            for tty = instrfind('port', opt.port)
                if ~isempty(tty)
                    disp(['serPort ' tty.serPort.port 'is in use.   Closing it.'])
                    fclose(tty);
                    delete(tty);
                end
            end

            if opt.verbose
                disp('Establishing connection to Roomba...');
            end

            % defaults at 57600, can change
            robot.serPort = serial(opt.port,'BaudRate', 57600);
            set(robot.serPort,'Terminator','LF')
            set(robot.serPort,'InputBufferSize',100)
            set(robot.serPort, 'Timeout', 1)
            set(robot.serPort, 'ByteOrder','bigEndian');
            set(robot.serPort, 'Tag', 'Roomba')

            robot.serPort

            if opt.verbose
                disp('Opening connection to Roomba...');
            end

            pause(0.5);
            try
            fopen(robot.serPort);
            catch me
                disp('open failed');
                me.message
                return
            end

            %% Confirm two way connumication
            disp('Setting Roomba to Control Mode...');

            % Start! and see if its alive
            robot.mode('passive');

            robot.mode('full');

            % light LEDS
            robot.write([139 25 0 128]);

            % set song 1, default beep
            robot.write([140 1 1 48 20], 0.05);

            % sing it
            robot.write([141 1])

            disp('I am alive if my two outboard lights came on')

            confirmation = robot.read(4, 0.1);
        end

        function msg = read(robot, nbytes, delay)
            msg = fread(robot.serPort, nbytes);
            if nargin > 2
                pause(delay);
            end
        end

        function v = fread(robot, n, type, delay)
            v = fread(robot.serPort, n, type);

            if nargin > 3
                pause(delay);
            else
                pause(robot.delay)
            end
        end

        function write(robot, msg, delay)
            fwrite(robot.serPort, msg);
            fprintf('write: '); disp(msg);

            if nargin > 2
                pause(delay);
            else
                pause(robot.delay)
            end
        end

        function fwrite(robot, val, type, delay)
            fwrite(robot.serPort, val, type);
            disp(['fwrite:' val '(' type ')']);

            if nargin > 3
                pause(delay);
            else
                pause(robot.delay)
            end
        end

        function flush(robot)
            %Flush Buffer
            N = robot.serPort.BytesAvailable();
            while(N ~= 0)
                fread(robot.serPort,N);
                N = robot.serPort.BytesAvailable();
            end
        end

        function delete(robot)
            tty = instrfind('port', robot.serPort.port);
            fclose(tty);
            pause(1)
            delete(tty);
            pause(1)
        end

        function s = char(robot)
            s = sprintf('iCreate robot on port %s', robot.serPort.port);
        end
    end % methods
end % classdef

author: petercorke

repo: robotics-toolbox-matlab


text:
function r = lt(p,q)
% MEAS/LT  Implement p < q for meas.

% find the difference between the two
diff = meas();
diff = p - q;

% see if p < q at the 1-sigma level
if (diff.value < diff.error)
    r = 1;
else
    r = 0;
end


text:
function data = detrenddata(data,T)

for n=1:length(T)
    ind = sum(T(1:n-1))+ (1:T(n));
    if isstruct(data)
        data.X(ind,:) = detrend(data.X(ind,:));
    else
        data(ind,:) = detrend(data(ind,:));
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function val = checkFuncs(obj, x, derivative_level)
    % Evaluates the external functions (or function handles) with a given
    % input argument.
    %
    % Parameters:
    % x: the value of dependent variables @type colvec
    % derivative_level: determines the level of derivatives to be exported
    %  (1, 2) @type double
    %
    % Return values:
    % val: a structure of computed values @type struct

    % test function
    if nargin < 2
        x = vertcat([obj.DepVariables.InitialValue]);
    end

    if nargin < 3
        derivative_level = 1;
    end

    if isempty(obj.AuxData)
        f_value = feval(obj.Funcs.Func, x);
        j_value = feval(obj.Funcs.Jac, x);
        j_struct = feval(obj.Funcs.JacStruct, 0);
        if derivative_level == 2
            h_value = feval(obj.Funcs.Hess, x, ones(obj.Dimension,1));
            h_struct = feval(obj.Funcs.Hess, 0, ones(obj.Dimension,1));
        else
            h_value = [];
            h_struct = [];
        end
    else
        f_value = feval(obj.Funcs.Func, x, obj.AuxData);
        j_value = feval(obj.Funcs.Jac, x, obj.AuxData);
        j_struct = feval(obj.Funcs.JacStruct, 0);
        if derivative_level == 2
            h_value = feval(obj.Funcs.Hess, x, ones(obj.Dimension,1), obj.AuxData);
            h_struct = feval(obj.Funcs.Hess, 0, ones(obj.Dimension,1));
        else
            h_value = [];
            h_struct = [];
            h_mat = [];
        end
    end

    % validate computed values
    assert(length(f_value) == obj.Dimension,...
        'The dimension of the output value does not match the dimension specified.');

    assert(~isempty(j_value),...
        'All-zero Jacobian entries found.');
    if derivative_level == 2
        if obj.Type == NlpFunction.LINEAR
            assert(isempty(h_value),...
                'Non-zero Hessian entries found for the Linear constraints.');
        elseif obj.Type == NlpFunction.NONLINEAR
            assert(~isempty(h_value),...
                'All-zero Hessian entries found for the Non-Linear constraints.');
        end
        h_mat = sparse(h_struct(:,1), h_struct(:,2), h_value, length(x), length(x));
    end

    j_mat = sparse(j_struct(:,1), j_struct(:,2), j_value, length(f_value), length(x));

    val = struct(...
        'f', f_value, ...
        'J', j_mat, ...
        'H', h_mat);

end

author: ayonga

repo: frost-dev

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


text:
function b = r83p_ml ( n, a_lu, x, job )

%*****************************************************************************80
%
%% R83P_ML computes A * x or x * A, where A has been factored by R83P_FA.
%
%  Discussion:
%
%    The R83P storage format stores a periodic tridiagonal matrix as
%    a 3 by N array, in which each row corresponds to a diagonal, and
%    column locations are preserved.  The matrix value
%    A(1,N) is stored as the array entry A(3,N), and the matrix value
%    A(N,1) is stored as the array entry A(1,1).
%
%  Example:
%
%    Here is how a R83P matrix of order 5 would be stored:
%
%      A51 A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54 A15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 3.
%
%    Input, real A_LU(3,N), the factors computed by R83P_FA.
%
%    Input, real X(N), the vector to be multiplied by the matrix.
%
%    Input, integer JOB, indicates what product should be computed.
%    0, compute A * x.
%    nonzero, compute A' * x.
%
%    Output, real B(N), the result of the multiplication.
%

%
%  Multiply A(1:N-1,1:N-1) and X(1:N-1).
%
  b(1:n-1) = r83_np_ml ( n-1, a_lu(1:3,1:n-1), x, job );
%
%  Add terms from the border.
%
  if ( job == 0 )
    b(1) = b(1) + a_lu(3,n) * x(n);
    b(n-1) = b(n-1) + a_lu(1,n) * x(n);
    b(n) = a_lu(1,1) * x(1) + a_lu(3,n-1) * x(n-1) + a_lu(2,n) * x(n);
  else
    b(1) = b(1) + a_lu(1,1) * x(n);
    b(n-1) = b(n-1) + a_lu(3,n-1) * x(n);
    b(n) = a_lu(3,n) * x(1) + a_lu(1,n) * x(n-1) + a_lu(2,n) * x(n);
  end

  return
end


text:
function [params, names] = simwhiteKernExtractParam(kern)

% SIMWHITEKERNEXTRACTPARAM Extract parameters from the SIM-WHITE kernel
% structure.
% FORMAT
% DESC extracts parameters from the SIM-White (Single Input Motif - White)
% kernel structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the SIM-White (Single
% Input Motif - White) kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If the
% field 'transforms' is not empty in the kernel structure, the parameters
% will be transformed before optimisation (for example positive only
% parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO simwhiteKernParamInit, simwhiteKernExpandParam, kernExtractParam,
% scg, conjgrad
%
% COPYRIGHT : David Luengo, 2009
%
% KERN

params = [kern.decay kern.variance kern.sensitivity];
if nargout > 1
  names = {'decay', 'variance', 'sensitivity'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function mrcnn_build()
if ~exist('liblinear_train')
  fprintf('Compiling liblinear version 1.94\n');
  fprintf('Source code page:\n');
  fprintf('   http://www.csie.ntu.edu.tw/~cjlin/liblinear/\n');
  mex -outdir bin ...
      CFLAGS="\$CFLAGS -std=c99 -O3 -fPIC" -largeArrayDims ...
      external/liblinear-1.94/matlab/train.c ...
      external/liblinear-1.94/matlab/linear_model_matlab.c ...
      external/liblinear-1.94/linear.cpp ...
      external/liblinear-1.94/tron.cpp ...
      "external/liblinear-1.94/blas/*.c" ...
      -output liblinear_train;
end

if ~exist('adaptive_region_pooling_mex')
  fprintf('Compiling adaptive_region_pooling_mex\n');

  mex -outdir bin ...
      -largeArrayDims ...
      code/adaptive_region_pooling/adaptive_region_pooling_mex.cpp ...
      -output adaptive_region_pooling_mex;
end

if ~exist('nms_mex')
  fprintf('Compiling nms_mex\n');

  mex -outdir bin ...
      -largeArrayDims ...
      code/postprocessing/nms_mex.cpp ...
      -output nms_mex;
end
end

author: gidariss

repo: mrcnn-object-detection