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

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

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

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

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

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

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

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

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

h = get(gca, 'children');

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

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

for i = 1:barnum
    if ~doerrbar_width
        errwidth = [-.03 .03];
    end

    if ~verLessThan('matlab', '8.4') % HG2
        error('this function won''t work properly versions later 2014b; please use bar_wani_2016.m instead');
    else
        hh = get(h(i), 'children');
        xdata = get(hh(2), 'xData');
    end

    for j = 1:(length(xdata)/9)
        xdata(9*(j-1)+4:9*(j-1)+5) = xdata(9*(j-1)+1:9*(j-1)+2)+errwidth ;
        xdata(9*(j-1)+7:9*(j-1)+8) = xdata(9*(j-1)+1:9*(j-1)+2)+errwidth ;
        if doast
            text(double(xdata(9*(j-1)+1)+ast_adj_x), double(ast_loc(j,barnum-i+1)), ast{j,barnum-i+1}, 'fontsize', 18, 'HorizontalAlignment','center'); % **ADJUST** font size
            if dotext
                text_y = ymax*1.15 + text_adj_y;
                text(xdata(9*(j-1)+1), text_y, num2str(text_bar(j, (barnum-i+1))), 'fontsize', 20, 'HorizontalAlignment','center'); % **ADJUST** font size
            end
        end

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

        xdata_all(i,j) = xdata(9*(j-1)+1);

        if grnum == 1
            break
        end

    end

    set(hh(2), 'xdata', xdata);
    set(h(i), 'lineWidth', 2);
end

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

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

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

xtickdata = sort(xdata_all(:));

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

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

xdata_all = flipud(xdata_all);

end

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

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

ast = cell(size(p));

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

end

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

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

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

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

change_axis = 0;
ymax = 0;

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

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

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

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

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

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

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

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

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

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

	handles.ax = gca;

	hold off
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function lagrange_basis_display_test ( )

%*****************************************************************************80
%
%% LAGRANGE_BASIS_DISPLAY_TEST tests the LAGRANGE_BASIS_DISPLAY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 August 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_BASIS_DISPLAY_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LAGRANGE_BASIS_DISPLAY_TEST library.\n' );
%
%  Plot the 11 standard "basis vectors".
%
  a = 0.0;
  b = 1.0;
  m = 10;

  standard_basis_display ( a, b, m )
%
%  Plot the 11 Lagrange "basis vectors" for evenly spaced nodes.
%
  a = 0.0;
  b = 1.0;
  m = 10;
  xd = linspace ( a, b, m + 1 )

  lagrange_basis_display ( a, b, m, xd )
%
%  Plot the 11 Lagrange "basis vectors" for Chebyshev-spaced nodes.
%
  a = -1.0;
  b = 1.0;
  m = 10;
  xd(1:m+1) = cos ( pi * ( 2 * ( m + 1 ) + 1 - 2 * (1:m+1) ) / ( 2 * ( m + 1 ) ) );

  lagrange_basis_display ( a, b, m, xd )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_BASIS_DISPLAY_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = double(t)
%DOUBLE Convert tenmat to double array.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

a = t.data;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [rc,fval,it] = ARhat_misd(ng,xg,rcinit,rc0,lag_max)

%ARHAT_MISD AR model from measurements with missing data
%  [rc,fval,it] = ARhat_misd(ng,xg,rcinit,rc0,lag_max) estimates
%  reflection coefficients from missing data (ng,xg). ng contains
%  measurement times and xg the corresponding measurements.
%
%  The starting value for the reflection coefficients is given by
%  rc_start = [1 rc0 rcinit]
%  An AR model is estimated using approximate Maximum Likelihood (ML)
%  estimation, where rc0 is fixed and the ML is sought over the last
%  reflection coefficients, starting in rcinit.
%
%  See also ARMLFIT, ARHAT_MISD.

nseg = length(ng);
n_obs = 0;
for seg = 1:nseg
   n_obs = n_obs+length(ng{seg});
end

opties = optimset('Display','iter','TolX',.001/sqrt(n_obs),'TolFun',.0001);

[rc_tan,fval,exitflag,output]= fminunc('ARMLfit',tan(.5*pi*rcinit),opties,ng,xg,rc0,lag_max);
rc = 2/pi*atan(rc_tan);
it = output.iterations;

author: Sable

repo: mcbench-benchmarks


text:
function DBExportGallery(name,query,varargin)

%DBExportGallery - Create figure gallery from current database.
%
% Select figures and create an HTML gallery (although not required, installing
% export_fig will yield better PNG images).
%
%  USAGE
%
%    DBExportGallery(name,query,<options>)
%
%    name           gallery name
%    query          optional figure list query (WHERE clause; see Example)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'size'        horizontal and vertical thumbnail size
%                   (default = [240 240])
%     'nColumns'    number of columns (default = 3)
%     'info'        display detailed info with thumbnails (default = 'off')
%     'code'        create generation code pages (default = 'off')
%    =========================================================================
%
%  EXAMPLES
%
%    % Export entire gallery with extensive information
%    DBExportGallery('Complete Gallery','info','on');
%
%    % Export only figures for the session named 'SESSION1'
%    DBExportGallery('Partial Gallery','eid="SESSION1"');
%
%    % Export figures for all sessions with names starting with 'SLEEP'
%    DBExportGallery('Sleep Gallery','eid like "SLEEP%"');

% Copyright (C) 2007-2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Make sure MyM is installed and functional
CheckMyM;

% Default values
hSize = 240;
vSize = 240;
nColumns = 3;
showInfo = 'off';
showCode = 'off';

% Check number of parameters
if nargin < 1,
	error('Incorrect number of parameters (type ''help DBExportGallery'' for details).');
end

% Optional query provided?
if nargin == 1,
	query = '';
elseif isstring_FMAT(query,'size','ncolumns','info','code'),
	varargin = {query varargin{:}};
	query = '';
end

% Edit query
query = strtrim(query);
query = regexprep(query,'^where','');
if ~isempty(query), query = [' where ' query]; end

% Parse parameter list
for j = 1:2:length(varargin),
	if ~ischar(varargin{j}),
		error(['Parameter ' num2str(j+7) ' is not a property (type ''help DBExportGallery'' for details).']);
	end
	switch(lower(varargin{j})),
		case 'size',
			sizes = varargin{j+1};
			if ~isivector(sizes,'>0','#2'),
				error('Incorrect value for ''size'' (type ''help DBExportGallery'' for details).');
			end
			hSize = sizes(1);
			vSize = sizes(2);

		case 'ncolumns',
			nColumns = varargin{j+1};
			if ~isiscalar(nColumns,'>0'),
				error('Incorrect value for ''ncolumns'' (type ''help DBExportGallery'' for details).');
			end

		case 'info',
			showInfo = lower(varargin{j+1});
			if ~isstring_FMAT(showInfo,'on','off'),
				error('Incorrect value for ''info'' (type ''help DBExportGallery'' for details).');
			end

		case 'code',
			showCode = lower(varargin{j+1});
			if ~isstring_FMAT(showCode,'on','off'),
				error('Incorrect value for ''code'' (type ''help DBExportGallery'' for details).');
			end

		otherwise,
			error(['Unknown property ''' num2str(varargin{j}) ''' (type ''help DBExportGallery'' for details).']);

	end
end

showInfo = strcmp(showInfo,'on');
showCode = strcmp(showCode,'on');

% Edit query
query = strtrim(query);
fullQuery = regexprep(query,'^where','');
if ~isempty(fullQuery), fullQuery = [' where ' fullQuery]; end

% Query database
f = mym(['select png,eid,name,parameters,comments,mfiles,code,date,user from figures' fullQuery]);
if isempty(f.eid),
	warning(['No figures match ''' query '''.']);
end

% Reformat code
code = {};
for i = 1:length(f.code),
	for j = 1:length(f.code{i}),
		code{i}{j} = char(f.code{i}{j})';
	end
end
f.code = code;

% Correct eids and names: replace / with _ to avoid conflicts with filesystem
nFigures = length(f.png);
for i = 1:nFigures,
	f.eid{i} = strrep(f.eid{i},'/','_');
        f.name{i} = strrep(f.name{i},'/','_');
end

% Create root directory
try
	mkdir(name);
	mkdir([name '/images']);
	mkdir([name '/thumbs']);
	if showCode,
		mkdir([name '/code']);
	end
catch
	error('Cannot create gallery (check file access permissions).');
end
[path,name] = fileparts(name);

% Export png images + thumbnails (and code pages if required)
for i = 1:nFigures,
	% png image
	figureName{i} = [f.eid{i} '-' f.name{i} '.png'];
	file = fopen([path name '/images/' figureName{i}],'wb');
	if file == -1,
		error(['Could not create figure ' figureName{i} '.']);
	end
	fwrite(file,f.png{i});
	fclose(file);

	% Thumbnail
	pngData = imread([path name '/images/' figureName{i}],'png');
   [height,width,nz] = size(pngData);
   hScale = width/hSize;
   vScale = height/vSize;
	x = floor(1:hScale:width);
	y = floor(1:vScale:height);
   thumbnail = pngData(y,x,:);
   imwrite(thumbnail,[path name '/thumbs/' figureName{i}],'png');

	% Code
	if showCode && ~isempty(f.code) && ~isempty(f.code{i}),
		file = fopen([path name '/code/' figureName{i}(1:end-4) '.m'],'w');
		if file == -1,
			error('Could not create code output.');
		end
		nFunctions = size(f.code{i},2);
		for j = 1:nFunctions,
			if j > 1, fprintf(file,'\n\n'); end
			fprintf(file,'%% =======================================================================================================\n');
			fprintf(file,['%%     ' f.mfiles{i}{j} '\n']);
			fprintf(file,'%% =======================================================================================================\n\n\n');
			fprintf(file,'%s',f.code{i}{j});
		end
%  		if isempty(j), j = 1; end
%  		fprintf(file,'%s',f.code{i}{j});
		fclose(file);
	end
end

% Output HTML code
% Header...
fid = fopen([path name '/index.html'],'w');
fprintf(fid,'<?xml version="1.0" encoding="UTF-8" ?>\n');
fprintf(fid,'<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "http://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd">\n');
fprintf(fid,'\n');
fprintf(fid,' \n');
fprintf(fid,'  <title>%s</title>\n',name);
fprintf(fid,'  <meta http-equiv="content-type" content="text/html; charset=UTF-8"/>\n');
fprintf(fid,'  <meta name="GENERATOR" content="FMAToolbox"/>\n');
fprintf(fid,'  \n');
fprintf(fid,' \n');
fprintf(fid,' \n');
fprintf(fid,'
%s
\n',name);
fprintf(fid,'
%d images (%s)

\n',nFigures,datestr(now));
fprintf(fid,'
\n');
fprintf(fid,'
\n');

% ... figure table...
nLines = ceil(nFigures/nColumns);
i = 1;
for line = 1:nLines,
	fprintf(fid,'
\n');
	for column = 1:nColumns,
		if i > nFigures, break; end
		fprintf(fid,'    \n');
		fprintf(fid,'     \n',figureName{i},figureName{i},hSize,vSize);
		fprintf(fid,'
\n');
		fprintf(fid,'
\n');
		fprintf(fid,'       \n');
		fprintf(fid,'        EID\n');
		fprintf(fid,'       	\n');
		fprintf(fid,'       \n');
		fprintf(fid,'
%s
\n',f.eid{i});
		fprintf(fid,'       	\n');
		fprintf(fid,'
\n');
		fprintf(fid,'
\n');
		fprintf(fid,'       \n');
		fprintf(fid,'        Name\n');
		fprintf(fid,'       	\n');
		fprintf(fid,'       \n');
		fprintf(fid,'
%s
\n',f.name{i});
		fprintf(fid,'       	\n');
		fprintf(fid,'
\n');
		if showInfo,
			fprintf(fid,'
\n');
			fprintf(fid,'       \n');
			fprintf(fid,'        Comments\n');
			fprintf(fid,'       	\n');
			fprintf(fid,'       \n');
			fprintf(fid,'
%s
\n',f.comments{i});
			fprintf(fid,'       	\n');
			fprintf(fid,'
\n');
			fprintf(fid,'
\n');
			fprintf(fid,'       \n');
			fprintf(fid,'        Parameters\n');
			fprintf(fid,'       	\n');
			fprintf(fid,'       \n');
			fprintf(fid,'
%s
\n',f.parameters{i});
			fprintf(fid,'       	\n');
			fprintf(fid,'
\n');
		end
		if showCode && ~isempty(f.code) && ~isempty(f.code{i}),
			codeName = [figureName{i}(1:end-4) '.m'];
			fprintf(fid,'
\n');
			fprintf(fid,'       \n');
			fprintf(fid,'        Code\n');
			fprintf(fid,'       	\n');
			fprintf(fid,'       \n');
			fprintf(fid,'        [M-File]\n',codeName);
			fprintf(fid,'       	\n');
			fprintf(fid,'
\n');
		end
		if showInfo,
			fprintf(fid,'
\n');
			fprintf(fid,'       \n');
			fprintf(fid,'        Date\n');
			fprintf(fid,'       	\n');
			fprintf(fid,'       \n');
			fprintf(fid,'
%s
\n',f.date{i});
			fprintf(fid,'       	\n');
			fprintf(fid,'
\n');
			fprintf(fid,'
\n');
			fprintf(fid,'       \n');
			fprintf(fid,'        User\n');
			fprintf(fid,'       	\n');
			fprintf(fid,'       \n');
			fprintf(fid,'
%s
\n',f.user{i});
			fprintf(fid,'       	\n');
			fprintf(fid,'
\n');
		end
		fprintf(fid,'
\n');
		fprintf(fid,'    	\n');
		i = i+1;
	end
	fprintf(fid,'
\n');
end

% ... and footer
fprintf(fid,'
\n');
fprintf(fid,' \n');
fprintf(fid,'\n');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [nav_e] = ins_visual(imu, gnss, visual, att_mode)
% ins_visual: loosely-coupled integrated navigation system.
%
% ins_visual integrates IMU and visual measurements by using an Extended Kalman filter.
%
% INPUT
%   imu, IMU data structure.
%         t: Ix1 time vector (seconds).
%        fb: Ix3 accelerations vector in body frame XYZ (m/s^2).
%        wb: Ix3 turn rates vector in body frame XYZ (radians/s).
%       arw: 1x3 angle random walks (rad/s/root-Hz).
%       vrw: 1x3 velocity random walks (m/s^2/root-Hz).
%     g_std: 1x3 gyros standard deviations (radians/s).
%     a_std: 1x3 accrs standard deviations (m/s^2).
%    gb_sta: 1x3 gyros static biases or turn-on biases (radians/s).
%    ab_sta: 1x3 accrs static biases or turn-on biases (m/s^2).
%    gb_dyn: 1x3 gyros dynamic biases or bias instabilities (radians/s).
%    ab_dyn: 1x3 accrs dynamic biases or bias instabilities (m/s^2).
%   gb_corr: 1x3 gyros correlation times (seconds).
%   ab_corr: 1x3 accrs correlation times (seconds).
%    gb_psd: 1x3 gyros dynamic biases root-PSD (rad/s/root-Hz).
%    ab_psd: 1x3 accrs dynamic biases root-PSD (m/s^2/root-Hz);
%      freq: 1x1 sampling frequency (Hz).
% ini_align: 1x3 initial attitude at t(1).
% ini_align_err: 1x3 initial attitude errors at t(1).
%
%	gnss, GNSS data structure.
%         t: Gx1 time vector (seconds).
%       lat: Gx1 latitude (radians).
%       lon: Gx1 longitude (radians).
%         h: Gx1 altitude (m).
%       vel: Gx3 NED velocities (m/s).
%       std: 1x3 position standard deviations (rad, rad, m).
%      stdm: 1x3 position standard deviations (m, m, m).
%      stdv: 1x3 velocity standard deviations (m/s).
%      larm: 3x1 lever arm from IMU to GNSS antenna (x-fwd, y-right, z-down) (m).
%      freq: 1x1 sampling frequency (Hz).
%   zupt_th: 1x1 ZUPT threshold (m/s).
%  zupt_win: 1x1 ZUPT time window (seconds).
%       eps: 1x1 time interval to compare current IMU time to current GNSS time vector (s).
%
%   visual, visual data structure
%         t: Vx1 time vector (seconds)
%       lat: Vx1 latitude (radians).
%       lon: Vx1 longitude (radians).
%         h: Vx1 altitude (m).
%       vel: Vx3 NED velocities (m/s).
%     covvm: Vx36 velocity and position covariance matrices (m^2 and m^2/s^2).
%      covv: Vx9 velocity covariance matrices (m^2/s^2).
%      larm: 3x1 lever arm (m).
%      freq: 1x1 sampling frequency (Hz).
%       eps: 1x1 time interval to compare current IMU time to current visual time vector (s).
%
%  att_mode: attitude mode string.
%      'quaternion': attitude updated in quaternion format. Default value.
%             'dcm': attitude updated in Direct Cosine Matrix format.
%
% OUTPUT
%   nav_e, INS/visual navigation estimates data structure.
%         t: Ix1 INS time vector (seconds).
%        tg: Vx1 VISUAL time vector, when Kalman filter was executed (seconds).
%      roll: Ix1 roll (radians).
%     pitch: Ix1 pitch (radians).
%       yaw: Ix1 yaw (radians).
%       vel: Ix3 NED velocities (m/s).
%       lat: Ix1 latitude (radians).
%       lon: Ix1 longitude (radians).
%         h: Ix1 altitude (m).
%        xi: Vxn Kalman filter a priori states.
%        xp: Vxn Kalman filter a posteriori states.
%         z: Vxr INS/visual measurements
%         v: Vxr Kalman filter innovations.
%         b: Vxr Kalman filter biases compensations, [gb_dyn ab_dyn].
%         A: Vxn^2 Kalman filter transition-state matrices, one matrix per
%            row ordered by columns.
%        Pp: Vxn^2 Kalman filter a posteriori covariance matrices, one
%         matrix per row ordered by columns.
%        Pi: Vxn^2 Kalman filter a priori covariance matrices, one matrix
%            per row ordered by columns.
%         K: Vx(n*r) Kalman gain matrices
%         S: Vxr^2 Innovation matrices
%        ob: Vx1 Number of observable states after each visual data arriving
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics, vol. 17,
% issue 2, pp. 110-120, 2015. Alg. 2.
%
%   Groves, P.D. (2013), Principles of GNSS, Inertial, and
% Multisensor Integrated Navigation Systems (2nd Ed.). Artech House.
%
%   ZUPT algothim based on Groves, Chapter 15, "INS Alignment, Zero Updates,
% and Motion Constraints".
%
%   ins_gps.m, ins_gnss function is based on that previous NaveGo function.
%
% Version: 002
% Date:    2022/07/19
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

if nargin < 3, att_mode  = 'quaternion'; end

%% ZUPT ALGORITHM

zupt_flag = false;

%% PREALLOCATION

% Kalman filter dimensions
n = 15; % number of states
r = 6;  % number of sensors
q = 12; % number of inputs

% Constant matrices
I = eye(3);
O = zeros(3);

% Length of INS time vector
LI = length(imu.t);

% Length of visual time vector
LV = length(visual.t);

% Preallocation of attitude vectors
roll_e  = zeros (LI, 1);
pitch_e = zeros (LI, 1);
yaw_e   = zeros (LI, 1);

% Preallocation of velocity vector
vel_e   = zeros (LI, 3);

% Preallocation of gravity vector
gn_e   = zeros (LI, 3);

% Preallocation of position vectors
lat_e    = zeros (LI, 1);
lon_e    = zeros (LI, 1);
h_e      = zeros (LI, 1);

% Preallocation of Kalman filter matrices for later performance analysis
xi = zeros(LV, n);      % Evolution of Kalman filter a priori states
xp = zeros(LV, n);      % Evolution of Kalman filter a posteriori states
z = zeros(LV, r);       % INS/visual measurements
v = zeros(LV, r);       % Kalman filter innovations

A  = zeros(LV, n^2);    % Transition-state matrices
Pi = zeros(LV, n^2);    % A priori covariance matrices
Pp = zeros(LV, n^2);    % A posteriori covariance matrices
K  = zeros(LV, n*r);    % Kalman gain matrices
S  = zeros(LV, r^2);    % Innovation matrices
ob = zeros(LV, 1);      % Number of observable states at each visual data arriving

b = zeros(LV, r);       % Biases compensantions after Kalman filter correction

%% INITIAL VALUES AT INS TIME = 1

% Initial attitude
roll_e(1)  = imu.ini_align(1);
pitch_e(1) = imu.ini_align(2);
yaw_e(1)   = imu.ini_align(3);
DCMnb = euler2dcm([roll_e(1); pitch_e(1); yaw_e(1);]);
DCMbn = DCMnb';
qua   = euler2qua([roll_e(1) pitch_e(1) yaw_e(1)]);

% Initial velocity
vel_e(1,:) = gnss.vel(1,:);

% Initial position
lat_e(1) = gnss.lat(1);
lon_e(1) = gnss.lon(1);
h_e(1)   = gnss.h(1);

% Initial dynamic biases
gb_dyn = imu.gb_dyn';
ab_dyn = imu.ab_dyn';

% Turn-rates update with both updated velocity and position
omega_ie_n = earth_rate(lat_e(1));
omega_en_n = transport_rate(lat_e(1), vel_e(1,1), vel_e(1,2), h_e(1));

% Gravity update
gn_e(1,:) = gravity(lat_e(1), h_e(1));

%% INITIALIZATION OF KALMAN FILTER MATRICES

% Prior estimates
kf.xi = [ zeros(1,9), imu.gb_dyn, imu.ab_dyn ]';  % Error vector state
kf.Pi = diag([imu.ini_align_err, gnss.stdv, gnss.std, imu.gb_dyn, imu.ab_dyn].^2);

kf.Q  = diag([imu.arw, imu.vrw, imu.gb_psd, imu.ab_psd].^2);

fb_corrected = imu.fb(1,:)' - ab_dyn - imu.ab_sta';
fn = DCMbn * fb_corrected;
wn_corrected = imu.wb(1,:)' - gb_dyn - imu.gb_sta';
wn = DCMbn * wn_corrected;

% Vector to update matrix F
upd = [gnss.vel(1,:) gnss.lat(1) gnss.h(1) fn' wn'];

% Update matrices F and G
[kf.F, kf.G] = F_update(upd, DCMbn, imu);

[RM,RN] = radius(gnss.lat(1));
Tpr = diag([(RM + gnss.h(1)), (RN + gnss.h(1)) * cos(gnss.lat(1)), -1]);  % radians-to-meters

% Update matrix H
kf.H = [ O I O O O ;
    O O Tpr O O ; ];
kf.R = diag([gnss.stdv gnss.stdm]).^2;
kf.z = [ gnss.stdv, gnss.stdm ]';

% Propagate prior estimates to get xp(1) and Pp(1)
kf = kf_update( kf );

% Initial matrices for Kalman filter performance analysis
xi(1,:) = kf.xi';
xp(1,:) = kf.xp';
Pi(1,:) = reshape(kf.Pi, 1, n^2);
Pp(1,:) = reshape(kf.Pp, 1, n^2);
K(1,:)  = reshape(kf.K, 1, n*r);
S(1,:)  = reshape(kf.S, 1, r^2);
v(1,:)  = kf.v';
z(1,:)  = kf.z';
b(1,:) = [gb_dyn', ab_dyn'];

%% INS (IMU) TIME IS THE MASTER CLOCK
for i = 2:LI

    %% INERTIAL NAVIGATION SYSTEM (INS)

    % Print a dot on console every 10,000 INS executions
    if (mod(i,10000) == 0), fprintf('. ');  end
    % Print a return on console every 200,000 INS executions
    if (mod(i,200000) == 0), fprintf('\n'); end

    % IMU sampling interval
    dti = imu.t(i) - imu.t(i-1);

    % Inertial sensors corrected with a posteriori KF biases estimation and
    % deterministic static biases
    wb_corrected = imu.wb(i,:)' - gb_dyn - imu.gb_sta';
    fb_corrected = imu.fb(i,:)' - ab_dyn - imu.ab_sta';
    fn = DCMbn * fb_corrected;
    wn = DCMbn * wb_corrected;

    % Velocity update
    vel = vel_update(fn, vel_e(i-1,:), omega_ie_n, omega_en_n, gn_e(i,:)', dti);
    vel_e (i,:) = vel;

    % Position update
    pos = pos_update([lat_e(i-1) lon_e(i-1) h_e(i-1)], vel_e(i,:), dti);
    lat_e(i) = pos(1);
    lon_e(i) = pos(2);
    h_e(i)   = pos(3);

    % Turn-rates update with both updated velocity and position
    omega_ie_n = earth_rate(lat_e(i));
    omega_en_n = transport_rate(lat_e(i), vel_e(i,1), vel_e(i,2), h_e(i));

    % Gravity update
    gn_e(i,:) = gravity(lat_e(i), h_e(i));

    % Attitude update
    [qua, DCMbn, euler] = att_update(wb_corrected, DCMbn, qua, ...
        omega_ie_n, omega_en_n, dti, att_mode);
    roll_e(i) = euler(1);
    pitch_e(i)= euler(2);
    yaw_e(i)  = euler(3);

    %% ZUPT DETECTION ALGORITHM
    idz = floor( gnss.zupt_win / dti ); % Index to set ZUPT window time

    if ( i > idz )

        % Mean velocity value for the ZUPT window time
        vel_m = mean (vel_e(i-idz:i , :));

        % If mean velocity value is under the ZUPT threshold velocity...
        if (abs(vel_m) < gnss.zupt_th)
            disp("ZUPT condition applied");

            % Current attitude is equal to the mean of previous attitudes
            % inside the ZUPT window time
            roll_e(i)  = mean (roll_e(i-idz:i , :));
            pitch_e(i) = mean (pitch_e(i-idz:i , :));
            yaw_e(i)   = mean (yaw_e(i-idz:i , :));

            % Current position is equal to the mean of previous positions
            % inside the ZUPT window time
            lat_e(i) = mean (lat_e(i-idz:i , :));
            lon_e(i) = mean (lon_e(i-idz:i , :));
            h_e(i)   = mean (h_e(i-idz:i , :));

            % Alternative attitude ZUPT correction
            % roll_e(i)  = (roll_e(i-idz , :));
            % pitch_e(i) = (pitch_e(i-idz , :));
            % yaw_e(i)   = (yaw_e(i-idz, :));
            % lat_e(i) = (lat_e(i-idz:i , :));
            % lon_e(i) = (lon_e(i-idz:i , :));
            % h_e(i)   = (h_e(i-idz:i , :));

            zupt_flag = true;

            %             fprintf(' z\n')       % DEBUG
        end
    end

    %% KALMAN FILTER UPDATE

    % Check if there is a new visual measurement to process at current INS time
    gdx =  find (visual.t >= (imu.t(i) - visual.eps) & visual.t < (imu.t(i) + visual.eps));

    if (~isempty(gdx) && gdx(1) > 1)
        for index = 1:size(gdx,1)
            % disp("updated with visual measurements");

            %% MEASUREMENTS

            % Meridian and normal radii of curvature update
            [RM,RN] = radius(lat_e(i));

            % Radians-to-meters matrix
            Tpr = diag([(RM + h_e(i)), (RN + h_e(i)) * cos(lat_e(i)), -1]);

            % Position innovations in meters with lever arm correction
            zp = Tpr * ([lat_e(i); lon_e(i); h_e(i);] - [visual.lat(gdx(index)); visual.lon(gdx(index)); visual.h(gdx(index));]) ...
                + (DCMbn * visual.larm);

            % Velocity innovations with lever arm correction
            zv = (vel_e(i,:) - visual.vel(gdx(index),:) - ((omega_ie_n + omega_en_n) * (DCMbn * visual.larm))' ...
                + (DCMbn * skewm(wb_corrected) * visual.larm )' )';

            %% KALMAN FILTER

            % visual sampling interval
            dtg = visual.t(gdx(index)) - visual.t(gdx(index)-1);

            % Vector to update matrix F
            upd = [vel_e(i,:) lat_e(i) h_e(i) fn' wn'];

            % Matrices F and G update
            [kf.F, kf.G] = F_update(upd, DCMbn, imu);

            % Matrix H update
            if(zupt_flag == false)
                kf.H = [ O I O O O ;
                    O O Tpr O O ; ];
                kf.R = reshape(visual.covvm(gdx(index),:),[6,6])';
                kf.z = [ zv' zp' ]';
            else
                kf.H = [ O I O O O ; ];
                kf.R = reshape(visual.covv(gdx(index),:),[3,3])';
                kf.z = zv;
            end

            % a posteriori states are forced to be zero (error-state approach)
            kf.xp = zeros(n , 1);
            % Execution of the extended Kalman filter
            kf = kalman(kf, dtg);

            %% OBSERVABILITY

            % Number the observable states at current GNSS time
            ob(gdx) = rank(obsv(kf.F, kf.H));

            %% INS/GNSS CORRECTIONS

            % Quaternion correction
            qua_skew = -skewm(qua(1:3));    % According to Crassidis, qua_skew should be
            % positive, but if positive NaveGo diverges.
            % Crassidis, Eq. A.174a
            Xi = [qua(4)*eye(3) + qua_skew; -qua(1:3)'];

            % Crassidis, Eq. 7.34
            qua = qua + 0.5 .* Xi * kf.xp(1:3);
            qua = qua / norm(qua);          % Brute-force normalization

            % DCM correction
            DCMbn = qua2dcm(qua);

            % Attitude correction, method 1
            %         euler = qua2euler(qua);
            %         roll_e(i) = euler(1);
            %         pitch_e(i)= euler(2);
            %         yaw_e(i)  = euler(3);

            % Attitude correction, method 2
            roll_e(i)  = roll_e(i)  - kf.xp(1);
            pitch_e(i) = pitch_e(i) - kf.xp(2);
            yaw_e(i)   = yaw_e(i)   - kf.xp(3);

            % Velocity correction
            vel_e(i,1) = vel_e(i,1) - kf.xp(4);
            vel_e(i,2) = vel_e(i,2) - kf.xp(5);
            vel_e(i,3) = vel_e(i,3) - kf.xp(6);

            % Position correction
            lat_e(i) = lat_e(i) - kf.xp(7);
            lon_e(i) = lon_e(i) - kf.xp(8);
            h_e(i)   = h_e(i)   - kf.xp(9);

            % Biases estimation
            gb_dyn   = -kf.xp(10:12);
            ab_dyn   = -kf.xp(13:15);

            % Matrices for later Kalman filter performance analysis
            xi(gdx(index),:) = kf.xi';
            xp(gdx(index),:) = kf.xp';
            b(gdx(index),:) = [gb_dyn', ab_dyn'];
            A(gdx(index),:)  = reshape(kf.A,  1, n^2);
            Pi(gdx(index),:) = reshape(kf.Pi, 1, n^2);
            Pp(gdx(index),:) = reshape(kf.Pp, 1, n^2);

            if(zupt_flag == false)
                v(gdx(index),:)  = kf.v';
                z(gdx(index),:)  = kf.z';
                K(gdx(index),:)  = reshape(kf.K, 1, n*r);
                S(gdx(index),:)  = reshape(kf.S, 1, r^2);
            else
                zupt_flag = false;
                z(gdx(index),:)  = [ kf.z' 0 0 0 ]';
                v(gdx(index),:)  = [ kf.v' 0 0 0 ]';
                K(gdx(index),1:n*3) = reshape(kf.K, 1, n*3);
                S(gdx(index),1:9)  = reshape(kf.S, 1, 3^2);
            end
        end
    end

    if imu.t(i) > max(visual.t)
        break;
    end
end

%% Summary from INS/visual integration

nav_e.t     = imu.t(1:i, :);    % INS time vector
nav_e.tg    = visual.t;         % visual time vector, which is the time vector when the Kalman filter was executed
nav_e.roll  = roll_e(1:i, :);   % Roll
nav_e.pitch = pitch_e(1:i, :);  % Pitch
nav_e.yaw   = yaw_e(1:i, :);    % Yaw
nav_e.vel   = vel_e(1:i, :);    % NED velocities
nav_e.lat   = lat_e(1:i, :);    % Latitude
nav_e.lon   = lon_e(1:i, :);    % Longitude
nav_e.h     = h_e(1:i, :);      % Altitude
nav_e.gn    = gn_e(1:i, :);     % Gravity estimation in the nav-frame.

nav_e.xi    = xi;       % A priori states
nav_e.xp    = xp;       % A posteriori states
nav_e.z     = z;        % INS/visual measurements
nav_e.v     = v;        % Kalman filter innovations
nav_e.b     = b;        % Biases compensations

nav_e.A     = A;        % Transition matrices
nav_e.Pi    = Pi;       % A priori covariance matrices
nav_e.Pp    = Pp;       % A posteriori covariance matrices
nav_e.K     = K;        % Kalman gain matrices
nav_e.S     = S;        % Innovation matrices
nav_e.ob    = ob;       % Number of observable states after each visual data arriving

fprintf('\n');
end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
clc
clear
close all
addpath(genpath(pwd))

%  basic plotting
tmp_ = 5;
t1 = 0:pi/20:8*pi;
t2 = 8*pi:pi/20:16*pi;
y1_ = exp(-t1/tmp_ );
y2_ = exp(-t1/tmp_ ).*sin(tmp_ *t1);
t = [t1, t2];
y1 = [y1_, fliplr(y1_)];
y2 = [y2_, fliplr(y2_)];

figure;
plot(t, y2, 'Color', 'r', 'LineStyle', '-', 'LineWidth', 1.5)
hold on
plot(t, y1, 'Color', 'b', 'LineStyle', ':', 'LineWidth', 1.5)
plot(t, -y1, 'Color', 'b', 'LineStyle', ':','LineWidth', 1.5)
xlim([min(t), max(t)])

% add 2 zoomed zones
zp = BaseZoom();
zp.plot;
zp.plot;

author: iqiukp

repo: ZoomPlot-MATLAB


text:
function x = l1_rand(mu,sigma,M)
error('this function was just for testing')
nu = 1;
u = gamrnd(nu*1,2/nu,M);
mean_invu = mean(1./u)
x = normrnd(mu,sigma.*u,M);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [beta_gibbs sigma_gibbs]=dogibbs(It,Bu,Bcap,phicap,Scap,alphacap,alphatop,n,k)

% start iterations
for ii=1:It-Bu

% draw B from a matrix-variate student distribution with location Bcap, scale Scap and phicap and degrees of freedom alphatop
B=bear.matrixtdraw(Bcap,Scap,phicap,alphatop,k,n);

% then draw sigma from an inverse Wishart distribution with scale matrix Scap and degrees of freedom alphacap (step 3)
sigma=bear.iwdraw(Scap,alphacap);

% record values before starting next iteration
beta_gibbs(:,ii)=B(:);
sigma_gibbs(:,ii)=sigma(:);

% go for next iteration
end

author: european-central-bank

repo: BEAR-toolbox


text:
function test_failed=test_fbwarped_framebounds
test_failed = 0;
Ls = 44100;
fs = 44100;
fmax = fs/2;
bins = 1;
fac = [1,7/8,3/4,5/8,1/2];

eigstol = 1e-4;
eigsmaxit = 100;

pcgmaxit = 150;
pcgtol = 1e-4;

warpfun = cell(4,1);
invfun = cell(4,1);

% ERBlet warping
 warpfun{3} = @freqtoerb;
 invfun{3} = @erbtofreq;
% constant-Q warping
 warpfun{4} = @(x) 10*log(x);
 invfun{4} = @(x) exp(x/10);
% sqrt-warping
 warpfun{2} = @(x) sign(x).*((1+abs(x)).^(1/2)-1);
 invfun{2} = @(x) sign(x).*((1+abs(x)).^2-1);
% Linear warping
 warpfun{1} = @(x) x/100;
 invfun{1} = @(x) 100*x;

 fmin = [0,0,0,50];

A = zeros(4,length(fac));
B = ones(4,length(fac));
red = A;

for jj = 1:4
    [g,a,fc,L]=warpedfilters(warpfun{jj},invfun{jj},fs,fmin(jj),fmax,bins,Ls,'bwmul',1.5,'fractional','complex');

    gf=filterbankresponse(g,a,Ls); framebound_ratio = max(gf)/min(gf);
    disp(['Painless system frame bound ratio: ', num2str(framebound_ratio)]);

    for kk = 1:length(fac)
        %[jj,kk]
        atemp = a;
        idx = [2:length(fc)/2,length(fc)/2+2:length(fc)];
        atemp(idx,2) = ceil(atemp(idx,2).*fac(kk));
        red(jj,kk) = sum(atemp(:,2)./atemp(:,1));

        [g,asan]=filterbankwin(g,atemp,L,'normal');
        gtemp=comp_filterbank_pre(g,asan,L,10);
        gtemp{1}.H = gtemp{1}.H.*sqrt(fac(kk));
        gtemp{length(fc)/2+1}.H = gtemp{length(fc)/2+1}.H.*sqrt(fac(kk));

        F = frame('filterbank',gtemp,asan,numel(gtemp));
        [A(jj,kk),B(jj,kk)] = framebounds(F,Ls,'tol',eigstol,'pcgtol',pcgtol,'maxit',eigsmaxit,'pcgmaxit',pcgmaxit);

     end
end

disp('This is a ratio B/A. Rows - warping sunction, Cols - redundancy compared to te minimal painless case')
B./A

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%TRNORM Normalize an SO(3) or SE(3) matrix
%
% TRNORM(R) is guaranteed to be a proper orthogonal matrix rotation
% matrix (3x3) which is "close" to the input matrix R (3x3). If R
% = [N,O,A] the O and A vectors are made unit length and the normal vector
% is formed from N = O x A, and then we ensure that O and A are orthogonal
% by O = A x N.
%
% TRNORM(T) as above but the rotational submatrix of the homogeneous
% transformation T (4x4) is normalised while the translational part is
% unchanged.
%
% If R (3x3xK) or T (4x4xK) representing a sequence then the normalisation
% is performed on each of the K planes.
%
% Notes::
% - Only the direction of A (the z-axis) is unchanged.
% - Used to prevent finite word length arithmetic causing transforms to
%   become `unnormalized'.
% - There is no Toolbox function for SO(2) or SE(2).
%
% See also OA2TR, SO3.trnorm, SE3.trnorm.

%## 3d homogeneous

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function TR = trnorm(T)

    assert(ishomog(T) || isrot(T), 'SMTB:trnorm:badarg', 'expecting 3x3xN or 4x4xN hom xform');

    if ndims(T) == 3
        % recurse for transform sequence
        n = size(T);
        TR = zeros(n);
        for i=1:n(3)
            TR(:,:,i) = trnorm(T(:,:,i));
        end
        return
    end

    o = T(1:3,2); a = T(1:3,3);
    n = cross(o, a);         % N = O x A
    o = cross(a, n);         % O = A x N
    R = [unit(n) unit(o) unit(a)];

    if ishomog(T)
        TR = rt2tr( R, T(1:3,4) );
    elseif isrot(T)
        TR = R;
    end

author: petercorke

repo: spatialmath-matlab


text:
function [M,N,G] = logfmap(I,L,H)
% [M,N] = logfmap(I,L,H)
%     Return a maxtrix for premultiplying spectrograms to map
%     the rows into a log frequency space.
%     Output map covers bins L to H of input
%     L must be larger than 1, since the lowest bin of the FFT
%     (corresponding to 0 Hz) cannot be represented on a
%     log frequency axis.  Including bins close to 1 makes
%     the number of output rows exponentially larger.
%     N returns the recovery matrix such that N*M is approximately I
%     (for dimensions L to H).
%
% 2004-05-21 dpwe@ee.columbia.edu

% Convert base-1 indexing to base-0
L = L-1;
H = H-1;

ratio = (H-1)/H;
opr = round(log(L/H)/log(ratio));
%ibin = H*exp((opr-[1:opr])*log((H-1)/H));
ibin = L*exp([0:(opr-1)]*-log(ratio));

M = zeros(opr,I);

for i = 1:opr
  % Where do we sample this output bin?
  % Idea is to make them 1:1 at top, and progressively denser below
  % i.e. i = max -> bin = topbin, i = max-1 -> bin = topbin-1,
  % but general form is bin = A exp (i/B)
%  M(i,round(ibin(i))) = 1;
  tt = pi*([0:(I-1)]-ibin(i));
  M(i,:) = (sin(tt)+eps)./(tt+eps);
end

% Normalize rows, but only if they are boosted by the operation
%G = 1./max(1,diag(M'*M))';
%% Fixup gain in bottom bins
%G(1:find(G==min(G))) = min(G);

G = ones(1,I);
G(1:(H+1)) = [0:H]./H;

% Inverse is just transpose plus scaling
N = (M.*repmat(G,opr,1))';

author: posenhuang

repo: deeplearningsourceseparation


text:
function d=dcirc(p,xc,yc,r)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

d=sqrt((p(:,1)-xc).^2+(p(:,2)-yc).^2)-r;

repo: ifem


text:
function node_num = grid_q4_node_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_Q4_NODE_NUM counts the nodes in a grid of Q4 elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer NODE_NUM, the number of nodes in the grid.
%
  node_num = ( nelemx + 1 ) * ( nelemy + 1 );

  return
end


text:
function [F] = mci_ramsay_fx (x,U,P,M)
% State equation for Ramsay model
% FORMAT [F] = mci_ramsay_fx (x,U,P,M)
%
% x     State vector
%       x(1) Voltage variable
%       x(2) Recovery variable
% U     inputs
% P     vector of model parameters - 2 params only
% M     model
%
% F     dx/dt
%
% J Ramsay et al (2007) Parameter estimation for differential equations:
% a generalised smoothing approach. J Roy Stat Soc B, 69(5):741-796.
%
% See also section 10 (page 26) and contribution by W.Penny on page 75 of:
%
% Girolami and Calderhead (2011) Riemann manifold Langevin and Hamiltonian
% Monte Carlo methods. J Roy Stat Soc B,73(2):123-214.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_ramsay_fx.m 6548 2015-09-11 12:39:47Z will $

c=3;

P(1)=exp(P(1));
P(2)=exp(P(2));

F=zeros(2,1);
F(1)=c*(x(1)-(1/3)*x(1)^3+x(2));
F(2)=-(1/c)*(x(1)-P(1)+P(2)*x(2));

author: spm


text:
function t=Tc(v,fc)

% t=Tc(v,fc)
%
% Return the channel coherence time for a mobile moving v m/s and communicating
% over carrier frequenc fc in Hz.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

error(nargchk(2,2,nargin));
error(chk_param(v,'v','scalar','real','>=',0));
error(chk_param(fc,'fc','scalar','real','>=',0));

dop_spread=2*doppler_shift(fc,v);
t=1/4/dop_spread;

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function pole_cplx = cplxpole(zta,wn,w,flag)
%
% Utility function: CPLXPOLE
%
% The purpose of this function is to compute the magnitude response of a
% complex pole or zero.

% Author: Craig Borghesani
% Date: 8/8/94
% Revised:
% Copyright (c) 1999, Prentice-Hall

if length(flag)==1,
 s = sqrt(-1)*w(:)';
 ht = s.^2 + 2*zta*wn*s + wn^2;
else
 T = flag(2);
 z = exp(sqrt(-1)*w(:)'*T);
 a = zta*wn; b = wn*sqrt(1-zta^2);
 ht = z - 2*exp(-a*T)*cos(b*T) + exp(-2*a*T)./z;
end
zero=find(abs(ht)==0);
if length(zero), ht(zero)=ones(1,length(zero))*eps; end
pole_cplx = ht.^flag(1);


text:
function plot_vector(X,f,plt,Xerr,c,w)
% Function to plot a frequency dependent vector X. If error bars are specified in Xerr,
% it also plots them. Xerr can either contain upper and lower confidence intervals
% on X, or simply a theoretical confidence level (for the coherence). Used
% to plot the spectrum and coherency.
% Usage: plot_vector(X,f,plt,Xerr,c)
% Inputs:
% X: input vector as a function of frequency (f), see third argument
% f: f axis grid for plot. Default. [1:length(X)]
% plt: 'l' for log, 'n' for no log.
% Xerr: lower and upper confidence intervals for X1: lower/upper x f. Or
%       simply a single number specifying an f-independent confidence
%       level.
% c: controls the color of the plot - input 'b','g','r' etc. Default 'b'
% w: controls the width of the lines - input 1, 1.5, 2 etc

if nargin < 1; error('Need data'); end;
N=length(X);
if nargin < 2 || isempty(f);
    f=1:N;
end;
if length(f)~=N; error('frequencies and data have incompatible lengths'); end;
if nargin < 3 || isempty(plt) ;
    plt='l';
end;
if nargin < 4 || isempty(Xerr);
    Xerr=[];
end;
if nargin < 5 || isempty(c)
    c='b';
end;
if nargin < 6 || isempty(w);
    w=1;
end;

if strcmp(plt,'l');
    X=10*log10(X);
    if nargin >=4 & ~isempty(Xerr); Xerr=10*log10(Xerr); end;
end;

if nargin < 4 || isempty(Xerr);
    plot(f,X,c,'Linewidth',w);
else
    if length(Xerr)==1;
       plot(f,X,c);
       line(get(gca,'xlim'),[Xerr,Xerr],'Color',c,'LineStyle','--','Linewidth',w);
    elseif ~isempty(Xerr);
       plot(f,X,c);
       hold on; plot(f,Xerr(1,:),[c '--'],'Linewidth',w); plot(f,Xerr(2,:),[c '--'],'Linewidth',w);
    end
end
xlabel('f');
if strcmp(plt,'l'); ylabel('10*log10(X)'); else ylabel('X'); end;


text:
function vox=tal2vox(tal,VOL)
% converts from talairach coordinate to voxel coordinate
% based on variables from SPM.M (passed here for
% faster operation)
%
% :Example:
% ::
%
%    foo = tal2vox([-30 28 -30], VOL)

if(isfield(VOL, 'M'))
    M = VOL.M;
elseif(isfield(VOL, 'mat'))
    M = VOL.mat;
else
    error('Error in %s: VOL does not have an "M" or a "mat" field.', mfilename);
end

vox=[0 0 0];
vox(1)=(tal(1)-M(1,4))/M(1,1);
vox(2)=(tal(2)-M(2,4))/M(2,2);
vox(3)=(tal(3)-M(3,4))/M(3,3);

return


text:
function [x_best, psi_best, out] = CSDF(mapp, x0, options)
% CSDF is a derivative-free algorithm for solving systems of nonlinear
% equations :math:`f(x) = 0`, x in :math:`R^m`
% using the nonlinear unconstrained minimization :math:`\textrm{min}\ \psi(x) = 1/2 ||f(x)||^2` s.t. `x` in :math:`R^m`.
%
% USAGE:
%
%    [x_best, psi_best, out] = CSDF(mapp, x0, options)
%
% INPUTS:
%    mapp:        function handle provides `f(x)` and gradient `f(x)`
%    x0:          initial point
%    options:     structure including the parameteres of scheme
%
%                   * .MaxNumIter - maximum number of iterations
%                   * .MaxNumMapEval - maximum number of function evaluations
%                   * .TimeLimit - maximum running time
%                   * .epsilon - accuracy parameter
%                   * .x_opt - optimizer
%                   * .psi_opt - optimum
%                   * .sigma - strong duplomonotone parameter
%                   * .l - Lipschitz continuity constant of `f`
%                   * .tauBar - a constant for determining the step-size
%                   * .flag_x_error - 1: saves :math:`x_{error}`, 0: do not saves :math:`x_{error}` (default)
%                   * .flag_psi_error - 1:saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                   * .flag_time - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                   * .Stopping_Crit - stopping criterion:
%
%                     1. stop if :math:`||nfxk|| \leq \epsilon`
%                     2. stop if `MaxNumIter` is reached
%                     3. stop if `MaxNumMapEval` is reached
%                     4. stop if `TimeLimit` is reached
%                     5. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
%
% OUTPUTS:
%    x_best:      the best approximation of the optimizer
%    psi_best:    the best approximation of the optimum
%    out:         structure including more output information
%
%                   * .T - running time
%                   * .Niter - total number of iterations
%                   * .Nmap - total number of mapping evaluations
%                   * .merit_func - array including all merit function values
%                   * .x_error - relative error :math:`norm(x_k(:)-x_{opt}(:))/norm(x_{opt})`
%                   * .psi_error - relative error :math:`(\psi_k-\psi_{opt})/(\psi_0-\psi_{opt}))`
%                   * .Status - reason of termination
%
% .. REFERENCE:
% .. Algorithm 2 of [1]: F.J. Aragon Artacho, R.M.T. Fleming, Globally convergent algorithms for finding zeros of duplomonotone mappings, Optimization Letter, 9, 569-584 (2015)
% .. Author: - Masoud Ahookhosh, System Biochemistry Group, Luxembourg Center for System Biomedicine, University of Luxembourg, Luxembourg
%            - Update July 2017 - M. Ahookhosh

format longG ;

% ================ Error messages for input and output =================
if nargin > 3
    error('The number of input arguments is more than what is needed');
elseif nargin < 3
    error('The number of input arguments is not enough');
end;

if isempty(mapp)
    error('the function handle mapp has to be defined');
elseif ~isa(mapp,'function_handle')
    error('mapp should be a function handle');
end

if isempty(x0)
    error('The starting point x0 has to be defined');
elseif ~isa(x0,'numeric')
    error('x0 should be a numeric vector');
end

% =================== initializing the parameters ======================
% ===== user has requested viewing the default values of "options" =====
[MaxNumIter,MaxNumMapEval,TimeLimit,epsilon,alpha, ...
         beta,sigma,l,tauBar,lambda_min,lambda_max,flag_x_error, ...
         flag_psi_error,flag_time,Stopping_Crit] = InitialDuplo(options);

if isfield(options,'x_opt')
    x_opt=options.x_opt;
elseif flag_x_error==1
    error('x_error requires to x_opt be specified');
end

if flag_x_error == 1
    Nxopt      = sqrt(sum(x_opt(:).^2));
    x_error(1) = sqrt(sum((x0(:)-x_opt(:)).^2))/Nxopt;
end

if flag_psi_error == 1
    psi_error(1) = 1;
end

if flag_time == 1
    Time(1) = 0;
end

xk         = x0;
Xk         = xk;
Niter      = 1;
fx0        = mapp(x0);
Nmap       = 1;
nfx0       = norm(fx0);
fxk        = fx0;
nfxk       = nfx0;
merit_func = 0.5*nfxk^2;
lambda     = min(sigma/l^2,tauBar);
StopFlag   = 0;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Main body of CSDF.m %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T0 = tic;

% ======================= start of the main loop =======================
while ~StopFlag

    xk    = xk-lambda*fxk;
    Niter = Niter+1
    fxk   = mapp(xk);
    Nmap  = Nmap+1;
    nfxk2 = norm(fxk)^2;

    % ================= Gathering output information ===================
    psik              = 0.5*nfxk2;
    merit_func(Niter) = psik;
    if flag_time == 1
        Time(Niter+1) = toc(T0);
    end

    if flag_x_error == 1
        Nx_opt = norm(x_opt);
        x_error(Niter+1) = sqrt(sum((xk(:)-x_opt(:)).^2))/Nx_opt;
    end

    if flag_psi_error == 1
        psi_error(Niter+1) = (psik-psi_opt)/(psi0-psi_opt);
    end

    % ================== checking stopping criteria ====================
    T = toc(T0);

    [StopFlag,Status] = StopCritDuplo(nfxk,Niter,Nmap,T, ...
              MaxNumIter,MaxNumMapEval,TimeLimit,epsilon,Stopping_Crit);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Status
x_best         = xk;
psi_best       = psik;
out.Xk         = Xk;
out.T          = T;
out.nhx        = nfxk;
out.merit_func = merit_func';
out.Niter      = Niter;
out.Nmap       = Nmap;
out.Status     = Status;

if flag_x_error == 1
    out.x_error = x_error;
end
if flag_psi_error == 1
    out.psi_error = psi_error;
end
if flag_time == 1
    out.Time = Time;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of CSDF.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function H = mrdivide(F, G)
%/   SPHEREFUNV right divide.
%   F/G divides each component of the SPHEREFUNV F by the DOUBLE or
%   SPHEREFUN G.
%
%   Only allowed to divide by a DOUBLE or a SPHEREFUN.
%
% See also SPHEREFUNV/MLDIVIDE.

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

if ( isempty(F) || isempty(G) )
   H = spherefunv;
   return
end

if ( ~isa(G,'double') && ~isa(G,'chebfun2') )
    error('SPHEREFUN:SPHEREFUNV:mrdivide:nonScalar', ...
        'Division must be scalar valued.');
end

% componentwise divide.
if ( isa(G,'double') )
    H = F;
    for j = 1:3
        H.components{j} = mrdivide(F.components{j}, G);
    end
else
    H = F;
    for j = 1:3
        H.components{j} = rdivide(F.components{j}, G);
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% MATLAB xUnit Test Framework
% Version 3.0.1 (R2010a) 16-Jun-2010
%
% Running Unit Tests
%   runtests                  - Run unit tests
%
% Writing Unit Tests
%   assertElementsAlmostEqual - Assert floating-point array elements almost equal
%   assertEqual               - Assert that inputs are equal
%   assertFilesEqual          - Assert that two files have the same content
%   assertExceptionThrown     - Assert that specified exception is thrown
%   assertFalse               - Assert that input condition is false
%   assertTrue                - Assert that input condition is true
%   assertVectorsAlmostEqual  - Assert floating-point vectors almost equal in norm sense
%   initTestSuite             - Utility script used for subfunction-based tests
%
% Framework Classes
%   CommandWindowTestRunDisplay - Print test suite results to command window
%   VerboseCommandWindowTestRunDisplay - Print verbose test suite results to command window
%   FunctionHandleTestCase    - Test case based on a function handle
%   TestCase                  - Class defining interface for test cases
%   TestCaseInDir             - Test case requiring temporary directory change
%   TestCaseWithAddPath       - Test case requiring temporary path modification
%   TestComponent             - Abstract base class for TestCase and TestSuite
%   TestComponentInDir        - Test component requiring temporary directory change
%   TestLogger                - Collect data (silently) from running test suite
%   TestRunMonitor            - Abstract base class for monitoring test suite
%   TestSuite                 - Collection of TestComponent objects
%   TestSuiteInDir            - Test suite requiring temporary directory change

% Steven L. Eddins
% Copyright 2008-2010 The MathWorks, Inc.

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [warnings] = addplotlystartup(startupPaths)
%[1]looks at statup.m files specified by the entries of startupPaths
%[2]appends the addplotly function to startup.m files (if not already present)
%[3]checks for other plotly addpath calls within any startup.m and outputs warning

%output warnings
warnings = cell(size(startupPaths));

for locs = 1:size(startupPaths,1);
    %addpath string for Plotly API
    addString = 'addpath(genpath(fullfile(matlabroot,''toolbox'',''plotly'')),''-end'');';
    %open current startup.m to read
    currentStartupID = fopen(startupPaths{locs},'r');
    if currentStartupID == -1, error('plotly:startupRead',...
            ['\n\nShoot! It looks like something went wrong reading the file: ' ...
            '\n' startupPaths{locs} '\n' ...
            '\nPlease contact your system admin. or chuck@plot.ly for more \ninformation. In the ' ...
            'mean time you can add the Plotly API \nto your search path manually whenever you need it! \n\n']);
    end
    %check for any instances of the addplotlyapi function
    startupScan = textscan(currentStartupID, '%s', 'delimiter', '\n', 'whitespace', '');
    startupLines = startupScan{1};
    Index = find(strcmp(startupLines,addString));
    otherPlotlyOccurrence = findstr(fileread(startupPaths{locs}),'plotly');
    %if addString is not in startup.m add it
    if(~any(Index));
        %reopen current startup.m with new permission
        currentStartupID = fopen(startupPaths{locs},'a+');
        if currentStartupID == -1, error('plotly:startupWrite',...
            ['\n\nShoot! It looks like something went wrong writing to the file: ' ...
            '\n\n' startupPaths{locs} '\n' ...
            '\nPlease contact your system admin. or chuck@plot.ly for more \ninformation. In the ' ...
            'mean time you can add the Plotly API \nto your search path manually whenever you need it! \n']);
        end
        fprintf(currentStartupID,['\n' addString]);
    end
    if(length(Index) ~= length(otherPlotlyOccurrence));
        warnings{locs} = ['\n[WARNING]: \n\nWe found an addpath specification for another version of Plotly at: ' ...
                                '\n\n' startupPaths{locs} '\n\nyou may be forcing MATLAB to look for an older version of Plotly!\n\n'];
    else
        warnings{locs} = '';
    end
    fclose(currentStartupID);
end
end

author: goGPS-Project


text:
function polygon_struct = Read_shapefile( finputname, polygon, bbox, ...
                                          h0, boubox, plot_on, shapefile_3d)
% Read_shapefile: Reads a shapefile or a NaN-delimited vector
% containing polygons and/or segments in the the desired region
% of interest. Classifies the vector data as either a
% mainland, island, or outer boundary. The program will automatically
% trim islands that are have an area smaller than 4*h0^2 and will in gaps
% in the vector that are larger than h0/2. This is necessary for the
% use of boundary rejection method in dpoly.

% INPUTS:
% finputname : file name(s) of the shapefile listed in a cell array
% polygon    : a NaN-delimited vector of polygons and/or segments.
% bbox       : the bounding box that we want to extract out
% h0         : minimum edge length in region of interest.
% plot_on    : plot the final polygon or not (=1 or =0)
%
% OUTPUTS:
% polygon_struct    : a structure containing the vector features identified as
%              either islands, mainland, or outer.
% Written by William Pringle and Keith Roberts, CHL,UND, 2017
% Edits by Keith Roberts, July 2018.
%% Loop over all the filenames and get the shapefile within bbox
SG = [];
if bbox(1,2) > 180 && bbox(1,1) < 180
    % bbox straddles 180/-180 meridian
    loop = 2; minus = 0;
elseif all(bbox(1,:) > 180)
    % beyond 180 in 0 to 360 format
    loop = 1; minus = 1;
else
    loop = 1; minus = 0;
end
if (size(finputname,1)~=0)
    for fname = finputname
        for nn = 1:loop
            bboxt = bbox';
            if loop == 2
                if nn == 1
                    bboxt(2,1) = 180;
                else
                    bboxt(1,1) = -180; bboxt(2,1) = bboxt(2,1) - 360;
                end
            end
            if minus
               bboxt(:,1) = bboxt(:,1) - 360;
            end
            % Read the structure
            try
                % The shaperead is faster if it is available
                S = shaperead(fname{1},'BoundingBox',bboxt);
                % Get rid of unwanted components;
                D = struct2cell(S);
                S = cell2struct(D(3:4,:)',{'X','Y'},2);
                disp('Read shapefile with shaperead')
                sr = 1;
            catch
                % If only m_shaperead is available or if some error occured
                % with shaperead (e.g., 3D shapefile, m_shaperead may work)
                disp('Reading shapefile with m_shaperead')
                % This uses m_map (slower but free)
                S = m_shaperead(fname{1},reshape(bboxt',4,1));
                % Let's just keep the x-y data
                D = S.ncst;
                if isfield(S,'dbf')  || isfield(S,'dbfdata')
                    code = S.dbfdata(:,1);
                    S = cell2struct([D code]',{'points' 'type'},1);
                else
                    S = cell2struct(D','points',1);
                end
                sr = 0;
            end
            if ~isempty(S)
                % Keep the following polygons
                SG = [SG; S];
            end
        end
    end
else
    sr = 1 ;
    % convert NaN-delimited vector to struct
    count = 1; j=1;
    for i = 1 : length(polygon)
        % the end of the segment

        if(isnan(polygon(i,1))==1)
            % put NaN at end
            SG(count,1).X(:,j) =NaN;
            SG(count,1).Y(:,j) =NaN;
            % reset
            j = 1 ; count = count + 1;

            continue
        else
            % keep going
            SG(count,1).X(:,j) = polygon(i,1);
            SG(count,1).Y(:,j) = polygon(i,2);
            j=j+1;
        end
    end
end
% If we don't have an outer polygon already then make it by bbox
polygon_struct.outer = boubox;
% Densify the outer polygon (fills gaps larger than half min edgelength).
[latout,lonout] = my_interpm(polygon_struct.outer(:,2),...
                             polygon_struct.outer(:,1),h0/2);
polygon_struct.outer = [];
polygon_struct.outer(:,1) = lonout;
polygon_struct.outer(:,2) = latout;

%% Find whether the polygons are wholey inside the bbox..
%% Set as islands or mainland
disp('Partitioning the boundary into islands, mainland, ocean')
polygon_struct.inner = [];
polygon_struct.mainland = [];
polygon_struct.innerb = [];
polygon_struct.mainlandb = [];
polygon_struct.innerb_type = [];
polygon_struct.mainlandb_type = [];
edges = Get_poly_edges( polygon_struct.outer );

if isempty(SG); return; end

if sr
    tmpM = [[SG.X]',[SG.Y]'] ; % MAT
    if bbox(1,2) > 180
        tmpM(tmpM(:,1) < 0,1) =  tmpM(tmpM(:,1) < 0,1) + 360;
    end
    for i = 1 : length(SG)
       dims(i) = length(SG(i).X) ;
    end
    tmpC = mat2cell(tmpM,dims); % TO CELL
else
    tmpC =  struct2cell(SG)';
    tmpCC = []; nn = 0;
    for ii = 1:size(tmpC,1)
       % may have NaNs inside
       isnan1 = find(isnan(tmpC{ii,1}(:,1)));
       if isempty(isnan1)
           isnan1 = length(tmpC{ii,1})+1;
       elseif isnan1(end) ~= length(tmpC{ii,1})
           isnan1(end+1) = length(tmpC{ii,1})+1;
       end
       is = 1;
       for jj = 1:length(isnan1)
           nn = nn + 1;
           ie = isnan1(jj)-1;
           tmpCC{nn,1} = tmpC{ii,1}(is:ie,:);
           tmpCC{nn,2} = tmpC{ii,2};
           is = isnan1(jj)+1;
       end
    end
    if ~isempty(tmpCC)
        tmpC = tmpCC;
    end
    tmpM =  cell2mat(tmpC(:,1));
    if size(tmpM,2) == 3
       tmpM = tmpM(:,1:2);
    end
end
% Get current polygon
% Check proportion of polygon that is within bbox
tmpIn = inpoly(tmpM,polygon_struct.outer, edges);
tmpInC = mat2cell(tmpIn,cellfun(@length,tmpC(:,1)));

j = 0 ; k = 0 ; height = []; new_islandb = []; new_mainb = [];
new_islandb_type = []; new_mainb_type = [];
for i = 1 : size(tmpC,1)
    if sr
        % using shaperead
        points = tmpC{i,1}(1:end-1,:) ;
        In     = tmpInC{i,1}(1:end-1) ;
    else
        % using m_shaperead
        points = tmpC{i,1}(1:end,:) ;
        if size(points,2) == 3
            points = points(:,1:2);
        end
        if shapefile_3d
            % if 3-D shapefile
            height = points(:,3);
            points = points(:,1:2);
            type   = tmpC{i,2};
            if strcmp(type,'BA040')
                type = 'ocean';
            elseif strcmp(type,'BH080')
                type = 'lake';
            elseif strcmp(type,'BH140')
                type = 'river';
            end
        else
            height = [];
        end
        In     = tmpInC{i,1}(1:end) ;
    end
    if bbox(1,2) > 180
       lond = abs(diff(points(:,1)));
       if any(lond > 350)
           points(points(:,1) > 180,1) = 0;
       end
    end
    % lets calculate the area of the
    % feature using the shoelace algorithm and decided whether to keep or
    % not based on area.
    if all(points(1,:) == points(end,:))
        area = shoelace(points(:,1),points(:,2));
    else
        area = 999; % not a polygon
    end
    if length(find(In == 1)) == length(points)
        % Wholey inside box
        if area < 4*h0^2 % too small, then don't consider it.
            continue;
        end
        % Set as island (with NaN delimiter)
        k = k + 1 ;
        new_island{k} = [points; NaN NaN];
        if ~isempty(height)
            new_islandb{k} = [points height; NaN NaN NaN];
            new_islandb_type{k} = type;
        end
    else
        %Partially inside box
        if area < 100*h0^2 % too small, then don't consider it.
            continue;
        end
        % Set as mainland
        j = j + 1 ;
        new_main{j} = [points; NaN NaN];
        if ~isempty(height)
            new_mainb{j} = [points height; NaN NaN NaN];
            new_mainb_type{j} = type;
        end
    end
end
if k > 0
    polygon_struct.inner = [polygon_struct.inner; cell2mat(new_island')];
    polygon_struct.innerb = cell2mat(new_islandb');
    polygon_struct.innerb_type = new_islandb_type;
end
if j > 0
    polygon_struct.mainland = [polygon_struct.mainland; cell2mat(new_main')];
    polygon_struct.mainlandb = cell2mat(new_mainb');
    polygon_struct.mainlandb_type = new_mainb_type;
end
% Merge overlapping mainland and inner
if exist('polybool','file') || exist('polyshape','file')
    if ~isempty(new_mainb) && k > 0
        polym = polygon_struct.mainland;
        mergei = false(k,1);
        for kk = 1:k
            polyi = new_island{kk};
            IA = find(ismembertol(polym,polyi,1e-5,'ByRows',true));
            if length(IA) > 2
               if exist('polyshape','file')
                  polyout = union(polyshape(polym),polyshape(polyi));
                  polym = polyout.Vertices;
               else
                  [x,y] = polybool('union',polym(:,1),polym(:,2),...
                                   polyi(:,1),polyi(:,2));
                  polym = [x,y];
               end
               mergei(kk) = 1;
            end
        end
        if ~isnan(polym(end,1)); polym(end+1,:) = NaN; end
        polygon_struct.mainland = polym;
        new_island(mergei) = [];
        polygon_struct.inner = cell2mat(new_island');
    end
else
    warning(['no polyshape or polybool available to merge possible ' ...
             'overlapping of mainland and inner'])
end
% Remove parts of inner and mainland overlapping with outer
polygon_struct.outer = [polygon_struct.outer; polygon_struct.mainland];
%% Plot the map
if plot_on >= 1 && ~isempty(polygon_struct)
    figure(1);
    hold on
    plot(polygon_struct.outer(:,1),polygon_struct.outer(:,2))
    if ~isempty(polygon_struct.inner)
        plot(polygon_struct.inner(:,1),polygon_struct.inner(:,2))
    end
    if ~isempty(polygon_struct.mainland)
        plot(polygon_struct.mainland(:,1),polygon_struct.mainland(:,2))
    end
end
%EOF
end

author: CHLNDDEV


text:
function gmt_master(sOutput)

hAxes = gca;

gmt_init(hAxes, sOutput);

vKids = get(hAxes, 'Children');

for nCnt = length(vKids):-1:1
  sType = lower(get(vKids(nCnt), 'Type'));
  switch sType
  case 'surface'
    gmt_pcolor(vKids(nCnt), sOutput, num2str(nCnt));
  case 'line'
    gmt_xy(vKids(nCnt), sOutput, num2str(nCnt));
  end
end

gmt_done(sOutput);

author: CelsoReyes


text:
function [reduced_pot,successful] = reduce_pot(pot,tailnodes)
% Executes the reduce operation defined in
% Stable Local Computation with Conditional Gaussian Distributions
% Steffen L. Lauritzen
% Frank Jensen
% September 1999
% The potential pot is reduced if B contains any zero columns
% The test are restricted to the positions in tailnodes.
% Any columns successfully deleted are entered in the array successful
if nargin < 2
    tailnodes = pot.ctaildom;
end

successful = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Keep track of remaining tailnodes %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rem_tailnodes = pot.ctaildom;
for i = tailnodes
    pos = find(i==rem_tailnodes);
    successful_red = [pos];
    red_scgcpot = cell(1,pot.dsize);
    j = 1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Test whether all components of pot.scgpotc can be reduced %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    while ((j <= pot.dsize) & ~isempty(successful_red))
        [cpot,successful_red] = reduce_pot(pot.scgpotc{j},pos);
        red_scgcpot{j} = cpot;
        j = j + 1;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % If i is a reducible tailnode, then reduce the potential %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(successful_red)
        successful = [successful i];
        pot.scgpotc = red_scgcpot;
        rem_tailnodes = mysetdiff(rem_tailnodes,i);
    end;
end

pot.ctaildom = rem_tailnodes;
positions = find_equiv_posns(rem_tailnodes,pot.ctaildom);
pot.ctailsizes = pot.ctailsizes(positions);
pot.ctailsize = sum(pot.ctailsizes);
pot.domain = mysetdiff(pot.domain,successful);
reduced_pot = pot;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function tP = calcTriplePoints(gB)
%
% Input
%  gB - @grainBoundary

% list of phaseIds ordered as grainIds
grainId = gB.grainId;
grainPhaseId = zeros(max(grainId(:)),1);
grainPhaseId(grainId(grainId>0)) = gB.phaseId((grainId>0));

% compute triple points
[i,~,f] = find(gB.F);
I_VF = sparse(f,i,1,size(gB.V,1),size(gB.F,1));
I_VG = (I_VF * gB.I_FG)==2;
% triple points are those with exactly 3 neigbouring grains and 3
% boundary segments
itP = full(sum(I_VG,2)==3 & sum(I_VF,2)==3);
[tpGrainId,~] = find(I_VG(itP,:).');
tpGrainId = reshape(tpGrainId,3,[]).';
tpPhaseId = full(grainPhaseId(tpGrainId));

% compute ebsdId
% first step: compute faces at the triple point
% clean up incidence matrix
%I_FD(~any(I_FD,2),:) = [];
% incidence matrix between triple points and voronoi cells
%I_TD = I_VF(itP,:) * I_FD;
[tPBoundaryId,~] = find(I_VF(itP,:).');
tPBoundaryId = reshape(tPBoundaryId,3,[]).';

% get the three end vertices
iV = reshape(gB.F(tPBoundaryId,:),[],6).';
% TODO: the repmat can be removed in new versions of Matlab
iV = reshape(iV(iV ~= repmat(find(itP).',size(iV,1),1)).',3,[]).';

tP = triplePointList(find(itP),gB.V,...
  tpGrainId,tPBoundaryId,tpPhaseId,iV,gB.phaseMap,gB.CSList);

end

author: mtex-toolbox

repo: mtex

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


text:
function options = drawFluxVariability(map,model,minFlux,maxFlux,options,varargin)
%drawFluxVariablity Overlay flux variability data on a metabolic map
%
% options = drawFluxVariability(map,model,minFlux,maxflux,options,varargin)
%
%INPUTS
% map           COBRA map structure
% model         COBRA model structure
% minFlux       Vector containing minimum flux values for each reaction
% maxFlux       Vector containing maximum flux values for each reaction
%
%OPTIONAL INPUTS
% options       options structure
% varargin      parameter, parameter value inputs for optional parameters
%
%OUTPUT
% options       options structure used to generate map
%
%Richard Que (03/2010)

if nargin<5, options = []; end
if ~isfield(options,'fluxVarColor'), options.fluxVarColor = []; end
if mod(length(varargin),2)==0
    for i=1:2:length(varargin)-1
        options = setMapOptions(options,map,model,varargin{i},varargin{i+1});
%         switch lower(varargin{i})
%             case 'bidircolor', options.fluxVarColor.biDirColor = varargin{i+1};
%             case 'unidirirrcolor', options.fluxVarColor.uniDirIrrColor = varargin{i+1};
%             case 'unidirfwdcolor', options.fluxVarColor.uniDirFwdColor = varargin{i+1};
%             case 'unidirrevcolor', options.fluxVarColor.uniDirRevColor = varargin{i+1};
%             case 'filename', options.fileName = varargin{i+1};
%             case 'rxndirflag', options.rxnDirFlag = varargin{i+1};
%             case 'rxndirmultiplier', options.rxnDirMultiplier = varargin{i+1};
%             otherwise
%                 options.(varargin{i}) = varargin{i+1};
%                 fprintf('Unrecognized parameter %s\n',varargin{i});
%         end
    end
else
    error('Invalid number of parameters/values');
end
color = zeros(length(model.rxns),3);
%set defaults
if ~isfield(options.fluxVarColor,'biDirColor')
    options.fluxVarColor.biDirColor = [0 255 0];
end
if ~isfield(options.fluxVarColor,'uniDirIrrColor')
    options.fluxVarColor.uniDirIrrColor = [0 0 255];
end
if ~isfield(options.fluxVarColor,'uniDirFwdColor')
    options.fluxVarColor.uniDirFwdColor = [255 0 255];
end
if ~isfield(options.fluxVarColor,'uniDirRevColor')
    options.fluxVarColor.uniDirRevColor = [0 255 255];
end
if isfield(options,'rxnDirMultiplier')
    options.rxnDirFlag = true;
end
if ~isfield(options,'rxnDirFlag')
    rxnDirFlag = false;
else
    rxnDirFlag = options.rxnDirFlag;
end

%assign colors
for i=1:length(model.rxns)
    if minFlux(i)*maxFlux(i)>=0
        %unidirectional
        if model.lb(i) < 0
            %reversible
            if minFlux(i)<0 || maxFlux(i)<0
                %reverse direction: Cyan 0 255 255
                color(i,:) = options.fluxVarColor.uniDirRevColor;
            else
                %forward direction: Magenta 255 0 255
                color(i,:) = options.fluxVarColor.uniDirFwdColor;
            end
        else
            %irreversible: Blue 0 0 255
            color(i,:) = options.fluxVarColor.uniDirIrrColor;
        end
    else
        %bidirectional: Green 0 255 0
        color(i,:) = options.fluxVarColor.biDirColor;
    end
end

if rxnDirFlag
    options.rxnDir = zeros(length(map.rxnIndex),1);
    options.rxnDir(ismember(map.connectionAbb,model.rxns((maxFlux>0) & (minFlux>=0)))) = 1;
    options.rxnDir(ismember(map.connectionAbb,model.rxns((maxFlux<=0) & (minFlux<0)))) = -1;
end

options = setMapOptions(options,map,model,'edgeColor',color);
drawCbMap(map,options);