text:
function varargout = setaxes(varargin)
%SETAXES   fixes common problems with MATLAB figures.
%
%   SETAXES(AX, ...) applies to the axes AX which, when not supplied, defaults to the current axes.
%
%
%    SETAXES('xoffset') [2-D view only] adjusts the axes position in the
%    horizontal direction to make up for any clipped (left or right) part
%    of the axes based on some default settings.
%
%    SETAXES('yoffset') [2-D view only] adjusts the axes position in the
%    vertical direction to make up for any clipped (lower or upper) part
%    of the axes based on some default settings.
%
%    SETAXES('axesmoveresize', [L B R T]) adds the array [L B -R -T] to the
%    'OuterPosition' property  of the axes.
%
%
%    SETAXES('xtick2text', OFFSET) [2-D view only] replaces the tick labels
%    on the X-Axis by text objects where OFFSET is the factor of the axis
%    height used as the vertical margin from the X-Axis.
%
%       When not supplied, OFFSET defaults to 0.03.
%
%    SETAXES('xtick2text',..., PROP1, VALUE1, PROP2, VALUE2, ...)
%       sets the values of the specified properties of the text objects.
%
%    HT = SETAXES('xtick2text',...) returns in HT the handles of the text objects.
%
%
%    SETAXES('ytick2text', OFFSET) [2-D view only] replaces the tick labels
%    on the Y-Axis by text objects where OFFSET is the factor of the axis
%    width used as the horizontal margin from the Y-Axis.
%
%       When not supplied, OFFSET defaults to 0.02.
%
%    SETAXES('ytick2text',..., PROP1, VALUE1, PROP2, VALUE2, ...)
%       sets the values of the specified properties of the text objects.
%
%    HT = SETAXES('ytick2text',...) returns in HT the handles of the text objects.
%
%
%    SETAXES('xlabelcorner') [2-D view only] places the XLABEL in an
%    appropriate corner. The most common corner is the right-bottom when
%    the X-Axis is to the bottom and the ticks are increasing from left to
%    right. The corner changes accordingly for other combinations of axes properties.
%
%    SETAXES('xlabelcorner', PROP1, VALUE1, PROP2, VALUE2, ...)
%    sets the values of the specified properties of the XLabel object.
%
%
%    SETAXES('ylabelcorner') [2-D view only] places the YLABEL in an
%    appropriate corner. The most common corner is the top-left when the
%    Y-Axis is to the left and the ticks are increasing upward. The corner
%    changes accordingly for other combinations of axes properties.
%
%    SETAXES('ylabelcorner', PROP1, VALUE1, PROP2, VALUE2, ...)
%    sets the values of the specified properties of the YLabel object.
%
%
%    SETAXES('axesarrows') [2-D view only] adds arrow annotation objects to
%    X- and Y-Axis. The default 'HeadLength' and 'HeadWidth' are both 5.
%    Line width and Color are inherited from the the corresponfing axis line.
%
%    SETAXES('axesarrows','xx') adds arrow annotation objects to X-Axis only.
%    SETAXES('axesarrows','yy') adds arrow annotation objects to Y-Axis only.
%
%    SETAXES('axesarrows',..., PROP1, VALUE1, PROP2, VALUE2, ...)
%    sets the values of the specified properties of the arrow objects.
%
%    HA = SETAXES('axesarrows',...) returns in HA the handles of the arrow objects.
%
%    Note that since SETAXES('axesarrows',...) changes the axes position property, any
%    annotation objects should be added after SETAXES('axesarrows',...).
%
%=======================================================================
%                   EXAMPLE trying to cover a few features:
%-----------------------------------------------------------------------
%
% figure('Units', 'centimeters', ...
%        'Position', [2 2 8 6], ...
%        'PaperUnits', 'centimeters', ...
%        'PaperSize', [8 6], ...
%        'PaperPositionMode', 'auto')
%
% ax = axes('FontSize', 8, ...
%            'LineWidth', 0.4, ...
%            'Box', 'off', ...
%            'TickDir', 'out' );
%
% hold all
%
% t = linspace(0,5,51);
% n = 1:49;
% rand('state', sum(100*clock));
% W = [0 cumsum(-log(rand(size(n)))./n)];
% hPlot(1) = stairs(W, [n 50], '-b');
% hPlot(2) = line(t, exp(t), 'LineStyle', '--', 'Color', 'r');
% axis([0 5 0 50])
%
% xlabel('$T=\mathrm{Exp}(\frac{1}{X})$', 'Interpreter','LaTeX');
% ylabel('$X$', 'Interpreter','LaTeX');
% legend(hPlot, {'sample $n$'; 'mean $\mu(T)$'}, 'Interpreter','LaTeX')
% set(ax, 'XTickLabel', {'$0$';'$\tau_1$';'$\tau_2$'; ...
%                          '$\tau_3$';'$\tau_4$';'$\tau_5$'}, ...
%         'YTickLabel', {'$0$';'$X_1$';'$X_2$';'$X_3$';'$X_4$';'$X_5$'} )
%
% setaxes('xtick2text')
% setaxes('ytick2text')
% setaxes('ylabelcorner')
% setaxes('yoffset')
% setaxes('axesarrows')
%
% [xfig,yfig] = dsxy2figxy([1 2], [45 45]);
% har = annotation('doublearrow', xfig , yfig);
% set(har, 'LineWidth', .25, 'Color', 'k', ...
%          'Head1Length', 3.5, 'Head2Length', 3.5, ...
%          'Head1Width', 3.5, 'Head2Width', 3.5, ...
%          'Head1Style', 'vback3', 'Head2Style', 'vback3')
%
% text(1.5, 48, '$\Delta \tau$', ...
%             'VerticalAlignment', 'bottom', ...
%             'HorizontalAlignment', 'center', 'Interpreter','LaTeX');
% print -dpdf setaxes_example.pdf
%
%
%=========================================================================

% Mukhtar Ullah
% mukhtar.ullah@informatik.uni-rostock.de
% First created: Sep 18, 2008
% Last modified: Oct 18, 2010
%
%=========================================================================

persistent SetAxes_CurrentAxes SetAxes_Functions SetAxes_ValidNames

[ax,argsmain,nIn] = axescheck(varargin{:});

error(nargchk(1, inf, nIn));

if isempty(ax)
    ax = gca;
else
    assert(~strcmp(get(ax,'Tag'),'legend'), 'Legend Axes not supported')
end

% Bring to focus
figure(get(ax, 'Parent'))
set(gcf, 'CurrentAxes', ax)

% first time run or a successive run for a different axes
if isempty(SetAxes_CurrentAxes) ||  ax~=SetAxes_CurrentAxes
    SetAxes_CurrentAxes = ax;
    SetAxes_Functions = {};
end

if isempty(SetAxes_Functions)
    oldaxunits = get(ax, 'Units');
    set(ax, 'Units', 'normalized')

    axprop = get(ax, {'Position', ...
        'OuterPosition', ...
        'LooseInset', ...
        'XTick', ...
        'YTick', ...
        'XLim', ...
        'YLim', ...
        'Xlabel', ...
        'Ylabel', ...
        'Title'} );
    set(ax, 'Units', oldaxunits)

    numvar = numel(axprop);

    axprop  = [axprop, get(ax, {'XAxisLocation', ...
        'YAxisLocation', ...
        'XDir', ...
        'YDir', ...
        'XScale', ...
        'YScale', ...
        'XTickLabelMode', ...
        'YTickLabelMode'} )];

    [axpos, outpos, looseinset, xtick, ytick, axlim(1:2), axlim(3:4), ...
        hxlab, hylab, htitle, xxloc, yyloc] = axprop{1:numvar+2};

    proptf = num2cell(strcmp(axprop(numvar+1:end), ...
        {'bottom', 'left', 'normal', 'normal', 'log', 'log', 'auto', 'auto'}));

    [isxxbottom, isyyleft, isxdirnormal, isydirnormal, ...
        isxxlog, isyylog, autoxticklab, autoyticklab] = proptf{:};

    [dxt, dyt, dxtfig, dytfig, xtickroom, ytickroom] = deal(0);
    isxlabelcorner = false;
    isylabelcorner = false;

    notxxloc = setdiff({'top', 'bottom'}, xxloc);
    notyyloc = setdiff({'left', 'right'}, yyloc);
    set([hxlab, hylab, htitle], {'Tag'}, {'XLabel'; 'YLabel'; 'Title'})
    xlabpos = get(hxlab, 'Position');
    ylabpos = get(hylab, 'Position');

    [axwidth, axheight, ...
        xlabhalign, ylabhalign, xlabvalign, ylabvalign] = falign();

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

        SetAxes_ValidNames = {'axesmoveresize', 'legendmoveresize'};

        [az,el]=view(ax);
        camUp = get(ax,'CameraUpVector');
        if (az==0) && (el==90) && isequal(abs(camUp),[0 1 0]) % Is2D
            SetAxes_ValidNames = [ SetAxes_ValidNames, ...
                { 'xoffset', 'yoffset', ...
                'xtick2text', 'ytick2text', ...
                'xlabelcorner', 'ylabelcorner', 'axesarrows' } ];
        end

        SetAxes_Functions = cellfun(@(s) eval(['@', s]), SetAxes_ValidNames, 'UniformOutput', false);
end

funstr = validatestring(argsmain{1}, SetAxes_ValidNames);

fh = SetAxes_Functions{strcmp(funstr, SetAxes_ValidNames)};
nOut = nargout;

switch sprintf('%d',[nIn-1 nOut]>0)
    case '00'
        fh();
    case '10'
        fh(argsmain{2:end});
    case '01'
        [varargout{1:nOut}] = fh();
    otherwise
        [varargout{1:nOut}] = fh(argsmain{2:end});
end

%==========================================================================
    function xoffset()
        if ~isylabelcorner && get(hylab, 'Rotation')==0
            set(hylab, 'HorizontalAlignment', ylabhalign)
        end
        tightinset = get(ax, 'TightInset');
        set(ax, 'LooseInset', max(looseinset, tightinset))
        restoreleftright(hylab)
    end
%--------------------------------------------------------------------------
    function yoffset()
        tightinset = get(ax, 'TightInset');
        set(ax, 'LooseInset', max(looseinset, tightinset))
        restoreupdown(hxlab)
        restoreupdown(htitle)
    end
%--------------------------------------------------------------------------
    function axesmoveresize(offset)
        if nargin == 0
            disp('Drag a rectangle in the axes to define OuterPosition:')
            k = waitforbuttonpress;
            if k == 1
                return
            end
            % fig = gcf; if the figure is already in focus
            fig = get(ax, 'Parent');
            oldunits = get(fig, 'Units');
            set(fig, 'Units', 'normalized')
            outpos = rbbox;
            set(fig, 'Units', oldunits)
        else
            outpos = outpos + [1 1 -1 -1].*offset;
        end
        set(ax, 'OuterPosition', outpos)
        drawnow
        axpos = get(ax, 'Position');
    end
 %--------------------------------------------------------------------------
    function hxt = xtick2text(varargin)
        [xtickoffset, pvpairs] = extractoffset(.03, varargin{:});
        hxt = [];
        xticklab = cellstr(get(ax, 'XTickLabel'));
        xticklab = matchticks(xtick, xticklab);
        if ~isempty(xticklab)
            set(ax, 'XTickLabel', []);

            if xor(isxxbottom,isydirnormal)
                yloc = axlim(4) + xtickoffset*axheight;
            else
                yloc = axlim(3) - xtickoffset*axheight;
            end

            yloc = yloc(ones(size(xtick)));

            if isyylog
                yloc = 10.^yloc;
                xticklab = append10(xticklab, autoxticklab, pvpairs);
            end

            hxt = text(xtick, yloc, xticklab);

            set(hxt(strncmp('$', xticklab, 1)), 'Interpreter', 'latex')
%             mathtick = ~cellfun(@isempty, ...
%                 regexp(xticklab, '\${1,2}[^\$]*\${1,2}', 'match', 'once'));
%             set(hxt(strncmp(mathtick), 'Interpreter', 'latex')

            set(hxt, 'HorizontalAlignment', 'center', ...
                     'VerticalALignment', xlabvalign, pvpairs{:});
            %--------------------------------------------------------------
            oldunits = get(hxt(1), 'Units');
            set(hxt, 'Units', 'normalized')
            xtickext = get(hxt, 'Extent');
            xtickroom = max(cellfun(@(x) x(4), xtickext));
            set(hxt, 'Units', oldunits)

            xtickext = get(hxt, 'Extent');
            xlabposnow = get(hxlab, 'Position');

            if isxxbottom
                xlabposnow(2) = ...
                    min(cellfun(@(x) x(2), xtickext)) + xtickoffset;
            else
                xlabposnow(2) = ...
                    max(cellfun(@(x) x(2)+x(4), xtickext)) + xtickoffset;
            end

            set(hxlab, 'Position', xlabposnow)

        end
    end
%--------------------------------------------------------------------------
    function hyt = ytick2text(varargin)
        [ytickoffset, pvpairs] = extractoffset(.02, varargin{:});
        hyt = [];
        yticklab = cellstr(get(ax, 'YTickLabel'));
        yticklab = matchticks(ytick, yticklab);
        if ~isempty(yticklab)
            set(ax, 'YTickLabel', []);

            if xor(isyyleft,isxdirnormal)
                xloc = axlim(2) + ytickoffset*axwidth;
            else
                xloc = axlim(1) - ytickoffset*axwidth;
            end

            xloc = xloc(ones(size(ytick)));

            if isxxlog
                xloc = 10.^xloc;
                yticklab = append10(yticklab, autoyticklab, pvpairs);
            end

            hyt = text(xloc, ytick, yticklab);

            set(hyt(strncmp('$', yticklab, 1)), 'Interpreter', 'latex')
%             mathtick = ~cellfun(@isempty, ...
%                 regexp(yticklab, '\${1,2}[^\$]*\${1,2}', 'match', 'once'));
%             set(hyt(strncmp(mathtick), 'Interpreter', 'latex')

            set(hyt, 'VerticalAlignment', 'middle', ...
                     'HorizontalALignment', ylabhalign, pvpairs{:})
            %--------------------------------------------------------------
            oldunits = get(hyt(1), 'Units');
            set(hyt, 'Units', 'normalized')
            ytickext = get(hyt, 'Extent');
            ytickroom = max(cellfun(@(x) x(3), ytickext));
            set(hyt, 'Units', oldunits)

            ytickext = get(hyt, 'Extent');
            ylabposnow = get(hylab, 'Position');

            if isyyleft
                ylabposnow(1) = min(cellfun(@(x) x(1), ytickext));
            else
                ylabposnow(1) = max(cellfun(@(x) x(1)+x(3), ytickext));
            end

            set(hylab, 'Position', ylabposnow)

        end
    end
%--------------------------------------------------------------------------
    function xlabelcorner(varargin)
        if isxxlog
            xlabpos(1:2) = 10.^xlabpos(1:2);
        end

        set(hxlab, ...
            'Position', xlabpos, ...
            'VerticalALignment', 'middle', ...
            'HorizontalALignment', xlabhalign, varargin{:});

        isxlabelcorner = true;
        axpos = get(ax, 'Position');
        restoreleftright(hxlab)
    end
%--------------------------------------------------------------------------
    function ylabelcorner(varargin)
        if isyylog
            ylabpos(1:2) = 10.^ylabpos(1:2);
        end

        set(hylab, ...
            'Position', ylabpos, ...
            'Rotation', 0, ...
            'VerticalALignment', ylabvalign, ...
            'HorizontalALignment', 'center', varargin{:});

        isylabelcorner = true;
        axpos = get(ax, 'Position');
        restoreupdown(hylab)
    end
%--------------------------------------------------------------------------
    function h = axesarrows(varargin)
        args = varargin;
        xxarrow = true;
        yyarrow = true;

        if nargin>0
            switch args{1}
                case 'xx'
                    yyarrow = false;
                    args(1) = [];
                case 'yy'
                    xxarrow = false;
                    args(1) = [];
                otherwise
                    error('Unknown option')
            end
        end

       [left, bottom, right, top] = deal(0);

        if xxarrow && ~isxlabelcorner
            if isyyleft
                right = dxtfig;
            else
                left = dxtfig;
            end
        end

        if yyarrow && ~isylabelcorner
            if isxxbottom
                top = dytfig;
            else
                bottom = dytfig;
            end
        end

        axesmoveresize([left, bottom, right, top])

        if isyyleft
            xxarx = axpos(1)+axpos(3)+[0 dxtfig];
            yyarx = axpos([1 1]);
        else
            xxarx = axpos(1)+[axpos(3) -dxtfig];
            yyarx = axpos([1 1]) + axpos([3 3]);
        end

        if isxxbottom
            yyary = axpos(2)+axpos(4)+[0 dytfig];
            xxary = axpos([2 2]);
        else
            yyary = axpos(2)+[axpos(4) -dytfig];
            xxary = axpos([2 2]) + axpos([4 4]);
        end

        if xxarrow
            h(1) = annotation('arrow', xxarx, xxary, ...
                                      'Color', get(ax, 'XColor'));
        end

        if yyarrow
            h(2) = annotation('arrow', yyarx, yyary, ...
                                      'Color', get(ax, 'YColor'));
        end

        set(nonzeros(h), 'HeadLength', 5, 'HeadWidth', 5, ...
                         'LineWidth', get(ax, 'LineWidth'), args{:})

    end

%==========================================================================
    function restoreupdown(lab)
        withlab = ~isempty(get(lab, 'String')) && ...
                            strcmp(get(lab, 'Visible'), 'on');
        if withlab
            oldunits = get(lab, 'Units');
            set(lab, 'Units', 'normalized')
            labext = get(lab, 'Extent');

            if strcmpi(get(lab, 'Tag'), 'YLabel')
                labext = labext.*axpos([3 4 4 3]);
            else
                labext = labext.*axpos([3 4 3 4]);
            end

            extroom = labext(4);

            if strcmpi(get(lab, 'Tag'), 'XLabel')
                set(lab, 'Visible', 'off')
                keyword = xxloc;
%                 extroom = 1.05*max(extroom, xtickroom*axpos(4));
                extroom = max(extroom, xtickroom*axpos(4));
            else
                keyword = notxxloc;
            end

            tightinset = get(ax, 'TightInset');
            set(lab, 'Visible', 'on', 'Units', oldunits)
            if strcmp(keyword, 'bottom')
                looseinset(2) = max(looseinset(2), tightinset(2) + extroom);
            else
                %                 looseinset(4) = max(looseinset(4), labext(2)+labext(4)-axpos(4));
                looseinset(4) = max(looseinset(4), tightinset(4) + extroom);
            end

            set(ax, 'LooseInset', looseinset)
            axpos = get(ax, 'Position');
        end
    end
%--------------------------------------------------------------------------
    function restoreleftright(lab)
        withlab = ~isempty(get(lab, 'String')) && ...
                                strcmp(get(lab, 'Visible'), 'on');
        if withlab
            oldunits = get(lab, 'Units');
            set(lab, 'Units', 'normalized')
            labext = get(lab, 'Extent');

            if strcmpi(get(lab, 'Tag'), 'YLabel')
                set(lab, 'Visible', 'off')
                labext = labext.*axpos([3 4 4 3]);
                if get(lab, 'Rotation')==90
                    extroom = labext(4);
                else
                    extroom = labext(3);
                end
                extroom = max(extroom, ytickroom*axpos(3));
                keyword = yyloc;
            else
                labext = labext.*axpos([3 4 3 4]);
                extroom = labext(3);
                keyword = notyyloc;
            end

%             extroom = 1.05*extroom;

            tightinset = get(ax, 'TightInset');
            set(lab, 'Visible', 'on', 'Units', oldunits)

            if strcmp(keyword, 'left')
                looseinset(1) = max(looseinset(1), tightinset(1) + extroom);
            else
                %                 looseinset(3) = max(looseinset(3), labext(1)+labext(3)-axpos(3));
                looseinset(3) = max(looseinset(3), tightinset(3) + extroom);
            end

            set(ax, 'LooseInset', looseinset)
            axpos = get(ax, 'Position');
        end
    end
%--------------------------------------------------------------------------
    function [ axwidth, axheight, xlabhalign, ylabhalign, ...
                                  xlabvalign, ylabvalign ] = falign()

        if isxxlog
            axlim(1:2) = log10(axlim(1:2));
            xlabpos(1) = log10(xlabpos(1));
        end

        if isyylog
            axlim(3:4) = log10(axlim(3:4));
            ylabpos(2) = log10(ylabpos(2));
        end

        axwidth = diff(axlim(1:2));
        axheight = diff(axlim(3:4));

        dxtfig = .02;
        dxt = max(0.05*axwidth, dxtfig*axwidth/axpos(3));
        dxtfig = dxt*axpos(3)/axwidth;

        dytfig = .02;
        dyt = max(0.05*axheight, dytfig*axheight/axpos(4));
        dytfig = dyt*axpos(4)/axheight;

        if isyyleft
            xlabhalign = 'left';
            ylabhalign = 'right';
            if isxdirnormal
                xlabpos(1) = axlim(2) + 1.5*dxt;
                ylabpos(1) = axlim(1);
            else
                xlabpos(1) = axlim(1) + 1.5*dxt;
                ylabpos(1) = axlim(2);
            end
        else
            xlabhalign = 'right';
            ylabhalign = 'left';
            if isxdirnormal
                xlabpos(1) = axlim(1) - 1.5*dxt;
                ylabpos(1) = axlim(2);
            else
                xlabpos(1) = axlim(2) - 1.5*dxt;
                ylabpos(1) = axlim(1);
            end
        end

        if isxxbottom
            xlabvalign = 'top';
            ylabvalign = 'bottom';
            if isydirnormal
                ylabpos(2) = axlim(4)  + 1.5*dyt;
                xlabpos(2) = axlim(3);
            else
                ylabpos(2) = axlim(3)  + 1.5*dyt;
                xlabpos(2) = axlim(4);
            end
        else
            xlabvalign = 'bottom';
            ylabvalign = 'top';
            if isydirnormal
                ylabpos(2) = axlim(3)  - 1.5*dyt;
                xlabpos(2) = axlim(4);
            else
                ylabpos(2) = axlim(4)  - 1.5*dyt;
                xlabpos(2) = axlim(3);
            end
        end

    end
%--------------------------------------------------------------------------
    function [tickoffset, pvpairs] = extractoffset(tickoffset, varargin)
        pvpairs = varargin;
        npv = nargin - 1;
        if npv>0 && isnumeric(pvpairs{1})
            tickoffset = pvpairs{1};
            pvpairs(1) = [];
            npv = npv - 1;
        end

%         validateattributes(npv, {'numeric'}, {'even'})
%
%         assert((rem(npv,2) == 0), 'Incorrect number of input arguments')

        textpropnames = {'FontAngle', 'FontName', ...
            'FontUnits', 'FontSize', 'FontWeight'};

        textpv = [textpropnames; get(ax, textpropnames)];
        pvpairs = [textpv(:)' pvpairs];
    end
%--------------------------------------------------------------------------
    function ticklab = matchticks(tick, ticklab)
        ntick = numel(tick);
        nticklab = numel(ticklab);
        q = floor(ntick/nticklab);
        r = rem(ntick,nticklab);
        ticklab = ticklab([repmat(1:nticklab, 1, q) 1:r]);
    end
%--------------------------------------------------------------------------
    function ticklab = append10(ticklab, autoticklab, pvpairs)
        if autoticklab
            interpreter = get(0, 'DefaultTextInterpreter');
            if ~isempty(pvpairs)
                tparams = pvpairs(1:2:end-1);
                tvals = pvpairs(2:2:end);
                k = strcmpi('interpreter', tparams);
                if any(k), interpreter = tvals{k}; end
            end
            if strcmpi(interpreter, 'tex')
                d = '';
            else
                d = '$';
            end
            ticklab = cellfun(@(x) [d '10^{' x '}' d], ticklab, 'UniformOutput', false);
        end
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function [g, gdata, gprior] = mlpgrad_weighted(net, x, t, eso_w)
%MLPGRAD Evaluate gradient of error function for 2-layer network.
%
%	Description
%	G = MLPGRAD(NET, X, T) takes a network data structure NET  together
%	with a matrix X of input vectors and a matrix T of target vectors,
%	and evaluates the gradient G of the error function with respect to
%	the network weights. The error funcion corresponds to the choice of
%	output unit activation function. Each row of X corresponds to one
%	input vector and each row of T corresponds to one target vector.
%
%	[G, GDATA, GPRIOR] = MLPGRAD(NET, X, T) also returns separately  the
%	data and prior contributions to the gradient. In the case of multiple
%	groups in the prior, GPRIOR is a matrix with a row for each group and
%	a column for each weight parameter.
%
%	See also
%	MLP, MLPPAK, MLPUNPAK, MLPFWD, MLPERR, MLPBKP
%

%	Copyright (c) Ian T Nabney (1996-9)

% Check arguments for consistency
errstring = consist(net, 'mlp', x, t);
if ~isempty(errstring);
  error(errstring);
end
[y, z] = mlpfwd(net, x);
temp = y - t;
ndata = size(x, 1);
for m=1:ndata,
      delout(m,:)=eso_w(m,1)*temp(m,:);
end
clear temp;
gdata = mlpbkp(net, x, z, delout);

[g, gdata, gprior] = gbayes(net, gdata);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [x,y,info] = sedumi(A,b,c,K,pars)
%      [x,y,info] = sedumi(A,b,c,K,pars)
%
% SEDUMI  Self-Dual-Minimization/ Optimization over self-dual homogeneous
%         cones.
%
% >  X = SEDUMI(A,b,c) yields an optimal solution to the linear program
%      MINIMIZE c'*x SUCH THAT A*x = b, x >= 0
%      x is a vector of decision variables.
%      If size(A,2)==length(b), then it solves the linear program
%      MINIMIZE c'*x SUCH THAT A'*x = b, x >= 0
%
% >  [X,Y,INFO] = SEDUMI(A,b,c) also yields a vector of dual multipliers Y,
%      and a structure INFO, with the fields INFO.pinf, INFO.dinf and
%      INFO.numerr.
%
%    (1) INFO.pinf=INFO.dinf=0: x is an optimal solution (as above)
%      and y certifies optimality, viz.\ b'*y = c'*x and c - A'*y >= 0.
%      Stated otherwise, y is an optimal solution to
%      MAXIMIZE b'*y SUCH THAT c-A'*y >= 0.
%      If size(A,2)==length(b), then y solves the linear program
%      MAXIMIZE b'*y SUCH THAT c-A*y >= 0.
%
%    (2) INFO.pinf=1: there cannot be x>=0 with A*x=b, and this is certified
%      by y, viz. b'*y > 0 and A'*y <= 0. Thus y is a Farkas solution.
%
%    (3) INFO.dinf=1: there cannot be y such that c-A'*y >= 0, and this is
%      certified by x, viz. c'*x <0, A*x = 0, x >= 0. Thus x is a Farkas
%      solution.
%
%    (I)   INFO.numerr = 0: desired accuracy achieved (see PARS.eps).
%    (II)  INFO.numerr = 1: numerical problems warning. Results are accurate
%          merely to the level of PARS.bigeps.
%    (III) INFO.numerr = 2: complete failure due to numerical problems.
%
%    INFO.feasratio is the final value of the feasibility indicator. This
%    indicator converges to 1 for problems with a complementary solution, and
%    to -1 for strongly infeasible problems. If feasratio in somewhere in
%    between, the problem may be nasty (e.g. the optimum is not attained),
%    if the problem is NOT purely linear (see below). Otherwise, the reason
%    must lie in numerical problems: try to rescale the problem.
%
% >  [X,Y,INFO] = SEDUMI(A,b,0) or SEDUMI(A,b) solves the feasibility problem
%    FIND x>=0 such that A*x = b
%
% >  [X,Y,INFO] = SEDUMI(A,0,c) or SEDUMI(A,c) solves the feasibility problem
%    FIND y such that A'*y <= c
%
% >  [X,Y,INFO] = SEDUMI(A,b,c,K) instead of the constraint "x>=0", this
%      restricts x to a self-dual homogeneous cone that you describe in the
%      structure K. Up to 5 fields can be used, called K.f, K.l, K.q, K.r and
%      K.s, for Free, Linear, Quadratic, Rotated quadratic and Semi-definite.
%      In addition, there are fields K.xcomplex, K.scomplex and K.ycomplex
%      for complex-variables.
%
%    (1) K.f is the number of FREE, i.e. UNRESTRICTED primal components.
%      The dual components are restricted to be zero. E.g. if
%      K.f = 2 then x(1:2) is unrestricted, and z(1:2)=0.
%      These are ALWAYS the first components in x.
%
%    (2) K.l is the number of NONNEGATIVE components. E.g. if K.f=2, K.l=8
%      then x(3:10) >=0.
%
%    (3) K.q lists the dimensions of LORENTZ (quadratic, second-order cone)
%      constraints. E.g. if K.l=10 and K.q = [3 7] then
%          x(11) >= norm(x(12:13)),
%          x(14) >= norm(x(15:20)).
%      These components ALWAYS immediately follow the K.l nonnegative ones.
%      If the entries in A and/or c are COMPLEX, then the x-components in
%      "norm(x(#,#))" take complex-values, whenever that is beneficial.
%       Use K.ycomplex to impose constraints on the imaginary part of A*x.
%
%    (4) K.r lists the dimensions of Rotated LORENTZ
%      constraints. E.g. if K.l=10, K.q = [3 7] and K.r = [4 6], then
%          2*x(21)x(22) >= norm(x(23:24))^2,
%          2*x(25)x(26) >= norm(x(27:30))^2.
%      These components ALWAYS immediately follow the K.q ones.
%      Just as for the K.q-variables, the variables in "norm(x(#,#))" are
%      allowed to be complex, if you provide complex data. Use K.ycomplex
%      to impose constraints on the imaginary part of A*x.
%
%    (5) K.s lists the dimensions of POSITIVE SEMI-DEFINITE (PSD) constraints
%      E.g. if K.l=10, K.q = [3 7] and K.s = [4 3], then
%          mat( x(21:36),4 ) is PSD,
%          mat( x(37:45),3 ) is PSD.
%      These components are ALWAYS the last entries in x.
%
%    (a) K.xcomplex lists the components in f,l,q,r blocks that are allowed
%     to have nonzero imaginary part in the primal. For f,l blocks, these
%    (b) K.scomplex lists the PSD blocks that are Hermitian rather than
%      real symmetric.
%    (c) Use K.ycomplex to impose constraints on the imaginary part of A*x.
%
%    The dual multipliers y have analogous meaning as in the "x>=0" case,
%    except that instead of "c-A'*y>=0" resp. "-A'*y>=0", one should read that
%    c-A'*y resp. -A'*y are in the cone that is described by K.l, K.q and K.s.
%    In the above example, if z = c-A'*y and mat(z(21:36),4) is not symmetric/
%    Hermitian, then positive semi-definiteness reflects the symmetric/
%    Hermitian parts, i.e. Z + Z' is PSD.
%
%    If the model contains COMPLEX data, then you may provide a list
%    K.ycomplex, with the following meaning:
%      y(i) is complex if ismember(i,K.ycomplex)
%      y(i) is real otherwise
%    The equality constraints in the primal are then as follows:
%          A(i,:)*x = b(i)      if imag(b(i)) ~= 0 or ismember(i,K.ycomplex)
%          real(A(i,:)*x) = b(i)  otherwise.
%    Thus, equality constraints on both the real and imaginary part
%    of A(i,:)*x should be listed in the field K.ycomplex.
%
%    You may use EIGK(x,K) and EIGK(c-A'*y,K) to check that x and c-A'*y
%    are in the cone K.
%
% >  [X,Y,INFO] = SEDUMI(A,b,c,K,pars) allows you to override the default
%      parameter settings, using fields in the structure `pars'.
%
%    (1) pars.fid     By default, fid=1. If fid=0, then SeDuMi runs quietly,
%      i.e. no screen output. In general, output is written to a device or
%      file whose handle is fid. Use fopen to assign a fid to a file.
%
%    (2) pars.alg     By default, alg=2. If alg=0, then a first-order wide
%      region algorithm is used, not recommended. If alg=1, then SeDuMi uses
%      the centering-predictor-corrector algorithm with v-linearization.
%      If alg=2, then xz-linearization is used in the corrector, similar
%      to Mehrotra's algorithm. The wide-region centering-predictor-corrector
%      algorithm was proposed in Chapter 7 of
%        J.F. Sturm, Primal-Dual Interior Point Approach to Semidefinite Pro-
%        gramming, TIR 156, Thesis Publishers Amsterdam, 1997.
%
%    (3) pars.theta, pars.beta   By default, theta=0.25 and beta=0.5. These
%      are the wide region and neighborhood parameters. Valid choices are
%      0 < theta <= 1 and 0 < beta < 1. Setting theta=1 restricts the iterates
%      to follow the central path in an N_2(beta)-neighborhood.
%
%    (4) pars.stepdif, pars.w. By default, stepdif = 2 and w=[1 1].
%       This implements an adaptive heuristic to control ste-differentiation.
%       You can enable primal/dual step length differentiation by setting stepdif=1 or 0.
%      If so, it weights the rel. primal, dual and gap residuals as
%      w(1):w(2):1 in order to find the optimal step differentiation.
%
%    (5) pars.eps     The desired accuracy. Setting pars.eps=0 lets SeDuMi run
%      as long as it can make progress. By default eps=1e-8.
%
%    (6) pars.bigeps  In case the desired accuracy pars.eps cannot be achieved,
%     the solution is tagged as info.numerr=1 if it is accurate to pars.bigeps,
%     otherwise it yields info.numerr=2.
%
%    (7) pars.maxiter Maximum number of iterations, before termination.
%
%    (8) pars.stopat  SeDuMi enters debugging mode at the iterations specified in this vector.
%
%    (9) pars.cg      Various parameters for controling the Preconditioned conjugate
%     gradient method (CG), which is only used if results from Cholesky are inaccurate.
%    (a) cg.maxiter   Maximum number of CG-iterates (per solve). Theoretically needed
%          is |add|+2*|skip|, the number of added and skipped pivots in Cholesky.
%          (Default 49.)
%    (b) cg.restol    Terminates if residual is a "cg.restol" fraction of duality gap.
%          Should be smaller than 1 in order to make progress (default 5E-3).
%    (c) cg.refine    Number of refinement loops that are allowed. The maximum number
%          of actual CG-steps will thus be 1+(1+cg.refine)*cg.maxiter. (default 1)
%    (d) cg.stagtol  Terminates if relative function progress less than stagtol (5E-14).
%    (e) cg.qprec    Stores cg-iterates in quadruple precision if qprec=1 (default 0).
%
%    (10) pars.chol   Various parameters for controling the Cholesky solve.
%     Subfields of the structure pars.chol are:
%    (a) chol.canceltol: Rel. tolerance for detecting cancelation during Cholesky (1E-12)
%    (b) chol.maxu:   Adds to diagonal if max(abs(L(:,j))) > chol.maxu otherwise (5E5).
%    (c) chol.abstol: Skips pivots falling below abstol (1e-20).
%    (d) chol.maxuden: pivots in dense-column factorization so that these factors
%      satisfy max(abs(Lk)) <= maxuden (default 5E2).
%
%    (11) pars.vplot  If this field is 1, then SeDuMi produces a fancy
%      v-plot, for research purposes. Default: vplot = 0.
%
%    (12) pars.errors  If this field is 1 then SeDuMi outputs some error
%    measures as defined in the Seventh DIMACS Challenge. For more details
%    see the User Guide.
%
% Bug reports can be submitted at http://sedumi.mcmaster.ca.
%
% See also mat, vec, cellK, eyeK, eigK

% This file is part of SeDuMi 1.21 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2009 Imre Polik, Lehigh University, Bethlehem, PA, USA (since 1.21)
%
% Copyright (C) 2006 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

% J.F. Sturm, "Using SeDuMi 1.02, a MATLAB toolbox for optimization over
% symmetric cones," Optimization Methods and Software 11-12 (1999) 625-653.
% http://sedumi.mcmaster.ca

cputime0=cputime;
% ************************************************************
% INITIALIZATION
% ************************************************************
% ----------------------------------------
% Check input
% ----------------------------------------
if (nargin < 5)
    pars.fid = 1;
    if nargin < 3
        if nargin < 2
            error('Should have at least (A,b) or (A,c) arguments')
        end
        if length(b) == max(size(A))
            c = b; b = 0;           % given (A,c): LP feasibility problem
        else
            c = 0;                  % given (A,b): LP feasibility problem
        end
    end
    if isstruct(c)
        if nargin == 4
            pars = K;      % given (A,b,K,pars) or (A,c,K,pars)
        end
        K = c;          % given (A,b,K) or (A,c,K): cone feasibility problem
        if length(b) == max(size(A))
            c = b; b = 0;
        else
            c = 0;
        end
    elseif nargin < 4
        K.l = max(size(A));    % given (A,b,c): default to LP problem.
    end
end
% ----------------------------------------
% Bring data in (strict) internal SeDuMi form, check parameters,
% and print welcome.
% ----------------------------------------
[A,b,c,K,prep,origcoeff] = pretransfo(A,b,c,K,pars);
[N,m]=size(A);
if issparse(A)
    if sprank([A;b'])>sprank(A)
        info.pinf=1;
        x=[];
        %A dual improving direction could be computed easily
        y=[A;b']\[zeros(N,1);1];

        warning('SeDuMi:PrimalInfeasible','The problem is primal infeasible. there is no x such that Ax=b.')
        return
    elseif sprank(A)<m
        error('The coefficient matrix is not full row rank.')
    end
else
    if rank([A;b'])>rank(A)
        info.pinf=1;
        x=[];
        y=[A;b']\[zeros(N,1);1];
        warning('SeDuMi:PrimalInfeasible', 'The problem is primal infeasible, there is no x such that Ax=b.')
        return
    elseif rank(A)<m
        error('The coefficient matrix is not full row rank.')
    end
end
if prep.cpx.dim>0
    origcoeff=[];      % No error measures for complex problems.
end
lponly = (K.l == length(c));
pars = checkpars(pars,lponly);
% ----------------------------------------
% Print welcome and statistics of cone-problem
% ----------------------------------------
my_fprintf(pars.fid,'SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.\n');
switch pars.alg
    case 0
        my_fprintf(pars.fid,'Alg = 0: No corrector, ');
    case  1
        my_fprintf(pars.fid,'Alg = 1: v-corrector, ');
    case 2
        my_fprintf(pars.fid,'Alg = 2: xz-corrector, ');
end
switch pars.stepdif
    case 0
    case 1
        my_fprintf(pars.fid,'Step-Differentiation, ');
    case 2
        my_fprintf(pars.fid,'Adaptive Step-Differentiation, ');
end
my_fprintf(pars.fid,'theta = %5.3f, beta = %5.3f\n',pars.theta,pars.beta);
% --------------------------------------------------
% Print preprocessing information
% --------------------------------------------------
if pars.prep==1
    if isfield(prep,'sdp')
        blockcount=0;
        varcount=0;
        for sdpind=1:length(prep.sdp)
            if prep.sdp{sdpind}(1)==1
                blockcount=blockcount+1;
                varcount=varcount+prep.sdp{sdpind}(2);
            end
        end
        if blockcount>0
            my_fprintf(pars.fid,'Detected %i diagonal SDP block(s) with %i linear variables\n',blockcount,varcount);
        end
    end
    if isfield(prep,'freeblock1') && ~isempty(prep.freeblock1)
        my_fprintf(pars.fid,'Detected %i free variables in the linear part\n',length(prep.freeblock1));
    end
    if isfield(prep,'Kf') && prep.Kf>0
        switch pars.free
            case 0
                my_fprintf(pars.fid,'Split %i free variables\n',prep.Kf);
            case 1
                my_fprintf(pars.fid,'Put %i free variables in a quadratic cone\n',prep.Kf);
        end
    end
end
% --------------------------------------------------
% Remove dense columns (if any)
% --------------------------------------------------
Ablkjc = partitA(A,K.mainblks);
[dense,DAt.denq] = getdense(A,Ablkjc,K,pars);
if ~isempty(dense.cols)
    dense.A = A(dense.cols,:)';
    A(dense.cols,:) = 0.0;
    Ablkjc = partitA(A,K.mainblks);
else
    dense.A = sparse(length(b),0);
end
% ----------------------------------------
% Order constraints from sparse to dense, and find corresponding
% incremental nonzero pattern "Aord.dz" of At*dy in PSD part.
% ----------------------------------------
Aord.lqperm = sortnnz(A,[],Ablkjc(:,3));        % Sparse LP+Lorentz
DAt.q = findblks(A,Ablkjc,2,3,K.qblkstart);     % Lorentz ddotA-part
if ~isempty(DAt.q)
    DAt.q(dense.q,:) = 0.0;
    DAt.q = DAt.q + spones(extractA(A,Ablkjc,1,2,K.mainblks(1),K.mainblks(2)));
    Aord.qperm = sortnnz(DAt.q,[],[]);
else
    Aord.qperm = (1:length(b))';
end
[Aord.sperm, Aord.dz] = incorder(A,Ablkjc(:,3),K.mainblks(3));  % PSD
% ----------------------------------------
% Get nz-pattern of ADA.
% ----------------------------------------
ADA = getsymbada(A,Ablkjc,DAt,K.sblkstart);
% ----------------------------------------
% Ordering and symbolic factorization of ADA.
% ----------------------------------------
L = symbchol(ADA);
% --------------------------------------------------
% Symbolic fwsolve dense cols: L\[dense.A, dense.blkq],
% sparse ordering for dense column factorization
% --------------------------------------------------
symLden = symbcholden(L,dense,DAt);
% ----------------------------------------
% Initial solution
% ----------------------------------------
[d, v,vfrm, y,y0, R] = sdinit(A,b,c,dense,K,pars);
n = length(vfrm.lab);                         % order of K
merit = (sum(R.w) + max(R.sd,0))^2 * y0 / R.b0;  % Merit function
my_fprintf(pars.fid,'eqs m = %g, order n = %g, dim = %g, blocks = %g\n',...
    length(b),n,length(c),1 + length(K.q) + length(K.s));
my_fprintf(pars.fid,'nnz(A) = %d + %d, nnz(ADA) = %d, nnz(L) = %d\n',nnz(A),nnz(dense.A), nnz(ADA), nnz(L.L));
if ~isempty(dense.cols)
    my_fprintf(pars.fid,'Handling %d + %d dense columns.\n',...
        length(dense.cols),length(dense.q));
end
my_fprintf(pars.fid,' it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec\n');
my_fprintf(pars.fid,'  0 :            %8.2E %5.3f\n',merit,0);
% ----------------------------------------
% Initialize iterative statistics
% ----------------------------------------
STOP = 0;
err.maxb = 0.0;                 % Estimates the need to recompute residuals
iter = 0;
if pars.vplot == 1
    vlist = [vfrm.lab];
    ratelist = [];
end
wr.delta = 0.0;
wr.desc = 1;
%Seems unnecessary
%rate = 1.0;
feasratio = 0.0;
cputime1 = cputime;
% ************************************************************
% MAIN PREDICTOR-CORRECTOR LOOP
% ************************************************************
while STOP == 0
    iter = iter+1;
    if any(iter == pars.stopat)
        keyboard
    end

    if pars.stepdif==2 && ...
            (iter>20 || (iter>1 && (err.kcg + Lsd.kcg>3)) || ...
            (iter>5 && abs(1-feasratio)<0.05) )
        pars.stepdif=1;
    end
    % --------------------------------------------------
    % Compute ADA
    % --------------------------------------------------
    DAt = getDAtm(A, Ablkjc, dense, DAt.denq, d, K);
    ADA = getada1(ADA, A, Ablkjc(:,3), Aord.lqperm, d, K.qblkstart);
    ADA = getada2(ADA, DAt, Aord, K);
    [ADA,absd] = getada3(ADA, A, Ablkjc(:,3), Aord, invcholfac(d.u, K, d.perm), K);
    % ------------------------------------------------------------
    % Block Sparse Cholesky: ADA(L.perm,L.perm) = L.L*diag(L.d)*L.L'
    % ------------------------------------------------------------
    [L.L,L.d,L.skip,L.add] = blkchol(L,ADA,pars.chol,absd);
    % ------------------------------------------------------------
    % Factor dense columns
    % ------------------------------------------------------------
    [Lden, L.d] = deninfac(symLden, L,dense,DAt,d,absd, K.qblkstart,pars.chol);
    % ----------------------------------------
    % FACTORIZATION of self-dual embedding
    % ----------------------------------------
    Lsd = sdfactor(L,Lden, dense,DAt, d,v,y, A,c,K,R, y0,pars);
    % ------------------------------------------------------------
    % Compute and take IPM-step
    % from (v,y,v, y0) --> (xscl,y,zscl,y0)
    % ------------------------------------------------------------
    y0Old = y0;
    [xscl,yNxt,zscl,y0Nxt, w,relt, dxmdz,err, wr] = ...
        wregion(L,Lden,Lsd,...
        d,v,vfrm,A,DAt,dense, R,K,y,y0,b, pars, wr);
    % ------------------------------------------------------------
    % Evaluate the computed step.
    % ------------------------------------------------------------
    if y0Nxt > 0
        R.b = R.b + err.b / y0Nxt;
        R.sd = R.sd + err.g / y0Nxt;
        R.b0 = R.b0 + err.db0 / y0Nxt;
        y0 = y0Nxt;
    else
        R.b = (y0Nxt * R.b + err.b) / y0Old;      % In fact, we should have y0=0.
        R.sd = (y0Nxt * R.sd + err.g) / y0Old;
        R.b0 = (y0Nxt * R.b0 + err.db0) / y0Old;
        R.w(2) = abs(y0Nxt/y0Old) * R.w(2);        %=0: dual feasible
        R.c = (y0Nxt/y0Old) * R.c;
        R.maxRc = norm(R.c,inf);
        y0 = y0Old;
    end
    R.maxRb = norm(R.b,inf);                % Primal residual
    R.w(1) = 2 * pars.w(1) * R.maxRb / (1+R.maxb);
    meritOld = merit;
    merit = (sum(R.w) + max(R.sd,0))^2 * y0 / R.b0;
    rate = merit / meritOld;
    if (rate >= 0.9999) && (wr.desc == 1)
        % ------------------------------------------------------------
        % STOP = -1  --> Stop due to numerical problems
        % ------------------------------------------------------------
        STOP = -1;                  % insuf. progress in descent direction.
        iter = iter - 1;
        y0 = y0Old; %#ok
        my_fprintf(pars.fid,'Run into numerical problems.\n');
        break
    end
    feasratio = dxmdz(1) / v(1);            % (deltax0/x0) - (deltaz0/z0)
    % --------------------------------------------------
    % Primal-Dual transformation
    % --------------------------------------------------
    y = yNxt;
    by = full(sum(b.*y));
    [d,vfrm] = updtransfo(xscl,zscl,w,d,K);
    v = frameit(vfrm.lab,vfrm.q,vfrm.s,K);
    x0 = sqrt(d.l(1)) * v(1);
    % ----------------------------------------
    % SHOW ITERATION STATISTICS
    % ----------------------------------------
    my_fprintf(pars.fid,' %2.0f : %10.2E %8.2E %5.3f %6.4f %6.4f %6.4f %6.2f %2d %2d  ',...
        iter,by/x0,merit,wr.delta,rate,relt.p,relt.d,feasratio,err.kcg, Lsd.kcg);
    if pars.vplot == 1
        vlist = [vlist vfrm.lab/sqrt((R.b0*y0)/n)];
        ratelist = [ratelist rate];
    end
    % ----------------------------------------
    % If we get in superlinear region of LP,
    % try to guess optimal solution:
    % ----------------------------------------
    if lponly && (rate < 0.05)
        [xsol,ysol] = optstep(A,b,c, y0,y,d,v,dxmdz, ...
            K,L,symLden,dense, Ablkjc,Aord,ADA,DAt, feasratio, R,pars);
        if ~isempty(xsol)
            STOP = 2;                   % Means that we guessed right !!
            feasratio = 1 - 2*(xsol(1)==0);
            break
        end
    elseif (by > 0) && (abs(1+feasratio) < 0.05) && (R.b0*y0 < 0.5)
        if max(eigK(full(qreshape(Amul(A,dense,y,1),1,K)),K)) <= pars.eps * by
            STOP = 3;                   % Means Farkas solution found !
            break
        end
    end
    % --------------------------------------------------
    % OPTIMALITY CHECK: stop if y0*resid < eps * (x0+z0).
    % For feas. probs, we should divide the residual by x0, otherwise by z0.
    % Before stopping, recompute R.norm, since it may have changed due to
    % residual updates (the change should be small though).
    % --------------------------------------------------
    r0 = sum(R.w);
    cx = by + y0*R.sd - x0 / d.l(1);
    rgap = max(cx-by,0) / max([abs(cx),abs(by),1e-3 * x0]);
    precision1=y0*r0/(1+x0);
    precision2=(y0 * r0 + rgap)/x0;
    my_fprintf(pars.fid,'%1.1E\n',max(precision1,precision2));
    if precision1 < pars.eps       % P/D residuals small
        if precision2 < pars.eps    %Approx feasible and optimal
            STOP = 1;
            break
        elseif y0 * R.maxRb + x0 * R.maxb < -pars.eps * cx   % Approx Farkas
            STOP = 1;
            break
        elseif y0 * R.maxRc + x0 * R.maxc < pars.eps * by    % Approx Farkas
            STOP = 1;
            break;
        end
    end
    if iter >= pars.maxiter
        my_fprintf(pars.fid,'Maximum number of iterations reached.\n');
        STOP = -1;
    end
end % while STOP == 0.
my_fprintf(pars.fid,'\n');
clear ADA
nnzLadd=nnz(L.add);
nnzLskip=nnz(L.skip);
normLL=full(max(max(abs(L.L))));
clear L
% ************************************************************
% FINAL TASKS:
% ************************************************************
cputime2=cputime;
info.iter = iter;
info.feasratio = feasratio;
info.pinf = 0; info.dinf = 0;
info.numerr = 0;
% ------------------------------------------------------------
% Create x = D*v.
% ------------------------------------------------------------
if STOP == 2                % Exact optimal solution found (LP)
    x = xsol;
    y = ysol;
elseif STOP == 3            % Farkas solution y found (in early stage)
    x = zeros(length(c),1);
else
    x = [sqrt(d.l).*v(1:K.l); asmDxq(d,v,K); psdscale(d,v,K,1)];
end
% --------------------------------------------------
% Compute cx, Ax, etc.
% --------------------------------------------------
x0 = x(1);
cx = full(sum(c.*x));
abscx = sum(abs(c).*abs(x));
by = full(sum(b.*y));
Ax = Amul(A,dense,x,0);
Ay = full(Amul(A,dense,y,1));      % "full" since y may be scalar.
normy = norm(y);
normx = norm(x(2:end));
clear A

% ------------------------------------------------------------
% Determine infeasibility
% ------------------------------------------------------------

pinf = norm(x0*b-Ax);
z = qreshape(Ay-x0*c,1,K);
dinf = max(eigK(z,K));
if x0 > 0
    relinf = max(pinf / (1+R.maxb), dinf / (1+R.maxc)) / x0;
    % ------------------------------------------------------------
    % If infeasibility larger than epsilon, evaluate Farkas-infeasibility
    % ------------------------------------------------------------
    if relinf > pars.eps
        pdirinf = norm(Ax);
        ddirinf = max(eigK(qreshape(Ay,1,K),K));
        if cx < 0.0
            reldirinf = pdirinf / (-cx);
        else
            reldirinf = inf;
        end
        if by > 0.0
            reldirinf = min(reldirinf, ddirinf / by);
        end
        % ------------------------------------------------------------
        % If the quality of the Farkas solution is good and better than
        % the approx. feasible soln, set x0=0: Farkas solution found.
        % ------------------------------------------------------------
        if (reldirinf < pars.eps) || (relinf > max(pars.bigeps, reldirinf))
            x0 = 0.0;
            pinf = pdirinf;
            dinf = ddirinf;
        end
    end % relinf too large
end % x0 > 0
% ------------------------------------------------------------
% Interpret the solution as feasible:
% ------------------------------------------------------------
info.r0 = Inf;
if x0 > 0
    x = x / x0;
    y = y / x0;
    pinf = pinf /x0;
    dinf = dinf / x0;
    cx = cx/ x0;
    by = by / x0;
    normx = normx / x0;
    normy = normy / x0;
    if cx <= by                % zero or negative duality gap
        r0 = 0;
    elseif cx == 0.0           % Dual feasibility problem
        r0 = -by/(R.maxb*normy +1E-10 * x0);
    elseif by == 0.0           % Primal feasibility problem
        r0 = cx / (R.maxc*normx +1E-10 * x0);
    else                       % Optimization problem
        r0 = (cx-by)/(abs(by) + 1E-5 * (x0+abscx));
    end
    if r0 == 0,
        sigdig = Inf;
    else
        sigdig = -log10(r0);
    end
    my_fprintf(pars.fid,...
        'iter seconds digits       c*x               b*y\n');
    my_fprintf(pars.fid,'%3d %8.1f %5.1f %- 17.10e %- 17.10e\n',...
        iter,cputime2-cputime1,sigdig,cx,by);
    my_fprintf(pars.fid,'|Ax-b| = %9.1e, [Ay-c]_+ = %9.1E, |x|=%9.1e, |y|=%9.1e\n',...
        pinf,dinf,normx,normy);
    % ------------------------------------------------------------
    % Determine level of numerical problems with x0>0 (feasible)
    % ------------------------------------------------------------
    info.r0 = max([r0 ; pinf;dinf] ./ [1; 1+R.maxb+(1E-3)*R.maxRb;...
            1+R.maxc+(1E-3)*R.maxRc]);
    if STOP == -1
        if info.r0 > pars.bigeps
            my_fprintf(pars.fid, 'No sensible solution found.\n');
            info.numerr = 2;                          % serious numerical error
        elseif info.r0 > pars.eps
            info.numerr = 1;                          % moderate numerical error
        else
            info.numerr = 0;                          % achieved desired accuracy
        end
    else
        info.r0 = min( info.r0, pars.eps );
    end
else  % (if x0>0)
    % --------------------------------------------------
    % Infeasible problems: pinf==norm(Ax), dinf==max(eigK(At*y,K)).
    % --------------------------------------------------
    if pinf < -pars.bigeps * cx
        info.r0 = abs(pinf/cx);
        info.dinf = 1;
        abscx = -cx;
        pinf = pinf / abscx;
        normx = normx / abscx;
        x = x / abscx;
        my_fprintf(pars.fid, 'Dual infeasible, primal improving direction found.\n');
    end
    if dinf < pars.bigeps * by
        info.r0 = abs(dinf/by);
        info.pinf = 1;
        dinf = dinf / by;
        normy = normy / by;
        y = y / by;
        my_fprintf(pars.fid, 'Primal infeasible, dual improving direction found.\n');
    end
    my_fprintf(pars.fid,'iter seconds  |Ax|    [Ay]_+     |x|       |y|\n');
    my_fprintf(pars.fid,'%3d %8.1f %9.1e %9.1e %9.1e %9.1e\n',...
        iter,cputime2-cputime1,pinf,dinf,normx,normy);
    % --------------------------------------------------
    % Guess infeasible, but stopped due to numerical problems
    % --------------------------------------------------
    if info.pinf + info.dinf == 0
        my_fprintf(pars.fid, 'Failed: no sensible solution/direction found.\n');
        info.numerr = 2;
    elseif STOP == -1
        if (pinf > -pars.eps * cx) && (dinf > pars.eps * by)
            info.numerr = 1;
        else
            info.numerr = 0;
        end
    end
end
% ----------------------------------------
% - Bring xsol into the complex format of original (At,c),
% - transform q-variables into r-variables (Lorentz),
% - bring ysol into complex format, indicated by K.ycomplex.
% - at 0's in ysol where rows where removed.
% ----------------------------------------'
[x,y,K] = posttransfo(x,y,prep,K,pars); %#ok
% Detailed timing
%Preprocessing+IPM+Postprocessing
info.timing=[cputime1-cputime0 cputime2-cputime1 cputime-cputime2];
% Total time (for backward compatibility)
info.cpusec=sum(info.timing);
my_fprintf(pars.fid,'\nDetailed timing (sec)\n')
my_fprintf(pars.fid,'   Pre          IPM          Post\n')
my_fprintf(pars.fid,'%1.3E    ',info.timing)
my_fprintf(pars.fid,'\n')

% ----------------------------------------
% Make a fancy v-plot if desired
% ----------------------------------------
if pars.vplot == 1
    subplot(2,1,1)
    plot(0:iter,vlist,'o',[0 iter],[1 1],'b',...
        [0 iter],[pars.theta pars.theta],'g')
    title('Wide region v-plot')
    xlabel('iterations')
    ylabel('normalized v-values')
    subplot(2,1,2)
    plot(1:iter,ratelist)
    axis([0 iter 0 1])
    title('Reduction rates')
    xlabel('iterations')
    ylabel('reduction rate')
end
my_fprintf(pars.fid,'Max-norms: ||b||=%d, ||c|| = %d,\n',R.maxb,R.maxc);
my_fprintf(pars.fid,'Cholesky |add|=%d, |skip| = %d, ||L.L|| = %g.\n',...
    nnzLadd, nnzLskip, normLL);

% ----------------------------------------
% Compute error measures if needed
% ----------------------------------------
if ~isempty(origcoeff)
    %Reload the original coefficients
    s=(origcoeff.c)-(origcoeff.At)*sparse(y);       %To make s sparse
    cx=sum((origcoeff.c).*x);        %faster than c'*x
    by=sum((origcoeff.b).*y);
    xs=sum(x.*s);
    normb=norm(origcoeff.b,1);
    normc=norm(origcoeff.c,1);
    info.err=zeros(1,6);
    %     Error measures.
    %     Primal infeasibility
    info.err(1)=norm(x'*(origcoeff.At)-(origcoeff.b)',2)/(1+normb);
    %Let us get rid of the K.f part, since the free variables don't make
    %any difference in the cone infeasibility.
    %origcoeff.K.f=0;
    %     Primal cone infeasibility
    info.err(2)=max(0,-min(eigK(full(x(origcoeff.K.f+1:end)),origcoeff.K))/(1+normb));
    %     Dual infeasibility
    %info.err(3)=0.0; %s is not maintained explicitely
    %     Dual cone infeasibility
    info.err(4)=max(0,-min(eigK(full(s(origcoeff.K.f+1:end)),origcoeff.K))/(1+normc));
    %     Relative duality gap
    info.err(5)=(cx-by)/(1+abs(cx)+abs(by));
    %     Relative complementarity
    info.err(6)=xs/(1+abs(cx)+abs(by));

    my_fprintf(pars.fid,'\nDIMACS error measures\n')
    my_fprintf(pars.fid,'  PInf     PConInf     DInf    DConInf    RelGap    RelComp\n')
    my_fprintf(pars.fid,'%2.2E  ',info.err)
    my_fprintf(pars.fid,'\n')
end

repo: REST


text:
function [r, iclust, lam] = drawClusters(ops, r, mPix, mLam, Ly, Lx)

icell = size(mPix,2);
iclust = zeros(Ly, Lx);

lam      = zeros(Ly, Lx);
for i = icell:-1:1
    ipos = find(mPix(:,i)>0);
    ipix = lam(mPix(ipos,i))+1e-4 < mLam(ipos,i);
    iclust(mPix(ipos(ipix),i)) = i;
    lam(mPix(ipos(ipix),i))      = mLam(ipos(ipix),i);
end

% if nargin>5
%     lam = I;
% end

r = cat(1, r, rand(icell+1-numel(r), 1));

if ops.fig
%     figure(2)
    clf
    Sat = ones(Ly, Lx);
    H = zeros(Ly, Lx);

    H(iclust>0) = r(iclust(iclust>0));

    V = max(0, min(.75 * reshape(lam, Ly, Lx)/mean(lam(lam>1e-10)), 1));
    H = reshape(H, Ly, Lx);
    rgb_image = hsv2rgb(cat(3, H, Sat, V));
    imagesc(rgb_image)
    axis off
end

drawnow

author: cortex-lab


text:
function t = numel(obj)
% Number of simple file arrays involved.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

% Should be this, but it causes problems when accessing
% obj as a structure.
%t = prod(size(obj));

t  = numel(struct(obj));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ebsd = loadEBSD_Oxfordcsv(fname,varargin)

activePassive = {'passive','active'};

ebsd = EBSD;

try

  % read header
  fid = efopen(fname);
  s = fgetl(fid);
  assert(strcmpi(s,'Phase,ID'));

  % read phase names
  cs = {};
  s = fgetl(fid);
  hl = 0;
  while ~isempty(s)
    hl = hl + 1;
    phase = textscan(s,'%s%d','WhiteSpace',',');
    if phase{2}
      cs{hl} = crystalSymmetry('cubic','mineral',phase{1}{1});     %#ok<AGROW>
    else
      cs{hl} = phase{1}{1};     %#ok<AGROW>
    end
    s = fgetl(fid);
  end

  while ~strncmp(s,'Point',5)
    s = fgetl(fid);

    % extract active passive flag
    flag = sscanf(s,'Bunge Euler Convention, %f');
    if ~isempty(flag), activePassive = activePassive{1+flag}; end
    hl = hl + 1;
  end
  colNames = s;
  fclose(fid);

  % passive is default
  if iscell(activePassive), activePassive = activePassive{1}; end

  if check_option(varargin,'check'); return;end

  % extract column names
  colNames = textscan(colNames,'%s','delimiter',',');
  colNames = strrep(colNames{1},'Crystal ID','phase');

  % read data via generic interface
  ebsd = loadEBSD_generic(fname,'CS',cs,'header',hl+3,'delimiter',',',...
    'ColumnNames',colNames,'bunge',activePassive,varargin{:},'keepNaN');

catch
  interfaceError(fname)
end

author: mtex-toolbox

repo: mtex

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


text:
function percent = parfor_progress(N)
%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor.
%   PARFOR_PROGRESS works by creating a file called parfor_progress.txt in
%   your working directory, and then keeping track of the parfor loop's
%   progress within that file. This workaround is necessary because parfor
%   workers cannot communicate with one another so there is no simple way
%   to know which iterations have finished and which haven't.
%
%   PARFOR_PROGRESS(N) initializes the progress monitor for a set of N
%   upcoming calculations.
%
%   PARFOR_PROGRESS updates the progress inside your parfor loop and
%   displays an updated progress bar.
%
%   PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress
%   bar.
%
%   To suppress output from any of these functions, just ask for a return
%   variable from the function calls, like PERCENT = PARFOR_PROGRESS which
%   returns the percentage of completion.
%
%   Example:
%
%      N = 100;
%      parfor_progress(N);
%      parfor i=1:N
%         pause(rand); % Replace with real code
%         parfor_progress;
%      end
%      parfor_progress(0);
%
%   See also PARFOR.

% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/
% Copyright (c) 2011, Jeremy Scheff
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

error(nargchk(0, 1, nargin, 'struct'));

if nargin < 1
    N = -1;
end

percent = 0;
w = 50; % Width of progress bar

if N > 0
    f = fopen('parfor_progress.txt', 'w');
    if f<0
        error('Do you have write permissions for %s?', pwd);
    end
    fprintf(f, '%d\n', N); % Save N at the top of progress.txt
    fclose(f);

    if nargout == 0
        disp(['  0%[>', repmat(' ', 1, w), ']']);
    end
elseif N == 0
    delete('parfor_progress.txt');
    percent = 100;

    if nargout == 0
        disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']);
    end
else
    if ~exist('parfor_progress.txt', 'file')
        error('parfor_progress.txt not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize parfor_progress.txt.');
    end

    f = fopen('parfor_progress.txt', 'a');
    fprintf(f, '1\n');
    fclose(f);

    f = fopen('parfor_progress.txt', 'r');
    progress = fscanf(f, '%d');
    fclose(f);
    percent = (length(progress)-1)/progress(1)*100;

    if nargout == 0
        perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage
        disp([repmat(char(8), 1, (w+9)), char(10), perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']);
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% std_topo() - uses topoplot() to get the interpolated Cartesian grid of the
%               specified component topo maps. The topo map grids are saved
%               into a (.icatopo) file and a pointer to the file is stored
%               in the EEG structure. If such a file already exists,
%               loads the information from it.
%
%               Returns the topo map grids of all the requested components. Also
%               returns the EEG sub-structure etc (i.e EEG.etc), which is modified
%               with a pointer to the float file and some information about the file.
% Usage:
%               >> X = std_topo(EEG, components, option);
%
%                  % Returns the ICA topo map grid for a dataset.
%                  % Updates the EEG structure in the Matlab environment and re-saves
% Inputs:
%   EEG        - an EEG dataset structure.
%   components - [numeric vector] components in the EEG structure to compute topo maps
%                      {default|[] -> all}
%   option     - ['gradient'|'laplacian'|'none'] compute gradient or laplacian of
%                the scale topography. This does not acffect the saved file which is
%                always 'none' {default is 'none' = the interpolated topo map}
% Outputs:
%   X          - the topo map grid of the requested ICA components, each grid is
%                     one ROW of X.
%
% File output: [dataset_name].icatopo
%
% Authors:  Hilit Serby, Arnaud Delorme, SCCN, INC, UCSD, January, 2005
%
% See also  topoplot(), std_erp(), std_ersp(), std_spec(), std_preclust()

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

function [X] = std_topo(EEG, comps, option, varargin)

if nargin < 1
    help std_topo;
    return;
end;
if isfield(EEG,'icaweights')
   numc = size(EEG.icaweights,1);
else
   error('EEG.icaweights not found');
end
if nargin < 2
   comps = 1:numc;
elseif isempty(comps)
   comps = 1:numc;
end

if nargin < 3
    option = 'none';
end;

g = finputcheck( varargin, { 'recompute'   'string'   { 'on' 'off' }   'off' }, 'std_topo');
if isstr(g), error(g); end;

% figure; toporeplot(grid,'style', 'both','plotrad', 0.5, 'intrad', 0.5, 'xsurface' ,Xi, 'ysurface',Yi );

% Topo information found in dataset
% ---------------------------------
if exist(fullfile(EEG.filepath, [ EEG.filename(1:end-3) 'icatopo' ])) & strcmpi(g.recompute, 'off')
    for k = 1:length(comps)
        tmp = std_readtopo( EEG, 1, comps(k));
        if strcmpi(option, 'gradient')
            [tmpx, tmpy]  = gradient(tmp); %Gradient
            tmp = [tmpx(:); tmpy(:)]';
        elseif strcmpi(option, 'laplacian')
            tmp = del2(tmp); %Laplacian
            tmp = tmp(:)';
        else
            tmp = tmp(:)';
        end;

        tmp = tmp(find(~isnan(tmp)));
        if k == 1
            X = zeros(length(comps),length(tmp)) ;
        end
        X(k,:) =  tmp;
    end
    return
end

all_topos = [];
for k = 1:numc

    % compute topo map grid (topoimage)
    % ---------------------------------
    [hfig grid plotrad Xi Yi] = topoplot( EEG.icawinv(:,k), EEG.chanlocs(EEG.icachansind), ...
                                          'verbose', 'off',...
                                           'electrodes', 'on' ,'style','both',...
                                           'plotrad',0.55,'intrad',0.55,...
                                           'noplot', 'on', 'chaninfo', EEG.chaninfo);

    all_topos = setfield(all_topos, [ 'comp' int2str(k) '_grid' ], grid);
    all_topos = setfield(all_topos, [ 'comp' int2str(k) '_x' ]   , Xi(:,1));
    all_topos = setfield(all_topos, [ 'comp' int2str(k) '_y' ]   , Yi(:,1));

end

% Save topos in file
% ------------------
all_topos.datatype = 'TOPO';
tmpfile = fullfile( EEG.filepath, [ EEG.filename(1:end-3) 'icatopo' ]);
std_savedat(tmpfile, all_topos);

for k = 1:length(comps)
    tmp =  getfield(all_topos, [ 'comp' int2str(comps(k)) '_grid' ]);

    if strcmpi(option, 'gradient')
        [tmpx, tmpy]  = gradient(tmp); % Gradient
        tmp = [tmpx(:); tmpy(:)]';
    elseif strcmpi(option, 'laplacian')
        tmp = del2(tmp); % Laplacian
        tmp = tmp(:)';
    else
        tmp = tmp(:)';
    end;

    tmp = tmp(find(~isnan(tmp)));
    if k == 1
        X = zeros(length(comps),length(tmp)) ;
    end
    X(k,:) =  tmp;
end

author: PatternRecognition

repo: OpenBMI


text:
function this = read_matrix_from_workspace(this, inputMatrix)
% Reads in matrix from workspace, updates dimInfo according to data
% dimensions
%
%   Y = MrDataNd()
%   Y.read_matrix_from_workspace()
%
% This is a method of class MrDataNd.
%
% IN
%
% OUT
%
% EXAMPLE
%   read_matrix_from_workspace
%
%   See also MrDataNd

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-10-12
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% check whether valid dimInfo now

% TODO: update dimInfo, but keeping information that is unaltered by
% changing data dimensions...
% e.g. via dimInfo.merge
hasDimInfo = isa(this.dimInfo, 'MrDimInfo');

this.data = inputMatrix;

% remove singleton 2nd dimension kept by size command
nSamples = size(this.data);
if numel(nSamples) == 2
    nSamples(nSamples==1) = [];
end
resolutions = ones(1, numel(nSamples));

% set dimInfo or update according to actual number of samples
if ~hasDimInfo
    this.dimInfo = MrDimInfo('nSamples', nSamples, ...
        'resolutions', resolutions);
else
    if any(nSamples) % only update dimInfo, if any samples loaded
        if (numel(nSamples) ~= this.dimInfo.nDims)
            % only display the warning of an non-empty dimInfo (i.e. nDims
            % ~=0) has been given
            if (this.dimInfo.nDims ~=0)
                warning('Number of dimensions in dimInfo (%d) does not match dimensions in data (%d), resetting dimInfo', ...
                    this.dimInfo.nDims, numel(nSamples));
            end
            this.dimInfo = MrDimInfo('nSamples', nSamples, ...
                'resolutions', resolutions);
        elseif ~isequal(this.dimInfo.nSamples, nSamples)
            % if nSamples are correct already, leave it at that, otherwise:

            currentResolution = this.dimInfo.resolutions;

            isValidResolution = ~any(isnan(currentResolution)) || ...
                ~any(isinf(currentResolution)) && ...
                numel(currentResolution) == numel(nSamples);

            if isValidResolution
                % use update of nSamples to keep existing offset of samplingPoints
                this.dimInfo.nSamples = nSamples;
            else % update with default resolutions = 1
                this.dimInfo.set_dims(1:this.dimInfo.nDims, 'nSamples', ...
                    nSamples, 'resolutions', resolutions);
            end

        end
    end
end

author: translationalneuromodeling

repo: tapas


text:
function figureca(command,~)
    % This function is a callback for makemenus.
    % It should not be called from the command line.
    % figureca (previously internally known as FigureCallback)
    %
    % Keith Rogers 11/30/93

    report_this_filefun(mfilename('fullpath'));

    if (command == 1)
        if(strcmp(get(gcm2,'Label'),'Landscape'))
            set(gcf,'PaperOrientation','landscape');
            set(gcm2,'Label','Portrait');
        else
            set(gcf,'PaperOrientation','portrait');
            set(gcm2,'Label','Landscape');
        end
    end

author: CelsoReyes


text:
function [ num_int, pi ] = plane_imp_triangle_int_3d ( a, b, c, d, t )

%*****************************************************************************80
%
%% PLANE_IMP_TRIANGLE_INT_3D: intersection ( implicit plane, triangle ) in 3D.
%
%  Discussion:
%
%    The implicit form of a plane in 3D is:
%
%      A * X + B * Y + C * Z + D = 0
%
%    There may be 0, 1, 2 or 3 points of intersection returned.
%
%    If two intersection points are returned, then the entire line
%    between them comprises points of intersection.
%
%    If three intersection points are returned, then all points of
%    the triangle intersect the plane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(3,3), the vertices of the triangle.
%
%    Input, real A, B, C, D, the implicit plane parameters.
%
%    Output, integer NUM_INT, the number of intersection points returned.
%
%    Output, real PI(3,3), the intersection points.
%
  dim_num = 3;
  num_int = 0;
  pi = [];
%
%  Compute the signed distances between the vertices and the plane.
%
  dist1 = a * t(1,1) + b * t(2,1) + c * t(3,1) + d;
  dist2 = a * t(1,2) + b * t(2,2) + c * t(3,2) + d;
  dist3 = a * t(1,3) + b * t(2,3) + c * t(3,3) + d;
%
%  Consider any zero distances.
%
  if ( dist1 == 0.0 )
    num_int = num_int + 1;
    pi(1:dim_num,num_int) = t(1:dim_num,1);
  end

  if ( dist2 == 0.0 )
    num_int = num_int + 1;
    pi(1:dim_num,num_int) = t(1:dim_num,2);
  end

  if ( dist3 == 0.0 )
    num_int = num_int + 1;
    pi(1:dim_num,num_int) = t(1:dim_num,3);
  end
%
%  If 2 or 3 of the nodes intersect, we're already done.
%
  if ( 2 <= num_int )
    return
  end
%
%  If one node intersects, then we're done unless the other two
%  are of opposite signs.
%
  if ( num_int == 1 )

    if ( dist1 == 0.0 )

      [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,2), t(1:dim_num,3), ...
        dist2, dist3, num_int, pi );

    elseif ( dist2 == 0.0 )

      [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), t(1:dim_num,3), ...
        dist1, dist3, num_int, pi );

    elseif ( dist3 == 0.0 )

      [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), t(1:dim_num,2), ...
        dist1, dist2, num_int, pi );

    end

    return

  end
%
%  All nodal distances are nonzero, and there is at least one
%  positive and one negative.
%
  if ( dist1 * dist2 < 0.0 & dist1 * dist3 < 0.0 )

    [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), t(1:dim_num,2), ...
      dist1, dist2, num_int, pi );

    [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), t(1:dim_num,3), ...
      dist1, dist3, num_int, pi );

  elseif ( dist2 * dist1 < 0.0 & dist2 * dist3 < 0.0 )

    [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,2), t(1:dim_num,1), ...
      dist2, dist1, num_int, pi );

    [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,2), t(1:dim_num,3), ...
      dist2, dist3, num_int, pi );

  elseif ( dist3 * dist1 < 0.0 & dist3 * dist2 < 0.0 )

    [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,3), t(1:dim_num,1), ...
      dist3, dist1, num_int, pi );

    [ num_int, pi ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,3), t(1:dim_num,2), ...
      dist3, dist2, num_int, pi );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pd1d_steady_test ( )

%*****************************************************************************80
%
%% PD1D_STEADY_TEST runs all the peridynamics problems.
%
%  Discussion:
%
%    This function calls both convergence tests for all the problems.
%    You probably don't want to do it this way, since the tests take some
%    time, and generate graphic files that may get overwritten by the
%    next test.  However, this shows you how to use either test with any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 February 2012
%
%  Author:
%
%    Marta D'Elia
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PD1D_STEADY_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the PERIDYNAMICS_1D_STEADY library.\n' );
  fprintf ( 1, '  Run all the peridynamics test cases.\n' );
%
%  Test 1 decreases EPSILON and H in a fixed ratio.
%  Test 2 uses a fixed EPSILON and decreases H.
%
  pd1d_steady_test01 ( @problem1 );
  pd1d_steady_test02 ( @problem1 );

  pd1d_steady_test01 ( @problem2 );
  pd1d_steady_test02 ( @problem2 );

  pd1d_steady_test01 ( @problem3 );
  pd1d_steady_test02 ( @problem3 );

  pd1d_steady_test01 ( @problem4 );
  pd1d_steady_test02 ( @problem4 );

  pd1d_steady_test01 ( @problem5 );
  pd1d_steady_test02 ( @problem5 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PD1D_STEADY_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [res_map, scale_ind] = cal_scale(res_maps, scale_weights)

% this file selects the most confident scale
% Input:
%       res_maps      -  response maps
%       scale_weights - penalty weight of each scale
% Output:
%       res_map       - response map of the scale
%       scale_ind     - calculated scale index
% By Xin Li, 2017-12-2

% if num_scale == 1
%    res_map = res_maps;
%    scale_ind = 1;
%    return;
% end

num_scale = numel(scale_weights);
res_scale = bsxfun(@times, res_maps, reshape(scale_weights, 1,1, num_scale));

res = reshape(res_scale, [], num_scale);
res_max = max(res);

[~, scale_ind] = max(res_max);
res_map = res_maps(:,:, scale_ind);

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function cnn_save_model(epoch_dir, net_config_input, exp_info)

if ~isempty(net_config_input)
    net_config=gen_saved_net_config(net_config_input);
    net_filename=net_config.ref.cache_filename;
    net_config.ref.cache_filename=[];
    net_config_file=fullfile(epoch_dir, net_filename);
    fprintf('save net_config to: %s\n', net_config_file);
    try
        save(net_config_file, 'net_config', '-v7.3') ;
    catch err_info
        disp(err_info.getReport);
        dbstack;
        keyboard;
    end
end

if ~isempty(exp_info)
    info_filename=exp_info.cache_filename;
    info_file=fullfile(epoch_dir, info_filename);
    fprintf('save exp_info to: %s\n', info_file);
    try
        save(info_file, 'exp_info', '-v7.3') ;
    catch err_info
        disp(err_info.getReport);
        dbstack;
        keyboard;
    end
end

end

function net_config=gen_saved_net_config(net_config_input)

assert(isa(net_config_input, 'ref_obj'));

net_config=[];
make_ref_obj(net_config);
net_config.ref=net_config_input.ref;

for group_idx=1:length(net_config.ref.group_infos)
    group_info=net_config.ref.group_infos{group_idx};

    net_info=group_info.net_info;
    if ~isempty(net_info)
        new_net_info=[];
        make_ref_obj(new_net_info);
        new_net_info.ref=net_info.ref;
        do_clean_net_info(new_net_info);
        if new_net_info.ref.net_on_gpu
            my_move_net(new_net_info, 'cpu') ;
        end
        group_info.net_info=new_net_info;
    end

    if check_group_dag_net(group_info)

        fprintf('processing dag_net for saving, group_idx:%d, group_name: %s\n', group_info.group_idx, group_info.name);

        dag_net=group_info.dag_net;
        device_org=dag_net.device;
        saved_dag_net = dag_net.saveobj() ;

        if strcmp(device_org, 'gpu')
           dag_net.move('gpu');
        end

        group_info.dag_net=saved_dag_net;
    end

    net_config.ref.group_infos{group_idx}=group_info;
end

end

function do_clean_net_info(new_net_info)

new_net_info.ref.tmp_data=[];

for i=1:numel(new_net_info.ref.layers)

    if strcmp(new_net_info.ref.layers{i}.type,'conv')
          new_net_info.ref.layers{i}.filtersMomentum = [];
          new_net_info.ref.layers{i}.biasesMomentum = [];
    end

    if strcmp(new_net_info.ref.layers{i}.type,'my_custom')
        layer_update_fn=new_net_info.ref.layers{i}.layer_update_fn;
        if ~isempty(layer_update_fn)
            update_info=[];
            update_info.group_info=input_group_info;
            new_net_info.ref.layers{i}=new_net_info.ref.layers{i}.layer_update_fn(...
                'cache_clean', new_net_info.ref.layers{i}, update_info);
        end
    end
end

end

function use_dagnn=check_group_dag_net(group_info)

    use_dagnn=isfield(group_info, 'dag_net') && ~isempty(group_info.dag_net);
    if isfield(group_info, 'use_dagnn') && group_info.use_dagnn
        assert(use_dagnn==group_info.use_dagnn);
    end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function view = setBrightness(view,val);
%
% view = setBrightness(view,[val]);
%
% Set the brightness of a view, by
% changing the color maps in each
% view mode. Only those parts which
% are used for the anatomical underlay
% image (the first 1:numGrays entries)
% are brightened -- the overlays are
% unchanged.
%
% Val can range from 0 to 1. If omitted,
% it is read off of the view's brightness
% slider (so hopefully, in these circumstances,
% it has one).
%
%
% ras 01/05.
if notDefined('val')
    val = get(view.ui.brightness.sliderHandle,'Value');
end

setSlider(view,view.ui.brightness,val);

numGrays = view.ui.mapMode.numGrays;
cmap = gray(numGrays);

delta = 2*val - 1;
if delta ~= 0
    cmap = brighten(cmap,delta);
end

view.ui.anatMode.cmap = cmap;
view.ui.ampMode.cmap(1:numGrays,:) = cmap;
view.ui.phMode.cmap(1:numGrays,:) = cmap;
view.ui.coMode.cmap(1:numGrays,:) = cmap;
view.ui.mapMode.cmap(1:numGrays,:) = cmap;

view = refreshScreen(view, 1);

return


text:
% ANGLE_ERROR
% VARGPLVM

function [m,m_vec] = angle_error(Y1,Y2)

% usage [m,e] = angle_error(Y1,Y2);
% where Y1 and Y2 are matrices of same size
%
% e is a 54 vector corr to rms error in each angle across all data points
% m is the mean of e
%
% ERRORS GIVEN ARE IN DEGREES

[d,n] = size(Y1);

mat_error = acos(cos((pi/180)*(Y1-Y2)));
mat_error = (180/pi)*mat_error;
% cos and acos to handle wrap around effects in angle space

sq_mat_error = mat_error.*mat_error;

m_vec = zeros(d,1);

%m_vec = mean(mat_error,2);
%L1 norm
m_vec = sqrt(mean(sq_mat_error,2));

m = mean(m_vec);

author: SheffieldML

repo: vargplvm


text:
function error_frobenius = r8mat_is_eigen_left ( n, k, a, x, lambda )

%*****************************************************************************80
%
%% R8MAT_IS_EIGEN_LEFTT determines the error in a left eigensystem.
%
%  Discussion:
%
%    An R8MAT is a matrix of real values.
%
%    This routine computes the Frobenius norm of
%
%      X * A - LAMBDA * X
%
%    where
%
%      A is an N by N matrix,
%      X is an K by N matrix (each of K columns is an eigenvector)
%      LAMBDA is a K by K diagonal matrix of eigenvalues.
%
%    This routine assumes that A, X and LAMBDA are all real.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer K, the number of eigenvectors.
%    K is usually 1 or N.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real X(K,N), the K eigenvectors.
%
%    Input, real LAMBDA(K), the K eigenvalues.
%
%    Output, real ERROR_FROBENIUS, the Frobenius norm
%    of X*A-LAMBDA*X.
%
  c = x * a;

  for i = 1 : k
    c(i,1:n) = c(i,1:n) - lambda(i) * x(i,1:n);
  end

  error_frobenius = r8mat_norm_fro ( k, n, c );

  return
end


text:
function b = ne(f1,f2,varargin)
% ? sS1 ~= sS2

b = ~eq(f1,f2,varargin{:});


%LEARNALGOITML Wrapper class to the actual ITML code
classdef LearnAlgoITML < LearnAlgo

    properties
       p %parameters
       s %struct
       available
    end

    properties (Constant)
        type = 'itml'
    end

    methods
        function obj = LearnAlgoITML(p)
           if nargin < 1
              p = struct();
           end

           if ~isfield(p,'roccolor')
                p.roccolor = 'b';
           end
           obj.p  = p;
           check(obj);
        end

        function bool = check(obj)
           bool = exist('ItmlAlg') ~= 0;
           if ~bool
               fprintf('Sorry %s not available\n',obj.type);
           end
           obj.available = bool;
        end

        function s = learnPairwise(obj,X,idxa,idxb,matches)
            if ~obj.available
                s = struct();
                return;
            end

            tic;
            s.M = PairMetricLearning(@ItmlAlg, idxa', idxb', matches, X');
            s.t = toc;
            s.learnAlgo = obj;
            s.roccolor = obj.p.roccolor;
        end

        function s = learn(obj,X,y)
            if ~obj.available
                s = struct();
                return;
            end

            tic;
            s.M = MetricLearning(@ItmlAlg, y', X');
            s.t = toc;
            s.learnAlgo = obj;
            s.roccolor = obj.p.roccolor;
        end

        function d = dist(obj, s, X, idxa,idxb)
            d = cdistM(s.M,X,idxa,idxb);
        end
    end
end


text:
function [countV,ktrAll,veAll,kepAll,vpAll] = runDCETofts(planC)
% runDCETofts.m
% Function to generate and write out DCE parameter maps, compute and save statistics to file.
%
% AI 5/26/16
% AI 7/06/16  Updated to provide option to scale AIF using reference region time course.
%
% Note:
% Required parameters are read from user-input .txt file.
% See ReadDCEParameterFile.m for list of valid input parameters.

%% --- Get model parameters --------

indexS = planC{end};

% Read input parameter file
[inputS,filePath] = ReadDCEParameterFile();
outPath = fileparts(filePath);
[paramS,shiftS] = getInputs(inputS);

% Get DCE scan no.
scanListC = unique({planC{indexS.scan}.scanType});
if numel(scanListC)>1
    [series,ok] = listdlg('ListString',scanListC,'Name','Select series','ListSize',[300,250],...
        'PromptString','Select DCE series','SelectionMode','Single');
    if ~ok
        return
    end
    scanNumV = find(strcmp({planC{indexS.scan}.scanType},scanListC{series}));
else
    scanNumV = 1:length(planC{indexS.scan});
end
scanS = planC{indexS.scan}(scanNumV);

%Extract relevant parameters from DICOM header
headerS.manufacturer = scanS(1).scanInfo(1).DICOMHeaders.Manufacturer;
headerS.xSize = scanS(1).scanInfo(1).DICOMHeaders.Rows;
headerS.ySize = scanS(1).scanInfo(1).DICOMHeaders.Columns;
headerS.RepetitionTime = scanS(1).scanInfo(1).DICOMHeaders.RepetitionTime;
tDel = getTDel(planC);
headerS.tDel = double(tDel)/1000 ; % convert time step to seconds from ms
headerS.MagneticFieldStrength = scanS(1).scanInfo(1).DICOMHeaders.MagneticFieldStrength;
headerS.nSlices = size(scanS(1).scanArray,3);
FA = scanS(1).scanInfo(1).DICOMHeaders.FlipAngle;
headerS.FlipAngle = shiftS.FAF*FA; % correct flip angle by multiplicative factor FAF

switch(headerS.manufacturer)
    case 'GE MEDICAL SYSTEMS'
        totalImages = scanS(1).scanInfo(1).DICOMHeaders.ImagesinAcquisition;
        headerS.frames = totalImages/headerS.nSlices;
    case 'Philips Healthcare'
        headerS.frames = scanS(1).scanInfo(1).DICOMHeaders.NumberofTemporalPositions;
end
%Display:
fprintf(['\nVendor: %s\nNo. slices: %d\nNo. frames: %d\nTR(ms): %g\nDel_T(s): %g\n',...
    'Flip angle: %g\nImage Size(rows x cols): %d x %d\nB0: %g\n'],headerS.manufacturer,...
    headerS.nSlices,headerS.frames,headerS.RepetitionTime,headerS.tDel,...
    headerS.FlipAngle,headerS.xSize,headerS.ySize,headerS.MagneticFieldStrength);

%Time shift for ROI time course (user-input)
ROIPointShift = shiftS.TROIPointShift * strcmp(paramS.model,'T') + ...
    shiftS.ETROIPointShift * strcmp(paramS.model,'ET');
shiftS.ROIPointShift = ROIPointShift;

% r1  (Using relaxivities from Pintaske, et. al. Invest Radiol 2006) **
%     Gd-DTPA, 1.5T   =  3.9 L/(mmol x s)
%     Gd-DTPA, 3.0T   =  3.3 L/(mmol x s)
%     Gd-BOPTA (multihance) 1.5T  = 8.1 L/(mmol x s)
%     Gd-BOPTA (multihance) 3.0T = 6.3 L/(mmol x s)
r1 = 3.9*(headerS.MagneticFieldStrength < 2.0)*strcmp(paramS.cAgent,'gd')+ ...
     3.3*(headerS.MagneticFieldStrength > 2.0)*strcmp(paramS.cAgent,'gd') + ...
     8.1*(headerS.MagneticFieldStrength < 2.0)*strcmp(paramS.cAgent,'mh') + ...
     6.3*(headerS.MagneticFieldStrength > 2.0)*strcmp(paramS.cAgent,'mh');

%Vector of sampling times
ind = double(0:(headerS.frames-1));
tsec = ind*headerS.tDel;
tmin = tsec./60;         % ktrans, etc. are in /min.  Parker coefs are in  min or /min.
TR = headerS.RepetitionTime/1000;  % TR from DICOM header is in ms. Change to s.

% Parameter array for the concentrarion calculation:
%    pCalcConc(1) = R1 = relaxivity of contrast agent (units are mmolxs)
%    pCalcConc(2) = TR (s)
%    pCalcConc(3) = flip angle (degrees)
%    pCalcConc(4) = pre-contrast T1 for tissue (s)
pCalcConc = [r1, TR, headerS.FlipAngle, paramS.T10];

%Generate AIF
%-------------temp:--------
coefFile = fullfile(getCERRPath,'PlanAnalysis','DCE-MR analysis','sample_param_and_coeff_files','Pros_Iliac_17pts_Parker_coefs.txt');
%---------------------------
fprintf('\n Generating AIF from Parker Coefficient file...');
AIFP = genAIF(tsec,coefFile,headerS.frames,tmin,paramS,shiftS);
fprintf('\n AIF generation complete.\n');

%If using T1 maps, load T1 maps for all slices into a 3D array
if ~(paramS.T1Map==0)
    T13M = getScanArray(paramS.T1Map,planC);
else
    T13M = [];
end

%% --------- Get ROI mask ---------------
%Identify structures associated with selected scan
structNumV = 1:length(planC{indexS.structures});
assocScansV = getStructureAssociatedScan(structNumV,planC);
structNumV = structNumV(ismember(assocScansV,scanNumV));
structS =  planC{indexS.structures}(structNumV);
emptyStructV = cellfun(@isempty,{structS.rasterSegments}); %Discard empty structures
structS =  structS(~emptyStructV);

%Extract mask
ROImask3M = true(size(scanS(1).scanArray));
structListC = cellfun(@lower,{structS.structureName},'un',0);
%User-selection of structure:
if numel(structListC)>1
    [strSel,ok] = listdlg('ListString',structListC,'Name','Select structure','ListSize',[300,250],...
        'PromptString','Select ROI','SelectionMode','Single');
    if ~ok
        return
    end
else
    strSel = 1;
end
strNum = structNumV(strSel);
rasterSegments = getRasterSegments(strNum,planC);
[temp,maskSlicesV] = rasterToMask(rasterSegments,1, planC);
ROImask3M(:,:,maskSlicesV) = temp;

%% ---- Concatenate DCE images ----------
sliceC = cell(1,maskSlicesV);
for k = 1:numel(maskSlicesV)
    timeSlices = arrayfun(@(x) x.scanArray(:,:,maskSlicesV(k)),scanS,'un',0);
    sliceC{k} = cat(3,timeSlices{:});
end
sliceTimeCourse4M = cat(4,sliceC{:});
%sliceTimeCourse4M = uint16(sliceTimeCourse4M);   %To match Kristen's code

%% ----- Scale AIF using reference region time course if reqd ----
% Get reference region (muscle) mask
useRefRegion = questdlg('Use muscle reference to scale AIF?','Reference region','Yes','No','No') ;
if strcmp(useRefRegion, 'No')
    AIFPScaled = AIFP;
else
    structListC = {structS.structureName};
    [strNum,ok] = listdlg('ListString',structListC,'SelectionMode','Single','Name','Select reference structure');
    if ~ok
        return
    end
    muscleRasterSegS = getRasterSegments(strNum,planC);
    [muscleMaskM,muscleROISlice] = rasterToMask(muscleRasterSegS,1, planC); %Assumed single-slice
    % Scale AIF using reference region time course
    muscleSliceC = arrayfun(@(x) x.scanArray(:,:,muscleROISlice),scanS,'un',0);
    muscleSliceTimeCourse3M  = cat(3,muscleSliceC{:});
    muscleTimeCourse3M = ROIMaskedSlices(muscleSliceTimeCourse3M,muscleMaskM);
    AIFPScaled = refRegionScale(AIFP,inputS,shiftS,filePath,headerS,muscleTimeCourse3M,muscleMaskM);
end

%% ---- Get parameter maps ----------
getVals = 'Yes';
nRuns = 1;
while(strcmp(getVals,'Yes'))
    if nRuns~=1                   %Runs with value from parameter file if checkShift='n' the first time
        shiftS.checkShift = 'y';  %Set checkShift to 'y' if user chooses to run again with different shift
        [countV,ktrAll,veAll,kepAll,vpAll,dceRsqAll,UserInShift,UserInBase] = getDCEToftsParams(paramS,shiftS,headerS.xSize,headerS.ySize,maskSlicesV,sliceTimeCourse4M,pCalcConc,...
            ROImask3M,T13M,tmin,AIFPScaled,outFolder);
    else
        [countV,ktrAll,veAll,kepAll,vpAll,dceRsqAll,UserInShift,UserInBase,outFolder] = getDCEToftsParams(paramS,shiftS,headerS,maskSlicesV,sliceTimeCourse4M,pCalcConc,...
            ROImask3M,T13M,tmin,AIFPScaled,outPath);
    end
    getVals = questdlg(sprintf('Current shift = %d.\n Run again with different shift?',UserInShift),...
        'Change shift?','Yes','No','No');
    nRuns = nRuns + 1;
end

%% ---Write user-input shifts to param file-------
fid = fopen(filePath);
inputC = textscan(fid, '%s %s %n','endofline','\r\n');
fclose(fid);
fieldsC = fieldnames(paramS);
if strcmp(paramS.model,'T') %Tofts
    shiftIdx = strcmp(fieldsC,'TROIPointShift');
else  %Extended tofts
    shiftIdx = strcmp(fieldsC,'ETROIPointShift');
end
baseIdx = strcmp(fieldsC,'basepts');
inputC{3}(shiftIdx) = UserInShift;
inputC{3}(baseIdx) = UserInBase;
fid = fopen(filePath,'w+');
fmt = '\r\n%s\t%s\t%d';
for lineNum = 1:size(inputC{1},1)
    col1 = inputC{1}(lineNum);
    col2 = inputC{2}(lineNum);
    col3 = inputC{3}(lineNum);
    fprintf(fid,fmt,col1{1},col2{1},col3);
end
fclose(fid);

%Write statistics
writeStats(outFolder,outfileBase,paramS.model,paramS.cutOff,countV,ktrAll,veAll,kepAll,vpAll,dceRsqAll);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function daefun(dh,x,z,u,p)

  conf = dh.userdata;

  ddp = - 1/conf.m * z.lambda*x.p - [0;9.81] + [u.F;0];

  dh.setODE('p',x.v);
  dh.setODE('v',ddp);
  dh.setODE('time', 1);

  % The algebraic equation constraints the pendulum's mass to be on a circular
  % path if the initial conditions are satisfied.
  dh.setAlgEquation(dot(ddp,x.p)+x.v(1)^2+x.v(2)^2);
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [image] = LoadImage(directory, timestamp, LUT)

% LoadImage - load a rectified image from disk
%
% [image] = LoadImage(directory, timestamp, LUT)
%
% eg.
% timestamps = dlmread('<dataset_root>/stereo.timestamps');
% [ ~, ~, ~, ~, ~, LUT] = ...
%     ReadCameraModel('<models_dir>/stereo_wide_left_undistortion.bin');
% image = LoadImage('<dataset_root>/stereo/left', timestamps(100,1), LUT);
%
% INPUTS:
%   directory: directory containing images named <timestamp>.png
%   timestamp: timestamp of image to load
%   LUT (optional): lookup table for image rectification, as returned from
%     ReadCameraModel. If not supplied, original distorted image will be
%     returned.
%     See ReadCameraModel and UndistortImage
%
% OUTPUTS:
%   image: image at the given timestamp. If an undistortion lookup table is
%     supplied, image will be undistorted.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2016 University of Oxford
% Authors:
%  Geoff Pascoe (gmp@robots.ox.ac.uk)
%  Will Maddern (wm@robots.ox.ac.uk)
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if directory(end) ~= '/'
    directory = [directory '/'];
  end

  path = [directory num2str(timestamp) '.png'];
  if ~exist(path, 'file')
    image = false;
    return;
  end

  if regexp(directory, 'stereo')
    bayer_pattern = 'gbrg';
  else
    bayer_pattern = 'rggb';
  end

  image = demosaic(imread(path), bayer_pattern);

  if exist('LUT', 'var')
    image = UndistortImage(image, LUT);
  end

end

author: ori-mrg

repo: robotcar-dataset-sdk


text:
function [selectedSlices, nrows, ncols] = montageDims(vw);
%
% [selectedSlices, nrows, ncols] = montageDims(vw);
%
% Return the dimensions of a montage view (for inplanes, flat
% levels). selectedSlices is a vector of the slices currently
% being displayed in the montage, nrows and ncols are the rows
% and columns in the image. Takes into account an optional
% 'columnMajor' setting, so montages can run across rows
% or columns (though the interface to set this isn't written
% yet).
%
% ras, 05/05
viewType = viewGet(vw,'viewType');
ui = viewGet(vw,'ui');
switch viewType
    case 'Inplane',
        firstSlice = viewGet(vw, 'curSlice');
        nSlices = get(ui.montageSize.sliderHandle,'Value');
        selectedSlices = firstSlice:firstSlice+nSlices-1;
    case 'Flat',
        selectedSlices = getFlatLevelSlices(vw);
        nSlices = length(selectedSlices);
    otherwise,
        error('drawROIsMontage: no support for this vw type.');
end

selectedSlices = selectedSlices(selectedSlices <= viewGet(vw, 'numSlices'));
nSlices = length(selectedSlices);

if isfield(vw.ui,'settings') && isfield(vw.ui.settings,'columnMajor') ...
    && vw.ui.settings.columnMajor==1
    % montage is column-major
	ncols = ceil(sqrt(nSlices));
	nrows = ceil(nSlices/ncols);
else
    % montage is row-major
	nrows = ceil(sqrt(length(selectedSlices)));
	ncols = ceil(nSlices/nrows);
end

return


text:
function [V]=ellipseCoord3(e,t)

V=[e.radii(1).*cos(t(:)) e.radii(2).*sin(t(:)) zeros(numel(t),1)];
V=(e.axes*V')';
V=V+e.centre(ones(numel(t),1),:);

%%
% _*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