text:
function [handle, xdata_all] = bar_wani_2016(y, e, bar_width, varargin)
% Draw a bar plot with error bars with some additional useful features
% (doesn't work with older than matlab 2015. If you're using matlab 2014
%  or older, please use bar_wani, not this).
%
% :Usage:
% ::
%
%    h = bar_wani(y, e, bar_width, varargin)
%
% :Inputs:
%
%   **y:**
%        y values for bars (row: bar grouping, column: bar values
%        within a bar group) (e.g., if there are m bar groups and
%        n bars for each group, y should be a m x n matrix)
%
%   **e:**
%        error bars (m x n matrix)
%
%   **bar_width:**
%        value for bar width between 0 and 1. This will determine
%        the intervals between bars.
%
% :Optional Inputs: Enter keyword followed by variable with values
%
%   **'ylim':**
%        y axis range, (e.g., 'ylim', [-1 1])
%
%   **'ytick':**
%        y tick values (e.g., 'ytick', -.002:.001:.002)
%
%   **'errbar_width':**
%        the horizontal width of error bars (e.g., 'errbar_width', 0)
%
%   **'colors':**
%        bar colors: each row determines the color for each bar in order (n x 3 matrix)
%
%   **'ast':**
%        put asterisks according to p values, which should be
%        given. (e.g., 'ast', p [m x n]) *p<.05, **p<.01, ***p<.001
%
%   **'btwlines':**
%        this option puts lines between bar groups. This is
%        followed by the line style (e.g., 'btwlines', '--');
%
%   **'dosave':**
%        followed by savename (e.g., 'dosave', savename)
%
% :Some advanced options:
%
%   **'scatter':**
%        show individual data points, which should be in cell array
%
%   **'text':**
%        this will put a number for each bar (e.g., 'text', round(y) [m x n])
%
%   **'ast_adj_y_pos':**
%        When the asterisk locations (on y axis) for bars with positive
%        values are off, you can adjust it using this option
%
%   **'ast_adj_y_neg':**
%        When the asterisk locations (on y axis) for bars with negative
%        values are off, you can adjust it using this option
%
%   **'ast_adj_x':**
%        When the asterisk locations (on x axis) are off, you
%        can adjust it using this option
%
%   **'bar_edgecol':**
%        You can use different colors for bar edges (col [n x 3 matrix])
%
%   **'bar_edgewidth':**
%        You can change linewidths for bar edges
%
% :Output:
%
%   **h:**
%        graphic handles for a bar plot
%
% :Examples: you can see the output in
% http://wagerlab.colorado.edu/wiki/doku.php/help/core/figure_gallery
% :Examples:
% ::
%
%    % data
%    y = [-0.6518   -0.6934   -0.5417   -0.6496   -0.5946   -0.3839
%        1.1511    0.9090    1.1681    1.2892    0.9346    1.1383];
%    e = [0.3226    0.2936    0.3080    0.3203    0.3368    0.3167
%        0.4026    0.4088    0.4012    0.5586    0.3734    0.4257];
%    p = [0.0433    0.0182    0.0785    0.0426    0.0775    0.2255
%        0.0042    0.0262    0.0036    0.0210    0.0123    0.0075];
%
%    col =  [0    0.1157    0.2686
%            0.1157    0.2765    0.4725
%            0.4843    0.1157    0.1078
%            0.3667    0.4765    0.1353
%            0.2765    0.1902    0.3824
%            0.0922    0.4216    0.5118
%            0.7941    0.3235   0];
%
%    % draw
%    bar_wani(y, e, .8, 'colors', col, 'errbar_width', 0, 'ast', p, 'ylim', [-2.5 2.5], 'ytick', -2:2, 'ast_adj_x', 0, 'ast_adj_y_neg', .15);
%    set(gca, 'ytickLabel', num2str(get(gca, 'ytick')'));
%    set(gcf, 'position', [1   531   399   169]);
%
%    savename = 'example_barwani.pdf';
%
%    try
%        pagesetup(gcf);
%        saveas(gcf, savename);
%    catch
%        pagesetup(gcf);
%        saveas(gcf, savename);
%    end
%
% ..
%    Copyright (C) 2014  Wani Woo
%
%    Programmers' notes:
%    'yline'
% ..

dosave = 0;
doman_ylim = 0;
docolor = 0;
doast = 0;
doerrbar_width = 0;
doytick = 0;
dobtwlines = 0;
dotext = 0;
doscatter = 0;
ast_adj_y_pos = .2;
ast_adj_y_neg = .5;
ast_adj_x = 0;
text_adj_y = 0;
bar_edgewidth = 1.8;
doyline = 0;
use_samefig = false;

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % functional commands
            case {'dosave', 'save'}
                dosave = 1;
                savename = varargin{i+1};
            case {'ylim'}
                doman_ylim = 1;
                ylim = varargin{i+1};
            case {'colors', 'color'}
                docolor = 1; colors = varargin{i+1};
            case {'ast'}
                doast = 1; p = varargin{i+1};
            case {'errbar_width'}
                doerrbar_width = 1; errwidth = varargin{i+1};
            case {'ytick'}
                doytick = 1; ytick = varargin{i+1};
            case {'btwlines'}
                dobtwlines = 1; btwlnstyle = varargin{i+1};
            case {'yline'}
                doyline = 1; yln_y = varargin{i+1};
            case {'text'}
                dotext = 1; text_bar = varargin{i+1};
            case {'scatter'}
                doscatter = 1; sc_data = fliplr(varargin{i+1});
            case {'ast_adj_y_pos'}
                ast_adj_y_pos = ast_adj_y_pos + varargin{i+1};
            case {'ast_adj_y_neg'}
                ast_adj_y_neg = ast_adj_y_neg + varargin{i+1};
            case {'ast_adj_x'}
                ast_adj_x = ast_adj_x + varargin{i+1};
            case {'text_adj_y'}
                text_adj_y = varargin{i+1};
            case {'bar_edgecol'}
                bar_edgecol = flipud(varargin{i+1});
            case {'bar_edgewidth'}
                bar_edgewidth = varargin{i+1};
            case {'use_samefig'}
                use_samefig = true;
        end
    end
end

if ~use_samefig
    handle.main = create_figure('Bar plot');
    set(gcf, 'Position', [1   450   300.*size(y,1)   250]);
end

barnum = size(y,2);
grnum = size(y,1);

if ~doman_ylim
    step = (max(max(y+e))-min(min(y+e)))*.2;
    ymax = max(max(y+e)) + step;
    ymin = min(min(y-e)) - step;
else
    ymin = ylim(1);
    ymax = ylim(2);
end

if doyline
    line([0.55 .45+size(y,1)], [yln_y yln_y], 'linewidth', 1, 'linestyle', '--', 'color', [.4 .4 .4]);
end

handle.bar = barweb(y, e, bar_width, [], [], [], [], [], [], [], [], []);
handle.bar = barweb(y, e, bar_width, [], [], [], [], [], [], [], [], []);
set(gcf, 'Color', 'w')
set(gca, 'ylim', [ymin ymax], 'XLim', [0.55 .45+size(y,1)], 'fontsize', 20, 'linewidth', 1.8); % **ADJUST**: adjust basic setting for axis

if doytick
    set(gca, 'ytick', ytick);
end

h = get(gca, 'children');

if docolor
    for i = (barnum+1):2*barnum
        try
            colors_flip = fliplr(colors);
            set(h(i), 'FaceColor', colors_flip{i-barnum});
        catch
            colors_flip = flipud(colors);
            set(h(i), 'FaceColor', colors_flip(i-barnum,:));
        end
    end
end

if doast
    ast = p_asterisk(p);
    for i = 1:barnum
        for j = 1:grnum
            if y(j,i) > 0
                ast_loc(j,i) = y(j,i) + e(j,i) + ast_adj_y_pos * max(max(e));
            else
                ast_loc(j,i) = y(j,i) - e(j,i) - ast_adj_y_neg * max(max(e));
            end
        end
    end
end

for i = 1:barnum
    if ~doerrbar_width
        errwidth = 8;
    end

    handle.bar.errors(i).CapSize = errwidth;
    xdata_all(:,i) = handle.bar.errors(i).XData';

    for j = 1:grnum %numel(handle.bar.errors(i).XData)
        if doast
            text(double(handle.bar.errors(i).XData(j)+ast_adj_x), double(ast_loc(j,i)), ast{j,i}, 'fontsize', 18, 'HorizontalAlignment','center'); % **ADJUST** font size
            if dotext
                text_y = ymax*1.15 + text_adj_y;
                text(handle.bar.errors(i).XData(j), text_y, num2str(text_bar(j,i)), 'fontsize', 20, 'HorizontalAlignment','center'); % **ADJUST** font size
            end
        end

        if doscatter
            hold on; scatter(double((handle.bar.errors(i).XData(j))+.02).*ones(size(sc_data{j,i})), sc_data{j,i}, 50, [.3 .3 .3], 'filled');
            if ymin > min(sc_data{j,i}), ymin = min(sc_data{j,i}); set(gca, 'ylim', [ymin ymax]); end
            if ymax < max(sc_data{j,i}), ymax = max(sc_data{j,i}); set(gca, 'ylim', [ymin ymax]); end
        end
    end

end

j = 1;
if ~exist('bar_edgecol', 'var'), bar_edgecol = repmat([0 0 0], barnum, 1); end
for i = barnum+1:2*barnum, set(h(i), 'lineWidth', bar_edgewidth, 'edgecolor', bar_edgecol(j,:)); j = j+1; end % **ADJUST** bar linewidth and edge color

% btwlines: drawing lines between bar groups
if dobtwlines
    btwx = 1.5:1:grnum;
    btwy = get(gca, 'ylim');
    btwx = repmat(btwx', 1, 2);
    btwy = repmat(btwy, size(btwx,1),1);
    btwlnwidth = 1;
    btwcol = [.3 .3 .3];
    for i = 1:size(btwx,1)
        line(btwx(i,:), btwy(i,:), 'linestyle', btwlnstyle, 'linewidth', btwlnwidth, 'color', btwcol)
    end
end

set(gca, 'tickdir', 'out', 'ticklength', [.01 .01]);

xtickdata = sort(xdata_all(:));

if size(y,1)==1
    set(gca, 'xtick', xtickdata(~isnan(y)));
else
    set(gca, 'xtick', xtickdata(~any(isnan(y(:)),2)));
end

if dosave
    try
        pagesetup(handle);
        saveas(handle, savename);
    catch
        pagesetup(handle);
        saveas(handle, savename);
    end
end

xdata_all = flipud(xdata_all);

end

function ast = p_asterisk(p)
% function ast = p_asterisk(p)
% change p values into asterisk
%

p = double(p < .1) + double(p < .05) + double(p < .01) + double(p < .001);

ast = cell(size(p));

for i = 1:size(p,1)
    for j = 1:size(p,2)
        if p(i,j) > 1
            ast{i,j} = repmat('*', 1, p(i,j)-1);
        elseif p(i,j) == 1
            ast{i,j} = '';
        else
            ast{i,j} = '';
        end
    end
end

end

function handles = barweb(barvalues, errors, width, groupnames, bw_title, bw_xlabel, bw_ylabel, bw_colormap, gridstatus, bw_legend, error_sides, legend_type)

% This function is from http://www.mathworks.com/matlabcentral/fileexchange/10803-barweb--bargraph-with-error-bars-

% Usage: handles = barweb(barvalues, errors, width, groupnames, bw_title, bw_xlabel, bw_ylabel, bw_colormap, gridstatus, bw_legend, error_sides, legend_type)
%
% Ex: handles = barweb(my_barvalues, my_errors, [], [], [], [], [], bone, [], bw_legend, 1, 'axis')
%
% barweb is the m-by-n matrix of barvalues to be plotted.
% barweb calls the MATLAB bar function and plots m groups of n bars using the width and bw_colormap parameters.
% If you want all the bars to be the same color, then set bw_colormap equal to the RBG matrix value ie. (bw_colormap = [1 0 0] for all red bars)
% barweb then calls the MATLAB errorbar function to draw barvalues with error bars of length error.
% groupnames is an m-length cellstr vector of groupnames (i.e. groupnames = {'group 1'; 'group 2'}).  For no groupnames, enter [] or {}
% The errors matrix is of the same form of the barvalues matrix, namely m group of n errors.
% Gridstatus is either 'x','xy', 'y', or 'none' for no grid.
% No legend will be shown if the legend paramter is not provided
% 'error_sides = 2' plots +/- std while 'error_sides = 1' plots just + std
% legend_type = 'axis' produces the legend along the x-axis while legend_type = 'plot' produces the standard legend.  See figure for more details
%
% The following default values are used if parameters are left out or skipped by using [].
% width = 1 (0 < width < 1; widths greater than 1 will produce overlapping bars)
% groupnames = '1', '2', ... number_of_groups
% bw_title, bw_xlabel, bw_ylabel = []
% bw_color_map = jet
% gridstatus = 'none'
% bw_legend = []
% error_sides = 2;
% legend_type = 'plot';
%
% A list of handles are returned so that the user can change the properties of the plot
% handles.ax: handle to current axis
% handles.bars: handle to bar plot
% handles.errors: a vector of handles to the error plots, with each handle corresponding to a column in the error matrix
% handles.legend: handle to legend
%
% See the MATLAB functions bar and errorbar for more information
%
% Author: Bolu Ajiboye
% Created: October 18, 2005 (ver 1.0)
% Updated: Dec 07, 2006 (ver 2.1)
% Updated: July 21, 2008 (ver 2.3)

% Get function arguments
if nargin < 2
	error('Must have at least the first two arguments:  barweb(barvalues, errors, width, groupnames, bw_title, bw_xlabel, bw_ylabel, bw_colormap, gridstatus, bw_legend, barwebtype)');
elseif nargin == 2
	width = 1;
	groupnames = 1:size(barvalues,1);
	bw_title = [];
	bw_xlabel = [];
	bw_ylabel = [];
	bw_colormap = jet;
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 3
	groupnames = 1:size(barvalues,1);
	bw_title = [];
	bw_xlabel = [];
	bw_ylabel = [];
	bw_colormap = jet;
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 4
	bw_title = [];
	bw_xlabel = [];
	bw_ylabel = [];
	bw_colormap = jet;
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 5
	bw_xlabel = [];
	bw_ylabel = [];
	bw_colormap = jet;
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 6
	bw_ylabel = [];
	bw_colormap = jet;
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 7
	bw_colormap = jet;
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 8
	gridstatus = 'none';
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 9
	bw_legend = [];
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 10
	error_sides = 2;
	legend_type = 'plot';
elseif nargin == 11
	legend_type = 'plot';
end

change_axis = 0;
ymax = 0;

if size(barvalues,1) ~= size(errors,1) || size(barvalues,2) ~= size(errors,2)
	error('barvalues and errors matrix must be of same dimension');
else
	if size(barvalues,2) == 1
		barvalues = barvalues';
		errors = errors';
	end
	if size(barvalues,1) == 1
		barvalues = [barvalues; zeros(1,length(barvalues))];
		errors = [errors; zeros(1,size(barvalues,2))];
		change_axis = 1;
	end
	numgroups = size(barvalues, 1); % number of groups
	numbars = size(barvalues, 2); % number of bars in a group
	if isempty(width)
		width = 1;
	end

	% Plot bars
	handles.bars = bar(barvalues, width,'edgecolor','k', 'linewidth', 2);
    h_temp = handles.bars(any(isnan(barvalues)));
    set(h_temp, 'linestyle', 'none')

	hold on
	if ~isempty(bw_colormap)
		colormap(bw_colormap);
	else
		colormap(jet);
	end
	if ~isempty(bw_legend) && ~strcmp(legend_type, 'axis')
		handles.legend = legend(bw_legend, 'location', 'best', 'fontsize',12);
		legend boxoff;
	else
		handles.legend = [];
	end

	% Plot erros
    % 	for i = 1:numbars
    % 		x =get(get(handles.bars(i),'children'), 'xdata');
    % 		x = mean(x([1 3],:));
    % 		handles.errors(i) = errorbar(x, barvalues(:,i), errors(:,i), 'k', 'linestyle', 'none', 'linewidth', 2);
    % 		ymax = max([ymax; barvalues(:,i)+errors(:,i)]);
    %         ymin = min([ymax; barvalues(:,i)+errors(:,i)]); % wani added
    % 	end
    for i = 1:numbars
        if ~verLessThan('matlab', '8.4') % HG2
            x =  handles.bars(i).XData + handles.bars(i).XOffset;
        else
            x =get(get(handles.bars(i),'children'), 'xdata');
            x = mean(x([1 3],:));
        end
        handles.errors(i) = errorbar(x, barvalues(:,i), errors(:,i), 'k', 'linestyle', 'none', 'linewidth', 2);
        ymax = max([ymax; barvalues(:,i)+errors(:,i)]);
        ymin = min([ymax; barvalues(:,i)+errors(:,i)]); % wani added
    end

    if error_sides == 1
        set(gca,'children', flipud(get(gca,'children')));
    end

    ylim([ymin*1.1 ymax*1.1]); % wani changed

	xlim([0.5 numgroups-change_axis+0.5]);

	if strcmp(legend_type, 'axis')
		for i = 1:numbars
			xdata = get(handles.errors(i),'xdata');
			for j = 1:length(xdata)
				text(xdata(j),  -0.03*ymax*1.1, bw_legend(i), 'Rotation', 60, 'fontsize', 12, 'HorizontalAlignment', 'right');
			end
		end
		set(gca,'xaxislocation','top');
	end

	if ~isempty(bw_title)
		title(bw_title, 'fontsize',14);
	end
	if ~isempty(bw_xlabel)
		xlabel(bw_xlabel, 'fontsize',14);
	end
	if ~isempty(bw_ylabel)
		ylabel(bw_ylabel, 'fontsize',14);
	end

	set(gca, 'xticklabel', groupnames, 'box', 'off', 'ticklength', [0 0], 'fontsize', 12, 'xtick',1:numgroups, 'linewidth', 2,'xgrid','off','ygrid','off');
	if ~isempty(gridstatus) && any(gridstatus == 'x')
		set(gca,'xgrid','on');
	end
	if ~isempty(gridstatus) && any(gridstatus ==  'y')
		set(gca,'ygrid','on');
	end

	handles.ax = gca;

	hold off
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% original data folder
hrdir = '/scratch/mfr/DIV2K/DIV2K_train_HR/';
% configuration
p = 128;

% stride
stride = 16;

fnum = 1;
files = dir([hrdir '/' '*.png']);
savedir = 'train';

% 3 channels
h5create([savedir '/hr.h5'],  '/data', [p p 3 Inf], 'Datatype', 'single', 'ChunkSize', [p p 3 1]);

for k = 1:length(files)
    file = files(k).name;
    image = imread([hrdir '/' file]);
    image = im2double(image);

    [w h c] = size(image);
    i = 1;
    j = 1;

    while ((j + p <= h))
        hrpatch = image(i:i+p-1, j:j+p-1, :);

        h5write([savedir '/hr.h5'],  '/data', single(hrpatch),  [1 1 1 fnum], [p p 3 1]);

        fnum = fnum + 1;

        if (mod(fnum, 1000) == 0)
            sprintf('already generated %d patches\n', fnum)
        end

        j = j + stride;
        if (i + p > w)
            i = 1;
            j = j + stride;
        end
    end
end
sprintf('generated %d patches\n', fnum)

author: IVRL

repo: Kernel-Modeling-Super-Resolution


text:
function d=sqdist(a,b)
% SQDIST - computes squared Euclidean distance matrix
%          computes a rectangular matrix of pairwise distances
% between points in A (given in columns) and points in B

% NB: very fast implementation taken from Roland Bunschoten

aa = sum(a.*a,1); bb = sum(b.*b,1); ab = a'*b;
d = abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [MMN] = representativeMolecule(MMN)
% For each moiety, identify a set of representative molecules, based on
% various criteria

[nMoiety,~]=size(MMN.L);

minimalMassMetabolite = cell(nMoiety,1);
fractMinimalMass = zeros(nMoiety,1);
numMinimalMassMetabolites = zeros(nMoiety,1);
for i=1:nMoiety
    bool=(MMN.L(i,:)~=0)';
    minimumMass=min(MMN.mol.Mass(bool));
    bool2 = bool &  MMN.mol.Mass==minimumMass;
    ind = find(bool2);
    %take the first one as a representative minimal Mass
    minimalMassMetabolite{i} = MMN.mol.Mets{ind(1)};
    numMinimalMassMetabolites(i) = length(ind);
    fractMinimalMass(i)=minimumMass/MMN.mol.Mass(ind(1));
end

MMN.moi = addvars(MMN.moi,minimalMassMetabolite,fractMinimalMass,numMinimalMassMetabolites, 'NewVariableNames',{'MinimalMassMol','minimalMassFraction','NumMinimalMassMol'});

author: opencobra

repo: cobratoolbox


text:
function T2 = shrink_obs_dims_in_table(T1, dom, evidence)
% SHRINK_OBS_DIMS_IN_TABLE Set observed dimensions to size 1
% T2 = shrink_obs_dims_in_table(T1, dom, evidence)
%
% If 'T1' contains observed nodes, it will have 0s in the positions that are
% inconsistent with the evidence. We now remove these 0s and set the corresponding dimensions to
% size 1, to be consistent with the way most inference engines handle evidence, which is to
% shrink observed nodes before doing inference.

% This is used by pearl and enumerative inf. engines.

odom = dom(~isemptycell(evidence(dom)));
vals = cat(1,evidence{odom});
ndx = mk_multi_index(length(dom), find_equiv_posns(odom, dom), vals(:));
T2 = T1(ndx{:});

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [w,s]=fram2wav(x,tt,mode)
%FRAM2WAV  converts frame values to a continuous waveform [W]=(X,TT,MODE)
%  Inputs:
%          x(nf,p)      is the input signal: one row per frame
%	       tt(nf,3)     specifies the frames. Each row has the form [start_sample end_sample flag]
%                       where flag = 1 for the start of a new spurt.
%                       If tt(:,3) is omitted, a new spurt will be started whenever there is a gap
%                       of more than one between the end of one frame and the beginning or the next.
%                       A new spurt is automatically started if x() = NaN.
%          mode         consists of one or more of the following letters:
%                          z for zero-order hold interpolation (i.e. constant within each frame)
%                          l for linear interpolation within each spurt [default]
%
% Outputs:
%          w(n,p)       contains the interpolated waveforms. Their length is n = tt(nf,2)
%          s(ns,2)      gives the starting and ending sample numbers of each spurt (excluding NaN spurts)
%
%    This routine converts frame-based values to continuous waveforms by performing
%    a chosen method of interpolation. Interpolation is restarted at the beginning of each spurt.

%    Bugs/Suggestions
%      (1)   Additional mode option for cubic interpolation
%      (2)   Additional mode option for interpolation in log domain
%      (3)   Additional mode option for x values being
%            frame averages rather than mid-frame values.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: fram2wav.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
    mode='l';
end
[nf,m]=size(x);
n=round(tt(end,2));
w=repmat(NaN,n,m);
nt=size(tt,2);
ix1=ceil(tt(:,1)); % start of frame sample
ix2=floor(tt(:,2)); % end of frame sample

% determine the start and end of spurts

if nt>2
    ty=tt(:,3)>0;   % frame type set by user
else
    ty=zeros(nf,1);
    ty(2:end)=ix1(2:end)>ix2(1:end-1)+1;    % new spurt whenever a gap
end
ty(1)=1;           % first frame always starts a spurt
ty(isnan(x))=1;    % NaN always ends previous spurt
ty(1+find(isnan(x(1:end-1))))=1; % NaN always forces a new spurt
ty=double(ty);
ty(1:end-1)=ty(1:end-1)+2*ty(2:end);
ty(end)=ty(end)+2;   % last frame always ends a spurtw=repmat(NaN,n,m);  % initialize output to all NaN
nx=ix2-ix1+1;

if any(mode=='z')   % zero-order hold
    for i=1:nf
        if nx(i)
            w(ix1(i):ix2(i),:)=repmat(x(i,:),nx(i),1);
        end
    end
else   % linear interpolation is the default
    ttm=(tt(:,1)+tt(:,2))/2;    % mid point of frame
    ixm=floor(ttm); % end of first half of frame
    for i=1:nf
        if i==176
            i
        end
        if nx(i)
            tyi=ty(i);
            if tyi==3    % use a zero order hold
                w(ix1(i):ix2(i),:)=repmat(x(i,:),nx(i),1);
            else
                nxm=ixm(i)-ix1(i)+1;
                if nxm
                    if tyi==1
                        grad=(x(i+1,:)-x(i,:))/(ttm(i+1)-ttm(i));
                    else
                        grad=(x(i,:)-x(i-1,:))/(ttm(i)-ttm(i-1));
                    end
                    w(ix1(i):ixm(i),:)=repmat(x(i,:),nxm,1)+((ix1(i):ixm(i))'-ttm(i))*grad;
                end
                if nx(i)>nxm
                    if tyi==2
                        grad=(x(i,:)-x(i-1,:))/(ttm(i)-ttm(i-1));
                    else
                        grad=(x(i+1,:)-x(i,:))/(ttm(i+1)-ttm(i));
                    end
                    w(ixm(i)+1:ix2(i),:)=repmat(x(i,:),ix2(i)-ixm(i),1)+((ixm(i)+1:ix2(i))'-ttm(i))*grad;
                end
            end
        end
    end
end

% now sort out the start and end spurt positions

ty(isnan(x))=0;    % Don't count NaN spurts
s=repmat(ix1(bitand(ty,1)>0),1,2);
s(:,2)=ix2(bitand(ty,2)>0);
if ~nargout
    tw=(1:n)';
    for i=size(s,1):-1:2
        j=s(i,1);   % start of new spurt
        tw=[tw(1:j-1); tw(j); tw(j:end)];
        w=[w(1:j-1); NaN; w(j:end)];        % insert a NaN to force a plotting break
    end
    plot(tt(:,1:2)',repmat(x(:)',2,1),'r-+',tw,w,'b-');
end

author: decouples


text:
function lbls = transformLabels(labels)

% TRANSFORMLABELS A small utility which transform the label vector: if labels are in 1-of-K encoding it transforms it into a vector of
% integers and vice versa.
% SEEALSO: utils_transformMultiLabel
%
% COPYRIGHT: Andreas C. Damianou, 2012

% VARGPLVM

if size(labels,1) == 1 || size(labels,2) ==1
    N = max(size(labels,1), size(labels,2));
    % The labels are not 1-K-encoding
    uniqueLabels = unique(labels);
    numClasses = length(uniqueLabels);
    lbls = zeros(N, numClasses);
    for c = 1:numClasses
        lbls(labels == uniqueLabels(c),c) = 1;
    end
else
    multilabel = false;
    N = size(labels,1);
    % This is in case we have e.g. 1 -1 -1 instead of 1 0 0
    labels(labels==-1)=0;
   % numClasses = size(labels,2);
    lbls = zeros(1,N);
    for n=1:N
        indx = find(labels(n,:));
        if length(indx) == 1
            lbls(1,n) = indx;
        else
            % WARNING!!!! When point n belongs to more than one classes,
            % this function will RANDOMLY select one!!!! This should be
            % fixed in future releases!
            permIndx = randperm(length(indx));
            lbls(1,n) = permIndx(1);
            multilabel = true;
        end
    end
    if multilabel
        warning('This is a multilabel dataset, transformLabels just kept one label per point randomly!')
    end
end

author: SheffieldML

repo: vargplvm


text:
function [outBoxes idsGood]= FilterBoxesWidth(inBoxes, minLen)
% [outBoxes idsGood]= FilterBoxesWidth(inBoxes, minLen)
%
% Filters out small boxes. Boxes have to have a width and height
% larger than minLen
%
% inBoxes:       M x 4 array of boxes
% minLen:   Minimum width and height of boxes
%
% outBoxes:      N x 4 array of boxes, N < M
% idsGood:       M x 1 logical array denoting boxes kept
%
%     Jasper Uijlings - 2013

[nr nc] = BoxSize(inBoxes);

idsGood = (nr >= minLen) & (nc >= minLen);
outBoxes = inBoxes(idsGood,:);

author: zouchuhang

repo: LayoutNet


text:
function cell2txtfile(varargin)

% function cell2txtfile(fileName,T,skipOpt,checkChar)
% ------------------------------------------------------------------------
%
% This function exports the content in the cell array T to
% the text file fileName. Each entry in the cell array will be a line in
% the txt file. Prior to text file creation the cell is converted to a
% column format. If the input skipOpt=1 cell entries which appear empty
% (after spaces are removed) will be skipped.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Change log:
% 2016/09/09: Updated for GIBBON
% 2016/09/09: Added conversion to char and associated warning for non-character content.
% 2018/02/14: Added whole cell conversion when skipOp==0
% 2018/02/14: Added varargin and defaults
%
% To do:
% * Proper handling of conversion of non-character entries (e.g. have user
% specify conversion format to use for numerical data).
%------------------------------------------------------------------------

%% parse input

switch nargin
    case 2
        fileName=varargin{1};
        T=varargin{2};
        skipOpt=0;
        checkChar=1;
    case 3
        fileName=varargin{1};
        T=varargin{2};
        skipOpt=varargin{3};
        checkChar=1;
    case 4
        fileName=varargin{1};
        T=varargin{2};
        skipOpt=varargin{3};
        checkChar=varargin{4};
end

% Force column if it isn't
if ~isvector(T)
    T=T(:);%Make column
end

% Convert non-character content to text
if checkChar==1
    indNotChar=find(~cellfun(@(x) ischar(x),T,'UniformOutput',1));
    if ~isempty(indNotChar)
        %warning('Non character entries detected');
        for q=indNotChar
            T{q}=vec2strIntDouble(T{q},'%6.7e');
        end
    end
end

%%
% Write to file

fid=fopen(fileName,'w');
switch skipOpt
    case 0 %Process in one go
        fprintf(fid,'%s\r\n', T{:});
    case 1 %Looping to check if lines are empty
        for q=1:numel(T)
            l=T{q}; %Get current cell entry
            if ~ischar(l) %If this isn't a char then attempt conversion
                l=sprintf('%u',l);
            end

            %Write if not empty
            if ~isempty(deblank(l))
                fprintf(fid,'%s\r\n',l);
            end
        end
end
fclose(fid);

%%
% _*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 [fn, FnVar, FnGrad, FnGradCov, constraint, ConstraintCov, ConstraintGrad, ConstraintGradCov] = Contamination(x, runlength, seed, ~)
% function [fn, FnVar, FnGrad, FnGradCov, constraint, ConstraintCov, ConstraintGrad, ConstraintGradCov] = Contamination(x, runlength, seed, other);
% x is a vector containing binary var for yes/no to prevention efforts done
% at the stage
% runlength is the number of independent generations of simulated time
% seed is the index of the substreams to use (integer >= 1)
% other is not used
% Returns cost of prevention treatment, constraint, and ConstraintCov
% If contraints not satisfied,
% prints comparison between probabilities that the contamination fractions are within
% the thresholds versus the probabilities (1-epsilon) that need to be exceeded for each stage.
%
%Note: RandStream.setGlobalStream(stream) can only be used for Matlab
%versions 2011 and later
%For earlier versions, use the method RandStream.setDefaultStream(stream)
%
%   *************************************************************
%   ***            Code written by Danielle Lertola           ***
%   ***          dcl96@cornell.edu    June 25th, 2012         ***
%   ***     Note for future update:  Example parameters need  ***
%   ***                              to be more reasonable.   ***
%   ***            Edited by Jennifer Shih                    ***
%   ***          jls493@cornell.edu    June 18th, 2014        ***
%   *************************************************************
%
% Last updated Jun 18, 2014

FnVar=NaN;
FnGrad = NaN;
FnGradCov = NaN;
ConstraintGrad = NaN;
ConstraintGradCov = NaN;

n=length(x);
if (length(x)~=n) ||(sum(x>ones(n,1))>0) || (sum(x<zeros(n,1))>0) || (runlength <= 0) || (seed <= 0) || (round(seed) ~= seed),
    fprintf('\nx has %u elements, elements of x are binary, \nrunlength should be positive and real, seed should be a positive integer.\n',n);
    fn = NaN;
    constraint = NaN;
    ConstraintCov = NaN;
else % main simulation
    %% *********************PARAMETERS*********************
    nGen=runlength;          %number of independent generations
    u=x;                     %prevention binary decision variable
    X=zeros(n,nGen);         %fraction contaminated at each stage for each generation
    epsilon=.05*ones(n,1);   %error probability
    p=.1*ones(n,1);          %proportion limit
    cost=ones(n,1);          %cost for prevention at stage i
    %Beta parameters for initial contamination, contamination rate,
    %restoration rate
    initialAlpha=1;
    initialBeta=30;
    contamAlpha=1;
    contamBeta=17/3;
    restoreAlpha=1;
    restoreBeta=3/7;

    %% GENERATE RANDOM NUMBER STREAMS
    % Generate new streams for
    [InitialStream, ContaminationStream, RestorationStream] = RandStream.create('mrg32k3a', 'NumStreams', 3);
    % Set the substream to the "seed"
    InitialStream.Substream = seed;
    ContaminationStream.Substream = seed;
    RestorationStream.Substream = seed;

    %% Generate initial fraction of contamination
    OldStream = RandStream.setGlobalStream(InitialStream); % Temporarily store old stream, for versions 2011 and later
    %OldStream = RandStream.setDefaultStream(InitialStream);%for versions 2010 and earlier
    % Generate initial fraction of contamination for stage 1 for each
    % generation
    initialX=betarnd(initialAlpha,initialBeta,1,nGen);

    %% Generate rates of contamination
    RandStream.setGlobalStream(ContaminationStream); %for Matlab versions 2011 and later
    %RandStream.setDefaultStream(ContaminationStream); % for versions 2010 and earlier

    % Generate rates of contamination for each stage and generation
    Lambda=betarnd(contamAlpha,contamBeta,n,nGen);

    %% Generate rates of restoration
    RandStream.setGlobalStream(RestorationStream); %for Matlab versions 2011 and later
    %RandStream.setDefaultStream(RestorationStream); %for Matlab versions 2010 and earlier
    % Generate rates of restoration for each stage and generation
    Gamma=betarnd(restoreAlpha,restoreBeta,n,nGen);

    RandStream.setGlobalStream(OldStream);                   % Restore old random number stream
    %RandStream.setDefaultStream(OldStream); %for versions 2010 and earlier
    %% Determinating fraction of contamination at each stage
    X(1,:)=Lambda(1,:)*(1-u(1)).*(1-initialX) + (1-Gamma(1,:)*u(1)).*initialX;
    for i= 2:n
        X(i,:)=Lambda(i,:)*(1-u(i)).*(1-X(i-1,:)) + (1-Gamma(i,:)*u(i)).*X(i-1,:);
    end

    %mu=mean(X,2);
    %sigma=std(X,0,2);
    limit=1-epsilon;
    %prob=normcdf(p,mu,sigma);
    %cost of contamination control
    fn=sum(cost.*u);
    %if sum(limit>=prob)==0,
    %    fprintf('\nGiven starting solution all contamination fractions are less than the threshold with probability > (1-epsilon)');
    %    fprintf('\nSuccessful with cost %4.2f.\n',sum(cost.*u));
    %else
    %    string=['\nGiven starting solution the probability that the contamination fractions \n'...
    %            'are less than the threshold is <= (1-epsilon) for at least one stage.\n' ...
    %            'Therefore the constraints are not satisfied.\n\n'...
    %            'Column 1 contains the probability that the contamination fractions are less than the threshold.\n' ...
    %            'Column 1 must be greater than column 2 (1-epsilon) for each stage in order for the constraint to be satsified.\n'];
    %    fprintf(string);
    %    results=[prob limit]
    %end
    %checking probability that Xi is <=pi is less than 1-eps
    con=zeros(n,runlength); %matrix of if Xi<pi for each trial and i
    for j=1:runlength
        con(:,j)= (X(:,j)<=p);
    end
    con=con';
    le=sum(sum(con,2)==n);
    constraint=zeros(1,n);
    for k=1:n
        constraint(k)=(sum(con(:,k))/runlength)-limit(k);
    end
    ConstraintCov=cov(con);

end
end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function l4mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% L4MAT_PRINT prints an L4MAT.
%
%  Discussion:
%
%    An L4MAT is an array of L4 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, logical A(M,N), the matrix.
%
%    Input, string TITLE, a title.
%
  l4mat_print_some ( m, n, a, 1, 1, m, n, title );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [k] = getObjCutMovie(inputMovie,inputImages,varargin)
	% Creates a movie cell array cut to a region around input cell images.
	% Biafra Ahanonu
	% started: 2014.01.03 [19:13:01]
	% inputs
		% inputMovie - [x y frames]
		% inputImages - [x y nSignals] - NOTE the dimensions, permute(inputImages,[3 1 2]) if you use [x y nSignals] convention
	% outputs
		%

	[k] = ciapkg.video.getObjCutMovie(inputMovie,inputImages,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function out = cfg_example_run_cumsum2(job)
% Example function that returns the cumulative sum of an vector given in
% job.a in out. The output is referenced as out(:), this is defined in
% cfg_example_vout_cumsum1.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

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

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

% The harvested cfg_repeat will return a cell array of numbers. These need
% to be cat'ed together to be useful as input to cumsum.
out.cs = cumsum(cat(1,job.a{:}));

author: spm


text:
function glscf = CO_glscf(y,alpha,beta,tau)
% CO_glscf  The generalized linear self-correlation function of a time series.
%
% This function was introduced in Queiros and Moyano in Physica A, Vol. 383, pp.
% 10--15 (2007) in the paper "Yet on statistical properties of traded volume:
% Correlation and mutual information at different value magnitudes"
% https://www.sciencedirect.com/science/article/pii/S0378437107004645
%
% The function considers magnitude correlations.
%
%---INPUTS:
% y, the input time series
% alpha and beta are real and nonzero parameters
% tau is the time-delay (can also be 'tau' to set to first zero-crossing of the ACF)
%
% When alpha = beta estimates how values of the same order of magnitude are
% related in time
% When alpha ~= beta, estimates correlations between different magnitudes of the
% time series.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check inputs and set defaults
% ------------------------------------------------------------------------------
if nargin < 4 || isempty(tau)
    tau = 'tau';
end

% Set tau to first zero-crossing of the autocorrelation function with the input 'tau'
if strcmp(tau,'tau')
    tau = CO_FirstCrossing(y,'ac',0,'discrete');
end

% Take magnitudes of time-delayed versions of the time series:
y1 = abs(y(1:end-tau));
y2 = abs(y(1+tau:end));

glscf = (mean((y1.^alpha).*(y2.^beta)) - mean(y1.^alpha)*mean(y2.^beta)) / ...
     		    (sqrt(mean(y1.^(2*alpha)) - mean(y1.^alpha)^2) ...
     		          * sqrt(mean(y2.^(2*beta)) - mean(y2.^beta)^2));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function CloseShapefile(fid,format)
%CLOSESHAPEFILE Writes footer to specified shapefile format
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////
%
% Function writes footer to Socet GRFX .grfx file, RemoteView .rvl file or
% Google Earth .kml file
%
%CLASSIFICATION:
%   Unclassified
%
%INPUTS:
%  fid      - required: File Handle of open Shape File
%  format   - required: 'grfx','rvl' or 'kml'
%
%OUTPUTS:
%
%VERSION:
%   1.0
%     - Tim Cox 20090324
%     - initial version

format = lower(format);

if (strcmp(format,'grfx') == 1)
    fprintf(fid,'</GrGmGraphicList>\n');
	fprintf(fid,'</GmLayer>\n');
	fprintf(fid,'</GmLayerList>\n');
	fprintf(fid,'</GmGraphicsLayout>\n');
	fprintf(fid,'</Graphics>\n');
    fclose(fid);
    return;
elseif (strcmp(format,'rvl') == 1)
    fprintf(fid,']\n');
	fprintf(fid,'>>\n');
	fprintf(fid,'>>\n');
    fclose(fid);
    return;
elseif (strcmp(format,'kml') == 1)
    fprintf(fid,'</Document>\n');
	fprintf(fid,'</kml>\n');
    fclose(fid);
    return;
else
    return;
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% cplxdual2D_plots
% DISPLAY 2D WAVELETS OF cplxdual2D.M

J = 4;
L = 3*2^(J+1);
N = L/2^J;
[Faf, Fsf] = FSfarras;
[af, sf] = dualfilt1;
x = zeros(2*L,6*L);
w = cplxdual2D(x, J, Faf, af);
w{J}{1}{2}{2}(N/2,N/2+0*N) = 1;
w{J}{1}{1}{3}(N/2,N/2+1*N) = 1;
w{J}{1}{2}{1}(N/2,N/2+2*N) = 1;
w{J}{1}{1}{1}(N/2,N/2+3*N) = 1;
w{J}{1}{2}{3}(N/2,N/2+4*N) = 1;
w{J}{1}{1}{2}(N/2,N/2+5*N) = 1;
w{J}{2}{2}{2}(N/2+N,N/2+0*N) = 1;
w{J}{2}{1}{3}(N/2+N,N/2+1*N) = 1;
w{J}{2}{2}{1}(N/2+N,N/2+2*N) = 1;
w{J}{2}{1}{1}(N/2+N,N/2+3*N) = 1;
w{J}{2}{2}{3}(N/2+N,N/2+4*N) = 1;
w{J}{2}{1}{2}(N/2+N,N/2+5*N) = 1;
y = icplxdual2D(w, J, Fsf, sf);
y = [y; sqrt(y(1:L,:).^2+y(L+[1:L],:).^2)];
figure(1)
clf
imagesc(y);
title('2D Dual-Tree Complex Wavelets')
axis image
axis off
colormap(gray(128))
print -djpeg95 cplxdual2D_plots

author: thomaskuestner


text:
%{
Test the trace-constrained least-squares problem with positive
semi-definite matrix variables

     minimize (1/2)*norm( A * X - b )^2
     with the constraint that X is positive semi-definite ( X >= 0 )
and optionally,
    the constraint trace(X) <= lambda

%}

% Try to load the problem from disk
fileName = fullfile('reference_solutions','traceLS_problem2_noisy');
if exist([fileName,'.mat'],'file')
    load(fileName);
    fprintf('Loaded problem from %s\n', fileName );
else
    % Generate a new problem
    randn('state',sum('Trace')+5);
    rand('state',sum('Trace2')+5);

    N = 30; R = 2;

    df = R*N;
    oversample = 5;
    Left  = randn(M,R);
    Right = Left;
    k = round(oversample*df);
    k = min( k, round(.8*N*N) );
    omega = randperm(N*N);
    omega = sort(omega(1:k)).';

    X_original = Left*Right';       % the "original" signal -- may not be optimal value
    b_original = X_original(omega);
    EPS = .001;        % noise level
    noise = EPS * randn(k,1);  % this destroys symmetry...
    b = b_original + noise;
    lambda       = trace( X_original );
    objective    = @(X) sum_square( X(omega) - b )/2;
    obj_original = objective(X_original);

    % get references via CVX
    cvx_begin
        cvx_precision best
        cvx_quiet true
        variable Xcvx(N,N)
        minimize objective(Xcvx)
        subject to
            Xcvx == semidefinite(N)
    cvx_end
    X_reference_noTraceConstraint = Xcvx;
        fprintf('Difference between convex solution (no trace constraint) and original signal: %.2e\n', ...
        norm( Xcvx - X_original, 'fro' ) );

    cvx_begin
        cvx_precision best
        cvx_quiet true
        variable Xcvx(N,N)
        minimize objective(Xcvx)
        subject to
            Xcvx == semidefinite(N)
            trace(Xcvx) <= lambda
    cvx_end
    X_reference = Xcvx;         % the trace norm minimizer
    obj_reference = objective(X_reference);
    fprintf('Difference between convex solution (with trace constraint) and original signal: %.2e\n', ...
        norm( Xcvx - X_original, 'fro' ) );
    save(fileName,'X_original','X_reference','X_reference_noTraceConstraint',...
        'omega','b','obj_original',...
        'Left','EPS','b_original','R','obj_reference','lambda');
    fprintf('Saved data to file %s\n', fileName);

end

[M,N]            = size(X_reference);
norm_X_reference = norm(X_reference,'fro');
er_reference     = @(x) norm(x-X_reference,'fro')/norm_X_reference;
norm_X_reference2= norm(X_reference_noTraceConstraint,'fro');
er_reference2    = @(x) norm(x-X_reference_noTraceConstraint,'fro')/norm_X_reference2;
objective        = @(X) norm( X(omega) - b )^2/2;

[omegaI,omegaJ] = ind2sub([M,N],omega);
mat = @(x) reshape(x,M,N);
vec = @(x) x(:);

k  = length(omega);
p  = k/(M*N);
df = R*(M+N-R);
fprintf('%d x %d rank %d matrix, observe %d = %.1f x df = %.1f%% entries\n',...
    M,N,R,k,k/df,p*100);
fprintf(' Trace norm solution and original matrix differ by %.2e\n',...
    norm(X_reference-X_original,'fro')/norm_X_reference );
%% Solve unconstrained version. No smoothing is necessary!
opts = struct('maxIts',500);
opts.errFcn{1} = @(f,x) er_reference2(x); % no trace constraint
opts.errFcn{2} = @(f,x) sum( abs(eig(x)) > 1e-5 ); % numerical rank
% tell it to use eigs instead of eig
% opts.largescale = true; % ( but not recommended)

% opts.symmetrize = true; % another option

A = sparse( omegaI,omegaJ,b,N,N);
[x,out] = solver_psdComp( A, opts );
% Note: we do not expect to get zero error because there might be more
%   than one optimal solution for this problem (since there are not
%   so many constraints)
% Check that we have a feasible solution
fprintf('Objective is %.2e, min eigenvalue is %.2e\n', objective(x),...
    min(eig(x)) );
%% Solve trace constrained version. No smoothing necessary!
opts = struct('maxIts',1500,'tol',1e-9);
opts.errFcn{1} = @(f,x) er_reference(x); % no trace constraint
opts.errFcn{2} = @(f,x) trace(x)-lambda;
opts.errFcn{3} = @(f,x) sum( abs(eig(x)) > 1e-5 ); % numerical rank
% we can also symmetrize "omega", but it makes little difference,
%   and gives slightly differen value than CVX, since CVX symmetrizes
%   differently:
% opts.symmetrize = true;
A = sparse( omegaI,omegaJ,b,N,N);
[x,out] = solver_psdCompConstrainedTrace( A,lambda, opts );

h=figure();
semilogy(out.err(:,1));

% Check that we are within allowable bounds
if out.err(end,1) < 1e-5
    disp('Everything is working');
else
    error('Failed the test');
end

%%
close(h)
% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
clear;clc;
util;
addpath(genpath('jsonlab-master'));
dataset = {};

%% pre-defined paramters
walllabel = 'wall';
floorlabel = 'floor';

%% load the scene list
scenedir = dir(house_folder);
for i = 3:length(scenedir)
    scenelist{i-2} = scenedir(i).name;
end

%% binary label vectors
labelset = load(labelsetname); % 'entirelabelset.mat' is a list of all fine-grained labels
labelset = labelset.labelset;

%% process each scene(house)
for i = 1:length(scenelist)
    scenename = scenelist{i};
    filename = [house_folder, filesep, scenename];
    filename = [filename, filesep, 'house.json'];
    % load the scene as a list of rooms
    try
        roomlist = parsehouse(filename);
        [ roomlist ] = query_obj_metadata( roomlist, scenename ); % fill in the object info (label, orientation, etc)
    catch
        fprintf([scenename,': ERROR in parsing house.json!\n']);
        continue;
    end

    for j = 1:length(roomlist)
        room = roomlist{j};
        roomname = room.roomname;
        try
            % only preserve the rooms meet our requirements
            [ iscontinue, floorobb, wallobblist, wallnum ] = selectroom( room, wcf_folder, roomtype );
            if(iscontinue)
                % transform the objects into [obb+label] format vector
                % and remove the objects with labels not in labelset
                room = convert_room_rep(room,labelset);

                % add wall objects to the room
                wallidx = strmatch(walllabel,labelset,'exact');
                wallonehot = zeros(length(labelset),1);
                wallonehot(wallidx) = 1;
                for v = 1:size(wallobblist,2)
                    room.labellist{length(room.labellist)+1} = walllabel;
                    room.obblist(1:size(room.obblist,1),size(room.obblist,2)+1) = [wallobblist(:,v);wallonehot];
                end

                % add floor to the room
                room.labellist{length(room.labellist)+1} = floorlabel;
                flooridx = strmatch(floorlabel,labelset,'exact');
                flooronehot = zeros(length(labelset),1);
                flooronehot(flooridx) = 1;
                room.obblist(1:size(room.obblist,1),size(room.obblist,2)+1) = [floorobb(:);flooronehot];

                % align the rooms, the floor is located at [0,0,0]
                room = alignroom( room );

                % extract the object-object relations, and relative
                % distance between objects
                [ Rmat, Dmat, sur ] = build_relationgraph( room.obblist', RELATIONS, length(labelset) );

                % add room to the dataset
                room.Rmat = Rmat;
                room.Dmat = Dmat;
                room.sur = sur;
                dataset{length(dataset)+1} = room;

                fprintf([num2str(length(dataset)), ':', scenename,', ',roomname, '\n']);
            end
        catch
            fprintf([scenename,',', roomname,': ERROR in parsing organizing room! \n']);
            continue;
        end
    end

end
save(output_filename, 'dataset');

repo: GRAINS


text:
function mr = mrReadMag(pth, t, slices);
% Read a P*.7.mag file from Gary's Recon code into
% an mr struct.
%
% Usage: mr = mrReadMag(pth, [time points], [slices]);
%
% ras, 07/05
if notDefined('slices'),	slices = [];			end
if notDefined('t'),			t = [];					end

% init mr struct
[p f ext] = fileparts(pth);
mr = mrCreateEmpty;
mr.path = pth;
mr.name = f;
mr.format = 'pmag';

% first read header file, determine if little endian
% format needs to be used (for Lucas 3T scanner @ Stanford):
[func, mr.hdr] = mrReadMagHeader(pth);
mr.info = mr.hdr;
mr.info.scanner = '';
mr.info.subject = mr.hdr.name;
mr.info.subjectSex = '';
mr.info.subjectAge = [];
[mr.info.date, mr.info.scanStart] = parsePmagDate(mr.hdr.date, mr.hdr.time);
if (mr.info.scanStart(1)<2005) & (mr.info.scanStart(2)<3)  % GE switched from big Endian to little Endian in March 2005
    littleEndian = 0;
else
    littleEndian = 1;
end
mr.info.effectiveResolution = func.effectiveResolution;
mr.info.coil = mr.hdr.coil;

% now read data, get other info:
mr.data = readMagFile(pth, slices, littleEndian);
if ~isempty(t)
	% sub-select time points
	mr.data = mr.data(:,:,:,t);
end
pixdim = mr.hdr.FOV/mr.hdr.equivMatSize;
mr.voxelSize = [pixdim pixdim mr.hdr.sliceThickness mr.hdr.tAcq/1000];
mr.dims = size(mr.data);
mr.dims(3) = size(mr.data, 3);  % in case of 2-D or 3-D data:
mr.dims(4) = size(mr.data, 4);  % enforce 1x4 dims vector
mr.extent = mr.voxelSize .* mr.dims;
mr.spaces = mrStandardSpaces(mr);
mr.dimUnits = {'mm' 'mm' 'mm' 'sec'};
mr.dataUnits = 'T2*-Weighted Intensity';
mr.dataRange = mrvMinmax(mr.data);

% build a space defining the scanner coords, where the
% xform maps from the pixCorners to the R|A|S coords of the
% the three corners from the header:
try
	mr.spaces(end+1).name = 'Scanner';
	mr.spaces(end).xform = mrReadMag_affineScannerXform(mr.hdr);
	mr.spaces(end).xform(1:3,4) = mr.extent(1:3) ./ 2;
	mr.spaces(end).dirLabels = {'L <--> R' 'P <--> A'  'I <--> S'};
	mr.spaces(end).sliceLabels =  {'Sagittal' 'Coronal' 'Axial'};
	mr.spaces(end).units = 'mm';
	mr.spaces(end).coords = [];
	mr.spaces(end).indices = [];

    % NOT YET IMPLEMENTED:
% 	% also update the direction labels on the pixel and L/R flipped spaces,
% 	% using the header info:
% 	dirs = mrEfileDirections(hdrPath);
% 	for i=1:3, mr.spaces(i).dirLabels = dirs; end
% 	mr.spaces(3).dirLabels{2} = dimFlip(mr.spaces(3).dirLabels{2});

catch
	if prefsVerboseCheck
		warning('Couldn''t read E-file header alignment fields.')
	end
end

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

% /-------------------------------------------------------------------/ %
function xform = mrReadMag_affineScannerXform(hdr);
%% given a Mag-file header, compute a 4x4 affine transform mapping from
%% pixel space into the scanner coordinate system.

% E-file headers give us the scanner coordinates of three points from
% each slice:
%	gw_point1 = upper left-hand corner  (x=1, y=1)
%	gw_point2 = lower left-hand corner  (x=1, y=nrows)
%	gw_point3 = upper right-hand corner (x=ncols, y=1)
% We know which points in pixel space coorrespond to these points, so we
% take the combined set of corner coordinates, and the corresponding
% scanner coords, and solve for the best xform to map between them.

% combine the scanner coordinates together
scanCoords = [hdr.gw_point1, hdr.gw_point2, hdr.gw_point3];

% get the corresponding pixel coordinates
[X Y Z] = meshgrid([1 hdr.imgsize], [1 hdr.imgsize], 1:hdr.slquant);
pixCoords = [Y(:) X(:) Z(:)]';

% now we have to clip and resort: every 4th coord here is the
% lower-right-hand corner, which we don't have, so remove these:
pixCoords(:,2:4:end) = [];

% do I have the order of points correct?
pixCoords = pixCoords([1 2 3],:);

% now re-sort to match the order of points in scanCoords:
% [ulhc across slices, llhc across slices, urhc across slices]:
pixCoords = [pixCoords(:,[1:3:end]), ...
			 pixCoords(:,[2:3:end]), ...
			 pixCoords(:,[3:3:end])];

% solve for the xform
xform = affineSolve(pixCoords, scanCoords);

% one last step: we need to account for the shift specified by the start
% location in the header:
xform(1:3,4) = hdr.start_loc;

xform = inv(xform);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = panel_sourcegrid(varargin)
% PANEL_SOURCEGRID: Options for the construction of volume source grid.
%
% USAGE:     bstPanel = panel_sourcegrid('CreatePanel', sProcess, sFiles)   : Called from the pipeline editor
%            bstPanel = panel_sourcegrid('CreatePanel', CortexFile)         : Called from the interactive interface
%                grid = panel_sourcegrid('GetPanelContents')                : When called from the interactive interface
%         GridOptions = panel_sourcegrid('GetPanelContents')                : When called from the pipeline editor
%                grid = panel_sourcegrid('GetGrid', GridOptions, CortexSurface)

% @=============================================================================
% 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, 2011-2021

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel(sProcess, sFiles)  %#ok<DEFNU>
    panelName = 'SourceGrid';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    % CALL: From GUI
    if (nargin == 1)
        isPreview = 1;
        CortexFile = sProcess;
        sProcess = [];
    % CALL: From Process
    elseif (nargin == 2)
        isPreview = 0;
        CortexFile = [];
    end
    TemplateFile = [];
    isShowGroup = 1;
    % Create main main panel
    jPanelNew = gui_river();

    % Prepare online computation of the grid and preview
    if isPreview
        GridOptions_name = 'GridOptions_headmodel';
        % Default options
        GridOptions = bst_get(GridOptions_name);
        % Find subject
        [sSubject, iSubject] = bst_get('SurfaceFile', CortexFile);
        % If this is the default subject: do not show the option "use group grid"
        if (iSubject == 0)
            isShowGroup = 0;
        % Else
        else
            % Get the default headmodel file from the group analysis subject
            [sSubjectGroup, iSubjectGroup] = bst_get('Subject', bst_get('NormalizedSubjectName'));
            if ~isempty(sSubjectGroup)
                sStudy = bst_get('DefaultStudy', iSubjectGroup);
                if ~isempty(sStudy.HeadModel) && ~isempty(sStudy.iHeadModel)
                    TemplateFile = sStudy.HeadModel(sStudy.iHeadModel).FileName;
                    GridOptions.FileName = TemplateFile;
                end
            end
        end
    % Calling non-interactive option panel from the pipeline editor
    else
        % Get type of process
        if strcmpi(func2str(sProcess.Function), 'process_ft_dipolefitting')
            GridOptions_name = 'GridOptions_dipfit';
        else
            GridOptions_name = 'GridOptions_headmodel';
        end
        % Get options from the process
        GridOptions = sProcess.options.volumegrid.Value;
        % Get default values if not defined yet
        if isempty(GridOptions)
            GridOptions = bst_get(GridOptions_name);
        end
    end

    % ===== GRID OPTIONS =====
    jPanelOpt = gui_river([4,5], [0,15,20,10], 'Grid options');
        jButtonGroup = ButtonGroup();
        jButtonGroupSurf = ButtonGroup();
        % RADIO: Generate grid
        jRadioGenerate = gui_component('radio', jPanelOpt, '', 'Generate from cortex surface (adaptive):', jButtonGroup, [], @(h,ev)UpdatePanel, []);
        % nLayers
        gui_component('label', jPanelOpt, 'br', '     ');
        jLabelLayers = gui_component('label', jPanelOpt, '', 'Number of layers: ', [], [], [], []);
        jTextLayers = gui_component('texttime', jPanelOpt, 'tab', num2str(GridOptions.nLayers, '%d'), [], [], [], []);
        java_setcb(jTextLayers, 'ActionPerformedCallback', @(h,ev)OptionsChanged_Callback, ...
                                'FocusLostCallback',       @(h,ev)OptionsChanged_Callback);
        jLabelLayersDef = gui_component('label', jPanelOpt, '', '  (default: 17)', [], [], [], []);
        % Reduction
        gui_component('label', jPanelOpt, 'br', '     ');
        jLabelReduction = gui_component('label', jPanelOpt, '', 'Downsampling factor: ', [], [], [], []);
        jTextReduction = gui_component('texttime', jPanelOpt, 'tab', num2str(GridOptions.Reduction, '%d'), [], [], [], []);
        java_setcb(jTextReduction, 'ActionPerformedCallback', @(h,ev)OptionsChanged_Callback, ...
                                   'FocusLostCallback',       @(h,ev)OptionsChanged_Callback);
        jLabelReductionDef = gui_component('label', jPanelOpt, '', '  (default: 3)', [], [], [], []);
        % nVerticesInit
        gui_component('label', jPanelOpt, 'br', '     ');
        jLabelVertInit = gui_component('label', jPanelOpt, '', 'Initial number of vertices: ', [], [], [], []);
        jTextVertInit = gui_component('texttime', jPanelOpt, 'tab', num2str(GridOptions.nVerticesInit, '%d'), [], [], [], []);
        java_setcb(jTextVertInit, 'ActionPerformedCallback', @(h,ev)OptionsChanged_Callback, ...
                                  'FocusLostCallback',       @(h,ev)OptionsChanged_Callback);
        jLabelVertInitDef = gui_component('label', jPanelOpt, '', '  (default: 4000)', [], [], [], []);
        gui_component('label', jPanelOpt, 'br', ' ');

        % RADIO: Isotropic
        jRadioIsotropic = gui_component('radio', jPanelOpt, 'br', 'Regular grid (isotropic):', jButtonGroup, [], @(h,ev)UpdatePanel, []);
        jRadioCortex = gui_component('radio', jPanelOpt, [], 'Cortex', jButtonGroupSurf, [], @(h,ev)UpdatePanel, []);
        jRadioSkull = gui_component('radio', jPanelOpt, [], 'Inner skull', jButtonGroupSurf, [], @(h,ev)UpdatePanel, []);
        jRadioHead  = gui_component('radio', jPanelOpt, [], 'Head', jButtonGroupSurf, [], @(h,ev)UpdatePanel, []);
        % Grid resolution
        gui_component('label', jPanelOpt, 'br', '     ');
        jLabelResolution = gui_component('label', jPanelOpt, '', 'Grid resolution: ', [], [], [], []);
        jTextResolution = gui_component('texttime', jPanelOpt, 'tab', num2str(GridOptions.Resolution * 1000), [], [], [], []);
        java_setcb(jTextResolution, 'ActionPerformedCallback', @(h,ev)OptionsChanged_Callback, ...
                                    'FocusLostCallback',       @(h,ev)OptionsChanged_Callback);
        jLabelResUnits = gui_component('label', jPanelOpt, '', ' mm');

        % RADIO: Load from file
        jRadioFile = gui_component('radio', jPanelOpt, 'br', 'Load from file [Nx3 double]', jButtonGroup, [], @(h,ev)UpdatePanel, []);
        % Filename
        gui_component('label', jPanelOpt, 'br', '     ');
        jTextFile = gui_component('text', jPanelOpt, 'hfill', '', [], [], [], []);
        jTextFile.setEditable(0);
        jButtonFile = gui_component('button', jPanelOpt, '', '...', [], [], @(h,ev)SelectFileVar, []);

        % RADIO: Load from variable
        jRadioVar = gui_component('radio', jPanelOpt, 'br', 'Load from Matlab variable [Nx3 double]', jButtonGroup, [], @(h,ev)UpdatePanel, []);
        gui_component('label', jPanelOpt, 'br', '     ');
        jTextVar = gui_component('text', jPanelOpt, 'hfill', '', [], [], [], []);
        java_setcb(jTextVar, 'ActionPerformedCallback', @(h,ev)UpdatePanel);
        java_setcb(jTextVar, 'FocusLostCallback', @(h,ev)UpdatePanel);

        % RADIO: Use template grid for group analysis
        if isShowGroup
            jRadioGroup = gui_component('radio', jPanelOpt, 'br', 'Use template grid for group analysis', jButtonGroup, [], @(h,ev)UpdatePanel, []);
            if ~isempty(TemplateFile) || ~isPreview
                gui_component('label', jPanelOpt, 'br', '       Loads the grid from the default head model in "Group analysis"');
            else
                gui_component('label', jPanelOpt, 'br', '       No head model available in "Group analysis"');
                jRadioGroup.setEnabled(0);
            end
        else
            jRadioGroup = [];
        end

        % Estimated number of vertices
        if isPreview
            gui_component('label', jPanelOpt, 'br', '     ');
            gui_component('label', jPanelOpt, 'br', 'Estimated number of grid points: ', [], [], [], []);
            jLabelPoints = gui_component('label', jPanelOpt, '', '2000', [], [], [], []);
        end
    jPanelNew.add('br hfill', jPanelOpt);
    % Select the appropriate radio button
    switch (GridOptions.Method)
        case 'adaptive'
            jRadioGenerate.setSelected(1);
            jRadioCortex.setSelected(1);
        case 'isotropic'
            jRadioIsotropic.setSelected(1);
            jRadioCortex.setSelected(1);
        case 'isoskull'
            jRadioIsotropic.setSelected(1);
            jRadioSkull.setSelected(1);
        case 'isohead'
            jRadioIsotropic.setSelected(1);
            jRadioHead.setSelected(1);
        case 'file'
            jRadioFile.setSelected(1);
            jTextFile.setText(GridOptions.FileName);
            jRadioCortex.setSelected(1);
        case 'var'
            jRadioVar.setSelected(1);
            jTextVar.setText(GridOptions.FileName);
            jRadioCortex.setSelected(1);
        case 'group'
            if ~isempty(jRadioGroup) && jRadioGroup.isEnabled()
                jRadioGroup.setSelected(1);
            else
                jRadioGenerate.setSelected(1);
            end
            jRadioCortex.setSelected(1);
        otherwise
            jRadioGenerate.setSelected(1);
    end

    % ===== VALIDATION BUTTONS =====
    if isPreview
        gui_component('button', jPanelNew, 'br right', 'Preview', [], [], @(h,ev)ShowGrid, []);
        constrNext = [];
    else
        constrNext = 'br right';
    end
    gui_component('button', jPanelNew, constrNext, 'Cancel', [], [], @ButtonCancel_Callback, []);
    gui_component('button', jPanelNew, [], 'OK', [], [], @ButtonOk_Callback, []);

    % ===== PANEL CREATION =====
    % Return a mutex to wait for panel close
    bst_mutex('create', panelName);
    % Controls list
    ctrl = struct('CortexFile',       CortexFile, ...
                  'isPreview',        isPreview, ...
                  'GridOptions_name', GridOptions_name, ...
                  'jRadioGenerate',   jRadioGenerate, ...
                  'jRadioIsotropic',  jRadioIsotropic, ...
                  'jRadioCortex',     jRadioCortex, ...
                  'jRadioSkull',      jRadioSkull, ...
                  'jRadioHead',       jRadioHead, ...
                  'jRadioFile',       jRadioFile, ...
                  'jRadioVar',        jRadioVar, ...
                  'jRadioGroup',      jRadioGroup, ...
                  'jTextLayers',      jTextLayers, ...
                  'jTextReduction',   jTextReduction, ...
                  'jTextVertInit',    jTextVertInit, ...
                  'jTextResolution',  jTextResolution, ...
                  'jTextFile',        jTextFile, ...
                  'jTextVar',         jTextVar);
    % Create the BstPanel object that is returned by the function
    bstPanelNew = BstPanel(panelName, jPanelNew, ctrl);
    % Update panel
    UpdatePanel();

%% =================================================================================
%  === INTERNAL CALLBACKS ==========================================================
%  =================================================================================
%% ===== CANCEL BUTTON =====
    function ButtonCancel_Callback(hObject, event)
        % Close preview window
        hFig = findobj(0, 'type', 'figure', 'tag', 'FigCheckGrid');
        if ~isempty(hFig)
            close(hFig);
        end
        % Close panel without saving (release mutex automatically)
        gui_hide(panelName);
    end

%% ===== OK BUTTON =====
    function ButtonOk_Callback(varargin)
        % Close preview window
        hFig = findobj(0, 'type', 'figure', 'tag', 'FigCheckGrid');
        if ~isempty(hFig)
            close(hFig);
        end
        % Release mutex and keep the panel opened
        bst_mutex('release', panelName);
    end

%% ===== OPTION CHANGED =====
    function OptionsChanged_Callback(varargin)
        % Get new options
        NewOptions.nLayers       = str2double(char(ctrl.jTextLayers.getText()));
        NewOptions.Reduction     = str2double(char(ctrl.jTextReduction.getText()));
        NewOptions.nVerticesInit = str2double(char(ctrl.jTextVertInit.getText()));
        NewOptions.Resolution    = str2double(char(ctrl.jTextResolution.getText())) / 1000;
        % If options changed: update panel
        if ~isequal(NewOptions.nLayers, GridOptions.nLayers) || ~isequal(NewOptions.Reduction, GridOptions.Reduction) || ~isequal(NewOptions.nVerticesInit, GridOptions.nVerticesInit) || ~isequal(NewOptions.Resolution, GridOptions.Resolution)
            GridOptions = NewOptions;
            UpdatePanel();
        end
    end

%% ===== UPDATE PANEL =====
    function UpdatePanel(varargin)
        global gGridLoc;
        % RADIO: Generate
        isGenerate = jRadioGenerate.isSelected();
        jTextLayers.setEnabled(isGenerate);
        jLabelLayers.setEnabled(isGenerate);
        jLabelLayersDef.setEnabled(isGenerate);
        jTextReduction.setEnabled(isGenerate);
        jLabelReduction.setEnabled(isGenerate);
        jLabelReductionDef.setEnabled(isGenerate);
        jTextVertInit.setEnabled(isGenerate);
        jLabelVertInit.setEnabled(isGenerate);
        jLabelVertInitDef.setEnabled(isGenerate);
        % RADIO: Isotropic
        isIsotropic = jRadioIsotropic.isSelected();
        jLabelResolution.setEnabled(isIsotropic);
        jTextResolution.setEnabled(isIsotropic);
        jRadioCortex.setEnabled(isIsotropic);
        jRadioSkull.setEnabled(isIsotropic);
        jRadioHead.setEnabled(isIsotropic);
        jLabelResUnits.setEnabled(isIsotropic);
        % RADIO: File
        isFile = jRadioFile.isSelected();
        jTextFile.setEnabled(isFile);
        jButtonFile.setEnabled(isFile);
        % RADIO: Variable
        isVar = jRadioVar.isSelected();
        jTextVar.setEnabled(isVar);
        % Preview grid
        if isPreview
            % Get the options
            GridOptions = GetOptions(ctrl);
            % Get grid of source points
            gGridLoc = bst_call(@GetGrid, GridOptions, ctrl.CortexFile);
            % Variable not found
            if strcmpi(GridOptions.Method, 'var') && isempty(gGridLoc)
                ctrl.jTextVar.setText('');
            end
            % Count number of points
            if ~isempty(gGridLoc)
                nTotal = length(gGridLoc);
            else
                nTotal = 0;
            end
            % Display new estimation
            jLabelPoints.setText(num2str(nTotal));
            % Get previous window: If it exists, update it
            if (nTotal ~= 0)
                hFig = findobj(0, 'type', 'figure', 'tag', 'FigCheckGrid');
                if ~isempty(hFig)
                    ShowGrid();
                end
            end
        end
    end

%% ===== SELECT FILE =====
    function SelectFileVar()
        % Get file
        filename = java_getfile( 'open', 'Import grid of points', '', 'single', 'files', ...
                                {{'*'}, 'Matlab or ASCII files (*.mat;*.*)', 'ALL'}, 1);
        % Update panel
        if ~isempty(filename)
            jTextFile.setText(filename);
        end
        UpdatePanel();
    end
end

%% =================================================================================
%  === EXTERNAL CALLBACKS ==========================================================
%  =================================================================================
%% ===== GET PANEL CONTENTS =====
function s = GetPanelContents()
    global gGridLoc;
    % Get panel controls
    ctrl = bst_get('PanelControls', 'SourceGrid');
    if isempty(ctrl)
        return;
    end
    % Get grid options
    GridOptions = GetOptions(ctrl);
    % Save the new options
    bst_set(ctrl.GridOptions_name, GridOptions);
    % For "group" option: Get the default template grid
    if isequal(GridOptions.Method, 'group') && isempty(GridOptions.FileName)
        GridOptions.FileName = GetDefaultGridFile();
    end
    % GUI: Return previewed grid
    if ctrl.isPreview
        s.GridLoc = gGridLoc;
        s.GridOptions = GridOptions;
        gGridLoc = [];
        clear gGridLoc;
    % Process: Return selected options
    else
        s = GridOptions;
    end
end

%% ===== GET OPTIONS =====
function GridOptions = GetOptions(ctrl)
    % Get default structure
    GridOptions = bst_get(ctrl.GridOptions_name);
    % Get computation mode
    if ctrl.jRadioGenerate.isSelected()
        GridOptions.Method = 'adaptive';
    elseif ctrl.jRadioIsotropic.isSelected()
        if ctrl.jRadioCortex.isSelected()
            GridOptions.Method = 'isotropic';
        elseif ctrl.jRadioSkull.isSelected()
            GridOptions.Method = 'isoskull';
        elseif ctrl.jRadioHead.isSelected()
            GridOptions.Method = 'isohead';
        end
    elseif ctrl.jRadioFile.isSelected()
        GridOptions.Method = 'file';
    elseif ctrl.jRadioVar.isSelected()
        GridOptions.Method = 'var';
    elseif isfield(ctrl, 'jRadioGroup') && ~isempty(ctrl.jRadioGroup) && ctrl.jRadioGroup.isSelected()
        GridOptions.Method = 'group';
    else
        GridOptions.Method = 'adaptive';
    end
    % Adaptive options
    GridOptions.nLayers       = str2double(char(ctrl.jTextLayers.getText()));
    GridOptions.Reduction     = str2double(char(ctrl.jTextReduction.getText()));
    GridOptions.nVerticesInit = str2double(char(ctrl.jTextVertInit.getText()));
    % Check for errors
    if isnan(GridOptions.nLayers) || isnan(GridOptions.Reduction) || isnan(GridOptions.nVerticesInit)
        bst_error('Invalid values.', 'Generate grid', 0);
        GridOptions = [];
        return
    end
    % Isotropic options
    GridOptions.Resolution = str2double(char(ctrl.jTextResolution.getText())) / 1000;
    % External inputs
    if strcmpi(GridOptions.Method, 'file')
        GridOptions.FileName = char(ctrl.jTextFile.getText());
    elseif strcmpi(GridOptions.Method, 'var')
        GridOptions.FileName = char(ctrl.jTextVar.getText());
    else
        GridOptions.FileName = [];
    end
end

%% ===== GET GRID =====
% USAGE:  grid = GetGrid(GridOptions, CortexFile)
function grid = GetGrid(GridOptions, CortexFile)
    grid = [];
    % Progress bar
    bst_progress('start', 'Volume grid', 'Creating grid...');
    % Get subject surfaces
    sSubject = bst_get('SurfaceFile', CortexFile);
    if ~isempty(sSubject.iInnerSkull)
        InnerSkullFile = sSubject.Surface(sSubject.iInnerSkull).FileName;
    else
        InnerSkullFile = [];
    end
    if ~isempty(sSubject.iScalp)
        HeadFile = sSubject.Surface(sSubject.iScalp).FileName;
    else
        HeadFile = [];
    end

    % Switch between methods
    switch (GridOptions.Method)
        case 'adaptive'
            % Compute grid
            grid = bst_sourcegrid(GridOptions, CortexFile, InnerSkullFile);

        case 'isotropic'
            grid = bst_sourcegrid(GridOptions, CortexFile, InnerSkullFile);

        case 'isoskull'
            if isempty(InnerSkullFile)
                error('No inner skull surface available.');
            end
            grid = bst_sourcegrid(GridOptions, InnerSkullFile);

        case 'isohead'
            if isempty(HeadFile)
                error('No head surface available.');
            end
            grid = bst_sourcegrid(GridOptions, HeadFile);

        case {'file', 'var'}
            % Read grid
            if ~isempty(GridOptions.FileName)
                grid = ReadGrid(GridOptions.FileName, GridOptions.Method);
            end

        case 'group'
            % === GET SUBJECT ===
            % Get subject using cortex surface
            [sSubject, iSubject] = bst_get('SurfaceFile', CortexFile);
            if isempty(sSubject.Anatomy) || isempty(sSubject.iAnatomy)
                return;
            end
            % Load subject MRI
            SubjectMriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
            sMriSubject = in_mri_bst(SubjectMriFile);

            % === GET TEMPLATE SUBJECT ===
            % For "group" option: Get the default template grid
            if isempty(GridOptions.FileName)
                GridOptions.FileName = GetDefaultGridFile();
                if isempty(GridOptions.FileName)
                    return;
                end
            end
            % Get template study using the grid file
            [sStudyGroup, iStudyGroup] = bst_get('HeadModelFile', GridOptions.FileName);
            if isempty(sStudyGroup)
                return;
            end
            % Get template subject
            [sSubjectGroup, iSubjectGroup] = bst_get('Subject', sStudyGroup.BrainStormSubject);
            if isempty(sSubjectGroup) || isempty(sSubjectGroup.Anatomy) || isempty(sSubjectGroup.iAnatomy)
                return;
            end
            % Load template MRI
            GroupMriFile = sSubjectGroup.Anatomy(sSubjectGroup.iAnatomy).FileName;
            sMriGroup = in_mri_bst(GroupMriFile);

            % === CONVERT TEMPLATE GRID TO SUBJECT SPACE ===
            % Load grid file
            HeadModelMat = in_bst_headmodel(GridOptions.FileName, 0, 'GridLoc');
            gridGroup = HeadModelMat.GridLoc;
            if isempty(gridGroup)
                return;
            end
            % Convert grid coordinates: Template SCS => MNI
            gridMni = cs_convert(sMriGroup, 'scs', 'mni', gridGroup);
            if isempty(gridMni)
                return;
            end
            % Convert grid coordinates: MNI => Subject SCS
            grid = cs_convert(sMriSubject, 'mni', 'scs', gridMni);
            if isempty(grid)
                disp('PROJECT> Error: The MNI transformation is not available for this subject.');
            end
    end
    % Close progress bar
    bst_progress('stop');
end

%% ===== GET DEFAULT GRID FILE =====
function HeadModelFile = GetDefaultGridFile()
    HeadModelFile = [];
    % Get the group analysis subject
    [sSubjectGroup, iSubjectGroup] = bst_get('Subject', bst_get('NormalizedSubjectName'));
    if isempty(sSubjectGroup) || isempty(sSubjectGroup.Anatomy) || isempty(sSubjectGroup.iAnatomy)
        return;
    end
    % Get default study of group analysis subject
    sStudy = bst_get('DefaultStudy', iSubjectGroup);
    % Exit if there is no head model in this subject
    if isempty(sStudy.HeadModel) || isempty(sStudy.iHeadModel)
        return;
    end
    % Return default headmodel file
    HeadModelFile = sStudy.HeadModel(sStudy.iHeadModel).FileName;
end

%% ===== READ GRID =====
function grid = ReadGrid(varname, Method)
    grid = [];
    % Read matrix from file
    if strcmpi(Method, 'file')
        FileName = varname;
        % File doesn't exist
        if ~file_exist(FileName)
            disp(['Error: File does not exist: ' 10 FileName]);
            return;
        end
        % Get extension
        [tmp, tmp, fExt] = bst_fileparts(FileName);
        % Matlab matrix
        if strcmpi(fExt, '.mat')
            % Other kinds of files
            try
                % Load file
                filemat = load(varname);
                % Get field "GridLoc" if it exists
                if isfield(filemat, 'GridLoc')
                    grid = filemat.GridLoc;
                else
                    allfields = fieldnames(filemat);
                    if (length(allfields) == 1)
                        grid = filemat.(allfields{1});
                    end
                end
                % Return an error if there is more than one matrix
                if isempty(grid)
                    disp('Error: The file must be a valid Brainstorm headmodel file, or a Matlab matrix with only one variable.');
                end
            catch
                disp('Error: Cannot read as .mat file.');
            end
        % ASCII file
        else
            % Other kinds of files
            try
                grid = load(varname, '-ascii');
            catch
                disp('Error: Cannot read as ASCII file.');
            end
        end
    % Read matrix from variable in base workspace
    elseif strcmpi(Method, 'var')
        grid = in_matlab_var(varname, 'numeric');
    else
        error('Unknown method.');
    end
    % Re-orient matrix
    if ~isempty(grid)
        if (size(grid,1) ~= 3) && (size(grid,2) ~= 3)
            disp('BST> Invalid grid format. Matrix must be [Nx3].');
        elseif (size(grid,1) == 3) && (size(grid,2) ~= 3)
            grid = grid';
        end
    end
end

%% ===== SHOW GRID =====
function ShowGrid()
    global gGridLoc;
    % Get panel controls
    ctrl = bst_get('PanelControls', 'SourceGrid');
    if isempty(ctrl)
        return;
    end
    % Get previous window
    hFig = findobj(0, 'type', 'figure', 'tag', 'FigCheckGrid');
    % If there is no grid to show
    if isempty(gGridLoc)
        % Close existing figure
        if ~isempty(hFig)
            close(hFig);
        end
        return;
    end
    % Create figure if it doesnt exist + show surface
    if isempty(hFig)
        hFig = view_surface(ctrl.CortexFile, .9, [.6 .6 .6], 'NewFigure', 0);
        set(hFig, 'Tag', 'FigCheckGrid');
    % Figure exists: remove previous points
    else
        delete(findobj(hFig, 'tag', 'ptCheckGrid'));
        % Focus on figure
        figure(hFig);
    end
    % No points to show: exit
    if isempty(gGridLoc)
        return;
    end
    % Get axes
    hAxes = findobj(hFig, 'Tag', 'Axes3D');
    % Show grid points
    line(gGridLoc(:,1), gGridLoc(:,2), gGridLoc(:,3), 'LineStyle', 'none', ...
         'Color', [0 1 0], 'MarkerSize', 2, 'Marker', '.', ...
         'Tag', 'ptCheckGrid', 'Parent', hAxes);
end

author: brainstorm-tools


text:
% DESCRIPTION:
%       subscript to display time steps and maximum supported frequency.
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 8th July 2014
%       last update - 8th July 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% display time step information
disp(['  dt: ' scaleSI(dt) 's, t_end: ' scaleSI(kgrid.t_array(end)) 's, time steps: ' num2str(length(kgrid.t_array))]);

% if using the elastic code, get the minimum sound speeds (not including
% zero if set for the shear speed)
if ~elastic_code
    c_min = min(medium.sound_speed(:));
else
    c_min_comp = min(medium.sound_speed_compression(:));
    c_min_shear = min(medium.sound_speed_shear(medium.sound_speed_shear ~= 0));
end

% get suitable scaling factor
grid_size_metric = [kgrid.x_size, kgrid.y_size, kgrid.z_size];
[x_sc, scale, prefix] = scaleSI( min(grid_size_metric(grid_size_metric ~= 0)) ); %#ok<*ASGLU>
clear grid_size_metric;

% display the grid size and maximum supported frequency
switch kgrid.dim
    case 1

        % display grid size
        disp(['  input grid size: ' num2str(kgrid.Nx) ' grid points (' scaleSI(kgrid.x_size) 'm)']);

        % display maximum supported frequency
        disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max * c_min / (2*pi) ) 'Hz']);

    case 2

        % display grid size
        disp(['  input grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Ny) ' grid points (' num2str(kgrid.x_size*scale) ' by ' num2str(kgrid.y_size*scale) prefix 'm)']);

        if ~elastic_code

            % display maximum supported frequency
            if kgrid.kx_max == kgrid.ky_max
                disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max * c_min / (2*pi) ) 'Hz']);
            else
                disp(['  maximum supported frequency: ' scaleSI( kgrid.kx_max * c_min / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * c_min / (2*pi) ) 'Hz']);
            end

        else

            % display the maximum supported frequency
            if kgrid.kx_max == kgrid.ky_max
                disp(['  maximum supported compressional frequency: ' scaleSI( kgrid.k_max * c_min_comp / (2*pi) ) 'Hz']);
                if isempty(c_min_shear)
                    disp('  maximum supported shear frequency: 0Hz');
                else
                    disp(['  maximum supported shear frequency: ' scaleSI( kgrid.k_max * c_min_shear / (2*pi) ) 'Hz']);
                end
            else
                disp(['  maximum supported compressional frequency: ' scaleSI( kgrid.kx_max * c_min_comp / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * c_min_comp / (2*pi) ) 'Hz']);
                if isempty(c_min_shear)
                    disp('  maximum supported shear frequency: 0Hz');
                else
                    disp(['  maximum supported compressional frequency: ' scaleSI( kgrid.kx_max * c_min_shear / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * c_min_shear / (2*pi) ) 'Hz']);
                end
            end

        end

    case 3

        % display grid size
        disp(['  input grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Ny) ' by ' num2str(kgrid.Nz) ' grid points (' num2str(kgrid.x_size*scale) ' by ' num2str(kgrid.y_size*scale) ' by ' num2str(kgrid.z_size*scale) prefix 'm)']);

        if ~elastic_code

            % display maximum supported frequency
            if (kgrid.kx_max == kgrid.kz_max) && (kgrid.kx_max == kgrid.ky_max)
                disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max * c_min / (2*pi) ) 'Hz']);
            else
                disp(['  maximum supported frequency: ' scaleSI( kgrid.kx_max * c_min / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * c_min / (2*pi) ) 'Hz by ' scaleSI( kgrid.kz_max * c_min / (2*pi) ) 'Hz']);
            end

        else

            % display the maximum supported frequency
            if (kgrid.kx_max == kgrid.kz_max) && (kgrid.kx_max == kgrid.ky_max)
                disp(['  maximum supported compressional frequency: ' scaleSI( kgrid.k_max * c_min_comp / (2*pi) ) 'Hz']);
                if isempty(c_min_shear)
                    disp('  maximum supported shear frequency: 0Hz');
                else
                    disp(['  maximum supported shear frequency: ' scaleSI( kgrid.k_max * c_min_shear / (2*pi) ) 'Hz']);
                end
            else
                disp(['  maximum supported frequency: ' scaleSI( kgrid.kx_max * c_min_comp / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * c_min_comp / (2*pi) ) 'Hz by ' scaleSI( kgrid.kz_max * c_min_comp / (2*pi) ) 'Hz']);
                if isempty(c_min_shear)
                    disp('  maximum supported shear frequency: 0Hz');
                else
                    disp(['  maximum supported frequency: ' scaleSI( kgrid.kx_max * c_min_shear / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * c_min_shear / (2*pi) ) 'Hz by ' scaleSI( kgrid.kz_max * c_min_shear / (2*pi) ) 'Hz']);
                end
            end

        end
end

% cleanup unused variables
clear c_min c_min_comp c_min_shear grid_size_metric

repo: Aero-matlab


text:
%compute the ycoordinate of the center in the tails coordinate system

function [data,units]=compute_ytailcentral_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
ytailcentral_mm=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    ytailcentral_mm{1,i}=trx(larva).ycentral_mm-trx(larva).ytail_mm;
end

units=parseunits('mm');
data=ytailcentral_mm;

author: kristinbranson

repo: JAABA


text:
% MANIPULATOR TRAJECTORY GENERATION
% Generates combined transform (rotation and translation) trajectories
% using custom time scaling from a separate trajectory.
%
% Copyright 2019 The MathWorks, Inc.

%% Setup
clear, clc, close all

% Define waypoint information
createWaypointData;

% Define IK
ik = inverseKinematics('RigidBodyTree',gen3);
ikWeights = [1 1 1 1 1 1];
ikInitGuess = gen3.homeConfiguration;

% Set up plot
plotMode = 2; % 0 = None, 1 = Trajectory, 2 = Coordinate Frames
show(gen3,gen3.homeConfiguration,'Frames','off','PreservePlot',false);
xlim([-1 1]), ylim([-1 1]), zlim([0 1.2])
hold on
if plotMode == 1
    hTraj = plot3(waypoints(1,1),waypoints(2,1),waypoints(3,1),'b.-');
end
plot3(waypoints(1,:),waypoints(2,:),waypoints(3,:),'ro','LineWidth',2);

%% Generate and follow trajectory
% Loop through segments one at a time
trajType = 'trap';
numWaypoints = size(waypoints,2);
for w = 1:numWaypoints-1

    % Get the initial and final transforms and times for the segment
    T0 = trvec2tform(waypoints(:,w)') * eul2tform(orientations(:,w)');
    Tf = trvec2tform(waypoints(:,w+1)') * eul2tform(orientations(:,w+1)');
    timeInterval = waypointTimes(w:w+1);
    trajTimes = timeInterval(1):ts:timeInterval(2);

    % Generate time scaling trajectory for the segment on the range [0 1]
    switch trajType
        case 'trap'
            [s,sd,sdd] = trapveltraj([0 1],numel(trajTimes), ...
                                     'EndTime',diff(timeInterval));
        case 'cubic'
            [s,sd,sdd] = cubicpolytraj([0 1],timeInterval,trajTimes);
        case 'quintic'
            [s,sd,sdd] = quinticpolytraj([0 1],timeInterval,trajTimes);
        otherwise
            error('Invalid trajectory type! Use ''trap'', ''cubic'', or ''quintic''');
    end

    % Find the transforms from trajectory generation
    [T,vel,acc] = transformtraj(T0,Tf,timeInterval,trajTimes, ...
                                'TimeScaling',[s;sd;sdd]);

    % Trajectory visualization for the segment
    if plotMode == 1
        eePos = tform2trvec(T);
        set(hTraj,'xdata',eePos(:,1),'ydata',eePos(:,2),'zdata',eePos(:,3));
    elseif plotMode == 2
        plotTransforms(tform2trvec(T),tform2quat(T),'FrameSize',0.05)
    end

    % Trajectory following for the segment
    for idx = 1:numel(trajTimes)
        % Solve IK
        tgtPose = T(:,:,idx);
        [config,info] = ik(eeName,tgtPose,ikWeights,ikInitGuess);
        ikInitGuess = config;

        % Show the robot
        show(gen3,config,'Frames','off','PreservePlot',false);
        title(['Trajectory at t = ' num2str(trajTimes(idx))])
        drawnow
    end

end

author: mathworks-robotics

repo: trajectory-planning-robot-manipulators


text:
function [x,y] = allPairs(x,y)
% all pairs of elements of x and y modulo permutation
%

if nargin == 1

  [x,y] = meshgrid(x,x);

  x = x(tril(ones(size(x)))>0);
  y = y(tril(ones(size(y)))>0);

else

  %
  x = x(:);
  y = y(:);
  iseq = bsxfun(@eq,x,y.');

  [ix,iy] = find(iseq);

  x = [x(ix);x(~any(iseq,2))];
  y = [y(iy);y(~any(iseq,1))];

  % all pairs
  [x,y] = meshgrid(x,y);

  % remove double pars
  A = zeros(size(x));
  A(1:length(ix),1:length(ix)) = 1;
  A = ~triu(A,1);

  x = x(A);
  y = y(A);

end

if nargout < 2, x = [x(:),y(:)]; end

author: mtex-toolbox

repo: mtex

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