text:
function o2 = canlab_results_fmridisplay(input_activation, varargin)
% :Usage:
% ::
%
%    canlab_results_fmridisplay(input_activation, [optional inputs])
%
% purpose:  This function display fmri results.
%
% :Input:
%
%   **input_activation:**
%        nii, img,
%
%        This image has the blobs you want to
%        display. You can also enter a cl "clusters" structure or
%        "region" object.
%
%        you can also get a thresholded image like the examples used here
%        from a number of places - by thresholding your results in SPM
%        and using "write filtered" to save the image, by creating masks
%        from meta-analysis or anatomical atlases, or by using
%        mediation_brain_results, robust_results_threshold,
%        robust_results_batch_script, threshold_imgs, or object
%        oriented tools including fmri_data and statistic_image objects.
%
% :Optional Inputs:
%
%   **'noblobs':**
%        do not display blobs
%
%   **'outline':**
%        display blob outlines
%
%   **'nooutline':**
%        do not display blob outlines (default)
%
%   **'addmontages':**
%        when entering existing fmridisplay obj, add new montages
%
%   **'noremove':**
%        do not remove current blobs when adding new ones
%
%   **'nofigure':**
%        do not create a new figure (for selected montage sets only)
%
%   **'outlinecolor:**
%        followed by new outline color
%
%   **'splitcolor':**
%        followed by 4-cell new split colormap colors (help fmridisplay or edit code for defaults as example)
%
%   **'montagetype':**
%
%        'full'            Axial, coronal, and saggital slices, 4 cortical surfaces
%        'compact'         Midline saggital and two rows of axial slices [the default]
%        'compact2'        A single row showing midline saggital and axial slices
%        'multirow'        A series of 'compact2' displays in one figure for comparing different images/maps side by side
%        'regioncenters'   A series of separate axes, each focused on one region
%
%        'full' for full montages of axial and sagg slices.
%
%        'full hcp' for full montage, but with surfaces and volumes from
%        HCP data
%
%        'compact' [default] for single-figure parasagittal and axials slices.
%
%        'compact2': like 'compact', but fewer axial slices.
%
%        'multirow': followed by number of rows
%           e.g., o2 = canlab_results_fmridisplay([], 'multirow', 2);
%
%        {'blobcenters', 'regioncenters'}: Slices for the center of each blob/region
%        Note: this creates a new figure, tagged
%        'fmridisplay_regioncenters', and is not compatible with 'nofigure'
%
%   **'noverbose':**
%        suppress verbose output, good for scripts/publish to html, etc.
%
%   **'overlay':**
%        specify anatomical image for montage (not surfaces), followed by
%        image name
%        e.g., o2 = canlab_results_fmridisplay([], 'overlay', 'icbm152_2009_symmetric_for_underlay.img')';
%
%         The default brain for overlays is based on Keuken et al. 2014
%         For legacy SPM8 single subject, enter as arguments:
%         'overlay', which('SPM8_colin27T1_seg.img')
%
% Other inputs to addblobs (fmridisplay method) are allowed, e.g., 'cmaprange', [-2 2], 'trans'
%
% See help fmridisplay
% e.g., 'color', [1 0 0]
%
% You can also input an existing fmridisplay object, and it will use the
% one you have created rather than setting up the canonical slices.
%
% Try "brighten(.4) to make the images brighter.
%
% :Example Script:
% ::
%
%    input_activation = 'Pick_Atlas_PAL_large.nii';
%
%    % set up the anatomical overlay and display blobs
%    % (see the code of this function and help fmridisplay for more examples)
%
%    o2 = canlab_results_fmridisplay(input_activation);
%
%    %% ========== remove those blobs and change the color ==========
%
%    cl = mask2clusters(input_activation);
%    removeblobs(o2);
%    o2 = addblobs(o2, cl, 'color', [0 0 1]);
%
%    %% ========== OR
%
%    r = region(input_activation);
%    o2 = removeblobs(o2);
%    o2 = addblobs(o2, r, 'color', [1 0 0]);
%
%    %% ========== Create empty fmridisplay object on which to add blobs:
%    o2 = canlab_results_fmridisplay
%    o2 = canlab_results_fmridisplay([], 'compact2', 'noverbose');
%
%    %% ========== If you want to start over with a new fmridisplay object,
%    % make sure to clear o2, because it uses lots of memory
%
%    % This image should be on your path in the "canlab_canonical_brains" subfolder:
%
%    input_activation = 'pain-emotion_2s_z_val_FDR_05.img';
%    clear o2
%    close all
%    o2 = canlab_results_fmridisplay(input_activation);
%
%    %% ========== save PNGs of your images to insert into powerpoint, etc.
%    % for your paper/presentation
%
%    scn_export_papersetup(400);
%    saveas(gcf, 'results_images/pain_meta_fmridisplay_example_sagittal.png');
%
%    scn_export_papersetup(350);
%    saveas(gcf, 'results_images/pain_meta_fmridisplay_example_sagittal.png');
%
%    Change colors, removing old blobs and replacing with new ones:
%    o2 = canlab_results_fmridisplay(d, o2, 'cmaprange', [.3 .45], 'splitcolor', {[0 0 1] [.3 0 .8] [.9 0 .5] [1 1 0]}, 'outlinecolor', [.5 0 .5]);
%
% ..
%    Tor Wager
%    1/27/2012
% ..

if ~which('fmridisplay.m')
    disp('fmridisplay is not on path.  it is in canlab tools, which must be on your path!')
    return
end

if nargin == 0
    o2 = canlab_results_fmridisplay(region(), 'noblobs', 'nooutline');
    return
end

if ischar(input_activation)

    if strcmp(input_activation, 'compact') || strcmp(input_activation, 'compact2') || strcmp(input_activation, 'full') ...
            || strcmp(input_activation, 'multirow') || strcmp(input_activation, 'coronal') || strcmp(input_activation, 'sagittal')

        % Entered no data map; intention is not to plot blobs, just create underlay
        varargin{end + 1} = 'noblobs';
        varargin{end + 1} = 'nooutline';
        varargin{end + 1} = input_activation; % so code below finds it later
        % do nothing else for now - this is not an input image

    else
        % assume it is an input image

        cl = region(fmri_data(input_activation));  % mask2clusters(input_activation);

    end

elseif isstruct(input_activation) || isa(input_activation, 'region')
    cl = input_activation;
    if ~isa(input_activation, 'region'), cl = cluster2region(cl); end

elseif isa(input_activation, 'image_vector')
    cl = region(input_activation);

elseif isempty(input_activation)
    % do nothing for now

else
    error('I don''t recognize the format of input_activation.  It should be a thresholded mask, clusters, or region object');
end

% process input arguments
% --------------------------------------------
doblobs = true;
dooutline = false;
doaddmontages = false;
doremove = true;
outlinecolor = [0 0 0];
splitcolor = {[0 0 1] [0 .8 .8] [1 .4 .5] [1 1 0]}; % {[0 0 1] [.3 .6 .9] [.8 .3 0] [1 1 0]};  % more straight orange to yellow: {[0 0 1] [0 1 1] [1 .5 0] [1 1 0]}
montagetype = 'compact';
doverbose = true;
%overlay='SPM8_colin27T1_seg.img';
overlay = 'keuken_2014_enhanced_for_underlay.img';
dofigure = true;

wh = strcmp(varargin, 'overlay');
if any(wh), wh = find(wh); overlay = varargin{wh(1) + 1};  varargin([wh wh+1]) = []; end

wh = strcmp(varargin, 'noblobs');
if any(wh), doblobs = false; varargin(wh) = []; end

wh = strcmp(varargin, 'nooutline');
if any(wh), dooutline = false; varargin(wh) = []; end

wh = strcmp(varargin, 'outline');
if any(wh), dooutline = true; varargin(wh) = []; end

wh = strcmp(varargin, 'addmontages');
if any(wh), doaddmontages = true; varargin(wh) = []; end

wh = strcmp(varargin, 'outlinecolor');
if any(wh), wh = find(wh); outlinecolor = varargin{wh(1) + 1}; end

wh = strcmp(varargin, 'splitcolor');
if any(wh), wh = find(wh); splitcolor = varargin{wh(1) + 1}; end

wh = strcmp(varargin, 'noremove');
if any(wh), doremove = false; varargin(wh) = []; end

wh = strcmp(varargin, 'full');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'multirow');
if any(wh), montagetype = varargin{find(wh)};
    nrows = varargin{find(wh) + 1};
    varargin{find(wh) + 1} = [];
        varargin(wh) = [];
end

wh = strcmp(varargin, 'blobcenters');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'regioncenters');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'full hcp');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'compact');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'compact2');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'coronal');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'saggital');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

wh = strcmp(varargin, 'allslices');
if any(wh), montagetype = varargin{find(wh)}; varargin(wh) = []; end

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

wh = strcmp(varargin, 'nofigure');
if any(wh), dofigure = false; varargin(wh) = []; end

wh = false(1, length(varargin));
for i = 1:length(varargin)
    wh(i) = isa(varargin{i}, 'fmridisplay');
    if wh(i), o2 = varargin{wh}; end
end
varargin(wh) = [];

xyz = [-20 -10 -6 -2 0 2 6 10 20]';
xyz(:, 2:3) = 0;

if isempty(input_activation)
    % we will skip the blobs, but process other optional input args
    doblobs = false;
    dooutline = false;
end

if ~exist('o2', 'var')

    % set up fmridisplay
    % --------------------------------------------
    % you only need to do this once
    % then you can add montages, add and remove blobs, add and remove points (for
    % meta-analysis), etc.

    if doverbose

        disp('Setting up fmridisplay objects');
        % disp('This takes a lot of memory, and can hang if you have too little.');

    end

    [opath, ofname, oext] = fileparts(overlay);
    if isempty(opath)
        % check for file in matlab path
        o2 = fmridisplay('overlay', which(overlay));
    else
        % complete path specified
        o2 = fmridisplay('overlay',overlay);
    end

    % You can customize these and run them from the command line

    switch montagetype

        case {'blobcenters', 'regioncenters'}
            %Make a series of montages at center of each region and add blobs to that:

            if ~exist('cl', 'var'), error('Must enter region object to use blobcenters option.'); end

            xyz = cat(1, cl.mm_center);

            %           onerowstr = [];
            %           if length(cl) < 20, onerowstr = 'onerow'; end

            orientation = 'axial';

            % Make a grid - determine subplots
            nr = floor(sqrt(length(cl)));
            nc = ceil(length(cl) ./ nr);

            [~, axh] = create_figure('fmridisplay_regioncenters', nr, nc, false, true);

            set(axh,'Visible','off'); % turn off axis grid for all axes

            for i = 1:length(cl)

                %axh(i) = subplot(nr, nc, i);
                %axes(axh(i))

                if i == 1
                    [o2, dat] = montage(o2, orientation, 'wh_slice', xyz(i, :), 'onerow', 'existing_axes', axh(i), 'existing_figure', 'noverbose');
                else
                    o2 = montage(o2, 'volume_data', dat, orientation, 'wh_slice', xyz(i, :), 'onerow', 'existing_axes', axh(i), 'existing_figure', 'noverbose');
                end

            end

            wh_montages = 1:length(cl);

        case 'compact'
            % The default

            % saggital
            axh1 = axes('Position', [-0.02 0.4 .17 .17]);
            [o2, dat] = montage(o2, 'saggital', 'wh_slice', [-4 0 0], 'onerow', 'noverbose', 'existing_axes', axh1);
            text(50, -50, 'left');
            drawnow

            axh2 = axes('Position', [-0.02 0.6 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'saggital', 'wh_slice', [4 0 0], 'onerow', 'noverbose', 'existing_axes', axh2);
            text(50, -50, 'right');
            drawnow

            % sagg center
            axh3 = axes('Position', [.08 0.5 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'saggital', 'wh_slice', [0 0 0], 'onerow', 'noverbose', 'existing_axes', axh3);
            drawnow;
            o2.montage{3}.slice_mm_coords;

            % axial bottom
            axh4 = axes('Position', [.10 0.4 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 12, 'noverbose', 'existing_axes', axh4);

            % axial top
            axh5 = axes('Position', [.12 0.54 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-46 50], 'onerow', 'spacing', 12, 'noverbose', 'existing_axes', axh5);

            wh_montages = [1 2 4 5];

            % Lines
            axes(o2.montage{3}.axis_handles)
            locs = [o2.montage{4}.slice_mm_coords; o2.montage{5}.slice_mm_coords];
            for i = 1:length(locs)
                %draw_vertical_line(locs(i));
                hh(i) = plot( [-105 65], [locs(i) locs(i)], 'b', 'LineWidth', 1);
            end

            % brighten(.5)
            sz = get(0, 'screensize');
            set(gcf, 'Color', 'w', 'Position', [sz(3).*.1 sz(4).*.9 sz(3).*.6 sz(4).*.6]);

        case 'multirow'

            % Notes: for some reason, at least in Matlab 2017a, when you
            % use the existing figure the slices scale with the fig
            % position in size. When using montage to create a new figure,
            % they don't seem to do this. could be enlarge_axes, or ???

            if dofigure
                slices_fig_h = figure('Color', 'w');
            end

            ss = get(0, 'ScreenSize');
            myheightdivisor = 1.5; % 3/nrows;  % controls figure aspect ratio
            set(gcf, 'Position', [round(ss(3)/20) round(ss(4)*.5) round(ss(3)*.9) round(ss(4)/myheightdivisor) ])

            %shiftvals = [0:.17:nrows]; % more than we need, but works
            %shiftvals = [0:.24:nrows]; % more than we need, but works
            shiftvals = repmat([0:.24:.75], 1, ceil(nrows/4));  % repeat positions every 4, for adding new figures

            for i = 1:nrows

                % Can only put 4 on one figure, so create additional
                % figures as needed
                if i > 4 && rem(i, 4) == 1

                    figure('Color', 'w');
                    set(gcf, 'Position', [round(i*ss(3)/20) round(ss(4)*.5) round(ss(3)*.9) round(ss(4)/myheightdivisor) ])

                end

                % saggital
                axh = axes('Position', [-0.02 .75-shiftvals(i) .17 .17]);  % [-0.02 0.15+shiftvals(i) .17 .17]);
                axh(2) = axes('Position', [.022 .854-shiftvals(i) .17 .17]);

                if i == 1
                    [o2, dat] = montage(o2, 'saggital', 'slice_range', [-2 2], 'spacing', 4, 'onerow', 'noverbose', 'existing_axes', axh);

                else
                    o2 = montage(o2, 'volume_data', dat, 'saggital', 'slice_range', [-2 2], 'spacing', 4, 'onerow', 'noverbose', 'existing_axes', axh);

                end
                drawnow

                % axial
                axh = axes('Position', [.022 0.8-shiftvals(i) .17 .17]);   % [.015 0.15+shiftvals(i) .17 .17]);
                o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-32 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);
                drawnow

            end

        case 'full'
            % saggital
            [o2, dat] = montage(o2, 'saggital', 'wh_slice', xyz, 'onerow', 'noverbose');
            shift_axes(-0.02, -0.04);

            % coronal
            axh = axes('Position', [-0.02 0.37 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'coronal', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            % axial
            axh = axes('Position', [-0.02 0.19 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            axh = axes('Position', [-0.02 0.01 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-44 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            allaxh = findobj(gcf, 'Type', 'axes');
            disp(length(allaxh));
            for i = 1:(length(allaxh)-36)
                pos1 = get(allaxh(i), 'Position');
                pos1(1) = pos1(1) - 0.03;
                set(allaxh(i), 'Position', pos1);
            end

            % surface
            o2 = surface(o2, 'axes', [0.1 0.74 .25 .25], 'direction', 'hires left', 'orientation', 'medial');
            o2 = surface(o2, 'axes', [0.3 0.74 .25 .25], 'direction', 'hires right', 'orientation', 'medial');
            o2 = surface(o2, 'axes', [0.5 0.74 .25 .25], 'direction', 'hires left', 'orientation', 'lateral');
            o2 = surface(o2, 'axes', [0.7 0.74 .25 .25], 'direction', 'hires right', 'orientation', 'lateral');

            wh_montages = [1 2 3 4];
            wh_surfaces = [1 2 3 4];

        case 'compact2'  % creates a new figure
            %subplot(2, 1, 1);
            [o2, dat] = montage(o2, 'axial', 'slice_range', [-32 50], 'onerow', 'spacing', 8, 'noverbose');

            % shift all axes down and right
            allaxh = o2.montage{1}.axis_handles;
            for i = 1:length(allaxh)
                pos1 = get(allaxh(i), 'Position');
                pos1(2) = pos1(2) - 0.08;
                pos1(1) = pos1(1) + 0.03;

                % enlarge a bit
                pos1(3:4) = pos1(3:4) + .02;

                set(allaxh(i), 'Position', pos1);
            end

            enlarge_axes(gcf, 1);
            %axh = axes('Position', [0.0 0.08 .15 1]);
            axh = axes('Position', [-0.02 .75-.3 .17 .17]);  % [-0.02 0.15+shiftvals(i) .17 .17]);
            axh(2) = axes('Position', [.022 .854-.3 .17 .17]);

            %o2 = montage(o2, 'saggital', 'wh_slice', [0 0 0], 'existing_axes', axh, 'noverbose');
            o2 = montage(o2, 'volume_data', dat, 'saggital', 'slice_range', [-2 2], 'spacing', 4, 'onerow', 'noverbose', 'existing_axes', axh);

            %ss = get(0, 'ScreenSize');
            %set(gcf, 'Position', [round(ss(3)/12) round(ss(4)*.9) round(ss(3)*.8) round(ss(4)/5.5) ]) % this line messes p the
            %images, makes it too big an overlapping
            wh_montages = [1 2];

            brighten(.4)

        case 'coronal'
            % coronal
            o2 = montage(o2, 'coronal', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose');
             wh_montages = 1;

        case  'saggital'
            o2 = montage(o2, 'saggital', 'wh_slice', xyz, 'onerow', 'noverbose');

             wh_montages = 1;

        case 'allslices'

            [o2, dat] = montage(o2, 'saggital', 'wh_slice', xyz, 'onerow', 'noverbose');
            shift_axes(-0.02, -0.04);

            axh = axes('Position', [-0.02 0.37 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'coronal', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            % axial
            axh = axes('Position', [-0.02 0.19 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            wh_montages = [1 2 3];

        case 'full hcp'
            % saggital
            [o2, dat] = montage(o2, 'saggital', 'wh_slice', xyz, 'onerow', 'noverbose');
            shift_axes(-0.02, -0.04);

            % coronal
            axh = axes('Position', [-0.02 0.37 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'coronal', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            % axial
            axh = axes('Position', [-0.02 0.19 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            axh = axes('Position', [-0.02 0.01 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-44 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            allaxh = findobj(gcf, 'Type', 'axes');
            disp(length(allaxh));
            for i = 1:(length(allaxh)-36)
                pos1 = get(allaxh(i), 'Position');
                pos1(1) = pos1(1) - 0.03;
                set(allaxh(i), 'Position', pos1);
            end

            % surface
            o2 = surface(o2, 'axes', [0.1 0.74 .25 .25], 'direction', 'surface left', 'orientation', 'medial');
            o2 = surface(o2, 'axes', [0.3 0.74 .25 .25], 'direction', 'surface right', 'orientation', 'medial');
            o2 = surface(o2, 'axes', [0.5 0.74 .25 .25], 'direction', 'surface left', 'orientation', 'lateral');
            o2 = surface(o2, 'axes', [0.7 0.74 .25 .25], 'direction', 'surface right', 'orientation', 'lateral');

            wh_montages = [1 2 3 4];
            wh_surfaces = [1:8];

        otherwise error('illegal montage type. choose full or compact.');
    end

    % wh_montages = [1 2];

else
    if doverbose, disp('Using existing fmridisplay object'); end

    % Other inputs will be passed into addblobs
    existingmons = length(o2.montage);
    existingsurs = length(o2.surface);

    if doaddmontages
        % use same o2, but add montages
        switch montagetype
            case 'full'
            % saggital
            [o2, dat] = montage(o2, 'saggital', 'wh_slice', xyz, 'onerow', 'noverbose');
            shift_axes(-0.02, -0.04);

            % coronal
            axh = axes('Position', [-0.02 0.37 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'coronal', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            % axial
            axh = axes('Position', [-0.02 0.19 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            axh = axes('Position', [-0.02 0.01 .17 .17]);
            o2 = montage(o2, 'volume_data', dat, 'axial', 'slice_range', [-44 50], 'onerow', 'spacing', 8, 'noverbose', 'existing_axes', axh);

            allaxh = findobj(gcf, 'Type', 'axes');
            disp(length(allaxh));
            for i = 1:(length(allaxh)-36)
                pos1 = get(allaxh(i), 'Position');
                pos1(1) = pos1(1) - 0.03;
                set(allaxh(i), 'Position', pos1);
            end

            % surface
            o2 = surface(o2, 'axes', [0.1 0.74 .25 .25], 'direction', 'hires left', 'orientation', 'medial');
            o2 = surface(o2, 'axes', [0.3 0.74 .25 .25], 'direction', 'hires right', 'orientation', 'medial');
            o2 = surface(o2, 'axes', [0.5 0.74 .25 .25], 'direction', 'hires left', 'orientation', 'lateral');
            o2 = surface(o2, 'axes', [0.7 0.74 .25 .25], 'direction', 'hires right', 'orientation', 'lateral');

            wh_surfaces = existingsurs + [1 2 3 4];

            case 'compact'
                [o2, dat] = montage(o2, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 6, 'noverbose');
                axh = axes('Position', [0.05 0.4 .1 .5]);
                o2 = montage(o2, 'volume_data', dat, 'saggital', 'wh_slice', [0 0 0], 'existing_axes', axh, 'noverbose');

            case 'compact2'
                [o2, dat] = montage(o2, 'axial', 'slice_range', [-40 50], 'onerow', 'spacing', 8, 'noverbose');
                enlarge_axes(gcf, 1);
                axh = axes('Position', [-0.03 0.15 .2 1]);
                o2 = montage(o2, 'volume_data', dat, 'saggital', 'wh_slice', [0 0 0], 'existing_axes', axh, 'noverbose');

                % shift all axes down and right
                shift_axes(+0.03, -0.10);

                %ss = get(0, 'ScreenSize');
                %set(gcf, 'Position', [round(ss(3)/12) round(ss(4)*.9) round(ss(3)*.8) round(ss(4)/5.5) ])

            otherwise error('illegal montage type. choose full or compact when adding to existing montage set.')
        end

        wh_montages = existingmons + [1 2];

    else
        if doremove
            o2 = removeblobs(o2);
        end

        wh_montages = 1:existingmons;
        wh_surfaces = 1:existingsurs;

    end

end

% Now we can add blobs
% --------------------------------------------

% they are added to all montages by default, but you can specify selected
% montages if you want to as well.

% it's easy to remove them as well:
% o2 = removeblobs(o2);

if doblobs
    if exist('wh_surfaces', 'var')
        o2 = addblobs(o2, cl, 'splitcolor', splitcolor, 'wh_montages', wh_montages, 'wh_surfaces', wh_surfaces, 'nolegend', varargin{:});
    else
        o2 = addblobs(o2, cl, 'splitcolor', splitcolor, 'wh_montages', wh_montages, varargin{:});
    end
end

if dooutline
    o2 = addblobs(o2, cl, 'color', outlinecolor, 'outline', 'wh_montages', wh_montages, 'no_surface');
end

% ------------------------------------------------------
% INLINE FUNCTIONS
% ------------------------------------------------------

    function shift_axes(x_offset, y_offset)

        % usage:   function shift_axes(-0.03, -0.18);
        % purpose: To shift axes according of the habdles.
        %
        % input:   Function parameters are (x, y) offset relative to current position

        % shift all axes according to offset values
        allaxh = findobj(gcf, 'Type', 'axes');
        for i = 1:length(allaxh)
            pos1 = get(allaxh(i), 'Position');
            pos1(2) = pos1(2) + y_offset;
            pos1(1) = pos1(1) + x_offset;
            set(allaxh(i), 'Position', pos1);
        end
    end % shift_axes

end  % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function a = forsythe ( alpha, beta, n )

%*****************************************************************************80
%
%% FORSYTHE returns the FORSYTHE matrix.
%
%  Discussion:
%
%    The Forsythe matrix represents a Jordan canonical matrix, perturbed
%    by a rank one update.
%
%  Formula:
%
%    If ( I = J )
%      A(I,J) = BETA
%    else if ( J = I+1 )
%      A(I,J) = 1
%    else if ( I = N and J = 1 ) then
%      A(I,J) = ALPHA
%    else
%      A(I,J) = 0
%
%  Example:
%
%    ALPHA = 2, BETA = 3, N = 5
%
%    3 1 0 0 0
%    0 3 1 0 0
%    0 0 3 1 0
%    0 0 0 3 1
%    2 0 0 0 3
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    The characteristic equation of A is
%
%      ( BETA - LAMBDA )^N - (-1)^N*ALPHA = 0
%
%    The eigenvalues of A are
%
%      LAMBDA(I) = BETA
%        + abs ( ALPHA )^1/N * exp ( 2 * I * PI * sqrt ( - 1 ) / N )
%
%    Gregory and Karney consider the special case where BETA is 0,
%    and ALPHA is a "small" value.  In that case, the characteristic
%    equation is LAMBDA^N - ALPHA = 0, and the eigenvalues are the
%    N-th root of ALPHA times the N roots of unity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 5.22,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 103,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, real ALPHA, BETA, define the matrix.  A typical
%    value of ALPHA is the square root of the machine precision; a typical
%    value of BETA is 0.0.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( j == i )
        a(i,j) = beta;
      elseif ( j == i + 1 )
        a(i,j) = 1.0;
      elseif ( i == n & j == 1 )
        a(i,j) = alpha;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = make_sin_out(T,Sdim,parm)
% Sinusoidal signal

if isfield(parm,'Tstart'),
	Tstart = parm.Tstart;
else
	Tstart = rand(1);
end

% period
Tmin = parm.Tmin;
Tmax = parm.Tmax;

if Sdim==1,
	Tperiod = Tmax;
elseif isfield(parm,'Tperiod')
	Tperiod = parm.Tperiod;
else
	Tstp = (Tmax - Tmin)/(Sdim-1);
	Tperiod = (0:(Sdim-1))*Tstp + Tmin;
	Tperiod = circshift(Tperiod, fix(Sdim/2));
%	Tperiod = ((Sdim-1):-1:0)*Tstp + Tmin;
%	Tperiod = Tmax*rand(1,Sdim) + Tmin;
end

% frequency
omega  = 2*pi./Tperiod(:);

% Sinusoidal signal
t = (1:T) + Tmax*Tstart;
%t = repadd( omega * t, - (2*pi) * rand(Sdim,1) );
x = sin( omega * t );

author: KamitaniLab

repo: GenericObjectDecoding


text:
% function h = plotclr(x,y,v,marker,vlim)
% plots the values of v colour coded
% at the positions specified by x and y.
% A colourbar is added on the right side of the figure.
%
% The colourbar strectches from the minimum value of v to its
% maximum.
%
% 'marker' is optional to define the marker being used. The
% default is a point. To use a different marker (such as circles, ...) send
% its symbol to the function (which must be enclosed in '; see example).
%
% 'vlim' is optional, to define the limits of the colourbar.
% v values outside vlim are not plotted
%
% modified by Stephanie Contardo, CSIRO, 2009
% from 'plotc' by Uli Theune, University of Alberta, 2004
%

function h = plotclr(x,y,v,marker,vlim)

if nargin <4
    marker='.';
end

map=colormap;
if nargin >4
    miv = vlim(1) ;
    mav = vlim(2) ;
else
    miv=min(v);
    mav=max(v);
end
clrstep = (mav-miv)/size(map,1) ;
% Plot the points
hold on
for nc=1:size(map,1)
    iv = find(v>miv+(nc-1)*clrstep & v<=miv+nc*clrstep) ;
    plot(x(iv),y(iv),marker,'color',map(nc,:),'markerfacecolor',map(nc,:))
end
hold off

% Re-format the colorbar
h=colorbar;

%set(h,'ylim',[1 length(map)]);
yal=linspace(1,length(map),10);
set(h,'ytick',yal);
% Create the yticklabels
ytl=linspace(miv,mav,10);
s=char(10,4);
for i=1:10
    if min(abs(ytl)) >= 0.001
        B=sprintf('%-4.3f',ytl(i));
    else
        B=sprintf('%-3.1E',ytl(i));
    end
    s(i,1:length(B))=B;
end
set(h,'yticklabel',s);
grid on
view(2)

author: Sable

repo: mcbench-benchmarks


text:
% Test file for @chebfun/join.m.

function pass = test_join(pref)

if ( nargin == 0  )
    pref = chebfunpref();
end

% Generate a few random points in [-1 1] to use as test values.
seedRNG(7681);
xr = 2 * rand(1000, 1) - 1;

% Check scalar CHEBFUNs.
f_op = @sin;
f1 = chebfun(f_op, [-1 -0.5 0]);
f2 = chebfun(f_op, [0 0.5 1]);
f = join(f1, f2);
pass(1) = norm(feval(f, xr) - f_op(xr), inf) < 10*vscale(f)*eps;

% Check array-valued CHEBFUNs.
f_op = @(x) [sin(x) cos(x)];
f1 = chebfun(f_op, [-1 -0.5 0]);
f2 = chebfun(f_op, [0 0.5 1]);
f = join(f1, f2);
err = feval(f, xr) - f_op(xr);
pass(2) = norm(err(:), inf) < 10*vscale(f)*eps;

% Check for quasimatrices.
f1q = quasimatrix(f1);
f2q = quasimatrix(f2);
fq = join(f1q, f2q);
err = feval(fq, xr) - f_op(xr);
pass(3) = norm(err(:), inf) < 10*vscale(fq)*eps;

% Check row CHEBFUNs.
ft = join(f1.', f2.');
err = feval(ft, xr) - f_op(xr).';
pass(4) = ft.isTransposed && (norm(err(:), inf) < 10*vscale(ft)*eps);

% Check row quasimatrices.
ftq = join(f1q.', f2q.');
err = feval(ftq, xr) - f_op(xr).';
pass(5) = ftq(1,:).isTransposed && (norm(err(:), inf) < 10*vscale(ftq)*eps);

% Check operation when the domains don't match.
f = chebfun(@sin, [-1 -0.5 0]);
g = chebfun(@cos, [1 1.5 2]);
h = join(f, g);
pass(6) = isequal(h.domain, [-1 -0.5 0 0.5 1]) && ...
    norm(feval(h, xr) - h_exact_fun(xr), inf) < 10*vscale(h)*eps;

% Check error conditions.
try
    f = join(f1, f2.');
    pass(7) = false;
catch ME
    pass(7) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:join:columnJoin:trans') ...
        || strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:join:dim');
end

% Test for singular function:
op1 = @(x) sin(x)./(x+1);
op2 = @(x) sin(x);
f = chebfun(op1, [-1 -0.5 0], 'exps', [-1 0 0]);
g = chebfun(op2, [1 2]);
h = join(f, g);
x = sort(xr);
h_vals = feval(h, x);
ind = ( x < 0 );
indComp = ~ind;
h_exact1 = feval(op1, x(ind));
h_exact2 = feval(op2, x(indComp)+1);
h_exact = [h_exact1; h_exact2];
err = h_exact - h_vals;
pass(8) = isequal(h.domain, [-1 -0.5 0 1]) && ...
    norm(err, inf) < 1e5*vscale(h)*eps;

% Test for function defined on unbounded domain:

% Set the domain:
domf = [-Inf -3];
domg = [1 Inf];
domCheck = [-1e2 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

opf = @(x) x.*exp(x);
f = chebfun(opf, domf);
opg = @(x) (1-exp(-x))./x;
g = chebfun(opg, domg);
h = join(f, g);
x = sort(x);
h_vals = feval(h, x);
ind = ( x < -3 );
indComp = ~ind;
h_exactf = feval(opf, x(ind));
h_exactg = feval(opg, x(indComp) + (domg(1) - domf(2)));
h_exact = [h_exactf; h_exactg];
err = h_exact - h_vals;
pass(9) = isequal(h.domain, [-Inf -3 Inf]) && ...
    norm(err, inf) < 1e1*vscale(h)*eps;

end

function y = h_exact_fun(x)
    y = zeros(size(x));

    xl_ind = x <= 0;
    xl = x(xl_ind);
    y(xl_ind) = sin(xl);

    xr_ind = x > 0;
    xr = x(xr_ind);
    y(xr_ind) = cos(xr + 1);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%SF_HEX_Q4 Triquartic conforming shape function for hexahedrons (Q4).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_HEX_Q4( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming triquartic Q4 shape functions on hexahedrons with values defined
%   in the nodes, edges, faces, and cell center. XI is [-1..1]^3 reference coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar:  3             Number of space dimensions
%       n_vert      scalar:  8             Number of vertices per cell
%       i_dof       scalar: 1-n_ldof       Local basis function to evaluate
%       xi          [n_sdim]               Local coordinates of evaluation point
%       aInvJac     [n,n_sdim*n_sdim]      Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SFLAG4, SF_HEX_Q1

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [map,textons] = computeTextons(fim,k)
% function [map,textons] = computeTextons(fim,k)

d = numel(fim);
n = numel(fim{1});
data = zeros(d,n);
for i = 1:d,
  data(i,:) = fim{i}(:)';
end

[map,textons] = kmeansML(k,data,'maxiter',30,'verbose',0);
[w,h] = size(fim{1});
map = reshape(map,w,h);

author: zouchuhang

repo: LayoutNet


text:
function rx = rxSetXform(rx, xform, aboutCenter)
%
% rx = rxSetXform(rx, xform, [aboutCenter]);
%
% Sets the xform field of a mrRx struct,
% and also sets UI controls of the GUI
% to agree with this transform.
%
% The aboutCenter flag is an optional flag
% specifying whether the rotations in the
% new xform are intended to rotate about the
% center of prescription or not. (E.g., mrVista
% alignments rotate about the corner, b/c the
% math is more straightforward). If 0 [default],
% this means that further adjustment is needed to
% ensure the rotations act about the center, so the
% xform is adjusted to compensate (the mapping of
% points is the same though).
%
% ras 03/05.
if notDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if notDefined('aboutCenter')
    aboutCenter = 0;
end

if aboutCenter==0
    % we center the rx at 0,0,0 to rotate about the
	% center -- this compensatory translation ensures
	% that zero settings return an unchanged matrix.
	% Because of this, the UI settings reflect a diff't
	% set of rotations from that in the mrVista alignment.
	% Compute the modified xform (after a similar set of
	% changes in rxRefresh, it will return the original
	% mrVista xform matrix):
	shift = [eye(3) -rx.rxDims([2 1 3])'./2; 0 0 0 1];
	xform = shift*xform/shift;
end

% set prev xform field
rx.prevXform = rx.xform;

% set xform field
rx.xform = xform;

% set ui controls to agree w/ the xform
[trans rot flip] = affineDecompose(rx.xform);
rot = rad2deg(rot+pi) - 180;
if ishandle(rx.ui.controlFig)
    rxSetSlider(rx.ui.corTrans, trans(1));
    rxSetSlider(rx.ui.axiTrans, trans(2));
    rxSetSlider(rx.ui.sagTrans, trans(3));
    rxSetSlider(rx.ui.corRot, rot(1));
    rxSetSlider(rx.ui.axiRot, rot(2));
    rxSetSlider(rx.ui.sagRot, rot(3));
    set(rx.ui.corFlip, 'Value', (flip(1)<0));
    set(rx.ui.axiFlip, 'Value', (flip(2)<0));
    set(rx.ui.sagFlip, 'Value', (flip(3)<0));

    % set the control fig w/ the new rx:
    set(rx.ui.controlFig, 'UserData', rx);
    rxRefresh(rx, 0);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ ap, ipvt, info ] = chpfa ( ap, n )

%*****************************************************************************80
%
%% CHPFA factors a complex hermitian packed matrix.
%
%  Discussion:
%
%    To solve A*X = B, follow CHPFA by CHPSL.
%
%    To compute inverse(A)*C, follow CHPFA by CHPSL.
%
%    To compute determinant(A), follow CHPFA by CHPDI.
%
%    To compute inertia(A), follow CHPFA by CHPDI.
%
%    To compute inverse(A), follow CHPFA by CHPDI.
%
%  Packed storage:
%
%    The following program segment will pack the upper
%    triangle of a hermitian matrix.
%
%      k = 0
%      do j = 1, n
%        do i = 1, j
%          k = k + 1
%          ap(k) = a(i,j)
%        end
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex AP(N*(N+1)/2); the packed form
%    of a hermitian matrix.  The columns of the upper triangle are
%    stored sequentially in a one-dimensional array.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex AP(N*(N+1)/2); a block diagonal matrix and the
%    multipliers which were used to obtain it stored in packed form.
%    The factorization can be written A = U*D*hermitian(U) where U is a
%    product of permutation and unit upper triangular matrices, hermitian(U)
%    is the conjugate transpose of U, and D is block diagonal with 1 by 1
%    and 2 by 2 blocks.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO.
%    0, normal value.
%    K, if the K-th pivot block is singular.  This is not an error condition
%    for this subroutine, but it does indicate that CHPSL or CHPDI may divide
%    by zero if called.
%

%
%  Initialize.
%
%  ALPHA is used in choosing pivot block size.
%
  alpha = ( 1.0 + sqrt ( 17.0 ) ) / 8.0;

  info = 0;
%
%  Main loop on K, which goes from N to 1.
%
  k = n;
  ik = ( n * ( n - 1 ) ) / 2;

  while ( 1 )
%
%  Leave the loop if K = 0 or K = 1.
%
    if ( k == 0 )
      break
    end

    if ( k == 1 )
      ipvt(1) = 1;
      if ( cabs1 ( ap(1) ) == 0.0 )
        info = 1;
      end
      break
    end
%
%  This section of code determines the kind of
%  elimination to be performed.  When it is completed,
%  KSTEP will be set to the size of the pivot block, and
%  SWAP will be set to TRUE if an interchange is
%  required.
%
    km1 = k - 1;
    kk = ik + k;
    absakk = cabs1 ( ap(kk) );
%
%  Determine the largest off-diagonal element in column K.
%
    imax = icamax ( k-1, ap(ik+1:ik+k-1), 1 );
    imk = ik + imax;
    colmax = cabs1 ( ap(imk) );

    if ( alpha * colmax <= absakk )

      kstep = 1;
      swap = 0;
%
%  Determine the largest off-diagonal element in row IMAX.
%
    else

      rowmax = 0.0;
      im = ( imax * ( imax - 1 ) ) / 2;
      imj = im + 2 * imax;

      for j = imax + 1 : k
        rowmax = max ( rowmax, cabs1 ( ap(imj) ) );
        imj = imj + j;
      end

      if ( imax ~= 1 )
        jmax = icamax ( imax-1, ap(im+1:im+imax-1), 1 );
        jmim = jmax + im;
        rowmax = max ( rowmax, cabs1 ( ap(jmim) ) );
      end

      imim = imax + im;

      if ( alpha * rowmax <= cabs1 ( ap(imim) ) )
        kstep = 1;
        swap = 1;
      elseif ( alpha * colmax * ( colmax / rowmax ) <= absakk )
        kstep = 1;
        swap = 0;
      else
        kstep = 2;
        swap = ( imax ~= km1 );
      end

    end
%
%  Column K is zero.  Set INFO and iterate the loop.
%
    if ( max ( absakk, colmax ) == 0.0 )
      ipvt(k) = k;
      info = k;
      ik = ik - ( k - 1 );
      if ( kstep == 2 )
        ik = ik - ( k - 2 );
      end
      k = k - kstep;
      continue
    end

    if ( kstep ~= 2 )
%
%  1 x 1 pivot block.
%
      if ( swap )

        temp             = ap(im+1:im+imax);
        ap(im+1:im+imax) = ap(ik+1:ik+imax);
        ap(ik+1:ik+imax) = temp;

        imj = ik + imax;

        for jj = imax : k
          j = k + imax - jj;
          jk = ik + j;
          t       = conj ( ap(jk) );
          ap(jk)  = conj ( ap(imj) );
          ap(imj) = t;
          imj = imj - ( j - 1 );
        end

      end
%
%  Perform the elimination.
%
      ij = ik - ( k - 1 );
      for jj = 1 : km1
        j = k - jj;
        jk = ik + j;
        mulk = -ap(jk) / ap(kk);
        t = conj ( mulk );
        ap(ij+1:ij+j) = ap(ij+1:ij+j) + t * ap(ik+1:ik+j);
        ijj = ij + j;
        ap(ijj) = real ( ap(ijj) );
        ap(jk) = mulk;
        ij = ij - ( j - 1 );
      end
%
%  Set the pivot array.
%
      if ( swap )
        ipvt(k) = imax;
      else
        ipvt(k) = k;
      end
%
%  2 x 2 pivot block.
%
    else

      km1k = ik + k - 1;
      ikm1 = ik - ( k - 1 );

      if ( swap )

        temp                 = ap(im+1:im+imax);
        ap(im+1:im+imax)     = ap(ikm1+1:ikm1+imax);
        ap(ikm1+1:ikm1+imax) = temp;

        imj = ikm1 + imax;

        for jj = imax : km1
          j = km1 + imax - jj;
          jkm1 = ikm1 + j;
          t        = conj ( ap(jkm1) );
          ap(jkm1) = conj ( ap(imj) );
          ap(imj)  = t;
          imj = imj - ( j - 1 );
        end

        t        = ap(km1k);
        ap(km1k) = ap(imk);
        ap(imk)  = t;

      end
%
%  Perform the elimination.
%
      km2 = k - 2;

      if ( km2 ~= 0 )

        ak = ap(kk) / ap(km1k);
        km1km1 = ikm1 + k - 1;
        akm1 = ap(km1km1) / conj ( ap(km1k) );
        denom = 1.0 - ak * akm1;
        ij = ik - ( k - 1 ) - ( k - 2 );

        for jj = 1 : km2
          j = km1 - jj;
          jk = ik + j;
          bk = ap(jk) / ap(km1k);
          jkm1 = ikm1 + j;
          bkm1 = ap(jkm1) / conj ( ap(km1k) );
          mulk = ( akm1 * bk - bkm1 ) / denom;
          mulkm1 = ( ak * bkm1 - bk ) / denom;
          t = conj ( mulk );
          ap(ij+1:ij+j) = ap(ij+1:ij+j) + t * ap(ik+1:ik+j);
          t = conj ( mulkm1 );
          ap(ij+1:ij+j) = ap(ij+1:ij+j) + t * ap(ikm1+1:ikm1+j);
          ap(jk) = mulk;
          ap(jkm1) = mulkm1;
          ijj = ij + j;
          ap(ijj) = real ( ap(ijj) );
          ij = ij - ( j - 1 );
        end

      end
%
%  Set the pivot array.
%
      if ( swap )
        ipvt(k) = -imax;
      else
        ipvt(k) = 1 - k;
      end

      ipvt(k-1) = ipvt(k);

    end

    ik = ik - ( k - 1 );
    if ( kstep == 2 )
      ik = ik - ( k - 2 );
    end

    k = k - kstep;

  end

  return
end


text:
function left = r8vec_bracket4 ( nt, t, ns, s )

%*****************************************************************************80
%
%% R8VEC_BRACKET4 finds the interval to each of a vector of values.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    The routine always returns the index LEFT of the sorted array
%    T with the property that either
%    *  T is contained in the interval [ T(LEFT), T(LEFT+1) ], or
%    *  T < T(LEFT) = T(1), or
%    *  T > T(LEFT+1) = T(N).
%
%    The routine is useful for interpolation problems, where
%    the abscissa must be located within an interval of data
%    abscissas for interpolation, or the "nearest" interval
%    to the (extreme) abscissa must be found so that extrapolation
%    can be carried out.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, length of the input array.
%
%    Input, real T(NT), an array that has been sorted
%    into ascending order.
%
%    Input, integer NS, the number of points to be bracketed.
%
%    Input, real S(NS), values to be bracketed by entries of T.
%
%    Output, integer LEFT(NS).
%    LEFT(I) is set so that the interval [ T(LEFT(I)), T(LEFT(I)+1) ]
%    is the closest to S(I); it either contains S(I), or else S(I)
%    lies outside the interval [ T(1), T(NT) ].
%

%
%  Check the input data.
%
  if ( nt < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_BRACKET4 - Fatal error!\n' );
    fprintf ( 1, '  NT must be at least 2.\n' );
    error ( 'R8VEC_BRACKET4 - Fatal error!' );
  end

  for i = 1 : ns

    left(i) = floor ( ( nt + 1 ) / 2 );
%
%  CASE 1: S < T(LEFT):
%  Search for S in [T(I), T(I+1)] for intervals I = 1 to LEFT-1.
%
    if ( s(i) < t(left(i)) )

      if ( left(i) == 1 )
        continue
      elseif ( left(i) == 2 )
        left(i) = 1;
        continue
      elseif ( t(left(i)-1) <= s(i) )
        left(i) = left(i) - 1;
        continue
      elseif ( s(i) <= t(2) )
        left(i) = 1;
        continue
      end
%
%  ...Binary search for S in [T(I), T(I+1)] for intervals I = 2 to LEFT-2.
%
      low = 2;
      high = left(i) - 2;

      while ( 1 )

        if ( low == high )
          left(i) = low;
          break
        end

        mid = floor ( ( low + high + 1 ) / 2 );

        if ( t(mid) <= s(i) )
          low = mid;
        else
          high = mid - 1;
        end

      end
%
%  CASE2: T(LEFT+1) < S:
%  Search for S in [T(I),T(I+1)] for intervals I = LEFT+1 to N-1.
%
    elseif ( t(left(i)+1) < s(i) )

      if ( left(i) == nt - 1 )
        continue
      elseif ( left(i) == nt - 2 )
        left(i) = left(i) + 1;
        continue
      elseif ( s(i) <= t(left(i)+2) )
        left(i) = left(i) + 1;
        continue
      elseif ( t(nt-1) <= s(i) )
        left(i) = nt - 1;
        continue
      end
%
%  ...Binary search for S in [T(I), T(I+1)] for intervals I = LEFT+2 to NT-2.
%
      low = left(i) + 2;
      high = nt - 2;

      while ( 1 )

        if ( low == high )
          left(i) = low;
          break
        end

        mid = floor ( ( low + high + 1 ) / 2 );

        if ( t(mid) <= s(i) )
          low = mid;
        else
          high = mid - 1;
        end

      end
%
%  CASE3: T(LEFT) <= S <= T(LEFT+1):
%  S is in [T(LEFT), T(LEFT+1)].
%
    else

    end

  end

  return
end


text:
function grid = reshapeGrid(~, grid)
%RESHAPEGRID    Add the repeated endpoints in the lambda- and theta-directions
%and extract half of the points in the theta-direction.

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

% Note: We use periodic discretizations in longitude (lambda) and latitude
% (theta). These dicretizations do not include the repeated endpoints. Before
% plotting data with SPINOPERATOR/INITIALIZEMOVIE, we add these repeated
% endpoints to the grid, and extract half of the points in theta, those that
% correspond to [0, pi]. (We use the DFS method with doubled-up [-pi, pi]
% theta-grid.)

% Get the points:
ll = grid{1};
tt = grid{2};

% Get the number of grid points:
N = length(ll);

% Add the endpoints in lambda and theta:
ll = [ll, 2*ll(:,end) - ll(:,end-1)];
ll = [ll; ll(1,:)];
tt = [tt; 2*tt(end,:) - tt(end-1,:)];
tt = [tt, tt(:,1)];

% Extract half of the points in the theta-direction:
ll = ll(N/2+1:end,:);
tt = tt(N/2+1:end,:);

% Output the new grid:
grid{1} = ll;
grid{2} = tt;

end


text:
function stats = log_iterates(problem, x, stats)
% function stats = log_iterates(problem, x, stats)
%
% This small helper function is used to log the sequence of iterates
% generated by the Riemannian Staircase

stats.Yvals = x';

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function [decentr_state, dgs_stats, final_increments, ...
    final_increment_times, opt_handle, final_accuracy_measurements] = ...
    wrapUpOptimization(...
    decentr_state, opt_handle, distributed_mapper_location, params, ...
    final_time)

t = tic;
% Collect data from still running optimization and launch last one.
[decentr_state, opt_handle, dgs_stats, opt_increment] = ...
    manageAsyncGaussSeidel(...
    decentr_state, opt_handle, distributed_mapper_location, true, ...
    params.opt_max_iters, true);
t0 = toc(t);
for i = 1:numel(dgs_stats)
    dgs_stats{i}.end_time = final_time + t0;
end
final_increments = ...
    {cat(3, opt_increment, zeros(size(opt_increment)), ...
    zeros(size(opt_increment)))};
final_increment_times = [final_time + t0];
[accuracy_measurement, decentr_state] = ...
    evalAccuracy(...
    decentr_state, opt_handle, final_time + t0);
final_accuracy_measurements = {accuracy_measurement};
plotDecentrState(decentr_state);
pause(0.0001);

% Collect last one (final launch will be ignored).
[decentr_state, ~, last_dgs_stats, opt_increment] = ...
    manageAsyncGaussSeidel(...
    decentr_state, opt_handle, distributed_mapper_location, true, ...
    params.opt_max_iters, false);
t1 = toc(t);
for i = 1:numel(last_dgs_stats)
    last_dgs_stats{i}.end_time = final_time + t1;
end
dgs_stats = [dgs_stats; last_dgs_stats];
final_increments = [final_increments;
    cat(3, opt_increment, zeros(size(opt_increment)), ...
    zeros(size(opt_increment)))];
final_increment_times = [final_increment_times; final_time + t1];
[accuracy_measurement, decentr_state] = ...
    evalAccuracy(decentr_state, opt_handle, final_time + t1);
final_accuracy_measurements = [final_accuracy_measurements; ...
    accuracy_measurement];
plotDecentrState(decentr_state);
pause(0.0001);

end

author: uzh-rpg


text:
function hh = rotateXLabels( ax, angle, varargin )
%rotateXLabels: rotate any xticklabels
%
%   hh = rotateXLabels(ax,angle) rotates all XLabels on axes AX by an angle
%   ANGLE (in degrees). Handles to the resulting text objects are returned
%   in HH.
%
%   hh = rotateXLabels(ax,angle,param,value,...) also allows one or more
%   optional parameters to be specified. Possible parameters are:
%     'MaxStringLength'   The maximum length of label to show (default inf)
%
%   Examples:
%   >> bar( hsv(5)+0.05 )
%   >> days = {'Monday','Tuesday','Wednesday','Thursday','Friday'};
%   >> set( gca(), 'XTickLabel', days )
%   >> rotateXLabels( gca(), 45 )
%
%   See also: GCA
%             BAR

%   Copyright 2006-2010 The MathWorks Ltd.
%   $Revision: 52 $
%   $Date: 2010-09-30 11:19:49 +0100 (Thu, 30 Sep 2010) $

error( nargchk( 2, inf, nargin ) );
if ~isnumeric( angle ) || ~isscalar( angle )
    error( 'RotateXLabels:BadAngle', 'Parameter ANGLE must be a scalar angle in degrees' )
end
angle = mod( angle, 360 );

[maxStringLength] = parseInputs( varargin{:} );

% Get the existing label texts and clear them
[vals, labels] = findAndClearExistingLabels( ax, maxStringLength );

% Create the new label texts
h = createNewLabels( ax, vals, labels, angle );

% Reposition the axes itself to leave space for the new labels

%repositionAxes( ax );

% If an X-label is present, move it too
%repositionXLabel( ax );

% Store angle
setappdata( ax, 'RotateXLabelsAngle', angle );

% Only send outputs if requested
if nargout
    hh = h;
end

%-------------------------------------------------------------------------%
    function [maxStringLength] = parseInputs( varargin )
        % Parse optional inputs
        maxStringLength = inf;
        if nargin > 0
            params = varargin(1:2:end);
            values = varargin(2:2:end);
            if numel( params ) ~= numel( values )
                error( 'RotateXLabels:BadSyntax', 'Optional arguments must be specified as parameter-value pairs.' );
            end
            if any( ~cellfun( 'isclass', params, 'char' ) )
                error( 'RotateXLabels:BadSyntax', 'Optional argument names must be specified as strings.' );
            end
            for pp=1:numel( params )
                switch upper( params{pp} )
                    case 'MAXSTRINGLENGTH'
                        maxStringLength = values{pp};

                    otherwise
                        error( 'RotateXLabels:BadParam', 'Optional parameter ''%s'' not recognised.', params{pp} );
                end
            end
        end
    end % parseInputs
%-------------------------------------------------------------------------%
    function [vals,labels] = findAndClearExistingLabels( ax, maxStringLength )
        % Get the current tick positions so that we can place our new labels
        vals = get( ax, 'XTick' );

        % Now determine the labels. We look first at for previously rotated labels
        % since if there are some the actual labels will be empty.
        ex = findall( ax, 'Tag', 'RotatedXTickLabel' );
        if isempty( ex )
            % Store the positions and labels
            labels = get( ax, 'XTickLabel' );
            if isempty( labels )
                % No labels!
                return
            else
                if ~iscell(labels)
                    labels = cellstr(labels);
                end
            end
            % Clear existing labels so that xlabel is in the right position
            set( ax, 'XTickLabel', {}, 'XTickMode', 'Manual' );
            setappdata( ax, 'OriginalXTickLabels', labels );
        else
            % Labels have already been rotated, so capture them
            labels = getappdata( ax, 'OriginalXTickLabels' );
            delete(ex);
        end
        % Limit the length, if requested
        if isfinite( maxStringLength )
            for ll=1:numel( labels )
                if length( labels{ll} ) > maxStringLength
                    labels{ll} = labels{ll}(1:maxStringLength);
                end
            end
        end

    end % findAndClearExistingLabels

%-------------------------------------------------------------------------%
    function textLabels = createNewLabels( ax, vals, labels, angle )
        % Work out the ticklabel positions
%         ylim = get(ax,'YLim');
%         y = ylim(1);
        zlim = get(ax,'ZLim');
        z = zlim(1);

        % We want to work
        % in normalised coords, but switch to normalised after setting
        % position causes positions to be rounded to the nearest pixel.
        % Instead we can do the calculation by hand.
        xlim = get( ax, 'XLim' );
        normVals = (vals-xlim(1))/(xlim(2)-xlim(1));
        y = 0;

        % Now create new text objects in similar positions.
        textLabels = -1*ones( numel( vals ), 1 );
        for ll=1:numel(vals)
            textLabels(ll) = text( ...
                'Units', 'Normalized', ...
                'Position', [normVals(ll), y, z], ...
                'String', labels{ll}, ...
                'Parent', ax, ...
                'Clipping', 'off', ...
                'Rotation', angle, ...
                'Tag', 'RotatedXTickLabel', ...
                'UserData', vals(ll) );
        end

        % Now copy font properties into the texts
        updateFont();

        % Depending on the angle, we may need to change the alignment. We change
        % alignments within 5 degrees of each 90 degree orientation.
        if 0 <= angle && angle < 5
            set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Top' );
        elseif 5 <= angle && angle < 85
            set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Top' );
        elseif 85 <= angle && angle < 95
            set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Middle' );
        elseif 95 <= angle && angle < 175
            set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Bottom' );
        elseif 175 <= angle && angle < 185
            set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Bottom' );
        elseif 185 <= angle && angle < 265
            set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom' );
        elseif 265 <= angle && angle < 275
            set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Middle' );
        elseif 275 <= angle && angle < 355
            set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Top' );
        else % 355-360
            set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Top' );
        end
    end % createNewLabels

%-------------------------------------------------------------------------%
    function repositionAxes( ax )
        % Reposition the axes so that there's room for the labels
        % Note that we only do this if the OuterPosition is the thing being
        % controlled
        if ~strcmpi( get( ax, 'ActivePositionProperty' ), 'OuterPosition' )
            return;
        end

        % Work out the maximum height required for the labels
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        maxHeight = 0;
        for ii=1:numel(vals)
            ext = get( textLabels(ii), 'Extent' );
            if ext(4) > maxHeight
                maxHeight = ext(4);
            end
        end

        % Remove listeners while we mess around with things, otherwise we'll
        % trigger redraws recursively
        removeListeners( ax );

        % Change to normalized units for the position calculation
        oldUnits = get( ax, 'Units' );
        set( ax, 'Units', 'Normalized' );

        % Not sure why, but the extent seems to be proportional to the height of the axes.
        % Correct that now.
        set( ax, 'ActivePositionProperty', 'Position' );
        pos = get( ax, 'Position' );
        axesHeight = pos(4);
        % Make sure we don't adjust away the axes entirely!
        heightAdjust = min( (axesHeight*0.9), maxHeight*axesHeight );

        % Move the axes
        if isappdata( ax, 'OriginalAxesPosition' )
            pos = getappdata( ax, 'OriginalAxesPosition' );
        else
            pos = get(ax,'Position');
            setappdata( ax, 'OriginalAxesPosition', pos );
        end
        setappdata( ax, 'PreviousYLim', get( ax, 'YLim' ) );
        set( ax, 'Position', pos+[0 heightAdjust 0 -heightAdjust] )
        set( ax, 'Units', oldUnits );
        set( ax, 'ActivePositionProperty', 'OuterPosition' );

        % Make sure we find out if axes properties are changed
        addListeners( ax );

    end % repositionAxes

%-------------------------------------------------------------------------%
    function repositionXLabel( ax )
        % Try to work out where to put the xlabel

        removeListeners( ax );
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        maxHeight = 0;
        for ll=1:numel(vals)
            ext = get( textLabels(ll), 'Extent' );
            if ext(4) > maxHeight
                maxHeight = ext(4);
            end
        end

        % Use the new max extent to move the xlabel
        xlab = get(ax,'XLabel');
        set( xlab, 'Units', 'Normalized', 'Position', [0.5 -maxHeight 0] );
        addListeners( ax );
    end % repositionXLabel

%-------------------------------------------------------------------------%
    function updateFont()
        % Update the rotated text fonts when the axes font changes
        properties = {
            'FontName'
            'FontSize'
            'FontAngle'
            'FontWeight'
            'FontUnits'
            };
        propertyValues = get( ax, properties );
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        set( textLabels, properties, propertyValues );
    end % updateFont

%-------------------------------------------------------------------------%
    function onAxesFontChanged()
        updateFont();
        repositionAxes( ax );
        repositionXLabel( ax );
    end % onAxesFontChanged

%-------------------------------------------------------------------------%
    function onAxesPositionChanged()
        % We need to accept the new position, so remove the appdata before
        % redrawing
        if isappdata( ax, 'OriginalAxesPosition' )
            rmappdata( ax, 'OriginalAxesPosition' );
        end
        if isappdata( ax, 'OriginalXLabelPosition' )
            rmappdata( ax, 'OriginalXLabelPosition' );
        end
        repositionAxes( ax );
        repositionXLabel( ax );
    end % onAxesPositionChanged

%-------------------------------------------------------------------------%
    function onAxesLimitsChanged()
        % The limits have moved, so make sure the labels are still ok
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        xlim = get( ax, 'XLim' );
        for tt=1:numel( textLabels )
            xval = get( textLabels(tt), 'UserData' );
            pos(1) = (xval-xlim(1)) / (xlim(2)-xlim(1));
            pos(2) = 0;
            % If the tick is off the edge, make it invisible
            if xval<xlim(1) || xval>xlim(2)
                set( textLabels(tt), 'Visible', 'off', 'Position', pos )
            elseif ~strcmpi( get( textLabels(tt), 'Visible' ), 'on' )
                set( textLabels(tt), 'Visible', 'on', 'Position', pos )
            else
                % Just set the position
                set( textLabels(tt), 'Position', pos );
            end
        end

        repositionXLabel( ax );

%         xlab = get(ax,'XLabel');
%         if ~strcmpi( get( xlab, 'Units' ), 'Data' )
%             set( xlab, 'Units', 'data' );
%         end
%         pos = get( xlab, 'Position' );
%         orig_yrange = diff( orig_ylim );
%         new_yrange = diff( ylim );
%         offset = (pos(2)-orig_ylim(1)) / orig_yrange;
%         pos(2) = offset * new_yrange + ylim(1);
%         pos(1) = mean( xlim );
%         set( xlab, 'Position', pos );
%         setappdata( ax, 'PreviousYLim', ylim );
    end % onAxesPositionChanged

%-------------------------------------------------------------------------%
    function addListeners( ax )
        % Create listeners. We store the array of listeners in the axes to make
        % sure that they have the same life-span as the axes they are listening to.
        axh = handle( ax );
        listeners = [
            handle.listener( axh, findprop( axh, 'FontName' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontSize' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontWeight' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontAngle' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontUnits' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'OuterPosition' ), 'PropertyPostSet', @onAxesPositionChanged )
            handle.listener( axh, findprop( axh, 'XLim' ), 'PropertyPostSet', @onAxesLimitsChanged )
            handle.listener( axh, findprop( axh, 'YLim' ), 'PropertyPostSet', @onAxesLimitsChanged )
            ];
        setappdata( ax, 'RotateXLabelsListeners', listeners );
    end % addListeners

%-------------------------------------------------------------------------%
    function removeListeners( ax )
        % Rempove any property listeners whilst we are fiddling with the axes
        if isappdata( ax, 'RotateXLabelsListeners' )
            delete( getappdata( ax, 'RotateXLabelsListeners' ) );
            rmappdata( ax, 'RotateXLabelsListeners' );
        end
    end % removeListeners

end % EOF

repo: MARS-evaluation


text:
function [res,err] = Dot_(varargin)
%DOT_         Dot products in K-fold precision
%
%   res = Dot_(A1,x1,...,An,xn,K)
%
%   [res,err] = Dot_(A1,x1,...,An,xn,K)
%
% res or res+err is sum(Ai*xi)
%     K optional, default K=2
%     K   >0  approximate result as if computed in K-fold precision
%         <0  inclusion of result as if computed in |K|-fold precision
%
% First factor Ai may be either matrix or scalar; in the latter case
%   res is of the size of xi.
% Second factor xi must be vector or scalar.
% All factors may be complex, but must be non-interval.
% Sizes of every product Ai*xi must be the same.
%
% Precision K must be at least 2.
%
% Interval input makes only sense, if all intervals are degenerated (point intervals);
%   therefore omitted.
%
% For speed and to fight interpretation overhead, 2 extra matrices are necessary; for
%   the product of to nxn matrices this would be n^3 memory, therefore this is not allowed.
%
% Examples are
%
%   res = Dot_(x',y);                  Approximate scalar product x'*y in quadruple precision
%   res = Dot_(A,x,-1,b);              Approximate residual A*x-b in quadruple precision
%   res = Dot_(A,x,-1,b,-2);           Inclusion of residual A*x-b in quadruple precision
%   res = Dot_(A,x1,A,x2,-1,b,3);      Approximate residual A*x1+A*x2-b in 3-fold precision
%
%For randomly generated full matrices of dimension n with b=a*randn(n,1) the following
%  table lists the computing times
%    t1   for ordinary residual  b - A*x
%    t2   with improved residual by lssresidual
%    t3   with quadruple precision residual by Dot_
%as well as the ratios to t1. All times are in seconds on a 800 MHz Pentium III Laptop.
%
%    n     t1      t2     t3    t2/t1  t3/t1     approximation of b-A*x
%--------------------------------------------
%   50   0.0002  0.001  0.003    3.1   14.6
%  100   0.0002  0.003  0.010   13.7   45.7
%  200   0.0005  0.010  0.045   20.5   89.5
%  500   0.0064  0.089  0.416   13.7   64.6
% 1000   0.0246  0.347  1.698   14.1   68.9
%
%The same table for verified inclusion of the residual b-A*x is as follows.
%
%    n     t1      t2     t3    t2/t1  t3/t1     inclusion of b-A*x
%--------------------------------------------
%   50   0.0012  0.001  0.004    0.9    3.0
%  100   0.0013  0.002  0.010    1.9    7.5
%  200   0.0020  0.012  0.048    6.1   23.8
%  500   0.0138  0.103  0.381    7.5   27.5
% 1000   0.0505  0.421  1.863    8.3   36.9
%
%Finally we display the achieved accuracy of the three methods for
%  matrices of dimension 100 and different condition numbers, a randomly
%  generated right hand side b and computed solution x=A\b.
%  We treat 100 test cases each and display the median and maximum componentwise
%  relative error of the results. Matrices are normed to 1.
%
%                A*x-b          lssresidual         Dot_
%condition  median maximum    median maximum    median maximum
%---------------------------------------------------------------
%    1e2    6.8e-16 1.0e+0   5.0e-21 8.9e-16   0.0e+0  4.7e-30
%    1e5    1.8e-01 1.0e+0   1.8e-18 3.0e-01   5.0e-29 1.7e-27
%   1e10    2.9e-01 1.0e+0   2.4e-06 6.8e-01   3.3e-24 1.8e-22
%   1e14    2.8e-01 1.0e+0   2.1e-06 1.0e+0    2.7e-20 2.6e-18
%
%The more we pay, the more we get.
%
% Implements algorithms Dot2 and DotK from
%   T. Ogita, S.M. Rump, S. Oishi: Accurate Sum and Dot Product,
%     SIAM Journal on Scientific Computing (SISC), 26(6):1955-1988, 2005
%

% written  11/08/03     S.M. Rump
% modified 11/30/03     S.M. Rump  complex data allowed
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 08/08/04     S.M. Rump  sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 12/12/05     S.M. Rump  apapted to paper, improved performance
% modified 05/30/07     S.M. Rump  typo
% modified 08/07/10     S.M. Rump  upper case Dot_
% modified 06/06/11     S.M. Rump  result in two parts
% modified 08/26/12     S.M. Rump  global variables removed, rounding
% modified 09/05/12     S.M. Rump  sparse complex (thanks to M. Lange)
% modified 10/17/12     S.M. Rump  comments
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  % detect parameters
  len = length(varargin);
  if even(len)                      % default precision
    K = 2;
  else
    K = varargin{len};
  end
  idim = 0;

  % check precision K
  if abs(K)<2
    error('invalid parameter K for Dot_')
  end

  % detect dimension of result
  if prod(size(varargin{1}))==1
    [m n] = size(varargin{2});
  else
    m = size(varargin{1},1);        % m or n must be 1
    n = size(varargin{1},2);
  end

  % Initialization
  factor = 134217729;               % splitting factor 2^27+1
  if abs(K)==2
    res = zeros(m,1);
    if K>0
      s = 0;
    else
      sinf = 0;
      ssup = 0;
    end
  else
    res = [];
    S = [];
  end

  % Compute array to sum up
  for i=1:floor(nargin/2)

    % compute next product A*x
    A = varargin{2*i-1};
    x = varargin{2*i};
    sA = size(A);
    sx = size(x);
    if sx(2)~=1
      error('second factor in Dot_ must be vector')
    end

    if ( prod(sA)==1 ) & ( sx(1)~=m )
      error('dimensions of summands in Dot_ do not match')
    end

    if ( prod(sA)==1 ) & ( abs(A)==1 )   % first factor A scalar +/- 1

      % TwoProduct(A,x)
      if A==1
        [res,q] = sumK([res x],2);
      else
        [res,q] = sumK([res -x],2);
      end
      if abs(K)==2
        if K>0
          s = s + q;
        else
          setround(-1)
          sinf = sinf + q;
          setround(1)
          ssup = ssup + q;
          setround(0)
        end
      else
        S = [S q];
      end
      idim = idim + 1;

    else                                    % first factor non-scalar or ~=  +/- 1

      % matrix-vector product or dot product
      if ( prod(sA)~=1 ) & ( sA(2)~=sx(1) )
        error('inner dimensions in Dot_ do not match')
      end
      if issparse(A)
        [iA,jA,ssA] = find(A);
      end
      C = factor*A;
      A1 = C - ( C - A );                       % upper part of A
      A2 = A - A1;                              % A = A1+A2 exact splitting

      C = factor*x;
      x1 = C - ( C - x );                       % upper part of x
      x2 = x - x1;                              % x = x1+x2 exact splitting

      I = sqrt(-1);
      if sA(1)~=1                               % m~=1: matrix-vector product
        if issparse(A)                          % factors s.t. A*x = sum( (A1+A2).*(x1+x2) , 2 )
          x1 = sparse(iA,jA,x1(jA),m,sA(2));
          x2 = sparse(iA,jA,x2(jA),m,sA(2));
        else
          x1 = x1(:,ones(1,m)).';
          x2 = x2(:,ones(1,m)).';
        end
        if isreal(A) | isreal(x)                % A or x real
          if issparse(A)
            h = sparse(iA,jA,ssA.*x(jA),m,sA(2));
          else
            h = A .* ( x(:,ones(1,m)).' );
          end
          idim = idim + sx(1);
        else                                    % both A and x complex
          if issparse(A)
            h = [ sparse(iA,jA,ssA.*real(x(jA)),m,sA(2)) , sparse(iA,jA,I*(ssA.*imag(x(jA))),m,sA(2)) ];
          else
            xx = x(:,ones(1,m)).';
            h = [ A.*real(xx) , A.*(I*imag(xx)) ];
          end
          idim = idim + 2*sx(1);
        end
      elseif sA(2)~=1
        if isreal(A) | isreal(x)                % A or x real
          h = A .* x.';
          x1 = x1.';
          x2 = x2.';
          idim = idim + sx(1);
        else                                    % both A and x complex
          h = [ A.*real(x).' , A.*(I*imag(x)).' ];
          x1 = x1.';
          x2 = x2.';
          idim = idim + 2*sx(1);
        end
      else
        if isreal(A) | isreal(x)                % A or x real
          h = A * x;
          idim = idim + 1;
        else                                    % both A and x complex
          h = [ A*real(x) , A*(I*imag(x)) ];
          idim = idim + 2;
        end
      end

      % TwoProduct(A,x): A*x = h + r
      if isreal(A) | isreal(x)                % A or x real
        r = A2.*x2 - ((( h - A1.*x1 ) - A2.*x1 ) - A1.*x2 );
      else                                    % both A and x complex
        r = [A2.*real(x2),A2.*(I*imag(x2))] - ((( h - [A1.*real(x1),A1.*(I*imag(x1))] ) - ...
          [A2.*real(x1),A2.*(I*imag(x1))] ) - [A1.*real(x2),A1.*(I*imag(x2))] );
      end

      if issparse(h)                    % compress h
        h = compress(h);
      end
      if issparse(r)                    % compress r
        r = compress(r);
      end

      % Summation vector and partial summation for |K|=2
      if abs(K)==2
        [res,q] = sumK([res h],2);      % error-free
        sq2 = size(q,2);
        sr2 = size(r,2);
        if sq2>sr2
          r = [r zeros(size(r,1),sq2-sr2)];
        elseif sq2>sr2
          q = [q zeros(size(q,1),sr2-sq2)];
        end
        if K>0
          s = s + sum(q+r,2);
        else
          setround(-1)
          sinf = sinf + sum(q+r,2);
          setround(1)
          ssup = ssup + sum(q+r,2);
          setround(0)
        end
      else
        [res,q] = sumK([res h],2);
        S = [S q r];
      end

    end

  end

  % Compute result in K-fold precision
  if abs(K)==2                          % quadruple precision
    if K>0
      res = res + s;
    else
      INTLAB_INTVAL_ETA = realmin*eps;  % smallest positive denormalized fl-pt
      setround(-1)
      resinf = sinf - 5*idim*INTLAB_INTVAL_ETA;
      setround(1)
      ressup = ssup + 5*idim*INTLAB_INTVAL_ETA;
      err = infsup(resinf,ressup);
      if nargout==1
        res = res + err;
      end
    end
  else                                  % |K|-fold precision
    if K>0
      res = Sum_([S res],K,2);
    else
      INTLAB_INTVAL_ETA = realmin*eps;  % smallest positive denormalized fl-pt
      [res,q] = sumK([S res],-K);       % error-free transformation
      setround(-1)
      resinf = sum(q,2) - 5*idim*INTLAB_INTVAL_ETA;
      setround(1)
      ressup = sum(q,2) + 5*idim*INTLAB_INTVAL_ETA;
      err = infsup(resinf,ressup);
      if nargout==1
        res = res + err;
      end
    end
  end

  if rndold
    setround(rndold)
  end

function [res,q] = sumK(a,K)
%SUM          Summation of "a" in K-fold precision along 2nd dimension
%             res is approximate sum, sum of errors in q
%             res + sum(q,2) = sum(a,2) in exact arithmetic

  % number of summands
  N = size(a,2);

  % summation (i.e. quadruple for K==2)
  for i=1:K-1
    pi = cumsum(a,2);
    z = diff(pi,1,2);
    q = ( pi(:,1:N-1) - ( pi(:,2:N) - z ) ) + ( a(:,2:N) - z );
    res = pi(:,N);
    if i~=K-1
      a = [q res];
    end
  end

function a = compress(a)
%COMPRESS     horizontal compress of array "a" such that sum(a,2) remains unchanged

  [i,j,s]=find(a.');
  rows = sum(spones(a),2);
  m = max(rows);
  if m==0
    a = sparse([],[],[],size(a,1),m);
    return
  end
  rows = nonzeros(rows);
  index = cumsum([1;rows(1:end-1)]);
  h = ones(size(i));
  h(index(2:end))=-rows(1:end-1)+1;
  a = sparse(j,cumsum(h),s,size(a,1),m);

repo: OpenMAS


text:
function [x, output] = findExtremePathway(fbaModel, obj)
% Finds an extreme ray
%
% USAGE:
%
%    [x, output] = findExtremePathway(fbaModel, obj)
%
% INPUT:
%    fbaModel:    FBA type model
%
% OPTIONAL INPUT:
%    obj:         default = random vector with size depending on fbaModel.S
%
% OUTPUTS:
%    x:           vector from `result`, where `result` is an output of `solveCobraLP` function
%    output:      `output.objval` contains `result.obj`
%
% .. $Revision: 0.1 $  $Date: 2011/05/01 $

A    = fbaModel.S;
[nmet,nrxn] = size(A);

% Convert model to conic form
if isfield(fbaModel, 'ub') && isfield(fbaModel, 'lb')
    revRxns = fbaModel.lb < 0 & fbaModel.ub > 0;
else
    error('missing fields: revRxns or ub and lb\n');
end

A = [A, -A(:,revRxns)];
[~, n] = size(A);

if nargin < 2
    obj = rand(n,1);
end

% Set required model components
LPProblem = struct();
LPProblem.A = sparse([A; ones(1,n)]);
LPProblem.c   = obj;
LPProblem.csense = repmat('E',nmet+1,1);
LPProblem.b   = [zeros(nmet,1); 1];
% Set optional model components
LPProblem.osense= -1;
LPProblem.lb = 0* ones(n,1);
LPProblem.ub = inf* ones(n,1);

result = solveCobraLP(LPProblem);  % Find extreme ray (be aware, that this can easily be a loop of a reversible reaction.

x = result.full(1:nrxn);
x(revRxns) = x(revRxns) - result.full(nrxn+1:end);

output.objval = result.obj;

author: opencobra

repo: cobratoolbox


text:
classdef UniversalElementSet < AbstractElementSet
    %KeplerianElementSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        c3(1,1) double = 0; %km^2/s^2
        rP(1,1) double = 1000 %km
        inc(1,1) double %rad
        raan(1,1) double %rad
        arg(1,1) double %rad
        tau(1,1) double %s

        optVar UniversalElementSetVariable
    end

    properties(Constant)
        typeEnum = ElementSetEnum.UniversalElements
    end

    methods
        function obj = UniversalElementSet(time, c3, rP, inc, raan, arg, tau, frame)
            if(nargin > 0)
                obj.time = time;
                obj.c3 = c3;
                obj.rP = rP;
                obj.inc = inc;
                obj.raan = raan;
                obj.arg = arg;
                obj.tau = tau;
                obj.frame = frame;
            end
        end

        %vectorized
        function cartElemSet = convertToCartesianElementSet(obj)
            cartElemSet = convertToCartesianElementSet(convertToKeplerianElementSet(obj));
        end

        %vectorized
        function kepElemSet = convertToKeplerianElementSet(obj)
%             gmu = obj.frame.getOriginBody().gm;
            gmu = NaN(length(obj), 1);
            for(i=1:length(obj))
                gmu(i) = obj(i).frame.getOriginBody().gm;
            end

            sma = -gmu./[obj.c3];
            ecc = abs(1 - [obj.rP]./sma);

            n = computeMeanMotion(sma, gmu);
            mean = [obj.tau] .* n;
            tru = computeTrueAnomFromMean(mean, ecc);

%             kepElemSet = KeplerianElementSet(obj.time, sma, ecc, obj.inc, obj.raan, obj.arg, tru, obj.frame);
            kepElemSet = repmat(KeplerianElementSet.getDefaultElements(), size(obj));
            for(i=1:length(obj))
                kepElemSet(i) = KeplerianElementSet(obj(i).time, sma(i), ecc(i), obj(i).inc, obj(i).raan, obj(i).arg, tru(i), obj(i).frame);
            end
        end

        %vectorized
        function geoElemSet = convertToGeographicElementSet(obj)
            geoElemSet = convertToGeographicElementSet(convertToCartesianElementSet(obj));
        end

        %vectorized
        function univElemSet = convertToUniversalElementSet(obj)
            univElemSet = obj;
        end

        function elemVect = getElementVector(obj)
            elemVect = [obj.c3,obj.rP,rad2deg(obj.inc),rad2deg(obj.raan),rad2deg(obj.arg),obj.tau];
        end
    end

    methods(Access=protected)
        function displayScalarObject(obj)
            fprintf('Universal State \n\tTime: %0.3f sec UT \n\tC3 Energy: %0.3f km^2/s^2 \n\tRp: %0.9f km \n\tInc: %0.3f deg \n\tRAAN: %0.3f deg \n\tArg Peri: %0.3f deg \n\tTime Past Peri.: %0.3f s \n\tFrame: %s\n', ...
                    obj.time, ...
                    obj.c3, ...
                    obj.rP, ...
                    rad2deg(obj.inc), ...
                    rad2deg(obj.raan), ...
                    rad2deg(obj.arg), ...
                    obj.tau, ...
                    obj.frame.getNameStr());
        end
    end

    methods(Static)
        function elemSet = getDefaultElements()
            elemSet = UniversalElementSet();
        end

        function errMsg = validateInputOrbit(errMsg, hC3, hRp, hInc, hRaan, hArg, hTau, bodyInfo, bndStr, checkElement)
            if(isempty(bndStr))
                bndStr = '';
            else
                bndStr = sprintf(' (%s Bound)', bndStr);
            end

            if(checkElement(2))
                Rp = str2double(get(hRp,'String'));
                enteredStr = get(hRp,'String');
                numberName = ['Radius of Periapsis', bndStr];
                lb = 0;
                ub = Inf;
                isInt = false;
                errMsg = validateNumber(Rp, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(1))
                c3 = str2double(get(hC3,'String'));
                enteredStr = get(hC3,'String');
                numberName = ['C3 Energy', bndStr];
                lb = -Inf;
                ub = Inf;
                isInt = false;
                errMsg = validateNumber(c3, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(3))
                inc = str2double(get(hInc,'String'));
                enteredStr = get(hInc,'String');
                numberName = ['Inclination', bndStr];
                lb = 0;
                ub = 180;
                isInt = false;
                errMsg = validateNumber(inc, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(4))
                raan = str2double(get(hRaan,'String'));
                enteredStr = get(hRaan,'String');
                numberName = ['Right Asc. of the Asc. Node', bndStr];
                lb = -360;
                ub = 360;
                isInt = false;
                errMsg = validateNumber(raan, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(5))
                arg = str2double(get(hArg,'String'));
                enteredStr = get(hArg,'String');
                numberName = ['Argument of Periapsis', bndStr];
                lb = -360;
                ub = 360;
                isInt = false;
                errMsg = validateNumber(arg, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(6))
                tru = str2double(get(hTau,'String'));
                enteredStr = get(hTau,'String');
                numberName = ['Time Past Periapsis', bndStr];
                lb = -Inf;
                ub = Inf;
                isInt = false;
                errMsg = validateNumber(tru, numberName, lb, ub, isInt, errMsg, enteredStr);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function optimParam = tuneParamForModel(model,varargin)
% Optimizes cplex parameters to make model resolution faster.
% Particularly interetsing for large-scale MILP models and repeated runs of
% optimisation.
% While, it won't optimize memory space nor model constraints for numerical
% infeasibilities, tuneParam will provide the optimal set of solver
% parameters for feasbile models. It requires IBM ILOG cplex (for now).
%
% USAGE:
%
%    optimalParameters = tuneParam(model,contFunctName,1000,1000,0);
%
% INPUT:
%         model:         A COBRA model struct.
%         contFunctName: Parameters structure containing the name and value.
%                        A set of routine parameters will be added by the solver
%                        but won't be reported.
%         timelimit:     default is 10000 second
%         nrepeat:       number of row/column permutation of the original
%                        problem, reports robust results.
%                        sets the CPX_PARAM_TUNINGREPEAT parameter
%                        High values of nrepeat would require consequent
%                        memory and swap.
%         printLevel:    0/1/2/3
%
% OUTPUT:
%         optimParam: structure of optimal parameter values directly usable as
%                     contFunctName argument in solveCobraLP function
%
% NOTE:
%       This is just a wrapper function that calls the tuneParam function
%       using a Cobra model structure converted to a LP problem.
%
% .. Author: Thomas Pfau Dec 2017

optimParam = tuneParam(buildLPproblemFromModel(model),varargin{:});


text:
function net = cnnff(net, x)
    n = numel(net.layers);
    net.layers{1}.a{1} = x;
    inputmaps = 1;

    for l = 2 : n   %  for each layer
        if strcmp(net.layers{l}.type, 'c')
            %  !!below can probably be handled by insane matrix operations
            for j = 1 : net.layers{l}.outputmaps   %  for each output map
                %  create temp output map
                z = zeros(size(net.layers{l - 1}.a{1}) - [net.layers{l}.kernelsize - 1 net.layers{l}.kernelsize - 1 0]);
                for i = 1 : inputmaps   %  for each input map
                    %  convolve with corresponding kernel and add to temp output map
                    z = z + convn(net.layers{l - 1}.a{i}, net.layers{l}.k{i}{j}, 'valid');
                end
                %  add bias, pass through nonlinearity
                net.layers{l}.a{j} = sigm(z + net.layers{l}.b{j});
            end
            %  set number of input maps to this layers number of outputmaps
            inputmaps = net.layers{l}.outputmaps;
        elseif strcmp(net.layers{l}.type, 's')
            %  downsample
            for j = 1 : inputmaps
                z = convn(net.layers{l - 1}.a{j}, ones(net.layers{l}.scale) / (net.layers{l}.scale ^ 2), 'valid');   %  !! replace with variable
                net.layers{l}.a{j} = z(1 : net.layers{l}.scale : end, 1 : net.layers{l}.scale : end, :);
            end
        end
    end

    %  concatenate all end layer feature maps into vector
    net.fv = [];
    for j = 1 : numel(net.layers{n}.a)
        sa = size(net.layers{n}.a{j});
        net.fv = [net.fv; reshape(net.layers{n}.a{j}, sa(1) * sa(2), sa(3))];
    end
    %  feedforward into output perceptrons
    net.o = sigm(net.ffW * net.fv + repmat(net.ffb, 1, size(net.fv, 2)));

end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function linplus_test14 ( )

%*****************************************************************************80
%
%% TEST14 tests R8BB_FA, R8BB_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n1 = 0;
  n2 = 10;
  n = n1 + n2;
  ml = 0;
  mu = 0;
  na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST14\n' );
  fprintf ( 1, '  For a border banded matrix:\n' );
  fprintf ( 1, '  R8BB_FA factors;\n' );
  fprintf ( 1, '  R8BB_SL solves.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N     = %d\n', n );
  fprintf ( 1, '  Matrix suborder N1 = %d\n', n1 );
  fprintf ( 1, '  Matrix suborder N2 = %d\n', n2 );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );
%
%  Set the matrix.
%
  [ a, seed ] = r8bb_random ( n1, n2, ml, mu, seed );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r8bb_mxv ( n1, n2, ml, mu, a, x );
%
%  Factor the matrix.
%
  [ a_lu, pivot, info ] = r8bb_fa ( n1, n2, ml, mu, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST14 - Fatal error!\n' );
    fprintf ( 1, '  R8BB_FA claims the matrix is singular.\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Solve the system.
%
  x = r8bb_sl ( n1, n2, ml, mu, a_lu, pivot, b );

  r8vec_print ( n, x, '  Solution:' );

  return
end


text:
%CALCOPTICALFLOWPYRLK  Calculates an optical flow for a sparse feature set using the iterative Lucas-Kanade method with pyramids
%
%     nextPts = cv.calcOpticalFlowPyrLK(prevImg, nextImg, prevPts)
%     [nextPts, status, err] = cv.calcOpticalFlowPyrLK(...)
%     [...] = cv.calcOpticalFlowPyrLK(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __prevImg__ First 8-bit single-channel or 3-channel input image, or
%   pyramid constructed by cv.buildOpticalFlowPyramid
% * __nextImg__ Second input image or pyramid of the same size and the same
%   type as `prevImg`.
% * __prevPts__ Vector of 2D points for which the flow needs to be found;
%   point coordinates must be single-precision floating-point numbers.
%   Cell-array of 2-element vectors is accepted: `{[x,y], ...}`
%
% ## Output
% * __nextPts__ Output vector of 2D points (with single-precision
%   floating-point coordinates) containing the calculated new positions of
%   input features in the second image. A cell array of 2-elements vectors of
%   the same size as `prevPts`.
% * __status__ Output status vector. Each element of the vector is set to 1 if
%   the flow for the corresponding features has been found. Otherwise, it is
%   set to 0.
% * __err__ Output vector of errors; each element of the vector is set to an
%   error for the corresponding feature, type of the error measure is
%   determined by `GetMinEigenvals` option; if the flow wasn't found then the
%   error is not defined (use `status` to find such cases).
%
% ## Options
% * __InitialFlow__ Vector of 2D points to be used for the initial estimate of
%   `nextPts`. If not specified, `prevPts` will be used as an initial
%   estimate. The vector must have the same size as in the input. Not set by
%   default.
% * __WinSize__ Size of the search window at each pyramid level. Default to
%   [21, 21].
% * __MaxLevel__ 0-based maximal pyramid level number. If set to 0, pyramids
%   are not used (single level). If set to 1, two levels are used, and so on;
%   if pyramids are passed to input then algorithm will use as many levels as
%   pyramids have but no more than `MaxLevel`. Default to 3.
% * __Criteria__ Parameter specifying the termination criteria of the
%   iterative search algorithm (after the specified maximum number of
%   iterations `Criteria.maxCount` or when the search window moves by less
%   than `Criteria.epsilon`. Struct with `{'type','maxCount','epsilon'}`
%   fields is accepted. The type field should have one of 'Count', 'EPS', or
%   'Count+EPS' to indicate which criteria to use. Default to
%   `struct('type','Count+EPS', 'maxCount',30, 'epsilon',0.01)`.
% * __GetMinEigenvals__ Use minimum eigen values as an error measure (see
%   `MinEigThreshold` description); if the flag is not set, then L1 distance
%   between patches around the original and a moved point, divided by number
%   of pixels in a window, is used as a error measure. Default to false.
% * __MinEigThreshold__ The algorithm calculates the minimum eigen value of a
%   2x2 normal matrix of optical flow equations (this matrix is called a
%   spatial gradient matrix in [Bouguet00]), divided by number of pixels in a
%   window; if this value is less than `MinEigThreshold`, then a corresponding
%   feature is filtered out and its flow is not processed, so it allows to
%   remove bad points and get a performance boost. Default to 1e-4.
%
% The function implements a sparse iterative version of the Lucas-Kanade
% optical flow in pyramids. See [Bouguet00]. The function is parallelized with
% the TBB library.
%
% ## References
% [Bouguet00]:
% > Jean-Yves Bouguet. "Pyramidal Implementation of the Lucas Kanade Feature
% > Tracker. Description of the algorithm", Intel Corporation, 2001.
%
% ## Example
%
%     prevIm = rgb2gray(imread('prev.jpg'));
%     nextIm = rgb2gray(imread('next.jpg'));
%     prevPts = cv.goodFeaturesToTrack(prevIm);
%     nextPts = cv.calcOpticalFlowPyrLK(prevIm, nextIm, prevPts);
%
% See also: cv.SparsePyrLKOpticalFlow, cv.calcOpticalFlowFarneback,
%  cv.goodFeaturesToTrack, vision.PointTracker
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Matlab Stuff for IE 753: Material Handling.
%
% This folder has IE753-specific Matlab stuff for the current semester.
%


text:
function [f] = elec_ellipse_fit_optim(r, X, Y, Z, xo, yo, zo)

% elec_ellipse_fit_optim - Optimization for elec_ellipse_fit.m
%
% Called from elec_ellipse_fit.m
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no implied or express warranties
% History:  02/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% r is a 3x1 vector of radius values for each x,y,z axis component of
% ellipse
%
% equation of ellipsoid with center (xo,yo,zo) and radius for each axis (x,y,z) = (a,b,c):
% (( x - xo )^2 / a^2) + (( y - yo )^2 / b^2) + (( z - zo )^2 / c^2) = 1
%
% This function below creates a scalar value to
% return to the fminsearch function in elec_ellipse_fit.

E = ( (X-xo).^2 )/r(1).^2 + ((Y-yo).^2)/r(2).^2 + ((Z-zo).^2)/r(3).^2  - 1;

f = sum( E .* E );  % sum of squares returned

return

author: PatternRecognition

repo: OpenBMI


text:
clear;
clc;
% This program solves the economic emission dispatch dispatch with Bmn coefficients by
% quadratic programming and equal incremental cost critetion
% the elddata matrix should have 5 columns of fuel cost coefficients and plant  limits.
% 1.a ($/MW^2) 2. b $/MW 3. c ($) 4.lower lomit(MW) 5.Upper limit(MW)
%no of rows denote the no of plants(n)
elddata=[0.15247	38.53973	756.79886	10	125
0.10587	46.15916	451.32513	10	150
0.02803	40.3965	1049.9977	35	225
0.03546	38.30553	1243.5311	35	210
0.02111	36.32782	1658.5596	130	325
0.01799	38.27041	1356.6592	125	315
];
% the emidata matrix should have 3 columns of fuel cost coefficients and plant  limits.
% 1.a (Kg/MW^2) 2. b Kg/MW 3. c (Kg)

emidata=[0.00419	0.32767	13.85932
0.00419	0.32767	13.85932
0.00683	-0.54551	40.2669
0.00683	-0.54551	40.2669
0.00461	-0.51116	42.89553
0.00461	-0.51116	42.8955];
% h1 1nd h2 are  is the weightage factor for economy and emission
h1=1; h2=44.788;
% Demand (MW)
Pd=700;
% Loss coefficients it should be squarematrix of size nXn where n is the no
% of plants
B=1e-4*[1.4	.17	.15	.19	.26	.22
.17	.6	.13	.16	.15	.2
.15	.13	.65	.17	.24	.19
.19	.16	.17	.71	.3	.25
.26	.15	.24	.3	.69	.32
.22	.2	.19	.25	.32	.85];
[P Fcost Emi Pl]=emield(elddata,emidata,h1,h2,B,Pd)


text:
function y = linef(lambda, fn, x, d, varargin)
%LINEF	Calculate function value along a line.
%
%	Description
%	LINEF(LAMBDA, FN, X, D) calculates the value of the function FN at
%	the point X+LAMBDA*D.  Here X is a row vector and LAMBDA is a scalar.
%
%	LINEF(LAMBDA, FN, X, D, P1, P2, ...) allows additional arguments to
%	be passed to FN().   This function is used for convenience in some of
%	the optimisation routines.
%
%	See also
%	GRADCHEK, LINEMIN
%

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

% Check function string
fn = fcnchk(fn, length(varargin));

y = feval(fn, x+lambda.*d, varargin{:});

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt