text:
function [ objects ] = buildObjStruct( list_path, list_img, list_event, list_event2 )
%BUILDOBJSTRUCT Builds the objects structure for storing all the
%   information.
%%%%

    %% Build struct for storing all objects found and their classes
    objects = struct('folder', [], 'imgName', [], 'idEvent', [], 'labelEvent', [], 'objects', struct('ULx', [], 'ULy', [], ...
        'BRx', [], 'BRy', [], 'objScore', [], 'eventAwareScore', [], 'features', [], 'label', 0));
    % label = 0 ---> not analyzed
    % label = 1 ---> no object

    %% Extract object candidates and objectness from each of them
    lenImgs = length(list_img);
    disp(['Starting extraction of objects from ' num2str(lenImgs) ' images.']);
    for i = 1:lenImgs
        % Initialize info for this image
        path = regexp(list_img{i}, '/', 'split');
        objects(i).folder = list_path{i};
        objects(i).imgName = path{length(path)};
        objects(i).idEvent = list_event(i);
        objects(i).labelEvent = list_event2(i);
    end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function N = find_steps_number(T,dt)
in = 10:16;
Ni = 2.^in;
[~,idx_N] = min(abs(bsxfun(@minus,Ni,round(T/dt))));
N = 2^in(idx_N);
end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function x = utl_purify_expression(x)
% Recursively remove impure elements from an expression.
% Output-Expression = utl_purify_expression(Input-Expression)
%
% In particular, x.tracking.expression is replaced by x. It is assumed that the head is not an impure expression.
%
%                                                     Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                                     2010-05-23

if isscalar(x) && all(isfield(x,{'head','parts'}))
    % x is a canonical expression: recurse
    for i=1:length(x.parts)
        x.parts{i} = utl_purify_expression(x.parts{i}); end
elseif isfield(x,'tracking') && isfield(x.tracking,'expression')
    % x is an impure expression: purify & recurse
    x = utl_purify_expression(x.tracking.expression);
end

repo: REST


text:
% -------------------------------------------------------------------------
function live_paths = incremental_linking(frames,iouth,costtype,jumpgap,threhgap)
% -------------------------------------------------------------------------

num_frames = length(frames);

%% online path building

live_paths = struct(); %% Stores live paths
dead_paths = struct(); %% Store the paths that has been terminated
dp_count = 0;
for  t = 1:num_frames
    num_box = size(frames(t).boxes,1);
    if t==1
        for b = 1 : num_box
            live_paths(b).boxes = frames(t).boxes(b,:);
            live_paths(b).scores = frames(t).scores(b);
            live_paths(b).allScores(t,:) = frames(t).allScores(b,:);
            live_paths(b).pathScore = frames(t).scores(b);
            live_paths(b).foundAT(t) = 1;
            live_paths(b).count = 1;
            live_paths(b).lastfound = 0; %less than 5 mean yes
        end
    else
        lp_count = getPathCount(live_paths);

        %         fprintf(' %d ', t);
        edge_scores = zeros(lp_count,num_box);

        for lp = 1 : lp_count
            edge_scores(lp,:) = score_of_edge(live_paths(lp),frames(t),iouth,costtype);
        end

        dead_count = 0 ;
        coverd_boxes = zeros(1,num_box);
        path_order_score = zeros(1,lp_count);
        for lp = 1 : lp_count
            if live_paths(lp).lastfound < jumpgap %less than 5 mean yes
                box_to_lp_score = edge_scores(lp,:);
                if sum(box_to_lp_score)>0 %%checking if atleast there is one match
                    [m_score,maxInd] = max(box_to_lp_score);
                    live_paths(lp).count = live_paths(lp).count + 1;
                    lpc = live_paths(lp).count;
                    live_paths(lp).boxes(lpc,:) = frames(t).boxes(maxInd,:);
                    live_paths(lp).scores(lpc) = frames(t).scores(maxInd);
                    live_paths(lp).allScores(lpc,:) = frames(t).allScores(maxInd,:);
                    live_paths(lp).pathScore = live_paths(lp).pathScore + m_score;
                    live_paths(lp).foundAT(lpc) = t;
                    live_paths(lp).lastfound = 0;
                    edge_scores(:,maxInd) = 0;
                    coverd_boxes(maxInd) = 1;
                else
                    live_paths(lp).lastfound = live_paths(lp).lastfound +1;
                end

                scores = sort(live_paths(lp).scores,'ascend');
                num_sc = length(scores);
                path_order_score(lp) = mean(scores(max(1,num_sc-jumpgap):num_sc));

            else
                dead_count = dead_count + 1;
            end
        end

        % Sort the path based on scoe of the boxes and terminate dead path

        [live_paths,dead_paths,dp_count] = sort_live_paths(live_paths,....
            path_order_score,dead_paths,dp_count,jumpgap);
        lp_count = getPathCount(live_paths);
        % start new paths using boxes that are not assigned
        if sum(coverd_boxes)<num_box
            for b = 1 : num_box
                if ~coverd_boxes(b)
                    lp_count = lp_count + 1;
                    live_paths(lp_count).boxes = frames(t).boxes(b,:);
                    live_paths(lp_count).scores = frames(t).scores(b);
                    live_paths(lp_count).allScores = frames(t).allScores(b,:);
                    live_paths(lp_count).pathScore = frames(t).scores(b);
                    live_paths(lp_count).foundAT = t;
                    live_paths(lp_count).count = 1;
                    live_paths(lp_count).lastfound = 0;
                end
            end
        end
    end
end

live_paths = fill_gaps(live_paths,threhgap);
dead_paths = fill_gaps(dead_paths,threhgap);
lp_count = getPathCount(live_paths);
lp = lp_count+1;
if isfield(dead_paths,'boxes')
    for dp = 1 : length(dead_paths)
        live_paths(lp).start = dead_paths(dp).start;
        live_paths(lp).end = dead_paths(dp).end;
        live_paths(lp).boxes = dead_paths(dp).boxes;
        live_paths(lp).scores = dead_paths(dp).scores;
        live_paths(lp).allScores = dead_paths(dp).allScores;
        live_paths(lp).pathScore = dead_paths(dp).pathScore;
        live_paths(lp).foundAT = dead_paths(dp).foundAT;
        live_paths(lp).count = dead_paths(dp).count;
        live_paths(lp).lastfound = dead_paths(dp).lastfound;
        lp = lp + 1;
    end
end

live_paths = sort_paths(live_paths);

% -------------------------------------------------------------------------
function sorted_live_paths = sort_paths(live_paths)
% -------------------------------------------------------------------------
sorted_live_paths = struct();

lp_count = getPathCount(live_paths);
if lp_count>0
    path_order_score = zeros(1,lp_count);

    for lp = 1 : length(live_paths)
        scores = sort(live_paths(lp).scores,'descend');
        num_sc = length(scores);
        path_order_score(lp) = mean(scores(1:min(20,num_sc)));
    end

    [~,ind] = sort(path_order_score,'descend');
    for lpc = 1 : length(live_paths)
        olp = ind(lpc);
        sorted_live_paths(lpc).start = live_paths(olp).start;
        sorted_live_paths(lpc).end = live_paths(olp).end;
        sorted_live_paths(lpc).boxes = live_paths(olp).boxes;
        sorted_live_paths(lpc).scores = live_paths(olp).scores;
        sorted_live_paths(lpc).allScores = live_paths(olp).allScores;
        sorted_live_paths(lpc).pathScore = live_paths(olp).pathScore;
        sorted_live_paths(lpc).foundAT = live_paths(olp).foundAT;
        sorted_live_paths(lpc).count = live_paths(olp).count;
        sorted_live_paths(lpc).lastfound = live_paths(olp).lastfound;
    end
end

% -------------------------------------------------------------------------
function gap_filled_paths = fill_gaps(paths,gap)
% -------------------------------------------------------------------------
gap_filled_paths = struct();
if isfield(paths,'boxes')
    g_count = 1;

    for lp = 1 : getPathCount(paths)
        if length(paths(lp).foundAT)>gap
            gap_filled_paths(g_count).start = paths(lp).foundAT(1);
            gap_filled_paths(g_count).end = paths(lp).foundAT(end);
            gap_filled_paths(g_count).pathScore = paths(lp).pathScore;
            gap_filled_paths(g_count).foundAT = paths(lp).foundAT;
            gap_filled_paths(g_count).count = paths(lp).count;
            gap_filled_paths(g_count).lastfound = paths(lp).lastfound;
            count = 1;
            i = 1;
            while i <= length(paths(lp).scores)
                diff_found = paths(lp).foundAT(i)-paths(lp).foundAT(max(i-1,1));
                if count == 1 || diff_found == 1
                    gap_filled_paths(g_count).boxes(count,:) = paths(lp).boxes(i,:);
                    gap_filled_paths(g_count).scores(count) = paths(lp).scores(i);
                    gap_filled_paths(g_count).allScores(count,:) = paths(lp).allScores(i,:);
                    i = i + 1;
                    count = count + 1;
                else
                    for d = 1 : diff_found
                        gap_filled_paths(g_count).boxes(count,:) = paths(lp).boxes(i,:);
                        gap_filled_paths(g_count).scores(count) = paths(lp).scores(i);
                        gap_filled_paths(g_count).allScores(count,:) = paths(lp).allScores(i,:);
                        count = count + 1;
                    end
                    i = i + 1;
                end
            end
            g_count = g_count + 1;
        end
    end
end

% -------------------------------------------------------------------------
function [sorted_live_paths,dead_paths,dp_count] = sort_live_paths(live_paths,...
    path_order_score,dead_paths,dp_count,gap)
% -------------------------------------------------------------------------

sorted_live_paths = struct();
[~,ind] = sort(path_order_score,'descend');
lpc = 0;
for lp = 1 : getPathCount(live_paths)
    olp = ind(lp);
    if live_paths(ind(lp)).lastfound < gap
        lpc = lpc + 1;
        sorted_live_paths(lpc).boxes = live_paths(olp).boxes;
        sorted_live_paths(lpc).scores = live_paths(olp).scores;
        sorted_live_paths(lpc).allScores = live_paths(olp).allScores;
        sorted_live_paths(lpc).pathScore = live_paths(olp).pathScore;
        sorted_live_paths(lpc).foundAT = live_paths(olp).foundAT;
        sorted_live_paths(lpc).count = live_paths(olp).count;
        sorted_live_paths(lpc).lastfound = live_paths(olp).lastfound;
    else
        dp_count = dp_count + 1;
        dead_paths(dp_count).boxes = live_paths(olp).boxes;
        dead_paths(dp_count).scores = live_paths(olp).scores;
        dead_paths(dp_count).allScores = live_paths(olp).allScores;
        dead_paths(dp_count).pathScore = live_paths(olp).pathScore;
        dead_paths(dp_count).foundAT = live_paths(olp).foundAT;
        dead_paths(dp_count).count = live_paths(olp).count;
        dead_paths(dp_count).lastfound = live_paths(olp).lastfound;

    end
end

% -------------------------------------------------------------------------
function score = score_of_edge(v1,v2,iouth,costtype)
% -------------------------------------------------------------------------

N2 = size(v2.boxes,1);
score = zeros(1,N2);

% try
bounds1 = [v1.boxes(end,1:2) v1.boxes(end,3:4)-v1.boxes(end,1:2)+1];
% catch
%     fprintf('catch here')
% end
bounds2 = [v2.boxes(:,1:2) v2.boxes(:,3:4)-v2.boxes(:,1:2)+1];
iou = inters_union(bounds1,bounds2);

for i = 1 : N2

    if iou(i)>=iouth

        scores2 = v2.scores(i);
        scores1 = v1.scores(end);
        score_similarity = sqrt(sum((v1.allScores(end,:) - v2.allScores(i,:)).^2));
        if strcmp(costtype, 'score')
            score(i) =  scores2;
        elseif strcmp(costtype, 'scrSim')
            score(i) = 1-score_similarity;
        elseif strcmp(costtype, 'scrMinusSim')
            score(i) = scores2 + (1 - score_similarity);
        end

    end

end

% -------------------------------------------------------------------------
function lp_count = getPathCount(live_paths)
% -------------------------------------------------------------------------

if isfield(live_paths,'boxes')
    lp_count = length(live_paths);
else
    lp_count = 0;
end

% -------------------------------------------------------------------------
function iou = inters_union(bounds1,bounds2)
% -------------------------------------------------------------------------

inters = rectint(bounds1,bounds2);
ar1 = bounds1(:,3).*bounds1(:,4);
ar2 = bounds2(:,3).*bounds2(:,4);
union = bsxfun(@plus,ar1,ar2')-inters;

iou = inters./(union+eps);

author: gurkirt

repo: realtime-action-detection


text:
function[]=inticks
%INTICKS  Sets the 'tickdir' property of the current axis to 'in'.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 J.M. Lilly --- type 'help jlab_license' for details
set(gca,'tickdir','in')

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function egrad = getPartialEuclideanGradient(problem, x, I, storedb, key)
% Computes the Euclidean gradient of a subset of terms in cost function.
%
% function egrad = getPartialEuclideanGradient(problem, x, I)
% function egrad = getPartialEuclideanGradient(problem, x, I, storedb)
% function egrad = getPartialEuclideanGradient(problem, x, I, storedb, key)
%
% Assume the cost function described in the problem structure is a sum of
% many terms, as
%
%    f(x) = sum_i f_i(x) for i = 1:d,

% where d is specified as d = problem.ncostterms.
%
% For a subset I of 1:d, getPartialEuclideanGradient obtains the Euclidean
% gradient of the partial cost function
%
%    f_I(x) = sum_i f_i(x) for i = I.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% See also: getGradient canGetPartialEuclidean Gradient getPartialGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 28, 2016
% Contributors:
% Change log:

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % Make sure I is a row vector, so that it is natural to loop over it
    % with " for i = I ".
    I = (I(:)).';

    if isfield(problem, 'partialegrad')
    %% Compute the partial Euclidean gradient using partialegrad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.partialegrad)
            case 2
                egrad = problem.partialegrad(x, I);
            case 3
                % Obtain, pass along, and save the store for x.
                store = storedb.getWithShared(key);
                [egrad, store] = problem.partialegrad(x, I, store);
                storedb.setWithShared(store, key);
            case 4
                % Pass along the whole storedb (by reference), with key.
                egrad = problem.partialegrad(x, I, storedb, key);
            otherwise
                up = MException('manopt:getPartialEuclideanGradient:badpartialegrad', ...
                    'partialegrad should accept 2, 3 or 4 inputs.');
                throw(up);
        end

    else
    %% Abandon computing the partial Euclidean gradient.

        up = MException('manopt:getPartialEuclideanGradient:fail', ...
            ['The problem description is not explicit enough to ' ...
             'compute the partial Euclidean gradient of the cost.']);
        throw(up);

    end

end

author: NicolasBoumal

repo: manopt


text:
function [THz] = MHz2THz(MHz)
% Convert frequency from megahertz to terahertz.
% Chad A. Greene 2012
THz = MHz*1e-6;

author: Sable

repo: mcbench-benchmarks


text:
%
% This code creates the plot of the D-value versus the b-value.
% The code is called from view_Dv.
%
%
[bvg5, ord] = sort(bvg(:,5));
bvg1 = bvg(ord,1);
bx = [bvg5(3):0.05: bvg5(end-1)]';
Dy = 2*bx;
clear ord;

figure_w_normalized_uicontrolunits('Numbertitle','off','Name','D versus b');
plot(bvg(:,5),bvg(:,1),'o');
hold on;
plot(bx, Dy, 'k-');
xlabel('b-value');
ylabel('D-value');
title('D-value versus b-value');

reg = [ones(size(bvg,1),1), bvg(:,5)];
[sl, cint, res, resint, stat] = regress(bvg(:,1), reg, 0.25);

sl
stat

rsl = [sl(2,1) sl(1,1)];
coef1 = [cint(2,1), cint(1,1)];
coef2 = [cint(2,2), cint(1,2)];
deltar = sl(2,1) - cint(2,1);
[line] = polyval(rsl,bx);
[line1] = polyval(coef1, bx);
[line2] = polyval(coef2, bx);

hold on;
plot(bx, line, 'r');
plot(bx, line1, 'g');
plot(bx, line2, 'g');

str2 = ['Correlation factor (D=xb): x =  ' sprintf('%.2f',sl(2,1)) '  +/- ' sprintf('%.2f', deltar)];
axes('pos',[0 0 1 1]); axis off; hold on;
te1 = text(0.15, 0.85, str2,'Fontweight','bold');

author: CelsoReyes


text:
function result=SideWindowBoxFilter(im, radius, iteration)
%papers: 1) Sub-window Box Filter, Y.Gong, B.Liu, X.Hou, G.Qiu, VCIP2018, Dec.09, Taiwan
%        2) Side Window Filtering, H.Yin, Y.Gong, G.Qiu. CVPR2019
%implemented by Yuanhao Gong

r = radius; %the radius of the side window
k = ones(2*r+1,1)/(2*r+1); %separable kernel
k_L=k; k_L(r+2:end)=0; k_L = k_L/sum(k_L); %half kernel
k_R=flipud(k_L);
m = size(im,1)+2*r; n = size(im,2)+2*r; total = m*n;
[row, col]=ndgrid(1:m,1:n);
offset = row + m*(col-1) - total;
im = single(im);
result = im;
d = zeros(m,n,8,'single');

for ch=1:size(im,3)
    U = padarray(im(:,:,ch),[r,r],'replicate');
    for i = 1:iteration
        %all projection distances
        d(:,:,1) = conv2(k_L, k_L, U,'same') - U;
        d(:,:,2) = conv2(k_L, k_R, U,'same') - U;
        d(:,:,3) = conv2(k_R, k_L, U,'same') - U;
        d(:,:,4) = conv2(k_R, k_R, U,'same') - U;
        d(:,:,5) = conv2(k_L, k, U,'same') - U;
        d(:,:,6) = conv2(k_R, k, U,'same') - U;
        d(:,:,7) = conv2(k, k_L, U,'same') - U;
        d(:,:,8) = conv2(k, k_R, U,'same') - U;

        %find the minimal signed distance
        tmp = abs(d);
        [~,ind] = min(tmp,[],3);
        index = offset+total*ind;
        dm = d(index); %signed minimal distance
        %update
        U = U + dm;
    end
    result(:,:,ch) = U(r+1:end-r,r+1:end-r);
end

author: YuanhaoGong

repo: SideWindowFilter


text:
function value = box_01_contains_point_nd ( dim_num, p )

%*****************************************************************************80
%
%% BOX_01_CONTAINS_POINT_ND determines if a point is inside a unit box in ND.
%
%  Discussion:
%
%    A unit box is assumed to be a rectangle with sides aligned on coordinate
%    axes.  It can be described by its low and high corner, P1 and P2:
%
%      0 <= P(1:DIM_NUM) <= 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real P(DIM_NUM), the point to be checked.
%
%    Output, logical VALUE, is TRUE if P is inside the box.
%
  value = 0;

  for i = 1 : dim_num
    if ( p(i) < 0.0 | 1.0 < p(i) )
      return
    end
  end

  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% function [ARSignal StatWT] = WaveletAnalysis(StatWT,L,wavename,iqr,SignalLength)
%
%
% Perform a wavelet motion correction of the dod data and computes the
% distribution of the wavelet coefficients. It sets the coefficient
% exceeding iqr times the interquartile range to zero, because these are probably due
% to motion artifacts. It applies the inverse discrete wavelet transform and
% reconstruct the signal.
%
% The algorithm follows in part the procedure described by
% Molavi et al.,Physiol Meas, 33, 259-270 (2012).
%
% INPUTS:
% StatWT:       matrix of wavelet coefficients (# of time points x # of
%               levels+1). The first column contains the approximation
%               coefficients, while the other all the details coefficients at
%               different levels
% L:            Lowest wavelet scale used in the analysis
% wavename:     name of the wavelet used in the reconstruction, should be the
%               same used in the previous decomposition
% iqr:          parameter used to compute the statistics (iqr = 1.5 means 1.5 times the
%               interquartile range and is usually used to detect outliers).
%               Increasing it, it will delete less coefficients.
% SignalLength: Length of the original signal before zero padding
%
%
% OUTPUTS:
% ARSIgnal:     signal reconstructed after the discrete inverse wavelet transform
%               and corrected for motion artifacts.
% StatWT:       matrix of wavelet coefficients corrected for motion
%               artifacts. Same size as StatWT input
%
% LOG:
% Script by Behnam Molavi bmolavi@ece.ubc.ca adapted for Homer2 by RJC
% modified 10/17/2012 by S. Brigadoi

function [ARSignal,StatWT]  = WaveletAnalysis(StatWT,L,wavename,iqr,SignalLength)

n=size(StatWT,1);       % Length of data vector with zero padding
N=log2(size(StatWT,1)); % Finest scale (original signal)
SignalLength_tmp = SignalLength;

for j=1:N-L-1
    SignalLength_tmp = fix(SignalLength_tmp/2);
    n_blocks = 2^j; % number of blocks in the level
    l_blocks = n/n_blocks; % length of the blocks in the level
    for b=0:(2^j-1)
        sr = StatWT(b*l_blocks+1:b*l_blocks+l_blocks,j+1);

        sr_temp = sr(1:SignalLength_tmp); % compute statistics only on original data
        quants = quantile(sr_temp,[.25 .50 .75]);  % compute quantiles
        IQR = quants(3)-quants(1);  % compute interquartile range
        prob1 = quants(3)+IQR*iqr;
        prob2 = quants(1)-IQR*iqr;
        outliers_1 = find(sr>prob1);
        outliers_2 = find(sr<prob2);
        outliers = [outliers_1' outliers_2'];
        sr(outliers) = 0;  % set outliers to 0
        StatWT(b*l_blocks+1:b*l_blocks+l_blocks,j+1) = sr;
    end
end
ARSignal=IWT_inv(StatWT,wavename);  % reconstruct the signal with the discrete inverse wavelet transform

author: BUNPC


text:
function hilite_virtual_blocks(modelName)
% This function highlight the virtual blocks that are not going to
% contribute to the simulation of a system.
%
% Useful to study the impact of optimization settings in the model
% configuration parameters.
%
% This function potentially can be used for studying model coverage.
%-----------------------------------------------------------------------
% Developed by: Neelakanda Bharathiraja
% Email: syslancer@gmail.com
%-----------------------------------------------------------------------

open_system(modelName);

% Set the start function call back which adds the listener to the blocks.
oldStartFcn = get_param(modelName, 'StartFcn');
newStartFcn = sprintf('%s\n%s ', oldStartFcn, 'callback_details');
set_param(modelName, 'StartFcn', newStartFcn);

% Just start and stop the model to have the callback run.
set_param(modelName,'SimulationCommand', 'start');
set_param(modelName, 'SimulationCommand', 'Stop');

% Cleanup the changes done in the model properties.
set_param(modelName, 'StartFcn', oldStartFcn);

uiwait(msgbox(['The virtual blocks are highlighted. Please review it.' char(10) ...
    'The model has unsaved changes, but they won''t affect the model behaviour.' char(10) ...
    'You can also close and open the model without saving it.'],'Highlight Virtual Blocks','modal'));

return;


text:
function stInfo=convertMatToStruct(allData, gtMaxFrame)

numCols=size(allData,2);
numLines=size(allData,1);

% quickly check format
assert(numCols>=7,'FORMAT ERROR: Each line must have at least 6 values');
% assert(all(allData(:,1)>0),'FORMAT ERROR: Frame numbers must be positive.');
% assert(all(allData(:,2)>0),'FORMAT ERROR: IDs must be positive.');

IDMap = unique(allData(:,2))';
nIDs = length(IDMap);
% fprintf('%d unique targets IDs discovered\n',nIDs);
% pause

% do we have bbox coordinates?
imCoord=0;
if any(any(allData(:,3:6)~= -1))
    imCoord=1;
end

% do we have world coordinates?
worldCoord=0;

% ... we do if they exist and are not -1 (default)
if size(allData,2)>=9
    if any(any(allData(:,[8 9])~= -1))
        worldCoord=1;
    end
end

% at least one representation must be defined
assert(imCoord || worldCoord, ...
    'FORMAT ERROR: Neither bounding boxes nor world coorsinates defined.');

% go through all lines
for l=numLines:-1:1
    lineData=allData(l,:);

    fr = lineData(1);   % frame number
    id = lineData(2);   % target id

    % map id to 1..nIDs
    id = find(IDMap==id);

    %%%% sanity checks
    % ignore non-positive frames and IDs
    if fr<1,
        continue;
    end

    % bounding box
    stInfo.W(fr,id) = lineData(5);
    stInfo.H(fr,id) = lineData(6);
    stInfo.Xi(fr,id) = lineData(3) + stInfo.W(fr,id)/2;
    stInfo.Yi(fr,id) = lineData(4) + stInfo.H(fr,id);

    % values should not be exactly 0
    if ~stInfo.W(fr,id)
        stInfo.W(fr,id) = stInfo.W(fr,id) + 0.0001;
    end
    if ~stInfo.H(fr,id)
        stInfo.H(fr,id) = stInfo.H(fr,id) + 0.0001;
    end
    if ~stInfo.Xi(fr,id)
        stInfo.Xi(fr,id) = stInfo.Xi(fr,id) + 0.0001;
    end
    if ~stInfo.Yi(fr,id)
        stInfo.Yi(fr,id) = stInfo.Yi(fr,id) + 0.0001;
    end
    % consider 3D coordinates
    if worldCoord
        stInfo.Xgp(fr,id) = lineData(7);
        stInfo.Ygp(fr,id) = lineData(8);

        % position should not be exactly 0
        if ~stInfo.Xgp(fr,id)
            stInfo.Xgp(fr,id)=stInfo.Xgp(fr,id)+0.0001;
        end
        if ~stInfo.Ygp(fr,id)
            stInfo.Ygp(fr,id)=stInfo.Ygp(fr,id)+0.0001;
        end
    end
end

% append empty frames?
if gtMaxFrame ~= -1

    Fgt=gtMaxFrame;
    F=size(stInfo.W,1);
    % if stateInfo shorter, pad with zeros
    if F<Fgt
        missingFrames = F+1:Fgt;
        stInfo.Xi(missingFrames,:)=0;
        stInfo.Yi(missingFrames,:)=0;
        stInfo.W(missingFrames,:)=0;
        stInfo.H(missingFrames,:)=0;
        if worldCoord
            stInfo.Xgp(missingFrames,:)=0;
            stInfo.Ygp(missingFrames,:)=0;
        end
    end
end

% set X,Y
stInfo.X=stInfo.Xi;stInfo.Y=stInfo.Yi;
if ~imCoord
    stInfo.X=stInfo.Xgp;stInfo.Y=stInfo.Ygp;

    % reset image coordinates to -1
    te = find(stInfo.X(:));
    stInfo.W(te)=-1;stInfo.H(te)=-1;
    stInfo.Xi(te)=-1;stInfo.Yi(te)=-1;
end

% remove empty target IDs
nzc=~~sum(abs(stInfo.Xi),1);

if isfield(stInfo,'X')
    stInfo.X=stInfo.X(:,nzc);
    stInfo.Y=stInfo.Y(:,nzc);
end

% nzc
% stInfo.Xgp'
if isfield(stInfo,'Xgp')
    stInfo.Xgp=stInfo.Xgp(:,nzc);
    stInfo.Ygp=stInfo.Ygp(:,nzc);
end

if isfield(stInfo,'Xi')
    stInfo.Xi=stInfo.Xi(:,nzc);
    stInfo.Yi=stInfo.Yi(:,nzc);
    stInfo.W=stInfo.W(:,nzc);
    stInfo.H=stInfo.H(:,nzc);
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function B = repelm(A,varargin)
%REPELM Replicate elements in an array.
%
%   REPELM(A,M,N) and REPELM(A,[M N]) replicates each element in the matrix
%   A so it becomes an M-by-N matrix.
%
%   REPELM(A,M,N,P,...) and REPELM(A,[M N P ...]) replicates each element in
%   A so it becomes an M-by-N-by-P-by-... array.  A can be N-D.
%
%   REPELM(A,M,N) when A is a matrix is the same as KRON(A,ONES(M,N)), but
%   the former is much faster since it does no multiplications.  The former
%   also works when A is not of class double.
%
%   Example:
%       repelm(magic(2),2,3)
%       repelm(NaN,2,3)
%
%   See also REPMAT, MESHGRID, NDGRID.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:21:10 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

if nargin < 2
   error( 'Not enough input arguments.' );
elseif nargin == 2
   if length(varargin{1}) == 1          % REPELM(A,M)
      siz = [ varargin{1} varargin{1} ];
   else                                 % REPELM(A,[M N P ...])
      siz = varargin{1};
   end
else                                    % REPELM(A,M,N,P,...)
   siz = [ varargin{:} ];
end

if length(A) == 1
   nelems = prod(siz);
   if nelems > 0
      % Since B doesn't exist, the first statement creates a B with
      % the right size and type.  Then use scalar expansion to
      % fill the array.  Finally reshape to the specified size.
      B(nelems) = A;
      B(:) = A;
      B = reshape(B,siz);
   else
      B = A(ones(siz));
   end
elseif ndims(A)==2 & length(siz)==2
   [m,n] = size(A);
   mind = 1:m;
   nind = 1:n;
   mind = mind(ones(1,siz(1)),:);
   nind = nind(ones(1,siz(2)),:);
   B = A(mind,nind);
else

%   Asiz = size(A);
%   Asiz = [Asiz ones(1,length(siz)-length(Asiz))];
%   siz = [siz ones(1,length(Asiz)-length(siz))];
%   for i = 1:length(Asiz)
%      ind = 1:Asiz(i);
%      subs{i} = ind(ones(1,siz(i)),:);
%   end
%   B = A(subs{:});

   Asiz = size(A);
   Adim = length(Asiz);
   Rdim = length(siz);
   Bdim = max(Adim, Rdim);
   Asiz  = [Asiz ones(1,Bdim-Adim)];
   subs = {':'};
   subs = subs(ones(1,Bdim));
   for i = 1:Rdim
      if siz(i) > 1
         ind = 1:Asiz(i);
         subs{i} = ind(ones(1,siz(i)),:);
      end
   end
   B = A(subs{:});

end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function test_suite=test_find_local_extrema
% tests for cosmo_find_local_extrema
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite

function test_find_local_maxima_basics
    % generate tiny dataset with 6 voxels
    ds=cosmo_synthetic_dataset('ntargets',1,'nchunks',1);
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
    % find local maxima within neighborhood of 1 voxel radius
    [feature_ids,scores]=cosmo_find_local_extrema(ds,nh);
    assertEqual(feature_ids,[1 5 3]);
    assertElementsAlmostEqual(scores,[2.0317 1.1908 -1.4437],...
                                        'absolute',1e-4);

    % only return two feature ids
    [feature_ids,scores]=cosmo_find_local_extrema(ds,nh,'count',2);
    assertEqual(feature_ids,[1 5]);
    assertElementsAlmostEqual(scores,[2.0317 1.1908],...
                                        'absolute',1e-4);

    % use another fitness function, namely local minima
    [feature_ids,scores]=cosmo_find_local_extrema(ds,nh,'fitness',@min);
    assertEqual(feature_ids,[3 4]);
    assertElementsAlmostEqual(scores,[-1.4437 -0.5177 ],...
                                        'absolute',1e-4);

    nh=cosmo_spherical_neighborhood(ds,'radius',2,'progress',false);
    [feature_ids,scores]=cosmo_find_local_extrema(ds,nh);
    assertEqual(feature_ids,[1 6]);
    assertElementsAlmostEqual(scores,[2.0317 -1.3265],...
                                        'absolute',1e-4);

    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
    ds.samples(:)=NaN;
    [feature_ids,scores]=cosmo_find_local_extrema(ds,nh);
    assertEqual(feature_ids,zeros(1,0));
    assertElementsAlmostEqual(scores,zeros(1,0));

function test_find_local_maxima_exceptions
    ds=cosmo_synthetic_dataset('ntargets',1,'nchunks',1);
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_find_local_extrema(varargin{:}),'');
    aet(cosmo_stack({ds,ds}),nh)
    ds.fa.i=ds.fa.i(end:-1:1);
    aet(ds,nh)

    ds.fa.i=ds.fa.i(end:-1:1);

    nh.a=struct();
    ds.samples=1;
    aet(ds,nh);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [gammaAngle, betaAngle, alphaAngle] = getAnglesFromInertialBodyAxes_Func(dcm, ut, rVect, vVect, bodyInfo, baseFrame)
    ce = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());
    [~, ~, ~, base_frame_2_inertial] = baseFrame.getOffsetsWrtInertialOrigin(ut, ce);

    angles = rotm2eulARH(base_frame_2_inertial' * dcm, 'zyx');

    angles = real(angles);

    gammaAngle = angles(3);
    betaAngle = angles(2);
    alphaAngle = angles(1);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% gray2RGBColorMap
% Below is a demonstration of the features of the |gray2RGBColorMap| function

%%
clear; close all; clc;

%% Syntax
% |[C]=gray2RGBColorMap(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function Ne = som_neighbors(sM,neigh)

% Ne = som_neighbors(sM,neigh)
%
% sM      (struct) map or data struct
%         (matrix) data matrix, size n x dim
% [neigh] (string) 'kNN' or 'Nk' (which is valid for a SOM only)
%                  for example '6NN' or 'N1'
%                  default is '10NN' for a data set and 'N1' for SOM
%
% Ne      (matrix) size n x n, a sparse matrix
%                  indicating the neighbors of each sample by value 1
%                  (note: the unit itself also has value 0)

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

if isstruct(sM),
  switch sM.type,
   case 'som_map',  M = sM.codebook;
   case 'som_data', M = sM.data; sM = [];
  end
else
  M = sM;
  sM = [];
end

n = size(M,1);

if nargin<2,
  if isempty(sM), neigh = '10NN'; else neigh = 'N1'; end
end

if strcmp(neigh(end-1:end),'NN'),
  k  = str2num(neigh(1:end-2));
  kmus = som_bmus(M,M,1:k+1);
  Ne = sparse(n,n);
  for i=1:n, Ne(i,kmus(i,:)) = 1; end
else
  if ~isstruct(sM), error('Prototypes must be in a map struct.'); end
  k  = str2num(neigh(2:end));
  N1 = som_unit_neighs(sM);
  Ne = sparse(som_neighborhood(N1,k)<=k);
end
Ne([0:n-1]*n+[1:n]) = 0; % remove self from neighbors

return;

author: ilarinieminen

repo: SOM-Toolbox


text:
function f = oka1( x )

	x1p = cos(pi./12.0).*x(:,1) - sin(pi./12.0).*x(:,2);
	x2p = sin(pi./12.0).*x(:,1) + cos(pi./12.0).*x(:,2);

	f(:,1) = x1p;
	f(:,2) = sqrt(2*pi) - sqrt(abs(x1p)) + 2 .* (abs(x2p-3.*cos(x1p)-3).^0.33333333);
end

author: Eric-Bradford

repo: TS-EMO


text:
function leap = halton_leap_get ( )

%*****************************************************************************80
%
%% HALTON_BASE_GET gets the leap vector of the leaped Halton subsequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer LEAP(1:DIM_NUM), the Halton leap vector.
%
  global halton_BASE
  global halton_LEAP
  global halton_DIM_NUM
  global halton_SEED
  global halton_STEP

  dim_num = halton_DIM_NUM;

  if ( ~halham_dim_num_check ( dim_num ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALTON_LEAP_GET - Fatal error!\n' );
    error ( 'HALTON_LEAP_GET - Fatal error!' );
  end

  if ( length ( halton_LEAP ) ~= dim_num )

    halton_LEAP = [];

    for i = 1 : dim_num
      halton_LEAP(i) = 1;
    end

  end

  leap(1:dim_num) = halton_LEAP(1:dim_num);

  return
end


text:
function params = retinoGetParams(vw, dt, scan)
% Get visual field mapping ("retinotopic", though it doesn't need
% to be) parameters for the selected scan.
%
% params = retinoGetParams(vw, <dt, scan>);
%
% dt and scan specify the data type and scans to use, and
% default to the view's current dt/scan if omitted.
%
% The parameters are stored in
%   dataTYPES(dt).blockedAnalysisParameters(scan).visualFieldMap.
%
% Use retinoSetParams to update the parameters.
%
%
% ras, 01/06: testing the waters if this code is needed. I see many
% other places where similar parameters are set, but none of them
% seem immediately useable to me.
if notDefined('vw'),    vw = getCurView;               end
if notDefined('dt'),    dt = viewGet(vw, 'curdt');     end
if notDefined('scan'),  scan = viewGet(vw, 'curscan'); end

mrGlobals;

if ischar(dt),    dt = existDataType(dt);               end

if checkfields(dataTYPES(dt), 'blockedAnalysisParams', 'visualFieldMap') && ...
    ~isempty(dataTYPES(dt).blockedAnalysisParams(scan).visualFieldMap)
        params = dataTYPES(dt).blockedAnalysisParams(scan).visualFieldMap;
else
    params = [];
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function z = cheby_t_poly_zero ( n )

%*****************************************************************************80
%
%% CHEBY_T_POLY_ZERO returns zeroes of Chebyshev polynomials T(n,x).
%
%  Discussion:
%
%    The I-th zero of T(n,x) is cos((2*I-1)*PI/(2*N)), I = 1 to N
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the polynomial.
%
%    Output, real Z(N), the zeroes of T(N)(X).
%
  for i = 1 : n
    angle = ( 2 * i - 1 ) * pi / real ( 2 * n );
    z(i) = cos ( angle );
  end

  return
end


text:
function A = y2a(y);

% ABCD = y2a(Y)
%
% Admittance to ABCD transformation
% only for 2x2xN matrices

% v1.1 - 03.02.2005 freq added

d = y(1,1,:).*y(2,2,:) - y(1,2,:).*y(2,1,:);

for i=1:size(y,3)
    while abs(y(2,1,i)) < 1e-8
        y(2,1,i) = y(2,1,i)*(1+rand*1e-8);
    end;
end;

A(1,1,:) = -y(2,2,:)./y(2,1,:);
A(1,2,:) = -1./y(2,1,:);
A(2,1,:) = -d./y(2,1,:);
A(2,2,:) = -y(1,1,:)./y(2,1,:);


text:
classdef ParadigmRWSRCSP < ParadigmBase
    % Riemannian Weighted Subjects Regularized CSP (RWSRCSP)
    %
    % This paradigm implements RWSRCSP [1], which is a generalization of the Common Spatial Patterns
    % algorithms to calibration data comprising multiple subjects (or recordings). To train a model
    % for a particular "goal" (or target) subject using auxiliary data from other subjects, this
    % algorithm blends the subject-specific covariance matrices used in standard CSP and LDA with
    % weighted averages of covariance matrices from other subjects. The weighting is determined
    % based on the riemannian distance between the target subject and each respective other subject.
    % Also, the method can optionally average predictions obtained with different values of the
    % regularization parameter instead of doing a (costly) parameter search as suggested in [1].
    %
    % References:
    % [1] Lotte, F.
    %     "Signal processing approaches to minimize or suppress calibration time in oscillatory activity-based Brain-Computer Interfaces",
    %     Proceedings of the IEEE, vol. 103, no. 6, pp. 871-890, 2015
    %
    % Name:
    %   Riemannian Weighted Subjects Regularized Common Spatial Patterns
    %
    %                            Christian Kothe, Syntrogi
    %                            2015-07-24

    methods

        function defaults = preprocessing_defaults(self)
            defaults = {'FIRFilter',{'Frequencies',[6 8 28 32],'Type','minimum-phase'}, 'EpochExtraction',[0.5 3.5], 'Resampling',100};
        end

        function defaults = machine_learning_defaults(self)
            % set up the default parameters for machine learning; not necessary
            defaults = {'lda'};
        end

        function model = calibrate(self,varargin)
            % calibrate an SSRCSP model from a corpus of training sets
            args = arg_define(varargin, ...
                arg_norep({'collection','Collection'}), ...
                arg_norep({'goal_identifier','GoalIdentifier'}), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'Number of CSP patterns (times two).'),...
                arg({'lambdas','Lambdas'},0.1:0.1:0.9,[0 1],'Covariance shrinkage. A range of shrinkage parameters to run over (classifiers will be averaged).'),...
                arg({'group_tasks_by','GroupTasksBy'},'subject',{'group','subject','day','montage','session','recording','block'},'Group tasks into. This allows to group the training data into tasks solved in a multi-task manner, e.g., such that data of a given subject forms one task. When hyper-parameters need to be optimized, this would usually be done using a basic blockwise cross-validation within each task.'), ...
                arg({'shrinkage_cov','ShrinkageCovariance','ShrinkageCov'},false,[],'Shrinkage covariance estimator. Whether to use shrinkage to estimate the covariance matrices.'), ...
                arg({'weight_bias','WeightedBias'}, false, [], 'Account for class priors in bias. If you do have unequal probabilities for the different classes, this should be enabled.'), ...
                arg({'weight_cov','WeightedCov'}, false, [], 'Account for class priors in covariance. If you do have unequal probabilities for the different classes, it makes sense to enable this.'), ...
                arg({'apply_to','ApplyTo'},'channels',{'channels','sources','components','full CSD'},'Apply classifier to. Allows to select the type of time series to apply this model to.'), ...
                arg_sub({'flt','SignalProcessing'}, self.preprocessing_defaults(), @flt_pipeline, 'Signal processing stages. These parameters control filter stages that run on the signal level; they can be enabled, disabled and configured for the given paradigm. The prediction operates on the outputs of this stage.'), ...
                arg_sub({'ml','MachineLearning'},{'Learner',self.machine_learning_defaults()},@ml_train,'Machine learning stage of the paradigm. Operates on the feature vectors that are produced by the feature-extraction stage.'),...
                arg({'arg_dialogsel','ConfigLayout'},self.dialog_layout_defaults(),[],'Parameters displayed in the config dialog. Cell array of parameter names to display (dot-notation allowed); blanks are translated into empty rows in the dialog. Referring to a structure argument lists all parameters of that struture, except if it is a switchable structure - in this case, a pulldown menu with switch options is displayed.','type','cellstr','shape','row'));

            % if this is run on a worker, we'll set the cache capacity to zero since no machine
            % has enough RAM to hold multiple workers' copies of the corpus in memory
            on_worker = hlp_iscaller('par_worker');
            if on_worker
                global tracking;
                tracking.cache.capacity = 0;
            end

            % split data into reference data (of goal subject) and remaining data
            [refsets,remaining] = utl_collection_closest(args.collection,args.goal_identifier);

            % recombine and move ref data to the beginning (because this paradigm will later
            % extract the weights learned for the first task) and turn into struct array for
            % convenience
            corpus = [refsets{:} remaining{:}];

            % determine group membership
            group_membership = {corpus.(args.group_tasks_by)};
            if iscellstr(group_membership)
                groups = unique(group_membership);
            else
                groups = num2cell(unique([group_membership{:}]));
            end

            % initialize per-class covariance matrices
            CovCSP = deal(cell(2,length(groups)));
            CovLDA = deal(cell(2,length(groups)));
            % for each group...
            for s=length(groups):-1:1
                matches = find(cellfun(@(x)isequal(x,groups{s}),group_membership));
                matchdata = cell(1,length(matches));

                % collect all matching data sets...
                for p=1:length(matches)
                    matchdata{p} = corpus(matches(p));
                    if length(matchdata{p}.streams) > 1
                    disp_once('Note: ParadigmRWCSP will use only the first data stream of a recording (no support for multi-modal data).'); end
                    matchdata{p} = matchdata{p}.streams{1};
                end

                % concatenate them into a single set
				procdata = set_concat(matchdata{:});
				% and preprocess the result (with further settings/overrides according to args.flt)
                procdata = flt_pipeline('Signal',procdata, args.flt); %#ok<*NODEF>
                if on_worker
                    % if we're running on a worker we don't cache the result due to memory
                    % constraints
                    procdata = exp_eval(procdata);
                else
                    procdata = exp_eval_optimized(procdata);
                end
                % extract data
                switch args.apply_to
                    case 'channels'
                        X = procdata.data;
                    case 'components'
                        X = reshape((procdata.icaweights*procdata.icasphere)*procdata.data(procdata.icachansind,:),[],procdata.pnts,procdata.trials);
                    case 'sources'
                        X = procdata.srcpot;
                    case 'full CSD'
                        X = procdata.srcpot_all;
                end
                X(~isfinite(X(:))) = 0;
                for k=1:2
                    trials{k} = exp_eval(set_picktrials(procdata,'rank',k));
                    % calculate the CSP class covariance matrix
                    if args.shrinkage_cov
                        CovCSP{k,s} = hlp_diskcache('featuremodels',@cov_shrink,reshape(trials{k}.data,size(trials{k}.data,1),[])');
                    else
                        CovCSP{k,s} = cov(reshape(trials{k}.data,size(trials{k}.data,1),[])');
                    end
                    CovCSP{k,s}(~isfinite(CovCSP{k,s})) = 0;
                end
                % solve CSP
                [V,D] = eig(CovCSP{1,s},CovCSP{1,s}+CovCSP{2,s}); %#ok<NASGU>
                featuremodel.filters = V(:,[1:args.patterns end-args.patterns+1:end]);
                features = self.feature_extract(procdata,featuremodel); % nT x nF
                % get the LDA covariance matrix per class
                %targets = [signal.epoch.target];
                for k=1:2
                    CovLDA{k,s} = cov(features([trials{k}.epoch.targnum],:)); end
            end

%             % check for and remove bad data
%             remove = [];
%             for s=1:length(features)
%                 if size(features{s},3) < 2
%                     fprintf('Encountered bad data at subject %s/%i.\n',corpus(s).streams{1}.parts{2}.parts{1:2});
%                     remove(end+1) = s;
%                 end
%             end
%             if ~isempty(remove)
%                 fprintf('Removing bad data...\n');
%                 scales(remove) = [];
%                 features(remove) = [];
%                 targets(remove) = [];
%                 transforms(remove) = [];
%                 % data_weights(remove) = [];
%             end

            % calculate Riemannian weights for each subject
            Covs = vertcat(CovCSP,CovLDA);
            % for each type of covariance matrix...
            for c=size(Covs,1):-1:1
                Ctarg = Covs{c,1};

                % calc distance to all other covs
                for s=size(Covs,2):-1:2
                    Cother = Covs{c,s};
                    [V,D] = eig(Ctarg\Cother);
                    dist(s) = sqrt(sum(log(diag(D)).^2));
                end

                % calc weighting
                sumdist = sum(dist(2:end));
                for s=2:size(Covs,2)
                    weights(s) = 1./(dist(s)/sumdist); end

                % use it to avg the covs
                Avg = zeros(length(Ctarg));
                for s=2:size(Covs,2)
                    Avg = Avg + weights(s)*Covs{c,s}; end
                AvgCovs{c} = Avg;
            end

            % for each reg. param, solve a classifier
            signal = procdata;
            targets = [signal.epoch.target];
            classes = unique(targets);
            for li = length(args.lambdas):-1:1
                lam = args.lambdas(li);
                for c=size(Covs,1):-1:1
                    BlendCovs{c} = lam*Covs{c,1} + (1-lam)*AvgCovs{c}; end
                [V,D] = eig(BlendCovs{1},BlendCovs{1}+BlendCovs{2}); P = inv(V); %#ok<NASGU>
                % train CSP part
                model.featuremodel.filters{li} = V(:,[1:args.patterns end-args.patterns+1:end]);
                model.featuremodel.patterns{li} = P([1:args.patterns end-args.patterns+1:end],:);
                % extract features
                trials = self.feature_extract(signal,model.featuremodel,li);
                % train LDA part
                for c = 1:2
                    X = trials(targets==classes(c),:);
                    n{c} = size(X,1);
                    mu{c} = mean(X,1);
                    sig{c} = BlendCovs{c+2};
                end
                ns = quickif(args.weight_cov,n,{1 1});
                nb = quickif(args.weight_bias,n,{1 1});
                % do the math
                mu_both = (mu{1}*nb{2} + mu{2}*nb{1}) / (nb{1}+nb{2});
                sig_both = (sig{1}*ns{1} + sig{2}*ns{2}) / (ns{1}+ns{2});
                w = (mu{2} - mu{1}) / sig_both;
                w = w / (mu{2}*w' - mu_both*w');
                model.predictivemodel.model{li} = struct('w',{w}, 'b',{mu_both*w'}, 'classes',{classes});
            end

            % set the filter graph based on the reference data
            model.tracking.filter_graph = signal;
            % also store channel locations for model visualization
            model.chanlocs = signal.chanlocs;
            model.classes = classes;
        end

        function predictions = predict(self,bundle,model)
            % for each lambda
            for m=length(model.featuremodel.filters):-1:1
                % extract features
                features = self.feature_extract(bundle.streams{1},model.featuremodel.filters{m});
                % apply classifier
                raw_preds(:,m) = features*model.predictivemodel.model{m}.w' - model.predictivemodel.model{m}.b;
            end
            % average all predictions
            raw_labels = mean(raw_preds,2);
            raw_labels = min(+1,max(-1,raw_labels));
            predictions = {'disc', [(1-raw_labels)/2 1-(1-raw_labels)/2], model.classes};
        end

        function features = feature_extract(self,signal,featuremodel,lam)
            if isstruct(featuremodel)
                featuremodel = featuremodel.filters; end
            if iscell(featuremodel)
                featuremodel = featuremodel{lam}; end
            % extract log-variance features from an epoched and preprocessed recording
            features = zeros(size(signal.data,3),size(featuremodel,2));
            for t=1:size(signal.data,3)
                features(t,:) = sum((signal.data(:,:,t)'*featuremodel).^2,1); end
            features = log(features/size(signal.data,2));
        end

        function visualize(self,varargin) %#ok<*INUSD>
            % visualize an mklCSP model
            args = arg_define(varargin, ...
                arg_norep({'model','Model'},[],[],'BCI Model to visualize.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'paper','PaperFigure'},false,[],'Use paper-style font sizes. Whether to generate a plot with font sizes etc. adjusted for paper.'));

            f = figure;
            % get number of pairs, and index of pattern per subplot
            np = size(args.model.featuremodel.patterns,1)/2;
            idx = [1:np 2*np:-1:np+1];
            % for each CSP pattern...
            for p=1:np*2
                subplot(2,np,p,'Parent',f);
                if args.patterns
                    topoplot(args.model.featuremodel.patterns(idx(p),:),args.model.featuremodel.chanlocs);
                else
                    topoplot(args.model.featuremodel.filters(:,idx(p)),args.model.featuremodel.chanlocs);
                end
                t = title(['CSP Pattern ' num2str(idx(p))]);
                if args.paper
                    set(t,'FontUnits','normalized');
                    set(t,'FontSize',0.1);
                end
            end
        end

        function layout = dialog_layout_defaults(self)
            % define the default configuration dialog layout
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.FIRFilter.Frequencies', ...
                'SignalProcessing.FIRFilter.Type', 'SignalProcessing.EpochExtraction', '', ...
                'PatternPairs', 'CovarianceShrinkage', '', 'MachineLearning.Learner'};
        end

    end
end

% (turn off a few editor warnings because some actual implementations are missing in this file)
%#ok<*INUSD,*STOUT,*MANU>


text:
%%
% Copyright (c) 2019-present, Mahmoud Afifi
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
%
% Please, cite the following paper if you use this code:
%
% Mahmoud Afifi and Michael S. Brown. What else can fool deep learning?
% Addressing color constancy errors on deep neural network performance.
% ICCV, 2019
%
% Email: mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%%

classdef synthWBmodel
% WB emulator model
    properties
        features % training features
        mappingFuncs9 % training mapping functions 9x3 poly
        K % K value for KNN
        encoder % PCA object
        wb_photo_finishing % WB & photo finishing (PF) styles
    end
    methods

        function feature = encode(obj,hist)
		% Generates a compacted feature of a given histogram tensor.
            feature =  obj.encoder.encode(hist);
        end

        function hist = RGB_UVhist(obj,I)
		% Generates an RGB-histogram tensor of a given image I.
            I = im2double(I);
            if size(I,1)*size(I,2) > 202500 % if input image > (450*450),
                factor = sqrt(202500/(size(I,1)*size(I,2))); % scale factor
                newH = floor(size(I,1)*factor); % new dimensions
                newW = floor(size(I,2)*factor);
                I = imresize(I,[newH,newW],'nearest'); % scale it down
            end

            h= sqrt(max(size(obj.encoder.weights,1),...
                size(obj.encoder.weights,2))/3); % histogram dimension
            eps= 6.4/h; % threshold
            I=(reshape(I,[],3)); % reshape input image
            A=[-3.2:eps:3.19]; % dummy vector
            hist=zeros(size(A,2),size(A,2),3); % histogram tensor will be stored here
            i_ind=I(:,1)~=0 & I(:,2)~=0 & I(:,3)~=0; % remove zero pixels
            I=I(i_ind,:);
            Iy=sqrt(I(:,1).^2+I(:,2).^2+I(:,3).^2); % intensity vector
            for i = 1 : 3 % for each layer in the histogram
                r = setdiff([1,2,3],i); % extract the current color channel
                Iu=log((I(:,i))./(I(:,r(1)))); % Iu vector
                Iv=log((I(:,i))./(I(:,r(2)))); % Iv vector
                diff_u=abs(Iu-A); % differences in u space
                diff_v=abs(Iv-A); % differences in v space
                diff_u=(reshape((reshape(diff_u,[],1)<=eps/2),...
                    [],size(A,2))); % set 1's for all pixels below the threshold
                diff_v=(reshape((reshape(diff_v,[],1)<=eps/2),...
                    [],size(A,2))); % the same in the v space
                hist(:,:,i)= ... % hist = Iy .* diff_u' * diff_v (.* element-wise mult)
                    (Iy.*double(diff_u))'*double(diff_v);
                hist(:,:,i)= ... %final hist is sqrt(hist/sum(hist))
                    sqrt(hist(:,:,i)/sum(sum(hist(:,:,i))));
            end
            hist = imresize(hist,[h h],'bilinear');
        end

        function [synthWBimages,wb_pf] = generate_wb_srgb (obj,I, ...
                outNum, feature, sigma)
		% Generates outNum new images from a given image I, where outNum should be <=10
            I = im2double(I); % convert to double
            if nargin == 2
                outNum = length(obj.wb_photo_finishing); % use all WB & PF styles (default)
                feature = obj.encode(obj.RGB_UVhist(I)); % encode histogram of I
                sigma = 0.25; % fall-off factor
            elseif nargin == 3
                feature = obj.encode(obj.RGB_UVhist(I));
                sigma = 0.25;
            elseif nargin == 4
                sigma = 0.25;
            end

            if outNum > length(obj.wb_photo_finishing) % if selected styles > the available WB & PF styles
                error('Error: number of new images should be <= %d',...
                    length(obj.wb_photo_finishing));
            end

            if outNum ~= length(obj.wb_photo_finishing) % if selected styles < the available WB & PF styles
                inds = randperm(length(obj.wb_photo_finishing)); % randomize from the available WB & PF styles
                wb_pf = obj.wb_photo_finishing(...
                    inds(1:outNum));
            else
                wb_pf = obj.wb_photo_finishing; % otherwise, use all available WB & PF styles
                inds = [1:length(wb_pf)];
            end

            synthWBimages = zeros(size(I,1),size(I,2),size(I,3),...
                length(wb_pf)); % new images will be stored here

            [dH,idH] = pdist2(obj.features,feature,...
                'euclidean','Smallest',obj.K); % K nearest neighbor
            weightsH = exp(-((dH).^2)/(2*sigma^2)); % compute blending weights
            weightsH = weightsH/sum(weightsH); % normalize weights
            count = 1;
            for i = inds(1:outNum) % for each WB & PF style, do
                mf = sum(weightsH .* ... % compute mapping funciton
                    obj.mappingFuncs9((idH-1)*10 + i,:),1);
                mf = reshape(mf,[9,3]); % reshape it to be 9*3
                synthWBimages(:,:,:,count) = obj.change_wb(I,mf); % apply mf
                count = count + 1;
            end

        end

        function out = change_wb(obj,input, m)
		% Applies a given mapping function m to input image input.
            sz=size(input);
            input=reshape(input,[],3); % reshape image to be n*3 (n total number of pixels)
            %higher degree (N-D)
            input=obj.kernelP9(input); % raise it to a higher-dim space
            out=input*m; % apply m
            out = obj.out_of_gamut_clipping(out); % clip out-of-gamut pixels
            out=reshape(out,[sz(1),sz(2),sz(3)]); % reshape the image back to its original shape
        end

        function O=kernelP9(obj,I)
		% Kernel function: \phi(r,g,b) -> (r,g,b,r2,g2,b2,rg,rb,gb)
            O=[I,... %r,g,b
                I.*I,... %r2,g2,b2
                I(:,1).*I(:,2),I(:,1).*I(:,3),I(:,2).*I(:,3),... %rg,rb,gb
                ];
        end

        function [I] = out_of_gamut_clipping(obj,I)
		% Clips out-of-gamut pixels.
            I(I>1)=1;
            I(I<0)=0;
        end

    end
end

author: mahmoudnafifi


text:
%% Intermediate Usage (for building complicated systems)
% This document builds upon <basic.html |Basic Usage|> and explains how to
% describe the problem in a more flexible way.

%% Example Code
% We examine a portion of the last code discussed in <basic.html |Basic Usage|>:
[E, H, obj_array, src_array] = maxwell_run(...
	'OSC', 1e-9, 1550, ...
	'DOM', {'vacuum', 'none', 1.0}, [-1100 1100; -1100 2600; 0 10], 10, BC.p, [100 100 0],...
	'OBJ', {'CRC/Ag', 'k'}, ...
		Box([-1100 -80; 0 1000; 0 10]), ...
		Box([80 1100; 0 1000; 0 10]), ...
	'SRCJ', PlaneSrc(Axis.y, -500, Axis.x), ...
	true);

%% Passing Variables as Arguments
% You can store any arugments of <../comp/maxwell_run.html |maxwell_run|> in
% variables and pass the variables instead of the arguments to |maxwell_run|.
% This increases the readability of the code.  For example, the meaning of the
% last logical argument |true| in the above code is hard to figure out from the
% code itself.  We can define a variable |inspect_only| and modify the code as:
inspect_only = true;
[E, H, obj_array, src_array] = maxwell_run(...
	{PARAMETER GROUPS}, ...
	inspect_only);

%%%
% Now, it is more obvious that |inspect_only = true| makes |maxwell_run| inpect
% the problem without solving it.
%
% Similarly, you can define a variable for the color code used for a material.
% Because the actual color of silver is closer to gray than black, let's define
% a RGB vector and substitute it for the color code 'k' used for the material
% |'CRC/Ag'|:
inspect_only = true;
gray = [0.5 0.5 0.5];  % [r g b]
[E, H, obj_array, src_array] = maxwell_run(...
	{OTHER PARAMETER GROUPS}, ...
	'OBJ', {'CRC/Ag', gray}, ...
		Box([-1100 -80; 0 1000; 0 10]), ...
		Box([80 1100; 0 1000; 0 10]), ...
	inspect_only);

%%%
% Upon execution, the above code generates the following figure describing the
% simulation domain:
%
% <<../img/intermediate_01.png>>
%
% Compare this with the last figure in <basic.html |Basic Usage|> and note that
% the two silver boxes are now drawn in gray.

%% Building Shapes Beforehand
% You can build the shapes of objects outside |maxwell_run|, store them in
% variables, and pass the variables to |maxwell_run|:
inspect_only = true;
shape1 = Box([-1100 -80; 0 1000; 0 10]);
shape2 = Box([80 1100; 0 1000; 0 10]);
gray = [0.5 0.5 0.5];  % [r g b]
[E, H, obj_array, src_array] = maxwell_run(...
	{OTHER PARAMETER GROUPS}, ...
	'OBJ', {'CRC/Ag', gray}, shape1, shape2, ...
	inspect_only);

%%%
% Defining shape variables before |maxwell_run| this way could be useful when,
% for example, you want to calculate the size of the shapes using complicated
% equations.
%
% You can even create an array of shapes and pass it as a single argument:
inspect_only = true;
shape1 = Box([-1100 -80; 0 1000; 0 10]);
shape2 = Box([80 1100; 0 1000; 0 10]);
shape_array = [shape1, shape2];
gray = [0.5 0.5 0.5];  % [r g b]
[E, H, obj_array, src_array] = maxwell_run(...
	{OTHER PARAMETER GROUPS}, ...
	'OBJ', {'CRC/Ag', gray}, shape_array, ...
	inspect_only);

%%%
% Passing an array of shapes as an argument can be useful when the number of
% shapes is not fixed.

%% Complete Code
% As shown above, |maxwell_run| allows flexible ways to construct arguments.  In
% fact, you can construct most arguments outside |maxwell_run|, store them in
% variables, and pass the variables to |maxwell_run|.  Below is a complete
% example that does so.
inspect_only = true;

% Geometric Parameters
L0 = 1e-9;  % length unit
wvlen = 1550;  % wavelength
dL = 10;  % grid size
Lpml = 10*dL;  % PML thickness
xn = -1000-Lpml; xp = 1000+Lpml;  % x boundaries
yn = -1000-Lpml; yp = 2500+Lpml;  % y boundaries
zn = 0; zp = dL;  % z boundaries
w = 160;  % slit width
h = 1000;  % metal thickness

% Materials and Shapes
vacuum = Material('vacuum', 'none', 1.0);
shape1 = Box([xn -w/2; 0 h; zn zp]);
shape2 = Box([w/2 xp; 0 h; zn zp]);
shape_array = [shape1, shape2];
gray = [0.5 0.5 0.5];  % [r g b]

% Source
y_src = -500;
polarization = Axis.x;
src = PlaneSrc(Axis.y, y_src, polarization);

% Solution
[E, H, obj_array, src_array] = maxwell_run(...
	'OSC', L0, wvlen, ...
	'DOM', vacuum, [xn xp; yn yp; zn zp], dL, BC.p, [Lpml Lpml 0],...
	'OBJ', {'CRC/Ag', gray}, shape_array, ...
	'SRCJ', src, ...
	inspect_only);

%%% See Also
% <basic.html |Basic Usage|>

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
filename='Throne_Tetrahedra_SYM';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';

nsteps = 10;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef GLMO < ALGORITHM
% <multi> <real/integer> <large/none>
% Grouped and linked mutation operator algorithm
% optimiser      --- 3 --- The optimisation method. 1 = SMPSO, 2 = NSGA-II, 3 = NSGA-III. Default = NSGA-III
% typeOfGroups   --- 2 --- Grouping method, 1 = linear, 2 = ordered, 3 = random. Default = ordered
% numberOfGroups --- 4 --- The number of varibale Groups. Default = 4

%------------------------------- Reference --------------------------------
% H. Zille, Large-scale Multi-objective Optimisation: New Approaches and a
% Classification of the State-of-the-Art, PhD Thesis, Otto von Guericke
% University Magdeburg, 2019.
% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille, Hisao Ishibuchi, Sanaz Mostaghim and Yusuke Nojima
%     "Mutation Operators Based on Variable Grouping for Multi-objective Large-scale Optimization"
%     IEEE Symposium Series on Computational Intelligence (SSCI), IEEE, Athens, Greece, December 2016
%     https://ieeexplore.ieee.org/document/7850214
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            [optimiser,typeOfGroups,numberOfGroups] = Algorithm.ParameterSet(3,2,4);
            if optimiser == 1
                GLMO_SMPSO(Algorithm,Problem,numberOfGroups,typeOfGroups);
            elseif optimiser == 2
                GLMO_NSGAII(Algorithm,Problem,numberOfGroups,typeOfGroups);
            else
                GLMO_NSGAIII(Algorithm,Problem,numberOfGroups,typeOfGroups);
            end
        end
        function GLMO_SMPSO(Algorithm,Problem,numberOfGroups,typeOfGroups)
            Population       = Problem.Initialization();
            Pbest            = Population;
            [Gbest,CrowdDis] = GLMO_UpdateGbest(Population,Problem.N);
            while Algorithm.NotTerminated(Gbest)
                Population       = GLMO_SMPSOOperator(Problem, Population, Pbest, Gbest(TournamentSelection(2,Problem.N,-CrowdDis)), numberOfGroups, typeOfGroups);
                [Gbest,CrowdDis] = GLMO_UpdateGbest([Gbest,Population],Problem.N);
                Pbest            = GLMO_UpdatePbest(Pbest,Population);
            end
        end
        function GLMO_NSGAII(Algorithm,Problem,numberOfGroups,typeOfGroups)
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = GLMO_NSGAIIEnvironmentalSelection(Population,Problem.N);
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                Offspring  = GLMO_GA(Problem, Population(MatingPool), numberOfGroups, typeOfGroups);
                [Population,FrontNo,CrowdDis] = GLMO_NSGAIIEnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
        function GLMO_NSGAIII(Algorithm,Problem,numberOfGroups,typeOfGroups)
            [Z,Problem.N] = UniformPoint(Problem.N,Problem.M);
            Population    = Problem.Initialization();
            Zmin          = min(Population(all(Population.cons<=0,2)).objs,[],1);
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,sum(max(0,Population.cons),2));
                Offspring  = GLMO_GA(Problem, Population(MatingPool), numberOfGroups, typeOfGroups);
                Zmin       = min([Zmin;Offspring(all(Offspring.cons<=0,2)).objs],[],1);
                Population = GLMO_NSGAIIIEnvironmentalSelection([Population,Offspring],Problem.N,Z,Zmin);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_matern52_to_ss(magnSigma2, lengthScale)
% CF_MATERN52_TO_SS - Matern covariance functions, nu=5/2, to state space
%
% Syntax:
%   [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_matern52_to_ss(magnSigma2, lengthScale)
%
% In:
%   magnSigma2  - Matern magnitude scale parameter (default: 1)
%   lengthScale - Matern distance scale parameter (default: 1)
%
% Out:
%   F           - Feedback matrix
%   L           - Noise effect matrix
%   Qc          - Spectral density of white noise process w(t)
%   H           - Observation model matrix
%   Pinf        - Covariance of the stationary process
%   dF          - Derivatives of F w.r.t. parameters
%   dQc         - Derivatives of Qc w.r.t. parameters
%   dPinf       - Derivatives of Pinf w.r.t. parameters
%   params      - Input and output parameter information
%
% Description:
%   This function converts one-dimensional covariance functions of
%   the Matern class to state space models. The covariance function
%   parametrization is as follows
%
%     k(tau) = magnSigma2 (1+sqrt(5) |tau|/lengthScale+5 |tau|^2/(3 lengthScale^2))
%               exp(-sqrt(3) |tau|/lengthScale),
%
%   where magnSigma2 is the magnitude scale parameter, lengthScale the
%   distance scale parameter, and tau the time difference between states,
%   tau = t-t'.
%     This function takes the covariance function parameters as inputs and
%   outputs the corresponding state space model matrices. The state space
%   model is given as follows in terms of a stochastic differential
%   equation
%
%      df(t)/dt = F f(t) + L w(t),
%
%   where w(t) is a white noise process with spectral denisty Qc. The
%   observation model for discrete observation y_k of f(t_k) at step k,
%   is as follows
%
%      y_k = H f(t_k) + r_k, r_k ~ N(0, R),
%
%   where r_k is the Gaussian measurement noise with covariance R.
%     Pinf is the stationary covariance, where the value of Pinf(i,j),
%   is defined as follows
%
%      Pinf(i,j) = E[(f_i(t)-E[f_i(t)])(f_j(t)-E[f_j(t)])],
%
%   where f_i(t) is component i of state vector f(t).
%     Derivatives: All have same form. For example, dF has the following
%   form:
%
%       dF(:,:,1) = dF/d(magnSigma2 = input parameter_1),
%       dF(:,:,i) = dF/d(input parameter_i).
%
% References:
%
%   [1] Simo Sarkka, Arno Solin, Jouni Hartikainen (2013).
%       Spatiotemporal learning via infinite-dimensional Bayesian
%       filtering and smoothing. IEEE Signal Processing Magazine,
%       30(4):51-61.
%
%   [2] Jouni Hartikainen and Simo Sarkka (2010). Kalman filtering and
%       smoothing solutions to temporal Gaussian process regression
%       models. Proceedings of IEEE International Workshop on Machine
%       Learning for Signal Processing (MLSP).
%
% See also:
%   COV_MATERN52, SPEC_MATERN52
%
% Copyright:
%   2012-2014   Arno Solin
%   2013-2014   Jukka Koskenranta
%
%  This software is distributed under the GNU General Public
%  License (version 3 or later); please refer to the file
%  License.txt, included with the software, for details.

%% Apply defaults

  % Check if magnSigm2 is given
  if nargin < 1 || isempty(magnSigma2), magnSigma2 = 1; end

  % Check if lengthScale is given
  if nargin < 2 || isempty(lengthScale), lengthScale = 1; end

%% Form state space model

  % Derived constants
  lambda = sqrt(5)/lengthScale;

  % Feedback matrix
  F = [ 0,          1,          0;
        0,          0,          1;
       -lambda^3, -3*lambda^2, -3*lambda];

  % Noise effect matrix
  L = [0; 0; 1];

  % Spectral density
  Qc = magnSigma2*400*sqrt(5)/3/lengthScale^5;

  % Observation model
  H = [1, 0, 0];

%% Stationary covariance

  % Calculate Pinf only if requested
  if nargout > 4,

    % Derived constant
    kappa = 5/3*magnSigma2/lengthScale^2;

    % Stationary covariance
    Pinf = [magnSigma2, 0,      -kappa;
            0,          kappa,  0;
            -kappa,     0,      25*magnSigma2/lengthScale^4];

  end

%% Calculate derivatives

  % Calculate derivatives only if requested
  if nargout > 5

    % Derivative of F w.r.t. parameter magnSigma2
    dFmagnSigma2    =  [0,  0,  0;
                        0,  0,  0;
                        0,  0,  0];

    % Derivative of F w.r.t parameter lengthScale
    dFlengthScale   =  [0,                          0,                  0;
                        0,                          0,                  0;
                        15*sqrt(5)/lengthScale^4,    30/lengthScale^3,   3*sqrt(5)/lengthScale^2];

    % Derivative of Qc w.r.t. parameter magnSigma2
    dQcmagnSigma2   =   400*sqrt(5)/3/lengthScale^5;

    % Derivative of Qc w.r.t. parameter lengthScale
    dQclengthScale  =   -magnSigma2*2000*sqrt(5)/3/lengthScale^6;

    % Derivative of Pinf w.r.t. parameter magnSigma2
    dPinfmagnSigma2 = Pinf/magnSigma2;

    % Derivative of Pinf w.r.t. parameter lengthScale
    kappa2 = -2*kappa/lengthScale;
    dPinflengthScale = [0,          0,       -kappa2;
                        0,          kappa2,  0;
                        -kappa2,    0,       -100*magnSigma2/lengthScale^5];

    % Stack all derivatives
    dF = zeros(3,3,2);
    dQc = zeros(1,1,2);
    dPinf = zeros(3,3,2);

    dF(:,:,1) = dFmagnSigma2;
    dF(:,:,2) = dFlengthScale;
    dQc(:,:,1) = dQcmagnSigma2;
    dQc(:,:,2) = dQclengthScale;
    dPinf(:,:,1) = dPinfmagnSigma2;
    dPinf(:,:,2) = dPinflengthScale;

  end

%% Return parameter names

  % Only return if requested
  if nargout > 8

    % Stationarity
    p.stationary = true;

    % Input parameter information
    p.in{1}.name = 'magnSigma2';  p.in{1}.default = 1; p.in{1}.opt = true;
    p.in{2}.name = 'lengthScale'; p.in{2}.default = 1; p.in{2}.opt = true;

    % Return parameter setup
    params = p;

  end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function output = minmod(a,b)
% minmod is the minmod function in (2.9b) in the paper
% m(a,b) = minmod(a,b) = (\frac{sgn(a)+sgn(b)}{2}) min(|a|,|b|)

output = (sign(a) + sign(b)) / 2 * min(abs(a), abs(b));

author: YimianDai

path: github-repos/MATLAB/YimianDai-Image-Processing-Codes-for-Easier-Understanding/Image-Processing-Codes-for-Easier-Understanding-874302799e48852624bc3760b58b46bd9360f238/src/(Physica D 1992) Nonlinear Total Variation based noise removal algorithms/version_1/support/minmod.m


text:
function dat = rebuild_volinfo_from_dat(dat, newdat)
% Will rebuild volInfo (the image space, or sometimes "mask") from a vectorized image.
% In other words, will rebuild dat.volInfo from newdat.
%
% Also resets all voxels to be significant, if a statistic image
%
% :Input:
%
%   **dat:**
%        an image_vector
%
%   **newdat:**
%        a vector that MUST be size of ENTIRE image (dat.volInfo.nvox)
%
% :Output:
%
%   **dat:**
%        dat.dat contains the non-zero values of newdat, and dat.volInfo is
%        correctly defining the image space

if length(newdat) ~= dat.volInfo.nvox
    error('newdat MUST be size of ENTIRE image (dat.volInfo.nvox).  Image_vector.reconstruct_image may be helpful');
end

wh_new = newdat ~=0 & ~isnan(newdat);

if isa(dat, 'statistic_image')
    % Rebuild fields specific to statistic_images

    dat = replace_empty(dat);
    k = size(dat.dat, 2);

    if isempty(dat.ste), dat.ste = single(zeros(size(dat.dat))); end
    if isempty(dat.sig), dat.sig = single(zeros(size(dat.dat))); end

    for i = 1:k

        p = ones(dat.volInfo.nvox, k);
        p(dat.volInfo.wh_inmask, i) = dat.p(:, i);

        ste = Inf .* ones(dat.volInfo.nvox, k);
        ste(dat.volInfo.wh_inmask, i) = dat.ste(:, i);

        sig = zeros(dat.volInfo.nvox, k);
        sig(dat.volInfo.wh_inmask, i) = dat.sig(:, i);

    end

    dat.p = p(wh_new, :);
    dat.ste = ste(wh_new, :);
    dat.sig = sig(wh_new, :);

end

% rebuild volInfo from newdat
dat.volInfo.image_indx = wh_new;
dat.volInfo.wh_inmask = find(dat.volInfo.image_indx);
dat.volInfo.n_inmask = length(dat.volInfo.wh_inmask);
[i, j, k] = ind2sub(dat.volInfo(1).dim(1:3), dat.volInfo(1).wh_inmask);
dat.volInfo(1).xyzlist = [i j k];

% rebuild dat.dat with new data
dat.dat = newdat(dat.volInfo.image_indx);

dat.removed_voxels = false(dat.volInfo.n_inmask, 1);

% rebuild volInfo.cluster
dat.volInfo.cluster = dat.volInfo.wh_inmask;
% this forces reparse_contiguous to rebuild the cluster correctly.
% maybe not ideal solution, but maybe yes -- don't know the code well enough.  Yoni 7/14

dat = reparse_contiguous(dat);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef Voigt2Tensor3Drepresentation < handle

    properties (Access = protected)
        indexTransformer
    end

    methods (Access = protected)

        function createIndexTransformer(obj)
            obj.indexTransformer = TensorVoigtIndexTransformer3D();
        end

    end

end


text:
function wn=dn(w,type);
w = w*length(w);
w = double(w);
D=sum(abs(w),2)+eps;

if type == 'ave'
    D=1./D;
    D=sparse(1:length(D),1:length(D),D);
    wn=D*w;
elseif type == 'gph'
    D=1./sqrt(D);
    D=sparse(1:length(D),1:length(D),D);
    wn=D*(w*D);
end

author: BatzoglouLabSU

repo: SIMLR


text:
function tests = test_fast_omp_chol
  tests = functiontests(localfunctions);
end

function [A, x, b, k] = problem_1()
    m = 100;
    n = 1000;
    k = 4;
    A = spx.dict.simple.gaussian_dict(m, n);
    gen = spx.data.synthetic.SparseSignalGenerator(n, k);
    % create a sparse vector
    x =  gen.biGaussian();
    b = A*x;
end

function [dict, reps, signals, k] = problem_2()
    m = 200;
    n = 1000;
    k = 10;
    s = 500;
    dict = spx.dict.simple.gaussian_dict(m, n);
    gen = spx.data.synthetic.SparseSignalGenerator(n, k, s);
    % create a sparse vector
    reps =  gen.biGaussian();
    signals = dict*reps;
end

function test_omp_chol_1(testCase)
    [A, x, b, k] = problem_1();
    A  = double(A);
    result = spx.fast.omp(A, b, k, 1e-12);
    cmpare = spx.commons.SparseSignalsComparison(x, result, k);
    cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

function test_omp_chol_2(testCase)
    [A, x, b, k] = problem_2();
    A  = double(A);
    result = spx.fast.omp(A, b, k, 1e-12);
    cmpare = spx.commons.SparseSignalsComparison(x, result, k);
    cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

author: indigits

repo: sparse-plex


text:
function ft_realtime_signalrecorder(cfg)

% FT_REALTIME_SIGNALRECORDER is an example realtime application for recording of data
% that is streaming to the buffer in real-time. It should work both for EEG and MEG.
%
% Use as
%   ft_realtime_signalrecorder(cfg)
% with the following configuration options
%   cfg.blocksize  = number, size of the blocks/chuncks that are processed (default = 1 second)
%   cfg.channel    = cell-array, see FT_CHANNELSELECTION (default = 'all')
%   cfg.bufferdata = whether to start on the 'first or 'last' data that is available (default = 'last')
%   cfg.jumptoeof  = whether to skip to the end of the stream/file at startup (default = 'yes')
%
% The source of the data, i.e. where it comes from, is configured as
%   cfg.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.datafile      = string
%   cfg.headerfile    = string
%   cfg.eventfile     = string
%   cfg.dataformat    = string, default is determined automatic
%   cfg.headerformat  = string, default is determined automatic
%   cfg.eventformat   = string, default is determined automatic
%
% The target for the data, i.e. where it goes to, is configured as
%   cfg.export.dataset    = string with the output file name
%   cfg.export.dataformat = string describing the output file format, see FT_WRITE_DATA
%
% Some notes about skipping data and catching up with the data stream:
%
% cfg.jumptoeof='yes' causes the realtime function to jump to the end
% when the function _starts_. It causes all data acquired prior to
% starting the realtime function to be skipped.
%
% cfg.bufferdata='last' causes the realtime function to jump to the last
% available data while _running_. If the realtime loop is not fast enough,
% it causes some data to be dropped.
%
% If you want to skip all data that was acquired before you start the
% RT function, but don't want to miss any data that was acquired while
% the realtime function is started, then you should use jumptoeof=yes and
% bufferdata='first'. If you want to analyze data from a file, then you
% should use cfg.jumptoeof='no' and cfg.bufferdata='first'.
%
% To stop this realtime function, you will have have to press Ctrl-C.

% Copyright (C) 2012, Robert Oostenveld
%
% $Id$

% set the default configuration options
if ~isfield(cfg, 'dataformat'),     cfg.dataformat = [];      end % default is detected automatically
if ~isfield(cfg, 'headerformat'),   cfg.headerformat = [];    end % default is detected automatically
if ~isfield(cfg, 'eventformat'),    cfg.eventformat = [];     end % default is detected automatically
if ~isfield(cfg, 'blocksize'),      cfg.blocksize = 1;        end % in seconds
if ~isfield(cfg, 'overlap'),        cfg.overlap = 0;          end % in seconds
if ~isfield(cfg, 'channel'),        cfg.channel = 'all';      end
if ~isfield(cfg, 'bufferdata'),     cfg.bufferdata = 'last';  end % first or last
if ~isfield(cfg, 'jumptoeof'),      cfg.jumptoeof = 'no';     end % jump to end of file at initialization

if ~isfield(cfg, 'dataset') && ~isfield(cfg, 'header') && ~isfield(cfg, 'datafile')
  cfg.dataset = 'buffer://localhost:1972';
end

% translate dataset into datafile+headerfile
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
cfg = ft_checkconfig(cfg, 'required', {'datafile' 'headerfile'});

% ensure that the persistent variables related to caching are cleared
clear ft_read_header

% start by reading the header from the realtime buffer
hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat, 'cache', true, 'retry', true);

% define a subset of channels for reading
cfg.channel = ft_channelselection(cfg.channel, hdr.label);
chanindx    = match_str(hdr.label, cfg.channel);
nchan       = length(chanindx);
if nchan==0
  ft_error('no channels were selected');
end

% make a copy of the header that will be passed to the writing function, update with the channel selection
writehdr          = hdr;
writehdr.nChans   = length(chanindx);
writehdr.label    = writehdr.label(chanindx);
writehdr.chantype = writehdr.chantype(chanindx);
writehdr.chanunit = writehdr.chanunit(chanindx);

% determine the size of blocks to process
blocksize = round(cfg.blocksize * hdr.Fs);
overlap   = round(cfg.overlap*hdr.Fs);

if strcmp(cfg.jumptoeof, 'yes')
  prevSample = hdr.nSamples * hdr.nTrials;
else
  prevSample = 0;
end
count = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while true

  % determine number of samples available in buffer
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat, 'cache', true);

  % see whether new samples are available
  newsamples = (hdr.nSamples*hdr.nTrials-prevSample);

  if newsamples>=blocksize

    % determine the samples to process
    if strcmp(cfg.bufferdata, 'last')
      begsample  = hdr.nSamples*hdr.nTrials - blocksize + 1;
      endsample  = hdr.nSamples*hdr.nTrials;
    elseif strcmp(cfg.bufferdata, 'first')
      begsample  = prevSample+1;
      endsample  = prevSample+blocksize;
    else
      ft_error('unsupported value for cfg.bufferdata');
    end

    % this allows overlapping data segments
    if overlap && (begsample>overlap)
      begsample = begsample - overlap;
      endsample = endsample - overlap;
    end

    % remember up to where the data was read
    prevSample  = endsample;
    count       = count + 1;
    fprintf('processing segment %d from sample %d to %d\n', count, begsample, endsample);

    % read data segment from buffer
    dat = ft_read_data(cfg.datafile, 'header', hdr, 'dataformat', cfg.dataformat, 'begsample', begsample, 'endsample', endsample, 'chanindx', chanindx, 'checkboundary', false);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % from here onward it is specific to saving of the data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if count==1
      ft_write_data(cfg.export.dataset, dat, 'header', writehdr, 'dataformat', cfg.export.dataformat, 'append', false);
    else
      ft_write_data(cfg.export.dataset, dat, 'header', writehdr, 'dataformat', cfg.export.dataformat, 'append', true);
    end

  end % if enough new samples
end % while true

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Interactive Perspective Transformation
% This program demonstrates Perspective Transformation.
%
% In this sample you will learn how to use the following OpenCV functions:
%
% * <matlab:doc('cv.getPerspectiveTransform') cv.getPerspectiveTransform>
% * <matlab:doc('cv.warpPerspective') cv.warpPerspective>
% * <matlab:doc('cv.perspectiveTransform') cv.perspectiveTransform>
%

function varargout = perspective_transform_gui(im)
    % load source image
    if nargin < 1
        img = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
    elseif ischar(im)
        img = imread(im);
    else
        img = im;
    end

    % create the UI
    h = buildGUI(img);
    if nargout > 0, varargout{1} = h; end
end

function onHelp(~,~)
    %ONHELP  Display usage help dialog

    helpdlg({
        'This program demonstrates Perspective Transformation.'
        ''
        'Drag the image corners using the mouse: Move the pointer over a'
        'vertex. The pointer changes to a circle. Click and drag the vertex'
        'to its new position.'
        ''
        'You can also drag the image itself: move the pointer inside the'
        'quadilateral. The pointer changes to a fleur shape. Click and drag'
        'the mouse to move the image.'
        ''
        'Note: you must not add or remove points from the polygon.'
        'The cv.getPerspectiveTransform function requires exactly'
        '4 points to estimate the homography.'
        'Also the function works best if no three points are collinear.'
    });
end

function onDrag(newpos, handles)
    %ONDRAG  Event handler for impoly

    % compute the perspective transform matrix from matching corners
    H = cv.getPerspectiveTransform(handles.pos, newpos);

    % trigger redraw with the new homography
    redraw(handles, H);
end

function redraw(handles, H)
    %REDRAW  Warp and repaint image

    % apply the perspective transformation on the source image
    img2 = cv.warpPerspective(handles.img, H);

    % display warped image and homography matrix
    set(handles.hImg, 'CData',img2);
    set(handles.hTxt, 'String',mat2str_latex(H));
    drawnow limitrate;
end

function showModalDialog(~,~,handles)
    %SHOWMODELDIALOG  Display dialog to edit matrix

    % prompt for matrix using a model dialog
    d = dialog('Position',[50 50 280 140], 'Resize','off', 'Name','Homography');
    movegui(d, 'center');
    uicontrol('Parent',d, 'Style','push', 'Position',[80 10 60 20], ...
        'String','Ok', 'Callback',{@onDialogClose,true});
    uicontrol('Parent',d, 'Style','push', 'Position',[140 10 60 20], ...
        'String','Cancel', 'Callback',{@onDialogClose,false});
    t = uitable(d, 'Position',[15 40 250 90], ...
        'Data',eye(3), 'ColumnWidth',{70}, 'ColumnEditable',true, ...
        'TooltipString','fill the 3x3 homography matrix');

    % wait for dialog to close before returning
    uiwait(d);

    function onDialogClose(~,~,flag)
        H = get(t, 'Data');  % get entered matrix
        delete(gcf);         % close dialog
        if flag && ~any(isnan(H(:)))
            % update impoly position
            newpos = cv.perspectiveTransform(handles.pos, H);
            setPosition(handles.hPoly, newpos);
            % transform image
            redraw(handles, H);
        end
    end
end

function str = mat2str_latex(M)
    %MAT2STR_LATEX  Convert numeric matrix to a latex table for display

    %str = mat2str(M,3);
    M = round(M, 9);  % nicely rounded numbers (for stuff like 1e-15)
    str = ['$$H = \left[\begin{array}{ccc}' ...
        sprintf('%.3g & %.3g & %.3g \\\\ ',M(1,:)) ...
        sprintf('%.3g & %.3g & %.3g \\\\ ',M(2,:)) ...
        sprintf('%.3g & %.3g & %.3g ',M(3,:)) ...
        '\end{array}\right]$$'];
end

function img = print_instructions(img)
    %PRINT_INSTRUCTIONS  Show help text on top of image

    if nargin < 1, img = zeros([512 512 3], 'uint8'); end
    opts = {'Color',[255 0 0], 'Thickness',3, 'FontScale',1.7};
    img = cv.putText(img, 'Drag the image', [50 200], opts{:});
    img = cv.putText(img, 'corners using', [50 300], opts{:});
    img = cv.putText(img, 'the mouse.', [50 400], opts{:});
end

function handles = buildGUI(img)
    %BUILDGUI  Creates the UI

    handles = struct();
    handles.img = print_instructions(img);

    % initial quadilateral (image corners), from top-left in clockwise order
    [h,w,~] = size(handles.img);
    handles.pos = [1 1; w 1; w h; 1 h];

    % display image and homography matrix
    handles.hImg = imshow(handles.img);
    handles.hTxt = text(10, 10, mat2str_latex(eye(3)), ...
        'Interpreter','latex', 'FontSize',20, 'Color','y', ...
        'HorizontalAlignment','left', 'VerticalAlignment','top');

    % create draggable polygon
    handles.hPoly = impoly(get(handles.hImg,'Parent'), handles.pos, 'Closed',true);
    setColor(handles.hPoly, 'y');
    if false
        % restrict polygon inside image limits
        setPositionConstraintFcn(hPoly, ...
            makeConstrainToRectFcn('impozy',[1 w], [1 h]));
    end

    uicontrol('Style','pushbutton', 'Position',[20 20 60 20], 'String','Help', ...
        'Callback',@onHelp);
    uicontrol('Style','pushbutton', 'Position',[80 20 60 20], 'String','Reset', ...
        'Callback',@(~,~) setPosition(handles.hPoly, handles.pos));
    uicontrol('Style','pushbutton', 'Position',[140 20 60 20], 'String','Preset 1', ...
        'Callback',@(~,~) setPosition(handles.hPoly, [100 100; w-100 100; w h; 1 h]));
    uicontrol('Style','pushbutton', 'Position',[200 20 60 20], 'String','Preset 2', ...
        'Callback',@(~,~) setPosition(handles.hPoly, [w+25 -25; 25 200; 175 h-10; w-100 h-50]));
    uicontrol('Style','pushbutton', 'Position',[260 20 60 20], 'String','Manual', ...
        'Callback',{@showModalDialog,handles});

    % set callback when dragging points
    addNewPositionCallback(handles.hPoly, @(p) onDrag(p, handles));
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function B = get_reorder(A, idx, dim)
%GET_REORDER    Returns A with one of its dimensions indexed.
%
%   B = GET_REORDER(A, IDX, DIM)
%
%   Returns A(:, ..., :, IDX, :, ..., :), where DIM determines
%   in which dimension to place the IDX.
%
%   See also SET_REORDER.

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

ndim = ndims(A);
s.type = '()';
s.subs = cell(1, ndim);
for k = 1:ndim
    if k == dim
        s.subs{k} = idx;
    else
        s.subs{k} = ':';
    end
end
B = subsref(A, s);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function Q = orth(A)
%ORTH   Array-valued CHEBFUN orthogonalization.
%   Q = ORTH(A) is an orthonormal basis for the range of the column CHEBFUN A.
%   That is, the columns of Q span the same space as the columns of A, Q'*Q = I,
%   and the number of columns of Q is the rank of A.
%
% See also NULL, SVD, RANK, QR.

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

if ( A(1).isTransposed )
	error('CHEBFUN:CHEBFUN:orth:row', ...
        'ORTH() only defined for column CHEBFUN objects.')
end

% Compute the SVD:
[U, S, V] = svd(A, 0);
s = diag(S);

% Choose a tolerance if none is given:
if ( nargin == 1 )
	tol = max(length(A)*eps(max(s)), vscale(A)*eps);
end

% Compute the rank:
r = sum(s > tol);

% Select these columns of U:
Q = extractColumns(U, 1:r);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pass = test_linearKDV( pref )
% Check that we can solve linear KDV equations.
% Alex Townsend, April 2013.

if ( nargin < 1 )
    pref = chebfunpref();
end
tol = 100*pref.cheb2Prefs.chebfun2eps;

% Simple example.
d = [-1 1 0 1];
exact = chebfun2(@(x,t) exp(-t).*exp(x),d);

N = chebop2(@(u) diffy(u) + diffx(u,3),d);
N.dbc = @(x) exp(x);
N.rbc = @(t,u) [u - exp(-t).*exp(1) ; diff(u)-exp(-t).*exp(1)];
N.lbc = @(t) exp(-t).*exp(-1);
u = N \ 0;

pass(1) = ( norm(u-exact) < tol);

% Another simple example.
d = [-1 1 0 1];
exact = chebfun2(@(x,t) exp(-t).*exp(x),d);

N = chebop2(@(u) diffy(u) + diffx(u,3),d);
N.dbc = @(x) exp(x);
N.rbc = @(t,u) [u - exp(-t).*exp(1) ; diff(u)-exp(-t).*exp(1)];
N.lbc = @(t,u) diff(u) - exp(-t).*exp(-1);
u = N \ 0;

pass(2) = ( norm(u-exact) < 2*tol);

% Different boundary conditions.
d = [-1 1 0 1];
exact = chebfun2(@(x,t) exp(-t).*exp(x),d);

N = chebop2(@(u) diffy(u) + diffx(u,3),d);
N.dbc = @(x) exp(x);
N.rbc = @(t,u) [u - exp(-t).*exp(1) ; diff(u,2)-exp(-t).*exp(1)];
N.lbc = @(t,u) diff(u) - exp(-t).*exp(-1);
u = N \ 0;

pass(3) = ( norm(u-exact) < 300*tol);

% Different boundary conditions.
d = [-1 1 0 1];
exact = chebfun2(@(x,t) exp(-t).*exp(x),d);

N = chebop2(@(u) diffy(u) + diffx(u,3),d);
N.dbc = @(x) exp(x);
N.rbc = @(t,u) [u - exp(-t).*exp(1) ; diff(u,2)-exp(-t).*exp(1)];
N.lbc = @(t,u) u - exp(-t).*exp(-1);
u = N \ 0;

pass(4) = ( norm(u-exact) < 100*tol);


text:
function [ data ] = sortWalls( data )
%SORTWALLS re-order the walls based on their orientation
% and delete the surround relations with walls

wallind = strmatch('wall',data.labellist,'exact');
wallfront = data.obblist([4,6],wallind);
wallidx = 1:size(wallfront,2);

% 1st floor
delta = wallfront - repmat([1;0],1,size(wallfront,2));
delta = delta.*delta;
delta = delta(1,:)+delta(2,:);
[~,I] = min(delta);
fid1 = wallidx(I);
ind = find([1:size(wallfront,2)]~=I);
wallfront = wallfront(:,ind);
wallidx = wallidx(ind);

% 2st floor
delta = wallfront - repmat([0;1],1,size(wallfront,2));
delta = delta.*delta;
delta = delta(1,:)+delta(2,:);
[~,I] = min(delta);
fid2 = wallidx(I);
ind = find([1:size(wallfront,2)]~=I);
wallfront = wallfront(:,ind);
wallidx = wallidx(ind);

% 3st floor
delta = wallfront - repmat([-1;0],1,size(wallfront,2));
delta = delta.*delta;
delta = delta(1,:)+delta(2,:);
[~,I] = min(delta);
fid3 = wallidx(I);
ind = find([1:size(wallfront,2)]~=I);
wallfront = wallfront(:,ind);
wallidx = wallidx(ind);

% 4st floor
delta = wallfront - repmat([0;-1],1,size(wallfront,2));
delta = delta.*delta;
delta = delta(1,:)+delta(2,:);
[~,I] = min(delta);
fid4 = wallidx(I);
ind = find([1:size(wallfront,2)]~=I);
wallfront = wallfront(:,ind);
wallidx = wallidx(ind);

wallorder = [fid1,fid2,fid3,fid4];
walllist = data.obblist(:,wallind);
data.walllist = walllist(:,wallorder);
data.obblist(:,wallind) = walllist(:,wallorder);

% update the Rmat and Dmat
wallRmat = data.Rmat(:,wallind);
wallRmat = wallRmat(:,wallorder);
data.Rmat(:,wallind) = wallRmat;
data.Rmat(wallind,:) = wallRmat';

wallDmat = data.Dmat(:,wallind);
wallDmat = wallDmat(:,wallorder);
data.Dmat(:,wallind) = wallDmat;
data.Dmat(wallind,:) = wallDmat';
%
% %% all the walls should be supported by the floor
% floorid = strmatch('floor',data.labellist,'exact');
% data.Rmat(wallind,floorid) = -2;
% data.Rmat(floorid,wallind) = -1;

%% delete the surround relations with walls
% nsur = {};
% for i = 1:length(data.sur)
%     list = data.sur{i};
%     ind = intersect(list,wallind);
%     if(length(ind)==0)
%         nsur{length(nsur)+1} = list;
%     end
% end
% data.sur = nsur;

%% floor is not supported by any object
% ind = find(data.Rmat(:,floorid)==-1|data.Rmat(:,floorid)==-3);
% data.Rmat(ind,floorid) = -2;
% data.Rmat(floorid,ind) = -1;

end

repo: GRAINS