text:
function [clusters, subcl] = extract_contrast_data(P, clusters, varargin)
% This function does not plot,  but separates subclusters using pca / cluster_princomp.m
% based on pattern across all conditions,  covariance (not correlation),
%
% :Usage:
% ::
%
%    function [clusters, subcl] = extract_contrast_data(P, clusters)
%
% :Inputs:
%
%   **P:**
%        cell array of strings containing image file names (data extracted from these)
%
%   **clusters:**
%
% :Optional Inputs:
%
%   **subclusters:**
%        to get sub-clustering based on pca and clustering of voxels
%        add the string 'subclusters' as an input argument
%
%   **cell array:** of strings with condition names
%
%   **split:**
%        value is 1: to split into 2 plots (first half and last half of P)
%
%        value is 2: to plot individual subjects over bars
%
%   **center:**
%        center parmeter values in plot (subtract row means)
%        this gives closer to correct "within subjects" error bars
%        and may be used when the overall parameter values have no meaning
%
%   **covs:**
%        followed by between-subject covariates (e.g.,  behavioral regressors)
%        plots remove these before plotting means and std. errors
%
%   **max:**
%        to make plots based on max z-values within region for each dataset P
%        not compatible with 'split' and 'center' (ignores these commands)
%        right now,  special for inhib - see also inhib2_cluster_barplot (good function)
%
%   **indiv:**
%        to threshold based on individual t-statistics or contrast values
%
%        FOLLOW with cell array of t-images or con images -- usually,  there will be one cell,  with images
%        for each subject in rows,  to define voxels for each ss.
%        BUT Tnames can be the same length as
%        contrast images,  one t-img per subject per contrast,  if
%        desired.
%
% :Outputs:
%
%   clusters struture,  with CONTRAST substructure added
%   substructure contains data extracted and image file names
%
% This program uses XYZmm millimeter coordinates in clusters to find voxels
% So clusters and data files may have different dimensions.
%
% :Examples:
% ::
%
%     cl = extract_contrast_data(EXPT.SNPM.P, cl, 'indiv', EXPT.FILES.Timgs{1});
%     cluster_barplot(EXPT.SNPM.P(7:12), clusters(2:3), {'ObjE' 'AttE' 'InteractE' 'ObjI' 'AttI' 'InteractI'}, 'split')
%     [clusters, subclusters] = cluster_barplot(EXPT.SNPM.P(17:24), clusters, 'subclusters', 'split')
%     RS2_8vs2_placeboCP = cluster_barplot(EXPT.SNPM.P([8 10 12 14 16]), RS2meta, 'indiv', T);
%
% also see mask2clusters.m,  a simpler version that just extracts clusters from a mask file.
%
% ..
%    by Tor Wager,  5/15/03
%    modified 3/19/04 by Tor to add individual subject plots
% ..

diary off

% -------------------------------------------------
% * Set up optional inputs
% -------------------------------------------------
splitat=0; dosubcl = 0; docenter = 0; covs = []; domax = 0; doindiv = 0;
for j = 1:length(P), cnames{j}=num2str(j); end

for i = 1:length(varargin),
    if iscell(varargin{i}),  cnames = varargin{i};  %varargin{i} = [];
    elseif strcmp('split', varargin{i}),  if varargin{i+1}==2,  splitat=2; else  splitat=round(length(P)./2);  end
    elseif strcmp('subclusters', varargin{i}),  dosubcl = 1;  disp('Finding sub-clusters')
    elseif strcmp('center', varargin{i}), docenter = 1;
    elseif strcmp('covs', varargin{i}), covs = varargin{i+1};
    elseif strcmp('max', varargin{i}), domax = 1;
    elseif strcmp('indiv', varargin{i}), doindiv = 1;
        Tnames = varargin{i+1}; varargin{i+1} = 'saved.';
    end
end

% convert to old format if necessary
if isa(clusters, 'region')
    clusters = region2struct(clusters);
end

if doindiv
    % replicate names,  in case only one set is entered
    if ~iscell(Tnames),
        Tnamestmp = Tnames; Tnames={};Tnames{1}=Tnamestmp;
    end
        if length(Tnames) < length(cnames),
            disp('Found only one set of T-images; replicating for each set of data extracted.');
        end
        while length(Tnames) < length(cnames), Tnames{end+1} = Tnames{1}; end
end

%if ~isempty(strmatch('split', varargin)),  splitat=round(length(P)./2);
%else  splitat=length(P); end
subcl = [];
if isempty(clusters)
    disp('Clusters is empty: No voxels for which to get data.');
    return
end

% -------------------------------------------------
% * load image files
% -------------------------------------------------
fprintf(1, 'Loading images.')
CLU = clusters2CLU(clusters);
dat = cell(1, length(clusters));
alldat = cell(1, length(clusters));

V = spm_vol(P{1}(1, :));
if any(any(V.mat - CLU.M)),
    disp('CLUSTERS and IMAGE DIMENSIONS DIFFERENT!  Resizing clusters.'),
    %CLU = transform_coordinates(CLU, V.mat);
    V.M = V.mat;
    for i = 1:length(clusters),
        clusters(i).XYZ = mm2voxel(clusters(i).XYZmm, V)'; clusters(i).M=V.mat;
        clusters(i).voxSize = diag(V.M(1:3, 1:3))';
        clusters(i).XYZmm = voxel2mm(clusters(i).XYZ, clusters(i).M);
        clusters(i).Z = ones(1, size(clusters(i).XYZ, 2));
    end
    CLU = clusters2CLU(clusters);
    for i = 1:length(clusters), clusters(i).u=CLU.u;clusters(i).VOX=CLU.VOX; end
    doslow = 1;
else
    doslow=0;
end

% check if any contiguous clusters we lose
%tmpCL = tor_extract_rois([], CLU, CLU);
%if length(tmpCL) < length(clusters),  doslow = 1;  end
%clear tmpCL

% reslice clusters if necessary.
if doslow
    clusters = cluster_interp(clusters, P{1}(1, :), 1);
    doslow = 0;
end

for i = 1:length(P)

    %if ~doslow
        % fast way: if all voxel sizes match
        cl{i} = tor_extract_rois(P{i}, clusters);  %CLU, CLU);
    %end

    % concatenated data for each cluster
    for j = 1:length(clusters)

        %if doslow
            % slow way: to keep clusters separate that otherwise touch
        %    cl{i}(j) = tor_extract_rois(P{i}, clusters(j), CLU);
        %end

        % alldat{cluster} is data x voxels,  data is subj within conditions strung
        % together
        alldat{j} = [alldat{j}; cl{i}(j).all_data];

        if length(cl{i}(j).timeseries) < size(dat{j}, 1),
            cl{i}(j).timeseries = [cl{i}(j).timeseries; NaN*zeros(size(dat{j}, 1)-length(cl{i}(j).timeseries), 1)];
        end

        dat{j} = [dat{j} cl{i}(j).timeseries];

    end

end % loop through P

%V{1}(1).M = V{1}(1).mat;
%imgdims = size(vols{1});
%mask = zeros(imgdims(1:3));

% -------------------------------------------------
% * subclusters
% -------------------------------------------------
for i = 1:length(clusters)

    if dosubcl && clusters(i).numVox > 4,
        % -------------------------------------------------
        % * subclusters
        % -------------------------------------------------

        [subcl{i}] = getsubcl(clusters(i), alldat{i}, length(P), size(P{1}, 1), cnames, splitat, i, docenter, covs);  % does subclusters within it
        if length(subcl{i}) > 1,
            subcluster_montage(subcl{i});    % now plot all subcluster locations
            saveas(gcf, ['cl' num2str(i) '_subcl_montage'], 'fig')
            saveas(gcf, ['cl' num2str(i) '_subcl_montage'], 'tif')
        end
    else

        % no subclusters

        sterr = nanstd(dat{i}) ./ sqrt(size(dat{i}, 1));

        if isfield(clusters, 'shorttitle'), mytit = clusters(i).shorttitle;
        else   mytit = clusters(i).title;
        end

        str = sprintf('%s %2.0f (%3.0f, %3.0f, %3.0f),  %3.0f voxels',  ...
        mytit, i, clusters(i).mm_center(1),  ...
        clusters(i).mm_center(2), clusters(i).mm_center(3),  clusters(i).numVox);

        if domax
            % set up
            len = size(clusters(1).CONTRAST.data, 1);
            tmpcl = []; tmpcl2 = [];
            for j = 1:round(size(alldat{1}, 1) ./ len)
                tmpcl{j}(i) = clusters(i); tmpcl{j}(i).all_data = alldat{i}(1+(j-1)*len:j*len, :);
                tmpcl{j}(i).timeseries = mean(tmpcl{j}(i).all_data, 2);

                [sterr, tmp] = ste(tmpcl{j}(i).all_data);
                clusters(i).Z(j, :) = spm_t2z(tmp, size(tmpcl{j}(i).all_data, 1)-1);

                wh = find(tmp == max(tmp)); wh = wh(1);
                clusters(i).CONTRAST.grp_peakdat(j, :) = tmpcl{j}(i).all_data(:, wh);

                %tmpcl2{j}(i) = cluster_ttest(tmpcl{j}(i), covs);  % add covariate here if necessary
            end
            % needs checking and/or debugging
            %inhib2_barplot(clusters(i), tmpcl2{1}(i), tmpcl2{2}(i), tmpcl2{3}(i), 'common');
            keyboard
        end

        %f = makefigure(dat{i}, str, cnames, splitat, docenter, covs);
        %saveas(gcf, ['cl' num2str(i) '_bar'], 'fig')
        %saveas(gcf, ['cl' num2str(i) '_bar'], 'tif')

        clusters(i).CONTRAST.files = P;
        clusters(i).CONTRAST.data = dat{i};
        clusters(i).CONTRAST.all_data = alldat{i};

        if dosubcl
            disp('Fewer than 5 voxels; skipping subclustering.')

            % make all field names match.
            % won't work if this is the first one!
            subcl{i} = subcl{1}(1);
            N = fieldnames(subcl{i});
            for NN = 1:length(N)
                eval(['subcl{i}(1).' N{NN} ' = [];'])
            end

            N = fieldnames(clusters(i));
            for NN = 1:length(N)
                eval(['subcl{i}(1).' N{NN} ' = clusters(i).' N{NN}])
            end

        end

        %close all

    end

    clusters(i).covs = covs;
end

fprintf(1, '\n')

% extract spatial peak data and locations
for i = 1:length(cl)
    cltmp = extract_ind_peak([], cl{i});
    for j = 1:length(cltmp)
        clusters(j).CONTRAST.peakdata{i} = cltmp(j).peakdata;
        clusters(j).CONTRAST.peakXYZ{i} = cltmp(j).peakXYZ;
        clusters(j).CONTRAST.peakXYZmm{i} = cltmp(j).peakXYZmm;
    end
end

% extract individual significant region center of mass
if doindiv
    disp(['getting individual significance regions'])
    disp(['Found ' num2str(size(Tnames{1}, 1))  ' T-images to get regions from'])
    indpeak = [];

    reverse_vals = input('Enter 0 to get most positive voxels for each image, or 1 to get most negative ones: ');

    for c = 1:size(Tnames, 2)    % for each contrast tested

        for i = 1:size(Tnames{1}, 1) % for each subject

            % pass in clusters with all data saved for contrast c,  Tnames
            % for contrast c; timeseries saves data from this

            [tmp] = cluster_tmask(cl{c}, Tnames{c}(i, :), i, alldat, reverse_vals);
            cl{c} = tmp;

        end

        % put stuff where it belongs
        for j = 1:length(clusters)

            clusters(j).CONTRAST.indiv_data(:, c) = cl{c}(j).timeseries;
            clusters(j).CONTRAST.indivXYZ{c} = cl{c}(j).INDIV.center;
            clusters(j).CONTRAST.indivXYZmm{c} = cl{c}(j).INDIV.mm_center;
            clusters(j).CONTRAST.indiv_sig{c} = cl{c}(j).INDIV.sigt;
            clusters(j).CONTRAST.indivXYZall{c} = cl{c}(j).INDIV.XYZ;
            clusters(j).CONTRAST.tname = cl{c}(j).INDIV.tname;

        end

    end

end

return

% -----------------------------------------------------------------------------
% * sub-functions
% -----------------------------------------------------------------------------

    function [subclusters] = getsubcl(clusters, adat, numimglists, numimgs, varargin)
    % always pass in only one cluster in clusters!! (vec of length 1)
    % if > 1 output,  subclusters is created and plotted here
    % if subclusters,  pass in cnames and splitat

    if length(varargin) > 0,  cnames = varargin{1};  splitat = varargin{2};
        clnum=varargin{3};  docenter = varargin{4};
    end
    if length(varargin) > 4,  covs = varargin{5}; else  covs = []; end

            % don't scale if this is across subjects,  because voxels
            % with more variation across ss SHOULD drive the pca
            %adat{j} = scale(adat{j});   % scale so high-noise voxels
            %                            % don't drive the pca

        % separate subclusters
        subclusters.CONTRAST = [];   % add for consistency of field names
        clusters.all_data = adat;

        [clusters, subclusters] = cluster_princomp(clusters, []);   %  , 0, 0, 1);

        subclusters = subclusters{1};

        % re-format averages within sub-regions for plotting

            for j = 1:length(subclusters)
                vind = 1;
                if size(subclusters(j).all_data, 2) > 1,
                    tmp = nanmean(subclusters(j).all_data')';
                else
                    tmp = subclusters(j).all_data;
                end

                for k = 1:numimglists
                    subclusters(j).CONTRAST.dat(:, vind) = ...
                        tmp((k-1)*numimgs+1:k*numimgs);
                        vind = vind + 1;
                end

                % now plot it
                sterr = nanstd(subclusters(j).CONTRAST.dat) ./ sqrt(size(subclusters(j).CONTRAST.dat, 1));
                str = sprintf('Cl %3.0f (%4.0f vox) subcl at (%3.0f, %3.0f, %3.0f),  %3.0f voxels',  ...
                clnum, clusters.numVox, subclusters(j).mm_center(1),  ...
                subclusters(j).mm_center(2), subclusters(j).mm_center(3),  subclusters(j).numVox);

                f = makefigure(subclusters(j).CONTRAST.dat, str, cnames, splitat, docenter, covs);
                saveas(gcf, ['cl' num2str(clnum) '_subc' num2str(j) '_bar'], 'fig')
                saveas(gcf, ['cl' num2str(clnum) '_subc' num2str(j) '_bar'], 'tif')
            end

    return

    function f = makefigure(dat, varargin)

    if length(varargin) > 1, cnames=varargin{2}; else  cnames = []; end
    if length(varargin) > 2, splitat=varargin{3}; else  splitat = 0; end
    if length(varargin) > 3, docenter=varargin{4}; else  docenter = 0; end
    if length(varargin) > 4, covs = varargin{5}; covs(:, end+1) = 1; else  covs = []; end
    if covs == 1,  covs = []; end    % if empty,  make it empty!

    dat1 = dat;     % save original data for ind subj plot

    if ~isempty(covs),
        covs(:, 1:end-1) = scale(covs(:, 1:end-1));
        for i = 1:size(dat, 2),
            b = pinv(covs) * dat(:, i);
            r = dat(:, i) - covs * b + b(end);
            dat(:, i) = r;
        end
    end

    if docenter,
        disp('Centering rows of data for plotting.')
        dat = dat - repmat(nanmean(dat')', 1, size(dat, 2));

    end

    doprintt = 1;

    if doprintt
        % print t values
        for i = 1:size(dat, 2)
            [h, p, ci, stats] = ttest(dat(:, i));
            nums{i} = sprintf('%3.2f', stats.tstat);
        end
    end

    sterr = nanstd(dat) ./ sqrt(size(dat, 1));
    dat = nanmean(dat);
    mysum = sign(dat) .* (sterr + abs(dat) + .15 .* abs(dat));

    if splitat > 2,
        f = figure('Color', 'w');

        h(1) = subplot(1, 2, 1); set(gca, 'FontSize', 16); hold on; grid on;
        bar(dat(1:splitat)); tor_bar_steplot(dat(1:splitat), sterr(1:splitat), {'b'});
        set(gca, 'XTick', 1:splitat)
        xlabel('Conditions', 'FontSize', 18), ylabel('fMRI Signal', 'FontSize', 18)
        if length(varargin) > 1,  set(gca, 'XTickLabel', varargin{2}(1:splitat)), end
        if length(varargin) > 0,  title(varargin{1}, 'FontSize', 20), end

        for i = 1:splitat,  text(i-.5, mysum(i), nums{i}, 'Color', 'k', 'FontWeight', 'b', 'FontSize', 14); end

        h(2) = subplot(1, 2, 2); set(gca, 'FontSize', 16); hold on; grid on;
        bar(dat(splitat+1:end)); tor_bar_steplot(dat(splitat+1:end), sterr(splitat+1:end), {'b'});
        set(gca, 'XTick', 1:length(dat)-splitat)
        xlabel('Conditions', 'FontSize', 18)
        if length(varargin) > 1,  set(gca, 'XTickLabel', varargin{2}(splitat+1:end)), end

        for i = 1:length(dat)-splitat,  text(i, mysum(splitat+i), nums{splitat+i}, 'Color', 'k', 'FontWeight', 'b', 'FontSize', 14); end

        equalize_axes(h);
        set(gcf, 'Position', [195         308        1259         674]),  drawnow

    elseif splitat == 2

        % ------------------------
        % this plot does individual subject estimates
        % ------------------------
        barplot_columns(dat1)

        % add reg lines,  if covs and sig
        for i = 1:size(dat1, 2)
            [B, dev, stat]=glmfit(covs(:, 1), dat1(:, i)); tmp = corrcoef(covs(:, 1), dat1(:, i));
            fprintf(1, 'Cond. %3.0f: Bo: t = %3.2f,  se = %3.4f,  p = %3.4f  B1: r = %3.2f,  t = %3.2f,  se = %3.4f,  p = %3.4f\n',  ...
                i, stat.t(1), stat.se(1), stat.p(1),  tmp(1, 2),  stat.t(2), stat.se(2), stat.p(2));

        end

        if ~isempty(cnames),  set(gca, 'XTickLabel', cnames), end

    else

        % ------------------------
        % standard bar plot
        % ------------------------

        f = figure('Color', 'w'); set(gca, 'FontSize', 16); %hold on; grid on;
        h = bar(dat); set(h, 'FaceColor', [.7 .7 .7])
        tor_bar_steplot(dat, sterr, {'k'});

        doinhib2 = 0;
        if doinhib2
            % special insert for inhib2 (triple inhibition)
            set(gca, 'XTickLabel', {'GNG' 'Flanker' 'SRC' 'Saccade'})
        end

        dointext = 0;
        if dointext
            % special insert for intext
            cla
            xp = dat; h = bar([xp(1:3);xp(4:6)], 'grouped'); cm = [0 1 0;1 0 0; 0 0 1];colormap(cm)
            xe = sterr;
            tor_bar_steplot([xp(1:3)],  xe(1:3), {'k'}, .55, .225)
            tor_bar_steplot([xp(4:6)],  xe(4:6), {'k'}, 1.55, .225)
            set(gca, 'FontSize', 16, 'XTickLabel', {'External' 'Internal'});legend(h, {'Object switching' 'Attribute Switching' 'Interaction'})
            ylabel('BOLD Contrast')
        end

        %set(gca, 'XTick', 1:size(dat, 2))
        %xlabel('Conditions', 'FontSize', 18), ylabel('fMRI Signal', 'FontSize', 18)
        %if length(varargin) > 1,  set(gca, 'XTickLabel', varargin{2}), end
        if length(varargin) > 0,  title(varargin{1}, 'FontSize', 20), end

        set(gcf, 'Position', [464   283   930   827]),  drawnow

    end

    return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [mm2] = yd22mm2(yd2)
% Convert area from square yards to square millimeters.
% Chad A. Greene 2012
mm2 = yd2*836127.36;

author: Sable

repo: mcbench-benchmarks


text:
function table = postproc_norm(table)

nrows = size(table, 1);

for nri = 1:nrows
  split = 'none';
  norm = 'none';
  descriptor = table{nri, 'descriptor'}{1};
  pos = strfind(descriptor, 'nsplit');
  if ~isempty(pos)
    split = descriptor(pos+7);
    norm = descriptor(pos+9:end);
    descriptor = descriptor(1:pos-2);
  end
  table.descriptor(nri) = {descriptor};
  table.norm_split(nri) = {split};
  table.norm_type(nri) = {norm};
end

end

author: hpatches

repo: hpatches-benchmark


text:
function [xpeak, ypeak, max_f] = findpeak(f,subpixel)
%FINDPEAK Find extremum of matrix.
%   [XPEAK,YPEAK,MAX_F] = FINDPEAK(F,SUBPIXEL) finds the extremum of F,
%   MAX_F, and its location (XPEAK, YPEAK). F is a matrix. MAX_F is the maximum
%   absolute value of F, or an estimate of the extremum if a subpixel
%   extremum is requested.
%
%   SUBPIXEL is a boolean that controls if FINDPEAK attempts to estimate the
%   extremum location to subpixel precision. If SUBPIXEL is false, FINDPEAK
%   returns the coordinates of the maximum absolute value of F and MAX_F is
%   max(abs(F(:))). If SUBPIXEL is true, FINDPEAK fits a 2nd order
%   polynomial to the 9 points surrounding the maximum absolute value of
%   F. In this case, MAX_F is the absolute value of the polynomial evaluated
%   at its extremum.
%
%   Note: Even if SUBPIXEL is true, there are some cases that result
%   in FINDPEAK returning the coordinates of the maximum absolute value
%   of F:
%   * When the maximum absolute value of F is on the edge of matrix F.
%   * When the coordinates of the estimated polynomial extremum would fall
%     outside the coordinates of the points used to constrain the estimate.

%   Copyright 1993-2004 The MathWorks, Inc.
%

% get absolute peak pixel
[max_f, imax] = max(abs(f(:)));
[ypeak, xpeak] = ind2sub(size(f),imax(1));

if ~subpixel || ...
    xpeak==1 || xpeak==size(f,2) || ypeak==1 || ypeak==size(f,1) % on edge
    return % return absolute peak

else
    % fit a 2nd order polynomial to 9 points
    % using 9 pixels centered on irow,jcol
    u = f(ypeak-1:ypeak+1, xpeak-1:xpeak+1);
    u = u(:);
    x = [-1 -1 -1  0  0  0  1  1  1]';
    y = [-1  0  1 -1  0  1 -1  0  1]';

    % u(x,y) = A(1) + A(2)*x + A(3)*y + A(4)*x*y + A(5)*x^2 + A(6)*y^2
    X = [ones(9,1),  x,  y,  x.*y,  x.^2,  y.^2];

    % u = X*A
    A = X\u;

    % get absolute maximum, where du/dx = du/dy = 0
    x_offset = (-A(3)*A(4)+2*A(6)*A(2)) / (A(4)^2-4*A(5)*A(6));
    y_offset = -1 / ( A(4)^2-4*A(5)*A(6))*(A(4)*A(2)-2*A(5)*A(3));

    if abs(x_offset)>1 || abs(y_offset)>1
        % adjusted peak falls outside set of 9 points fit,
        return % return absolute peak
    end

    % return only one-tenth of a pixel precision
    x_offset = round(10*x_offset)/10;
    y_offset = round(10*y_offset)/10;

    xpeak = xpeak + x_offset;
    ypeak = ypeak + y_offset;

    % Calculate extremum of fitted function
    max_f = [1 x_offset y_offset x_offset*y_offset x_offset^2 y_offset^2] * A;
    max_f = abs(max_f);

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
%MATCHLABLIST Match entries of lablist1 with lablist2
%
%    I = MATCHLABLIST(LABLIST1,LABLIST2)
%
% INPUT
%    LABLIST1  list of class names
%    LABLIST2  list of class names
%
% OUTPUT
%    I         indices for LABLIST1 appearing in LABLIST2
%
% DESCRIPTION
% Find the indices of places where the entries of LABLIST1 appear
% in LABLIST2, i.e. LABLIST1 = LABLIST2(I).
% Note that this operation is not symmetric, changing the order of
% LABLIST1 and LABLIST2 changes I!
% I(i) = 0 for labels appearing in LABLIST1 that are not in LABLIST2.
%
% SEE ALSO (PRTools Guide)
% DATASETS, RENUMLAB

% Copyright: D.M.J. Tax davidt@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function I = matchlablist(lablist1,lablist2)
	n = size(lablist1,1);
I = zeros(n,1); % the resulting vector

if isempty(lablist2) % nothing fits
	return
end

if iscell(lablist1) & iscell(lablist2)
	lablist1 = char(lablist1);
	lablist2 = char(lablist2);
end

for i=1:n
	if isstr(lablist1) & isstr(lablist2)
  	tmp = strmatch(deblank(lablist1(i,:)),lablist2,'exact');
	elseif ~isstr(lablist1) & ~isstr(lablist2)
		tmp = find(~sum((lablist2 ~= repmat(lablist1(i),size(lablist2,1),1)),2));
	else
		tmp = zeros(size(lablist1,1),1);
	end
  if ~isempty(tmp)
		I(i) = tmp(1);
	else
		I(i) = 0;
  end
end

return

author: marianux

repo: ecg-kit

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


text:
% TERNCOORDS calculate rectangular coordinates of fractions on a ternary plot
%   [X, Y] = TERNCOORDS(FA, FB) returns the rectangular X and Y coordinates
%   for the point with a fraction defined by FA and FB.  It is assumed that
%   FA and FB are sensible fractions.
%
%   [X, Y] = TERNCOORDS(FA, FB, FC) returns the same.  FC is assumed to be
%   the remainder when subtracting FA and FB from 1.

% Author: Carl Sandrock 20050211

% Modifications

% Modifiers

function [x, y] = terncoords(fA, fB, fC)
if nargin < 3
    fC = 1 - (fA + fB);
end

y = fB*sin(deg2rad(60));
x = fA + y*cot(deg2rad(60));


text:
% open binary files to write registered tiffs to
% takes ops, ops1, and red_binary (if 1 then write red channel to binary)
function [ops1, fid, fidRED, fidIntpol] = openBinFiles(ops, ops1, red_binary)
numPlanes = numel(ops.planesToProcess);
fid    = cell(numPlanes, size(ops1,2));
fidRED = cell(numPlanes, size(ops1,2));
fidIntpol = [];
if ops.interpolateAcrossPlanes == 1 && ~isempty(ops.RegFileBinLocation)
    fidIntpol = cell(numPlanes, size(ops1,2));
end

for i = 1:numPlanes
    for j = 1:size(ops1,2)
        ops1{i,j}.RegFile = fullfile(ops.RegFileRoot, ...
            sprintf('%s_%s_%s_plane%d.bin', ops.mouse_name, ops.date, ...
            ops.CharSubDirs, i + (j-1)*numPlanes));
        regdir = fileparts(ops1{i,j}.RegFile);
        if ~exist(regdir, 'dir')
            mkdir(regdir);
        end

        % open bin file for writing
        fid{i,j}              = fopen(ops1{i,j}.RegFile, 'w');

        if red_binary
            ops1{i,j}.RegFile2 = fullfile(ops.RegFileRoot, ...
                sprintf('%s_%s_%s_plane%d_RED.bin', ops.mouse_name, ops.date, ...
                ops.CharSubDirs, i + (j-1)*numPlanes));
            fidRED{i,j}              = fopen(ops1{i,j}.RegFile2, 'w');
        end

        if ops.interpolateAcrossPlanes && ~isempty(ops.RegFileBinLocation)
            % open separate files for result after averaging across
            % neighbouring planes
            str = sprintf('%d_',ops1{i,j}.expts);
            str(end) = [];
            folder = fullfile(ops1{i,j}.RegFileBinLocation, ops1{i,j}.mouse_name, ...
                ops1{i,j}.date, str, 'interpolated');
            if ~exist(folder, 'dir')
                mkdir(folder);
            end
            fidIntpol{i,j} = fopen(fullfile(folder, ...
                sprintf('%s_%s_%s_plane%d.bin', ops.mouse_name, ops.date, ...
                ops.CharSubDirs, i + (j-1)*numPlanes)), 'w');
        end
    end
end

author: cortex-lab


text:
% Resizes subs by resizVals.
%
% Resizes subs in subs/vals image representation by resizVals.
%
% This essentially replaces each sub by sub.*resizVals.  The only subtlety
% is that in images the leftmost sub value is .5, so for example when
% resizing by a factor of 2, the first pixel is replaced by 2 pixels and so
% location 1 in the original image goes to location 1.5 in the second
% image, NOT 2.  It may be necessary to round the values afterward.
%
% USAGE
%  subs = imSubsResize( subs, resizVals, [zeroPnt] )
%
% INPUTS
%  subs        - subscripts of point locations (n x d)
%  resizVals   - k element vector of shrinking factors
%  zeroPnt     - [.5] See comment above.
%
% OUTPUTS
%  subs        - transformed subscripts of point locations (n x d)
%
% EXAMPLE
%  subs = imSubsResize( [1 1; 2 2], [2 2] )
%
%
% See also IMSUBSTOARRAY

% Piotr's Image&Video Toolbox      Version NEW
% Written and maintained by Piotr Dollar    pdollar-at-cs.ucsd.edu
% Please email me if you find bugs, or have suggestions or questions!

function subs = imSubsResize( subs, resizVals, zeroPnt )

if( nargin<3 || isempty(zeroPnt) ); zeroPnt=.5;  end

[n d] = size(subs);
[resizVals,er] = checkNumArgs( resizVals, [1 d], -1, 2 ); error(er);

% transform subs
resizVals = repmat( resizVals, [n, 1] );
subs = (subs - zeroPnt) .* resizVals + zeroPnt;

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function sigVal = computeSigmoid(z)
% This function computes the value of the sigmoid of all of the numbers in
% an n x 1 vector, where n is the length of z.
%
% Input:
%   z: The n x 1 vector, where n is the length of z, of values whose
%   sigmoids need to be found
%
% Output:
%   sigVal: The value of the sigmoid

expz = exp(z);
sigVal = expz ./ (1 + expz);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [n0, n1] = versionstr2num(s)
n0 = 0;
if iscell(s)
    c = s;
else
	c = str2cell(s,'.');
end
n1 = zeros(1,length(c));
m = length(c);
b = 100;
for ii = 1:length(c)
    n0 = n0 + (str2num(c{ii}) * b^(m-ii));
    n1(1,ii) = str2num(c{ii});
end

author: BUNPC


text:
function [lambda,x,flag,its,x0] = eig_sshopm(A,varargin)
%EIG_SSHOPM Shifted power method for finding real eigenpair of real tensor.
%
%   [LAMBDA,X]=EIG_SSHOPM(A) finds an eigenvalue (LAMBDA) and eigenvector
%   (X) for the real tensor A such that Ax^{m-1} = lambda*x.
%
%   [LAMBDA,X]=EIG_SSHOPM(A,parameter,value,...) can specify additional
%   parameters as follows:
%
%     'Shift'    : Shift for eigenvalue calculation (Default: 'Adaptive')
%     'Margin'   : Margin for positive/negative definiteness in adaptive
%                  shift caluclation. (Default: 1e-6)
%     'MaxIts'   : Maximum power method iterations (Default: 1000)
%     'Start'    : Initial guess (Default: normal random vector)
%     'Tol'      : Tolerance on norm of change in |lambda| (Default: 1e-15)
%     'Concave'  : Treat the problem as concave rather than convex.
%                  (Default: true for negative shift; false otherwise.)
%     'Display'  : Display every n iterations (Default: -1 for no display)
%
%   [LAMBDA,X,FLAG]=EIG_SSHOPM(...) also returns a flag indicating convergence.
%
%      FLAG = 0  => Succesfully terminated
%      FLAG = -1 => Norm(X) = 0
%      FLAG = -2 => Maximum iterations exceeded
%
%   [LAMBDA,X,FLAG,IT]=EIG_SSHOPM(...) also returns the number of iterations.
%
%   [LAMBDA,X,FLAG,IT,X0]=EIG_SSHOPM(...) also returns the intial guess.
%
%   REFERENCES:
%   * T. G. Kolda and J. R. Mayo, Shifted Power Method for Computing Tensor
%     Eigenpairs, SIAM Journal on Matrix Analysis and Applications
%     32(4):1095-1124, October 2011, http://dx.doi/org/10.1137/100801482
%   * T. G. Kolda and J. R. Mayo, An Adaptive Shifted Power Method for
%     Computing Generalized Tensor Eigenpairs, SIAM Journal on Matrix
%     Analysis and Applications 35(4):1563-1582, December 2014,
%     http://dx.doi.org/0.1137/140951758
%
%   See also EIG_GEAP, EIG_SSHOPMC, TENSOR, SYMMETRIZE, ISSYMMETRIC.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Error checking on A
P = ndims(A);
N = size(A,1);

if ~issymmetric(A)
    error('Tensor must be symmetric.')
end

%% Check inputs
p = inputParser;
p.addParamValue('Shift', 'adaptive');
p.addParamValue('MaxIts', 1000, @(x) isscalar(x) && (x > 0));
p.addParamValue('Start', [], @(x) isequal(size(x),[N 1]));
p.addParamValue('Tol', 1.0e-15, @(x) isscalar(x) && (x > 0));
p.addParamValue('Display', -1, @isscalar);
p.addParamValue('Concave', false, @islogical);
p.addParamValue('Margin', 1e-6, @(x) isscalar(x) && (x > 0));
p.parse(varargin{:});

%% Copy inputs
maxits = p.Results.MaxIts;
x0 = p.Results.Start;
shift = p.Results.Shift;
tol = p.Results.Tol;
display = p.Results.Display;
concave = p.Results.Concave;
margin = p.Results.Margin;

%% Check shift
if ~isnumeric(shift)
    adaptive = true;
    shift = 0;
else
    adaptive = false;
end

%% Check starting vector
if isempty(x0)
    x0 = 2*rand(N,1)-1;
end

if norm(x0) < eps
    error('Zero starting vector');
end

%% Check concavity
if shift ~= 0
    concave = (shift < 0);
end

%% Execute power method
if (display >= 0)
    fprintf('TENSOR SHIFTED POWER METHOD: ');
    if concave
        fprintf('Concave ');
    else
        fprintf('Convex  ');
    end
    fprintf('\n');
    fprintf('----  --------- ----- ------------ -----\n');
    fprintf('Iter  Lambda    Diff  |newx-x|     Shift\n');
    fprintf('----  --------- ----- ------------ -----\n');
end

flag = -2;
x = x0 / norm(x0);
lambda = x'*ttsv(A,x,-1);
if adaptive
    shift = adapt_shift(A,x,margin,concave);
end

for its = 1:maxits

    newx = ttsv(A,x,-1) + shift * x;

    if (concave)
        newx = -newx;
    end

    nx = norm(newx);
    if nx < eps,
        flag = -1;
        break;
    end
    newx = newx / nx;

    newlambda = newx'* ttsv(A,newx,-1);

    if adaptive
        newshift = adapt_shift(A,newx,margin,concave);
    else
        newshift = shift;
    end

    if norm(abs(newlambda-lambda)) < tol
        flag = 0;
    end

    if (display > 0) && ((flag == 0) || (mod(its,display) == 0))
        fprintf('%4d  ', its);
        % Lambda
        fprintf('%9.6f ', newlambda);
        d = newlambda-lambda;
        if (d ~= 0)
            if (d < 0), c = '-'; else c = '+'; end
            fprintf('%ce%+03d ', c, round(log10(abs(d))));
        else
            fprintf('      ');
        end
        % Change in X
        fprintf('%8.6e ', norm(newx-x));

        % Shift
        fprintf('%5.2f', shift);

        % Line end
        fprintf('\n');
    end

    x = newx;
    lambda = newlambda;
    shift = newshift;

    if flag == 0
        break
    end
end

%% Check results
if (display >=0)
    switch(flag)
        case 0
            fprintf('Successful Convergence');
        case -1
            fprintf('Converged to Zero Vector');
        case -2
            fprintf('Exceeded Maximum Iterations');
        otherwise
            fprintf('Unrecognized Exit Flag');
    end
    fprintf('\n');
end

%% ----------------------------------------------------

function alpha = adapt_shift(A,x,tau,concave)

m = ndims(A);
Y = ttsv(A,x,-2);
e = eig(Y);

if concave
    if max(e) <= -tau/(m^2-m)
        alpha = 0;
    else
        alpha = -(tau/m) - ((m-1)*max(e));
    end
else
    if min(e) >= tau/(m^2-m)
        alpha = 0;
    else
        alpha = (tau/m) - ((m-1)*min(e));
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function r = mrdivide(p,q)
%RDIVIDE this function divides the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if isa(p, 'FaceVariable') && (~isa(q, 'FaceVariable'))
    r = (1/q)*p;
else
    error('FVMtool: Wrong use of mrdivide for a face variable. Try using ./ instead.');
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [elem,HB] = uniformcoarsenquad(elem)
%% UNIFORMCOARSENQUAD uniform coarsening of uniform quad mesh
%
% [elem,HB] = UNIFORMCOARSENQUAD(elem) remove grid points added by uniform
% refinement. See the illustration below:
%
%
%   See also: uniformcoarsen, coarsen, bisect, uniformcoarsen3, mg
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1);
N = max(elem(:));
d = size(elem,2);
if d == 3  % triangulation
    NT = NT/2;
end
HB = [];

%% Find ni and nj
b = N - NT + 1;
c = N;
delta = sqrt(b^2-4*c);
ni = (b+delta)/2;
nj = (b-delta)/2;
if ~isequal(mod(ni+1,2),0) || ~isequal(mod(nj+1,2),0)
    display('Not from refinement of uniform quad mesh');
    return
end
if d == 3 % triangulation
    if ~isequal(elem(1,1),ni+1)
        nj = ni;
        ni = elem(1,1) - elem(1,2);
    end
elseif d == 4
    if ~isequal(elem(1,2),ni+1)
        nj = ni;
        ni = elem(1,2) - elem(1,1);
    end
end

%% Coarse grids
nic = (ni+1)/2;
njc = (nj+1)/2;
Nc = nic*njc;
nodecidx = reshape(1:Nc,nic,njc);
t2nidxMap = nodecidx(1:nic-1,1:njc-1);
k = t2nidxMap(:);
elem = [k k+nic k+nic+1 k+1];

%% Record HB
HB = zeros(N,3);
nodeidx = reshape(1:N,ni,nj);
% fine nodes on vertical lines
i = 2:2:ni-1;
j = 1:2:nj;
k = nodeidx(i,j);
HB(k,1) = k(:);
HB(k,2) = k(:) - 1;
HB(k,3) = k(:) + 1;
% fine nodes on horizontal lines
i = 1:2:ni;
j = 2:2:nj-1;
k = nodeidx(i,j);
HB(k,1) = k(:);
HB(k,2) = k(:) - ni;
HB(k,3) = k(:) + ni;
% fine nodes in the center
i = 2:2:ni-1;
j = 2:2:nj-1;
k = nodeidx(i,j);
HB(k,1) = k(:);
HB(k,2) = k(:) - ni-1;
HB(k,3) = k(:) + ni+1;
HB(HB(:,2) == 0,:) = [];
% shift the index into coarse grid
i = 1:2:ni;
j = 1:2:nj;
k = nodeidx(i,j);
indexMap = zeros(N,1);
indexMap(k(:)) = 1:Nc;
HB(:,2:3) = indexMap(HB(:,2:3));

repo: ifem


text:
% DEMO_ZILOGGAUSSIAN  Regression problem demonstration for 2-input
%                     function with zero-inflated log-Gaussian likelihood
%
%  Description
%    The regression problem consist of a data with two input
%    variables and one output variable with zero-inflated log-Gaussian
%    likelihood. The inference is conducted with Laplace approximation for
%    the conditional posterior of latent function and (Laplace approximate)
%    MAP estimate for the hyperparameters.
%
% Copyright (c) 2016 Jarno Vanhatalo

% 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.

% Load the data
S = which('demo_regression1');
L = strrep(S,'demo_regression1.m','demodata/dat.1');
data=load(L);
x = [data(:,1) data(:,2)];
y = data(:,3);
[n, nin] = size(x);

% set zero observations and transfer the rest to exp scale
y(y<0) = 0;
y(y>0) = exp(y(y>0));

% Construct the model
lik = lik_ziloggaussian('sigma2_prior',prior_t('s2',0.1));
cfc = gpcf_constant;
gpcf = gpcf_sexp('lengthScale', [1.1 1.2], 'magnSigma2', 0.2^2)
gp = gp_set('lik', lik, 'cf', {cfc gpcf cfc gpcf},'comp_cf', {1:2 3:4});

disp(' MAP estimate for the parameters')

% Optimize with the scaled conjugate gradient method
opt=optimset('TolFun',1e-3,'TolX',1e-3,'Display','iter');
gp=gp_optim(gp,x,y,'opt',opt);

% Make predictions of the underlying function on a dense grid and plot it.
[xt1,xt2]=meshgrid(-1.8:0.1:1.8,-1.8:0.1:1.8);
xt=[xt1(:) xt2(:)];
nt = size(xt,1);
[Eft_map, Varft_map, Lpy_map, Ey_map, Vary_map] = gp_pred(gp, x, y, xt, 'yt', zeros(size(xt,1),1));

% Plot the prediction and data
figure(1)
clf
subplot(1,2,1)
mesh(xt1, xt2, reshape(1-exp(Lpy_map),37,37));
hold on
plot3(x(y>0,1), x(y>0,2), ones(sum(y>0),1), 'r*')
plot3(x(y==0,1), x(y==0,2), zeros(sum(y==0)), 'k*')
axis on;
title('The probability of being greater than zero');
subplot(1,2,2)
mesh(xt1, xt2, reshape(Ey_map,37,37));
hold on
plot3(x(:,1), x(:,2), y, '*')
axis on;
title('The posterior mean of y');

% %% Test the likelihood function
% % This block is a general purpose code for testing new likelihood
% % functions
% f = 2*randn(size(y,1)*2,1);
%
% lik.fh.ll(lik,y,f)
%
% % check llg with respect to latent
% fe = @(x) lik.fh.ll(lik,y,x');
% fg = @(x) lik.fh.llg(lik,y,x','latent')';
% gradcheck(randn(size(f')),fe,fg);
%
% % check llg with respect to param
% [w,s] = lik.fh.pak(lik)
% w=randn(size(w));
% fe = @(x) lik.fh.ll(lik.fh.unpak(lik,x),y,f);
% fg = @(x) lik.fh.llg(lik.fh.unpak(lik,x),y,f,'param');
% gradcheck(randn(size(w)),fe,fg);
%
% % check llg2 with respect to latent
% ind = unique( ceil(length(y)*rand(20,1)));
% h1 = lik.fh.llg2(lik,y(ind),f([ind ; ind+n]), 'latent');
% ny=length(ind);
% h1 = [diag(h1(1:ny,1)) diag(h1(ny+1:end,1)) ; diag(h1(1:ny,2)) diag(h1(ny+1:end,2))];
% fe = @(x) lik.fh.ll(lik,y(ind),x);
% h2 = hessian(fe,f([ind ; ind+n]));
% [min(min(h1)) max(max(h1))]
% [min(min(h2)) max(max(h2))]
% [min(min(h1-h2)) max(max(h1-h2))]
%
% gpla_e(gp_pak(gp),gp,x,y)
% gpla_g(gp_pak(gp),gp,x,y)
%
% gradcheck(randn(size(gp_pak(gp))),@gpla_e,@gpla_g,gp,x,y);

author: gpstuff-dev

repo: gpstuff

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


text:
function h = imgsc(I, varargin)
%IMGSC Alias for imagesc for images.
% Usage:
%   imgsc(I)
%   imgsc(I, ...)
%
% See also: imagesc, sc

% figure('KeyPressFcn',@KeyPressFcn_cb)
figclosekey
h = imagesc(I, varargin{:});
H = size(I,1); W = size(I,2);
if  max([H W]) / min([H W]) < 2 || all([H W] < 25)
    axis image
end
colorbar

if nargout < 1; clear h; end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
% EX_STOKES_DRIVEN_CAVITY_3D_TH: solve the Stokes problem in the driven cavity with generalized Taylor-Hood elements.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined from the aspect ratio using the NURBS toolbox
aspect_ratio=[1.5 0.3];
nrb_section = nrb4surf([0 0], [1 0], [0 aspect_ratio(1)], [1 aspect_ratio(1)]);
problem_data.geo_name = nrbextrude (nrb_section, [0 0 aspect_ratio(2)]);

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:6;
problem_data.nmnn_sides = [];

% Physical parameters
problem_data.viscosity = @(x, y, z) ones (size (x));

% Force and boundary terms
problem_data.f  = @(x, y, z) zeros ([3, size(x)]);
problem_data.h  = @test_stokes_3d_symdrivcav_h_drchlt;

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.element_name = 'th';       % Element type for discretization
method_data.degree       = [2  2  2];  % Degree of the splines (pressure space)
method_data.regularity   = [1  1  1];  % Regularity of the splines (pressure space)
method_data.nsub         = [5  5  5];  % Number of subdivisions
method_data.nquad        = [4  4  4];  % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = ...
                       solve_stokes (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file = 'Driven_cavity_3d_TH_Deg2_Reg1_Sub5';

fprintf ('The result is saved in the files %s \n and %s \n \n', ...
           [output_file '_vel'], [output_file '_press']);
vtk_pts = {linspace(0, 1, 15), linspace(0, 1, 15), linspace(0, 1, 15)};
sp_to_vtk (press, space_p, geometry, vtk_pts, [output_file '_press'], 'press')
sp_to_vtk (vel,   space_v, geometry, vtk_pts, [output_file '_vel'  ], {'velocity', 'divergence'}, {'value', 'divergence'})

%!test
%! aspect_ratio=[1.5 0.3];
%! nrb_section = nrb4surf([0 0], [1 0], [0 aspect_ratio(1)], [1 aspect_ratio(1)]);
%! problem_data.geo_name = nrbextrude (nrb_section, [0 0 aspect_ratio(2)]);
%! problem_data.drchlt_sides = 1:6;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y, z) ones (size (x));
%! problem_data.f  = @(x, y, z) zeros ([3, size(x)]);
%! problem_data.h  = @test_stokes_3d_symdrivcav_h_drchlt;
%! method_data.element_name = 'th';       % Element type for discretization
%! method_data.degree       = [2  2  2];  % Degree of the splines (pressure space)
%! method_data.regularity   = [1  1  1];  % Regularity of the splines (pressure space)
%! method_data.nsub         = [5  5  5];  % Number of subdivisions
%! method_data.nquad        = [4  4  4];  % Points for the Gaussian quadrature rule
%! [geometry, msh, space_v, vel, space_p, press] = ...
%!                        solve_stokes (problem_data, method_data);
%! assert (msh.nel, 125)
%! assert (space_p.ndof, 343)
%! assert (space_v.ndof, 5184)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
% MATLAB RRT*FN package for solving path/motion planning problems.
% Version 0.0.1 23-Aug-2013
%
% RRT and RRT* implementations are included as well
%
% FILES:
%   README.md                        - readme file, read it please.
%   rrt                              - Rapidly-Exploring Random Tree (RRT)
%   rrt_star                         - RRT*
%   rrt_star_fn                      - RRT*FN (Fixed Nodes)
%   FNRedundantManipulator           - Model of "n" DOF planar redundant
%                                      manipulator.
%   FNSimple2D                       - Model of 2D mobile robot
%   GridBased2Dimrotate              - Initial work on model of KUKA youbot
%                                      GPU usage employed (faster)
%   benchmark2D                      - benchmark script for 2D mobile robot
%   benchmarkRedundant               - benchmark script for nDOF redundant planar manipulator
%
% CONFIGURATION FILES:
%   configure_FNRedundantManipulator - Default configuration for Redundant Robotic Manipulator nDOF
%   configure_FNSimple2D             - Default configuration for 2D case
%   configure_GridBased2Dimrotate    - Default configuration for 2D case

author: olzhas


text:
%MDL_ONELINK Create model of a simple 1-link mechanism
%
% MDL_ONELINK is a script that creates the workspace variable tl which
% describes the kinematic and dynamic characteristics of a simple planar
% 1-link mechanism.
%
% Also defines the vector:
%   qz   corresponds to the zero joint angle configuration.
%
% Notes::
% - SI units are used.
% - It is a planar mechanism operating in the XY (horizontal) plane and is
%   therefore not affected by gravity.
% - Assume unit length links with all mass (unity) concentrated at the joints.
%
% References::
%  - Based on Fig 3-6 (p73) of Spong and Vidyasagar (1st edition).
%
% See also SerialLink, mdl_twolink, mdl_planar1.

% MODEL: generic, planar, 1DOF, standard_DH

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

a1 = 1;

onelink = SerialLink([
          Revolute('d', 0, 'a', a1, 'alpha', 0, 'standard')
    ], ...
    'name', 'one link');
qz = [0];
qn = [pi/6];

author: Allopart

repo: rbpf-gmapping


text:
function [N,Vn]=trinorm(F,V)

%% function [N,Vn]=trinorm(F,V)
% ------------------------------------------------------------------------
%
%
% 2021/09/14 Labelled as depricated
% ------------------------------------------------------------------------

%%

warning('The trinorm function is depricated. Use patchNormal instead.');
[N,Vn]=patchNormal(F,V);

%% Old

% %N.B. if F does not describe triangles this functions uses first three
% %vertices as a triangular description
%
%
% %Getting triangle surface normal (cross product of two edge vectors)
% vec1=[V(F(:,2),1)-V(F(:,1),1)  V(F(:,2),2)-V(F(:,1),2)  V(F(:,2),3)-V(F(:,1),3)];
% vec2=[V(F(:,3),1)-V(F(:,1),1)  V(F(:,3),2)-V(F(:,1),2)  V(F(:,3),3)-V(F(:,1),3)];
% N=cross(vec1,vec2,2);
%
% %Normalizing vector length
% N=N./(sqrt(sum(N.^2,2))*ones(1,size(N,2)));
%
% %Midface coordinates for normal vectors (mean of each face)
% X=V(:,1); Y=V(:,2); Z=V(:,3);
% Vn=[mean(X(F),2) mean(Y(F),2) mean(Z(F),2)];

end

%%
% _*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 [ds, S, p] = multivar_dist(X, varargin)
% multivariate normality checking and diagnostic plots
%
% :Usage:
% ::
%
%     [ds, S, p] = multivar_dist(X)
%
% :Input:
%
%     given matrix X with cases = rows, cols = variables
%
% :Optional input:
%
% 'noplot' : suppress plot
%
% :Outputs:
%
%   **ds:**
%        is matrix of squared distances, case numbers, and
%        expected chi2 values (in columns in this order)
%        rows are cases
%
%        NOTE: Sorted in order of ascending distance!
%
%   **S:**
%        estimated covariance matrix
%
%   **mv_distance:**
%        squared distances in original order of rows
%
%   **p:**
%        p-values in original order of rows
%
% ..
%    by Tor Wager
% ..

    % ..
    %    determine multivariate standard deviation matrix S
    % ..

    doplot = true;
    doverbose = true;

    if any(strcmp(varargin, 'noplot')), doplot = false; end
    if any(strcmp(varargin, 'noverbose')), doverbose = false; end

    % center
    Xs = X - repmat(mean(X), size(X, 1), 1);

    % covariance matrix S
    S = (Xs' * Xs) ./ (size(Xs, 1)-1);

    % -----------------------------------------------------
    % * get squared distance
    % -----------------------------------------------------

    % squared distance, Johnson & Wichern p. 201
    % (X-mean(X))' * inv(S) * (X - mean(X)) generalized to matrices
    d = Xs * inv(S) * Xs';
    d = diag(d);            % what do the off-diagonals signify?

    % -----------------------------------------------------
    % * compare with chi2 distribution
    % -----------------------------------------------------

    % calculate chi2 threshold
    % chi2 value compares the number of points within the ellipsoid
    % contour to those outside; roughly alpha of the squared distances
    % should be within the ellipsoid (for rough general test of normality).
    % outliers will have very high chi2 values

    t = chi2inv(.5, size(X, 2));  % for general test
    p50 = 100 * (sum(d > t) ./ length(d));

    if doverbose
        fprintf(1, 'Expected 50%% of points within 50%% normal ellipsoid, found %3.2f%%\n', p50);
    end

    t = chi2inv(.95, size(X, 2));  % for general test
    p95 = sum(d > t);

    if doverbose
        fprintf(1, 'Expected %3.2f outside 95%% ellipsoid, found %3.0f\n', .05*length(d), p95);
    end

    % -----------------------------------------------------
    % * get case numbers and sort by distance
    % -----------------------------------------------------
    d(:,2) = (1:length(d))';
    ds = sortrows(d, 1);

    % -----------------------------------------------------
    % * get chi2 quantiles for qd plot
    % -----------------------------------------------------
    q = (((1:size(ds, 1)) - .5) ./ size(ds, 1))';
    q = chi2inv(q, size(X, 2));
    ds(:,3) = q;

    % -----------------------------------------------------
    % * get distance^2 in original order and p-values
    % -----------------------------------------------------
    ds = sortrows(ds, 2);
    p = 1 - chi2cdf(ds(:,1), size(X, 2));

    % -----------------------------------------------------
    % * plot the results in a figure
    % -----------------------------------------------------
    if doplot

        figure('color', 'w');
        subplot(1, 3, 1); hold on; grid on
        plot([1:size(d, 1); 1:size(d, 1)], [zeros(size(d, 1), 1) d(:,1)]', 'b', 'LineWidth', 1.5);
        xlabel('Case number');
        ylabel('Squared stat. distance from origin');
        wh = (d(:,1) > t); d2 = d; d2(:,1) = d2(:,1) .* wh;  % zero out the non-"significant" chi2 cases
        plot([1:size(d2, 1); 1:size(d2, 1)], [zeros(size(d2, 1), 1) d2(:,1)]', 'r', 'LineWidth', 1.5);
        title('d^2, red cases outside 95% normal ellipsoid');
        plot([0 size(d, 1)], [t t], 'k');
        set(gca, 'YLim', [0 max(t+.5, max(d(:,1)))]);

        subplot(1, 3, 2); hold on; grid on
        plot(ds(:,3), ds(:,1), 'MarkerSize', 0.01, 'Color', 'w');
        for i = 1:size(ds, 1)
            text(ds(i,3), ds(i, 1), num2str(ds(i, 2)), 'Color', 'k');
        end
        xlabel('Expected chi2 value'), ylabel('Squared distance');
        plot([0 max([ds(:,3);ds(:,1)])], [0 max([ds(:,3);ds(:,1)])], 'k', 'LineWidth', 1.5);
        title('Line with slope = 1 is normal');

        subplot(1, 3, 3);
        imagesc(cov(X'));
        colorbar

    end

end

author: canlab


text:
function prob_test084 ( )

%*****************************************************************************80
%
%% TEST084 tests GOMPERTZ_CDF, GOMPERTZ_CDF_INV, GOMPERTZ_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST084\n' );
  fprintf ( 1, '  For the Gompertz PDF:\n' );
  fprintf ( 1, '  GOMPERTZ_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  GOMPERTZ_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  GOMPERTZ_PDF evaluates the PDF;\n' );

  a = 2.0;
  b = 3.0;

  check = gompertz_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST084 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =       %14f\n', a );
  fprintf ( 1, '  PDF parameter B =       %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = gompertz_sample ( a, b, seed );

    pdf = gompertz_pdf ( x, a, b );

    cdf = gompertz_cdf ( x, a, b );

    x2 = gompertz_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function V = construct_V_matrix(measurements)
%function V = construct_V_matrix(measurements)
%
% This function constructs and returns the translational data matrix V
% defined in equation (16) of the paper

% Copyright (C) 2016 by David M. Rosen

D = length(measurements.t{1}); % D = dimension of SE(d)
N = max(max(measurements.edges));  % N = number of nodes in the pose graph
M = size(measurements.edges,1); % M = number of edges in the pose graph

% The number of nonzero elements in V; there are D nonzero elements for
% each translational measurement, plus D*N slots to store the sum of the
% observations emanating from each node

NNZ = D*M + D*N;

rows = zeros(1, NNZ);
cols = zeros(1, NNZ);
vals = zeros(1, NNZ);

for e = 1:M

    %Extract measurement data
    i = measurements.edges(e, 1);  %The node that this edge *leaves*
    j = measurements.edges(e, 2);  %The node that this edge *enters*

    tij = measurements.t{e};
    tau_ij = measurements.tau{e};

    %Set V_ji = -tau_ij * t_ij'

    cmin = D*(e-1) + 1;
    cmax = D*(e-1) + D;

    rows(cmin:cmax) = j*ones(1, D);
    cols(cmin:cmax) = D*(i-1) + 1 : D*(i-1) + D;
    vals(cmin:cmax) = -tau_ij * tij';

    %Add this observation to the weighted sum of measurements emanating
    %from node i

    cmin = D*M + D*(i - 1) + 1;
    cmax = D*M + D*(i - 1) + D;
    vals(cmin : cmax) = vals(cmin : cmax) + tau_ij*tij';
end

% Fill in the indices for the running sums

for i = 1:N
    cmin = D*M + D*(i-1) + 1;
    cmax = D*M + D*(i-1) + D;

    rows(cmin:cmax) = i*ones(1,D);
    cols(cmin:cmax) = [D*(i-1) + 1 : D*(i - 1) + D];
end

V = sparse(rows, cols, vals, N, D*N);

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function facedemo ( n )

%*****************************************************************************80
%
% FACEDEMO: Example polygonal geometries for MESHFACES.
%
%  Author:
%
%    Darren Engwirda
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'The usual call is to FACEDEMO ( N ), with N = 1 or 2.\n' );
    fprintf ( 1, '  Since N was not specified, the value 1 will be chosen.\n' );
    n = 1;
  end

  switch n

    case 1

      node = [0.0, 0.0; 1.0,0.0; 1.0,1.0; 0.0,1.0; 1.01,0.0; 1.01,1.0; 3.0,0.0; 3.0,1.0];
      edge = [1,2; 2,3; 3,4; 4,1; 2,5; 5,6; 6,3; 5,7; 7,8; 8,6];
      face{1} = [1,2,3,4];
      face{2} = [5,6,7,2];
      face{3} = [8,9,10,6];

      meshfaces ( node, edge, face );

    case 2

      % Geometry
      dtheta = pi/36;
      theta = (-pi:dtheta:(pi-dtheta))';
      node1 = [cos(theta), sin(theta)];
      node2 = [-2.0,-2.0; 2.0,-2.0; 2.0,2.0; -2.0, 2.0];
      edge1 = [(1:size(node1,1))',[(2:size(node1,1))'; 1]];
      edge2 = [1,2; 2,3; 3,4; 4,1];

      edge = [edge1; edge2+size(node1,1)];
      node = [node1; node2];

      face{1} = 1:size(edge1,1);
      face{2} = 1:size(edge,1);

      meshfaces(node,edge,face);

    otherwise
      error('Invalid demo. N must be between 1-2');

  end

  return
end


text:
%  Figure 10.43      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_43.m is a script to generate Figure 10.43, the rootlocus for
% altitude control with inner-loop stabilization and altitude
% and altitude derivative feedback

clf;
ftq = [-0.0064    0.0263         0  -32.2000         0;
   -0.0941   -0.6240  761.1400 -196.2000         0;
   -0.0002   -0.0015   -4.4120  -12.4800         0;
         0         0    1.0000         0         0;
         0   -1.0000         0  830.0000         0];

g =[0;
  -32.7000;
   -2.0800;
         0;
         0];

h=[0 0 0 0 1];
hdot=ftq(5,1:5);
hf=hdot + 0.1*h;
j=0;

rlocus(ftq,g,-hf,j)
v=[-8, 14, -8, 8];
axis(v);
grid;
title('Fig. 10.43  Root locus for h and hdot feedback with inner-loop stab.')

author: Sable


text:
clear all;
close all;
clc;
check_yale_faces(4, @ssc_l1, 'ssc_l1');

author: indigits

repo: sparse-plex


text:
function [y Bitrate MSE Stepsize QNoise]=pcm(A,fm,fs,n)
%A=amplitute of cosine signal
%fm=frequency of cosine signal
%fs=sampling frequency
%n= number of bits per sample
%MSE=Mean Squar error, QNoise=Quantization Noise
%Example [y Bitrate MSE Stepsize QNoise]=pcm(2,3,20,3)
%If you have any problem or feedback please contact me @
%%===============================================
% NIKESH BAJAJ
% Asst. Prof., Lovely Professional University, India
% Almameter: Aligarh Muslim University, India
% +919915522564, bajaj.nikkey@gmail.com
%%===============================================

t=0:1/(100*fm):1;
x=A*cos(2*pi*fm*t);
%---Sampling-----
ts=0:1/fs:1;
xs=A*cos(2*pi*fm*ts);
%xs Sampled signal

%--Quantization---
x1=xs+A;
x1=x1/(2*A);
L=(-1+2^n); % Levels
x1=L*x1;
xq=round(x1);
r=xq/L;
r=2*A*r;
r=r-A;
%r quantized signal

%----Encoding---
y=[];
for i=1:length(xq)
    d=dec2bin(xq(i),n);
    y=[y double(d)-48];
end
%Calculations
MSE=sum((xs-r).^2)/length(x);
Bitrate=n*fs;
Stepsize=2*A/L;
QNoise=((Stepsize)^2)/12;

figure(1)
plot(t,x,'linewidth',2)
title('Sampling')
ylabel('Amplitute')
xlabel('Time t(in sec)')
hold on
stem(ts,xs,'r','linewidth',2)
hold off
legend('Original Signal','Sampled Signal');

figure(2)
stem(ts,x1,'linewidth',2)
title('Quantization')
ylabel('Levels L')
hold on
stem(ts,xq,'r','linewidth',2)
plot(ts,xq,'--r')
plot(t,(x+A)*L/(2*A),'--b')
grid
hold off
legend('Sampled Signal','Quantized Signal');

figure(3)
stairs([y y(length(y))],'linewidth',2)
title('Encoding')
ylabel('Binary Signal')
xlabel('bits')
axis([0 length(y) -1 2])
grid

repo: mcbench-benchmarks


text:
function varargout = process_channel_addloc( varargin )
% PROCESS_CHANNEL_ADDLOC: Import a raw file in the database.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015-2023

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % ===== EEG DEFAULTS =====
    strList = {''};
    % Get registered Brainstorm EEG defaults
    bstDefaults = bst_get('EegDefaults');
    % Build a list of strings representing all the defaults
    for iGroup = 1:length(bstDefaults)
        for iDef = 1:length(bstDefaults(iGroup).contents)
            strList{end+1} = [bstDefaults(iGroup).name ': ' bstDefaults(iGroup).contents(iDef).name];
        end
    end

    % ===== PROCESS =====
    % Description the process
    sProcess.Comment     = 'Add EEG positions';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Channel file'};
    sProcess.Index       = 36;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Epilepsy?highlight=%28Add+EEG+positions%29#Access_the_recordings';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw'};
    sProcess.OutputTypes = {'data', 'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Option: File to import
    sProcess.options.channelfile.Comment = 'File to import:';
    sProcess.options.channelfile.Type    = 'filename';
    sProcess.options.channelfile.Value   = {...
        '', ...                                % Filename
        '', ...                                % FileFormat
        'open', ...                            % Dialog type: {open,save}
        'Import channel file', ...             % Window title
        'ImportChannel', ...                   % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                          % Selection mode: {single,multiple}
        'files_and_dirs', ...                  % Selection mode: {files,dirs,files_and_dirs}
        bst_get('FileFilters', 'channel'), ... % Get all the available file formats
        'ChannelIn'};                          % DefaultFormats
    % Option: Default channel files
    sProcess.options.usedefault.Comment = 'Or use default:';
    sProcess.options.usedefault.Type    = 'combobox_label';
    sProcess.options.usedefault.Value   = {'', cat(1, strList, strList)};
    % Fix units
    sProcess.options.fixunits.Comment = 'Fix distance units automatically';
    sProcess.options.fixunits.Type    = 'checkbox';
    sProcess.options.fixunits.Value   = 1;
    % Apply vox2ras transformation
    sProcess.options.vox2ras.Comment    = 'Apply voxel=>subject transformation from the MRI';
    sProcess.options.vox2ras.Type       = 'checkbox';
    sProcess.options.vox2ras.Value      = 1;           % If value is set to 2 programatically, then it also removes the MRI coregistration (see process_import_bids)
    sProcess.options.vox2ras.Controller = 'Vox2ras';
    % File selection options
    SelectOptions = {...
        '', ...                            % Filename
        '', ...                            % FileFormat
        'open', ...                        % Dialog type: {open,save}
        'Reference MRI...', ...            % Window title
        'ImportAnat', ...                  % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                      % Selection mode: {single,multiple}
        'files', ...                       % Selection mode: {files,dirs,files_and_dirs}
        {{'_subjectimage'}, 'Volume from Brainstorm database (*subjectimage*.mat)', 'BST'}, ... % Accept only files from the database
        'MriIn'};                          % DefaultFormats: {ChannelIn,DataIn,DipolesIn,EventsIn,MriIn,NoiseCovIn,ResultsIn,SspIn,SurfaceIn,TimefreqIn
    % Option: MRI file
    sProcess.options.mrifile.Comment = 'Reference MRI:';
    sProcess.options.mrifile.Type    = 'filename';
    sProcess.options.mrifile.Value   = SelectOptions;
    sProcess.options.mrifile.Class   = 'Vox2ras';
    % Option: Fiducials
    sProcess.options.fiducials.Comment = 'Anatomical fiducials';
    sProcess.options.fiducials.Type    = 'label';
    sProcess.options.fiducials.Value   = [];
    sProcess.options.fiducials.Hidden  = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    ChannelMat = [];

    % ===== LOAD SELECTED FILE =====
    % Get filename to import
    ChannelFile = sProcess.options.channelfile.Value{1};
    FileFormat  = sProcess.options.channelfile.Value{2};

    % HANDLING ISSUE #591: https://github.com/brainstorm-tools/brainstorm3/issues/591
    % The list of default caps is changing between versions of Brainstorm, therefore the index of a "combobox" option can't be considered a reliable information.
    % On 6-Jan-2022: The option "usedefault" was changed from type "combobox" to "combobox_label" and the use of previous syntax is now an error.
    % Users with existing scripts will get an error and will be requested to update their scripts.
    if isempty(ChannelFile) && ~ischar(sProcess.options.usedefault.Value{1})
        bst_report('Error', sProcess, [], [...
            'On 6-Jan-2023, the option "usedefault" of process_channel_add loc was changed from type "combobox" to "combobox_label".' 10 ...
            'This parameter was previously an integer, indicating an index in a list that unfortunately changes across versions of Brainstorm.' 10 ...
            'The value now must be a string, which points at a specific default EEG cap with no amibiguity.' 10 10 ...
            'Scripts generated before 30-Jun-2022 and executed with a version of Brainstorm posterior to 30-Jun-2022' 10 ...
            'might have been selecting the wrong EEG cap, and should be fixed and executed again.' 10 10 ...
            'If you get this error, you must edit your processing script (or saved pipeline):' 10 ...
            'Use the pipeline editor to generate a new script to call the process "Add EEG position".' 10 10 ...
            'More information in GitHub issue #591: ' 10 ...
            'https://github.com/brainstorm-tools/brainstorm3/issues/591']);
        return
    end

    % Get other options
    if isfield(sProcess.options, 'fixunits') && isfield(sProcess.options.fixunits, 'Value')
        isFixUnits = sProcess.options.fixunits.Value;
    else
        isFixUnits = 1;
    end
    if isfield(sProcess.options, 'vox2ras') && isfield(sProcess.options.vox2ras, 'Value')
        isApplyVox2ras = sProcess.options.vox2ras.Value;
    else
        isApplyVox2ras = 1;
    end
    % Get reference MRI
    if isfield(sProcess.options, 'mrifile') && isfield(sProcess.options.mrifile, 'Value')
        MriFile = sProcess.options.mrifile.Value{1};
    else
        MriFile = [];
    end
    % Get fiducials
    if isfield(sProcess.options, 'fiducials') && isfield(sProcess.options.fiducials, 'Value') && ~isempty(sProcess.options.fiducials.Value)
        sFid = sProcess.options.fiducials.Value;
        isApplyVox2ras = 0;
    else
        sFid = [];
    end

    % Get channel studies
    [tmp, iChanStudies] = bst_get('ChannelForStudy', [sInputs.iStudy]);
    iChanStudies = unique(iChanStudies);
    % Load file
    if ~isempty(ChannelFile)
        ChannelMat = import_channel(iChanStudies, ChannelFile, FileFormat, [], [], 0, isFixUnits, isApplyVox2ras, MriFile);
        % Apply fiducials
        if ~isempty(sFid) && isfield(sFid, 'NAS') && (length(sFid.NAS)==3) && isfield(sFid, 'LPA') && (length(sFid.LPA)==3) && isfield(sFid, 'RPA') && (length(sFid.RPA)==3) && ~(isequal(sFid.NAS(:), [0;0;0]) && isequal(sFid.LPA(:), [0;0;0]) && isequal(sFid.RPA(:), [0;0;0]))
            ChannelMat.SCS.NAS = sFid.NAS;
            ChannelMat.SCS.LPA = sFid.LPA;
            ChannelMat.SCS.RPA = sFid.RPA;
            ChannelMat = channel_detect_type(ChannelMat, 1);
        end
    end

    % ===== USE DEFAULT =====
    if isempty(ChannelFile)
        % Get registered Brainstorm EEG defaults
        bstDefaults = bst_get('EegDefaults');
        % Get default channel file
        strDef = sProcess.options.usedefault.Value{1};
        % If there is something selected
        if ~isempty(strDef)
            % Format: "group: name"
            cDef   = strtrim(str_split(strDef, ':'));
            % Find the selected group in the list
            iGroup = find(strcmpi(cDef{1}, {bstDefaults.name}));
            % If group is found
            if ~isempty(iGroup)
                % Find the selected default in the list
                iDef = find(strcmpi(cDef{2}, {bstDefaults(iGroup).contents.name}));
                % If default was found
                if ~isempty(iDef)
                    ChannelFile = bstDefaults(iGroup).contents(iDef).fullpath;
                    isMni = strcmpi(bstDefaults(iGroup).name, 'ICBM152');
                end
            end
        end
        % If nothing was selected (file or default)
        if isempty(ChannelFile)
            bst_report('Error', sProcess, [], 'No channel file selected.');
            return
        end
        % Load channel file
        ChannelMat = in_bst_channel(ChannelFile);
    else
        isMni = 0;
    end

    % ===== ADD POSITIONS =====
    % Add channel positions
    if ~isempty(ChannelMat)
        channel_add_loc(iChanStudies, ChannelMat, 0, isMni);
    else
        bst_report('Warning', sProcess, [], 'No channel positions added.');
    end

    % Return input files
    OutputFiles = {sInputs.FileName};
end

author: brainstorm-tools


text:
function [objset,trainList]=readAnnotation(Name_batch,theConf)
MaxObjNum=1000000;

minArea=theConf.data.minArea;

fileID=fopen([theConf.data.catedir,'image_class_labels.txt'],'r');
idClassPair=textscan(fileID,'%d %d');
fclose(fileID);
if((length(Name_batch)==1)&&(Name_batch(1)-'0'>0)&&(Name_batch(1)-'0'<=9))
    logInd=idClassPair{2}==str2num(Name_batch);
    imgIds=idClassPair{1}(logInd);
else
    imgIds=idClassPair{1};
end

fileID=fopen([theConf.data.catedir,'train_test_split.txt'],'r');
idTrainTestPair=textscan(fileID,'%d %d');
fclose(fileID);
trainList=find(idTrainTestPair{2}==1);

fileID=fopen([theConf.data.catedir,'classes.txt'],'r');
batchClassnamePair=textscan(fileID,'%d %s');
fclose(fileID);
%logInd=batchClassnamePair{1}==str2num(Name_batch);
%classname=batchClassnamePair{2}(logInd);

fileID=fopen([theConf.data.catedir,'images.txt'],'r');
idNamePair=textscan(fileID,'%d %s');
fclose(fileID);
[~,idx,~]=intersect(idNamePair{1},imgIds);
imgnames=idNamePair{2}(idx);

fileID=fopen([theConf.data.catedir,'bounding_boxes.txt'],'r');
idBndboxPair=textscan(fileID,'%d %d %d %d %d');
fclose(fileID);
[~,idx,~]=intersect(idBndboxPair{1},imgIds);
x=idBndboxPair{2}(idx);
y=idBndboxPair{3}(idx);
width=idBndboxPair{4}(idx);
height=idBndboxPair{5}(idx);

clear logInd fileID idClassPair idNamePair idBndboxPair

objset(MaxObjNum).folder=[];
objset(MaxObjNum).filename=[];
objset(MaxObjNum).name=[];
objset(MaxObjNum).bndbox=[];
objset(MaxObjNum).ID=[];

j=0;
for i=1:length(imgnames)
    filename=imgnames(i);
    words=strsplit(filename{1},'.');
    classID=str2num(words{1});
    logInd=batchClassnamePair{1}==classID;
    classname=batchClassnamePair{2}(logInd);

    name=classname;
    bndbox.xmin=int2str(x(i));
    bndbox.xmax=int2str(x(i)+width(i));
    bndbox.ymin=int2str(y(i));
    bndbox.ymax=int2str(y(i)+height(i));
    if(~IsAreaValid(bndbox,minArea))
        continue;
    end
    j=j+1;
    objset(j).folder='.';
    objset(j).filename=filename{1};
    objset(j).name=name{1};
    objset(j).bndbox=bndbox;
    objset(j).ID=j;
    if(j>MaxObjNum)
        error('MaxObjNum is too small.');
    end
end
objset=objset(1:j);

%load([theConf.output.dir,Name_batch,'/sampling.mat'],'validObj');
%objset=objset(validObj);
for i=1:length(objset)
    objset(i).ID=i;
end
end

function pd=IsAreaValid(bndbox,minArea)
xmin=str2double(bndbox.xmin);
xmax=str2double(bndbox.xmax);
ymin=str2double(bndbox.ymin);
ymax=str2double(bndbox.ymax);
if((xmax-xmin+1)*(ymax-ymin+1)>=minArea)
    pd=true;
else
    pd=false;
end
end

repo: interpretableCNN


text:
%mmcfindls(lambda,mu,c)
%   This function finds the average system size
%   for an M/M/c queueing system.

function out = mmcfindls(lambda,mu,c)

pc = lambda/mu;
lq = mmcfindlq(lambda,mu,c);

ls = lq + pc;

out = ls;


text:
function overlayplot(c,scale,ord)

   % Private method. See ../plot for details.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % PREP PLOT
   figure('Color','w','Position',[50 50 850 200]);
   box on; hold on;

   % GET TEMP CORRELATION OBJECT
   c1 = c;
   c1 = subset(c1,ord);
   c1 = norm(c1);
   c1 = stack(c1);
   c1 = norm(c1);

   % GET MEAN TRACE AMPLITUDE FOR SCALING BELOW
   maxlist = [];
   for i = 1:length(c1.trig)
      maxlist(end+1) = max(abs( c1.traces(i).data ));
   end;
   normval = mean(maxlist);

   % LOOP THROUGH WAVEFORMS
   tmin =  999999;
   tmax = -999999;
   count = 0;
   for i = 1:length(c1.trig)-1
      count = count + 1;
      d = c1.traces(i).data;
      d = scale * d/normval;
      d = -1 * d; 				% because scale is reversed below
      wstartrel = c1.relativeStartTime(i);
      tr = wstartrel + [ 0:length(d)-1]'/ c1.traces(i).samplerate;
      plot(tr,d,'-','Color',[.4 .4 1],'LineWidth',1);
      % save min and max relative trace times
      if tr(1) < tmin
         tmin = tr(1);
      end;
      if tr(end) > tmax
         tmax = tr(end);
      end;

   end;

   % PLOT THE STACK OF TRACES
   d = c1.traces(end).data;
   d = scale * d/normval;
   d = -1 * d;
   wstartrel = c1.relativeStartTime(c1.ntraces);
   tr = wstartrel + [ 0:length(d)-1]' / c1.traces(end).samplerate;
   plot(tr,d,'k-','LineWidth',2);

   % adjust figure
   axis([tmin tmax -1.1 1.1]);
   set(gca,'YDir','reverse');
   set(gca,'YTick',[]);
   xlabel('Relative Time,(s)','FontSize',8);
   text(tmin,-0.8,[' (' num2str(length(c1.trig)-1) ' traces)'],'HorizontalAlignment','Left');

   %PRINT OUT FIGURE
   set(gcf, 'paperorientation', 'portrait');
   set(gcf, 'paperposition', [.25 5 8 2] );
   %print(gcf, '-depsc2', 'FIG_alignwfm.ps')

end

author: geoscience-community-codes

repo: GISMO


text:
%% Kernel Average Misorientation (KAM)
%
%%
% The kernel average misorientation (KAM) is a measure of local grain
% misorientation that is usually derived from EBSD data. For formaly
% defining the KAM we denote by $o_{i,j}$ the orientations at pixel
% position $(i,j)$ and by $N(i,j)$ the set of all neighboring pixels. Then
% the kernel average misorientation $\mathrm{kam}_{i,j}$ at pixel position
% $(i,j)$ is defined as
%
% $$\mathrm{KAM}_{i,j} = \frac{1}{|N(i,j)|}\sum_{(k,l) \in N(i,j)} \omega(o_{i,j}, o_{k,l}) $$
%
% Here $\lvert N(i,j) \rvert$ denotes the number of all neighboring pixels
% taking into account and $\omega(o_{i,j}, o_{k,l})$ the disorientation
% angle between the orientation $o_{ij}$ in the center and the neighbouring
% orientation $(o_{k,l})$. The specific choice of the set $N(i,j)$ of
% neighboring pixels is crucial for the compution of the KAM. Most commonly
% the following additional constrains are made
%
% * consider neighbors up to order $n$, e.g. $n=1,2,3,\ldots$
% * consider only neighbors belonging to the same grain
% * consider only neighbors with a misorientation angle smaller than a
% threshold angle $\delta$
%
% In the case of sqaure and hexagonal grids the order of neighbors is
% illustrated below

plotSquareNeighbours; nextAxis; plotHexNeighbours

%% A Deformed Ferrite Specimen
%
% Let us demonstrate the computation of the KAM at the example of a
% deformed Ferrite specimen. Lets import the data first and reconstruct the
% grain structure

mtexdata ferrite

[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));
% remove one-three pixel grains
ebsd(grains(grains.grainSize <= 3)) = [];
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));

grains = smooth(grains,5);

plot(ebsd('indexed'),ebsd('indexed').orientations)
hold on
plot(grains.boundary,'lineWidth',1.5)
hold off

%%
% Although MTEX allows the computation of the KAM from arbitrarily sampled
% EBSD maps the algorithms are much faster an memory efficient if the maps
% are measured on regular hexagonal or rectangular grid - as it is standard
% in most applications. The command <EBSD.gridify.html |gridify|> makes
% MTEX aware of such an underlying regular measurement grid.

ebsd = ebsd.gridify;

%%
%
% The kernel average misorientation is computed by the command
% <EBSD.KAM.html |ebsd.KAM|>. As all MTEX commands it return the mean
% disorientation angle in radiant. Hence, dividing by the constant |degree|
% gives the result in degree.

kam = ebsd.KAM / degree;

% lets plot it
plot(ebsd,kam,'micronbar','off')
caxis([0,15])
mtexColorbar
mtexColorMap LaboTeX
hold on
plot(grains.boundary,'lineWidth',1.5)
hold off

%%
% When computed with default parameters in MTEX neighbors up to order 1 are
% considered and no threshold angle $\delta$ is applied. If grains have
% been reconstructed and the property |ebsd.grainId| has been set (as we
% did above) only misorientations within the same grain are considered. As
% a consequence the resulting KAM map is dominated by the orientation
% gradients at the subgrain boundaries.
%
% Specifying a reasonable small theshold angle $\delta=2.5^{\circ}$ the
% subgrain boundaries can be effectively removed from the KAM.

plot(ebsd,ebsd.KAM('threshold',2.5*degree) ./ degree,'micronbar','off')
caxis([0,2])
mtexColorbar
mtexColorMap LaboTeX
hold on
plot(grains.boundary,'lineWidth',1.5)
hold off

%%
% Unfortunately, the remaining KAM becomes very sensitve to measurement
% errors and is often very noisy. The noise can be  reduced by considering
% heigher order neighbors

plot(ebsd,ebsd.KAM('threshold',2.5*degree,'order',5) ./ degree,'micronbar','off')
caxis([0,2])
mtexColorbar
mtexColorMap LaboTeX
hold on
plot(grains.boundary,'lineWidth',1.5)
hold off

%%
% Although this reduced noise it also smoothes away local dislocation
% structures. A much more effective way to reduce the effect of measurement
% errors to the kernel average misorientation is to denoise the EBSD map
% first and compute than the KAM from the first order neighbors.

% chose a denoising filter
F = halfQuadraticFilter;
F.alpha = 0.5;

% denoise the orientation map
ebsdS = smooth(ebsd,F,'fill',grains);

% plot the first order KAM
plot(ebsdS,ebsdS.KAM('threshold',2.5*degree) ./ degree,'micronbar','off')
caxis([0,2])
mtexColorbar
mtexColorMap LaboTeX
hold on
plot(grains.boundary,'lineWidth',1.5)
hold off

%%
% We observe that the KAM is not longer related to subgrain boundaries and
% nicely revalves local dislocation structures of the deformed material.
%
%% Some helper functions
%
% The functions below where only used to generate the neighborhood pictures
% of the first paragraph

function plotSquareNeighbours

N = [4 3 2 3 4;...
  3 2 1 2 3;...
  2 1 0 1 2;...
  3 2 1 2 3;...
  4 3 2 3 4];

cs = crystalSymmetry;
ebsd = EBSDsquare(rotation.nan(5,5),N,0:4,{cs,cs,cs,cs,cs},[10 10]);
plot(ebsd,'EdgeColor','black','micronbar','off','figSize','small')
legend off

text(ebsd,N)

end

function plotHexNeighbours

N = [3 2 2 2 3;...
  2 1 1 2 3;...
  2 1 0 1 2;...
  2 1 1 2 3;...
  3 2 2 2 3;...
  3 3 3 3 4];

cs = crystalSymmetry;
ebsd = EBSDhex(rotation.nan(6,5),N,0:4,{cs,cs,cs,cs,cs},10,1,1);
plot(ebsd,'edgecolor','k','micronbar','off','figSize','small')
legend off
text(ebsd,N)
axis off

end

author: mtex-toolbox

repo: mtex

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


text:
function [ objects ] = floorRectangle2hypothesis( xyzBox, rule, room_points )
%RECTANGLE2HYPOTHESIS Summary of this function goes here
%   Detailed explanation goes here
COUNTER_CLOCKWISE = true;
FITTINGRULES;
uni_room_points = zeros(8,3);
uni_room_points(POINTMAPPING{3},1) = room_points(:,1);
uni_room_points(POINTMAPPING{3},2) = room_points(:,2);
uni_room_points(POINTMAPPING{3},3) = room_points(:,3);

I = find(uni_room_points(1:4,1)<0 & uni_room_points(1:4,2)<0);
ID = rem([I I+1 I+2 I+3]-1,4) + 1;
uni_room_points = uni_room_points([ID ID+4],:);
wall_points = zeros(4,3,6);
wall_points(:,:,1) = uni_room_points([1 4 8 5],:);
wall_points(:,:,2) = uni_room_points([3 2 6 7],:);
wall_points(:,:,3) = uni_room_points([2 1 5 6],:);
wall_points(:,:,4) = uni_room_points([4 3 7 8],:);
wall_points(:,:,5) = uni_room_points([1 2 3 4],:);
wall_points(:,:,6) = uni_room_points([8 7 6 5],:);

LRDU_SURF = zeros(6,4);
LRDU_SURF(1,:) = [3 4 5 6];
LRDU_SURF(2,:) = [4 3 5 6];
LRDU_SURF(3,:) = [2 1 5 6];
LRDU_SURF(4,:) = [1 2 5 6];
LRDU_SURF(5,:) = [1 2 3 4];
LRDU_SURF(6,:) = [1 2 4 3];

normals = [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1];
normal_ids = [1 1 2 2 3 3];
main_angle = [pi/2 -pi/2 -pi 0 0];

objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
object = struct('out_points_w',[],'name','other','type',[],'points',[],'x_w',zeros(7,1));
xyzBox = reshape(xyzBox', 3, 4)';

%% if left two points in lhs wall, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,1)), xyzBox([1 4],:), 0 );
if all(inside)
    LID = LRDU_SURF(rule,1);
    wall_point = wall_points(:,:,LID);
    point5 = LineFaceIntersection(wall_point(1,:), normals(LID,:), [0 0 0], xyzBox(1,:));
    point4 = LineFaceIntersection(wall_point(1,:), normals(LID,:), [0 0 0], xyzBox(4,:));
    point6 = LineFaceIntersection(point5, normals(rule,:), [0 0 0], xyzBox(2,:));
    point1 = LineFaceIntersection(point4, normals(rule,:), [0 0 0], xyzBox(3,:));
    point3 = point4;    point2 = point1;
    point3(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point2(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 13;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point5-point4);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

%% if right two points in rhs wall, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,2)), xyzBox([2 3],:), 0 );
if all(inside)
    RID = LRDU_SURF(rule,2);
    wall_point = wall_points(:,:,RID);
    point6 = LineFaceIntersection(wall_point(1,:), normals(RID,:), [0 0 0], xyzBox(2,:));
    point1 = LineFaceIntersection(wall_point(1,:), normals(RID,:), [0 0 0], xyzBox(3,:));
    point5 = LineFaceIntersection(point6, normals(rule,:), [0 0 0], xyzBox(1,:));
    point4 = LineFaceIntersection(point1, normals(rule,:), [0 0 0], xyzBox(4,:));
    point3 = point4;    point2 = point1;
    point3(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point2(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 13;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point5-point4);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

%% if bottom two points in bottom wall, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,3)), xyzBox([4 3],:), 0 );
if all(inside)
    BID = LRDU_SURF(rule,3);
    wall_point = wall_points(:,:,BID);
    point4 = LineFaceIntersection(wall_point(1,:), normals(BID,:), [0 0 0], xyzBox(4,:));
    point1 = LineFaceIntersection(wall_point(1,:), normals(BID,:), [0 0 0], xyzBox(3,:));
    point5 = LineFaceIntersection(point4, normals(rule,:), [0 0 0], xyzBox(1,:));
    point6 = LineFaceIntersection(point1, normals(rule,:), [0 0 0], xyzBox(2,:));
    point3 = point4;    point2 = point1;
    point3(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point2(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 13;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point5-point4);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

%% if top two points in upper wall, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,4)), xyzBox([1 2],:), 0 );
if all(inside)
    UID = LRDU_SURF(rule,4);
    wall_point = wall_points(:,:,UID);
    point5 = LineFaceIntersection(wall_point(1,:), normals(UID,:), [0 0 0], xyzBox(1,:));
    point6 = LineFaceIntersection(wall_point(1,:), normals(UID,:), [0 0 0], xyzBox(2,:));
    point4 = LineFaceIntersection(point5, normals(rule,:), [0 0 0], xyzBox(4,:));
    point1 = LineFaceIntersection(point6, normals(rule,:), [0 0 0], xyzBox(3,:));
    point3 = point4;    point2 = point1;
    point3(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point2(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 13;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point5-point4);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

end

author: zouchuhang

repo: LayoutNet


text:
function gdout=filterbankdual(g,a,varargin)
%FILTERBANKDUAL  Dual filters
%   Usage:  gd=filterbankdual(g,a,L);
%           gd=filterbankdual(g,a);
%
%
%   `filterbankdual(g,a,L)` computes the canonical dual filters of *g* for a
%   channel subsampling rate of *a* (hop-size) and system length *L*.
%   *L* must be compatible with subsampling rate *a* as
%   `L==filterbanklength(L,a)`. This will create a dual frame valid for
%   signals of length *L*.
%
%   `filterabankrealdual(g,a)` does the same, but the filters must be FIR
%   filters, as the transform length is unspecified. *L* will be set to
%   next suitable length equal or bigger than the longest impulse response
%   such that `L=filterbanklength(gl_longest,a)`.
%
%   The input and output format of the filters *g* are described in the
%   help of |filterbank|.
%
%   In addition, the funtion recognizes a 'forcepainless' flag which
%   forces treating the filterbank *g* and *a* as a painless case
%   filterbank.
%
%   To actually invert the output of a filterbank, use the dual filters
%   together with the |ifilterbank| function.
%
%   REMARK: In general, perfect reconstruction can be obtained for signals
%   of length *L*. In some cases, using dual system calculated for shorter
%   *L* might work but check the reconstruction error.
%
%   See also: filterbank, ufilterbank, ifilterbank

complainif_notenoughargs(nargin,2,'FILTERBANKDUAL');

definput.import={'filterbankdual'};
definput.flags.outformat = {'fir','full','econ','asfreqfilter'};
definput.keyvals.efsuppthr = 10^(-5);

[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);

[g,asan,info]=filterbankwin(g,a,L,'normal');
if isempty(L)
    if info.isfir
        % Pick shortest possible length for FIR filterbank
        L = filterbanklength(info.longestfilter,asan);
    else
        % Just thow an error, nothing reasonable can be done without L
        error(['%s: L must be specified when not working with FIR ',...'
               'filterbanks.'], upper(mfilename));
    end
end
M=info.M;

% Force usage of the painless algorithm
if flags.do_forcepainless
    info.ispainless = 1;
end

% Check user defined L
if L~=filterbanklength(L,a)
     error(['%s: Specified length L is incompatible with the length of ' ...
            'the time shifts.'],upper(mfilename));
end;

% Prioritize painless over uniform algorithm if both are suitable
if info.isuniform && info.ispainless
    info.isuniform = 0;
end

% Factorization of frame operator to block-diagonal matrix
if info.isuniform
  % Uniform filterbank, use polyphase representation
  a=a(1);

  % Transfer functions of individual filters as cols
  G = filterbankfreqz(g,a,L);

  N=L/a;

  gd=zeros(M,N,class(G));

  for w=0:N-1
    idx = mod(w-(0:a-1)*N,L)+1;
    H = G(idx,:);

    H=pinv(H)';

    gd(:,idx)=H.';
  end;
  % gd was created transposed because the indexing gd(:,idx_a)
  % is much faster than gd(idx_a,:)
  gd =  gd.';
  thisclass = class(G);

switch flags.outformat
    case 'fir'
        gd=ifft(gd)*a;
        % Matrix cols to cell elements + cast
        gdout = cellfun(@(gdEl) cast(gdEl,thisclass), num2cell(gd,1),...
                  'UniformOutput',0);

  %      All filters in gdout will be treated as FIR of length L. Convert them
  %      to a struct with .h and .offset format.
        gdout = filterbankwin(gdout,a);
    case 'full'

        gdout = gd*a;
    case 'econ'
        Shorten filters to essential support
        gd = gd*a;
        gdout=economize_filters(gd,'efsuppthr',kv.efsuppthr);

    case 'asfreqfilter'
        gd = gd*a;
%        All filters in gdout will be treated as (numeric) freqfilter format.
%        Manually convert them to a struct with .H and .foff.
        template = struct('H',[],'foff',0,'realonly',0,'delay',0,'L',L);
        gdout = cell(1,M);
        gdout(:) = {template};

        [H,foff,~]=economize_filters(gd,'efsuppthr',kv.efsuppthr);
        for kk = 1:M
            gdout{kk} = setfield(gdout{kk},'H',H{kk});
            gdout{kk} = setfield(gdout{kk},'foff',foff(kk));
        end
    otherwise
        error('%s: Unknown filter format.', upper(mfilename));
end

elseif info.ispainless
   % Factorized frame operator is diagonal.
   gdout = comp_painlessfilterbank(g,asan,L,'dual',0);
else
        error(['%s: The canonical dual frame of this system is not a ' ...
               'filterbank. You must either call an iterative ' ...
               'method to perform the desired inverstion or transform ',...
               'or transform the filterbank to uniform one. Please see ' ...
               'FRANAITER or FRSYNITER for the former and ',...
               'NONU2UFILTERBANK for the latter case.'],upper(mfilename));

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function quadrule_test41 ( )

%*****************************************************************************80
%
%% QUADRULE_TEST41 tests RADAU_SET and SUM_SUB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2005
%
  order_max = 15;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  a = 0.0;
  b = 1.0;

  nsub = 1;
  xlo = -1.0;
  xhi = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST41\n' );
  fprintf ( 1, '  RADAU_SET sets up Radau quadrature;\n' );
  fprintf ( 1, '  SUM_SUB carries it out.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is [%f, %f]\n', a, b );
  fprintf ( 1, '  Number of subintervals is %d\n', nsub );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for n = 1 : order_max

      fprintf ( 1, '  %2d', n );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = radau_set ( n );

        result(i) = sum_sub ( @func, a, b, nsub, n, xlo, xhi, xtab, weight );

        fprintf ( 1, '  %12f', result(i) );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
classdef TestPyrUp
    %TestPyrUp

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_uint8_img
            img = cv.imread(TestPyrUp.im, 'ReduceScale',2);
            result = cv.pyrUp(img);
        end

        function test_double_img
            img = cv.imread(TestPyrUp.im, 'ReduceScale',2);
            img = double(img) ./ 255;
            result = cv.pyrUp(img);
        end

        function test_custom_size
            img = cv.imread(TestPyrUp.im, 'ReduceScale',2);
            img = img(1:200,1:256,:);
            [h,w,~] = size(img);
            sz = [w,h].*2;
            result = cv.pyrUp(img, 'DstSize',sz);
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef RegularizedExponentBestCoeffComputer < handle

   properties (Access = private)
        aValues
        bValues
        rValues
        qValues
        aValue
        bValue
        rValue
        qValue
        nA
        nB
        nR
        nQ
        nT
        error
    end

    properties (Access = private)
       errorComputer
    end

    methods (Access = public)

        function obj = RegularizedExponentBestCoeffComputer(cParams)
            obj.init(cParams);
            obj.initValues();
            obj.computeAxisValues();
        end

        function [a,b,r,q] = compute(obj)
            obj.evaluateAllSamples();
            [a,b,r,q] = obj.obtainBestParameterSample();
        end

        function [a,b,r,q] = obtainBestParameters(obj)
           a = 23.620689655172413;
           b = 1.444444444444444;
           r = 0.677777777777778;
           q = 11.322033898305085;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.errorComputer = cParams.errorComputer;
            obj.nA = 30;
            obj.nB = 10;
            obj.nR = 10;
            obj.nQ = 60;
            obj.nT = obj.nA*obj.nB*obj.nR*obj.nQ;
        end

        function initValues(obj)
            obj.aValue = zeros(obj.nT,1);
            obj.bValue = zeros(obj.nT,1);
            obj.rValue = zeros(obj.nT,1);
            obj.qValue = zeros(obj.nT,1);
            obj.error = zeros(obj.nT,1);
        end

        function computeAxisValues(obj)
            obj.aValues = linspace(22,25,obj.nA);
            obj.bValues = linspace(1.2,1.6,obj.nB);
            obj.rValues = linspace(0.65,0.68,obj.nR);
            obj.qValues = linspace(11,12,obj.nQ);
        end

        function evaluateAllSamples(obj)
            iter = 1;
            for iR = 1:length(obj.rValues)
                for iQ = 1:length(obj.qValues)
                    for iB = 1:length(obj.bValues)
                        for iA = 1:length(obj.aValues)
                            obj.obtainSampleValue(iter,iA,iB,iR,iQ);
                            obj.computeError(iter);
                            iter = iter + 1;
                        end
                    end
                end
            end
        end

        function obtainSampleValue(obj,iT,iA,iB,iR,iQ)
            obj.aValue(iT) = obj.aValues(iA);
            obj.bValue(iT) = obj.bValues(iB);
            obj.rValue(iT) = obj.rValues(iR);
            obj.qValue(iT) = obj.qValues(iQ);
        end

        function  computeError(obj,iT)
           a = obj.aValue(iT);
           b = obj.bValue(iT);
           r = obj.rValue(iT);
           q = obj.qValue(iT);
           e = obj.errorComputer(a,b,r,q);
           obj.error(iT) = e;
        end

        function [a,b,r,q] = obtainBestParameterSample(obj)
            [~,ind] = min(obj.error);
            a = obj.aValue(ind);
            b = obj.bValue(ind);
            r = obj.rValue(ind);
            q = obj.qValue(ind);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% The COBRAToolbox: testGapFind.m
%
% Purpose:
%     - testGapFind tests the validity of gapFind.
%
% Authors:
%     - Original file: Thomas Pfau Sept 2017
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testGapFind'));
cd(fileDir);

testModel = createToyModelForgapFind();
%test solver packages
solverPkgs = {'tomlab_cplex', 'ibm_cplex', 'gurobi', 'glpk'};

for k = 1:length(solverPkgs)

    solverOK  = changeCobraSolver(solverPkgs{k}, 'MILP', 0);

    if solverOK
        fprintf( 'Testing gapFind with solver %s ...\n',solverPkgs{k});
        %The Basic function should find exactly one Metabolite (H) as gapped
        %The Cycle goes undetected as it can carry flux and I is blocked but not
        %gapped.
        [allGaps, rootGaps, downstreamGaps] = gapFind(testModel);

        assert(isempty(downstreamGaps))
        assert(isequal(allGaps,{'H'}));
        assert(isequal(rootGaps,{'H'}));

        %Now, we change this to also find No Consumption gaps, It now detects I
        %alng H and does still not detect anything downstream.
        [allGaps, rootGaps, downstreamGaps] = gapFind(testModel,true);
        assert(isempty(downstreamGaps))
        assert(isequal(allGaps,{'H';'I'}));
        assert(isequal(rootGaps,{'H';'I'}));

        %Finally, we eliminate the cycle. This should yield E as additional Gap
        %(starting point of the disconnected part) and F and G as downstream Gaps.
        testModelWOCycle = removeRxns(testModel,'R6');
        [allGaps, rootGaps, downstreamGaps] = gapFind(testModelWOCycle);
        dsGaps = {'F';'G'};
        rGaps = {'E';'H'};
        aGaps = union(dsGaps,rGaps);
        assert(isequal(downstreamGaps,dsGaps))
        assert(isequal(allGaps,aGaps));
        assert(isequal(rootGaps,rGaps));

        %Only test for one solver (testing multiple solvers is a waste of
        %time)

        break
    end

end
if exist('MILPProblem.mat','file')
    delete('MILPProblem.mat')
end
if exist('MILPProblem.mat','file')
    delete('MILPProblem.mat')
end

author: opencobra

repo: cobratoolbox


text:
function result = is_orthogonal_matrix(P)
%
% Orthogonal Matrices
%
% is_orthogonal_matrix(P) determines if the matrix P is an orthogonal
% matrix. An error is returned if a matrix that is not square is attempted
% to be determined for orthogonality.
%
%  Function written by Anthony Russo, downloaded from MatlabCentral.
%

matrix_size = size(P);

m = matrix_size(1,1);
n = matrix_size(1,2);

tolerance = 10^-10;

if m ~= n
    error('Only square matrices can be orthogonal.');
else
    count = 0;

    identity_matrix = P*P';

    if det(P) ~= 0
        for i = 1:m
            if abs(identity_matrix(i,i) - 1) <= tolerance
                count = count + 1;
            else
                break
            end
        end
    end

    if count == m
        result = 1;
    else
        result = 0;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip