text:
function obj=draw(obj,do_redraw)
% draw Draws the plot in the current figure
%
% syntax: gramm_object.draw()
% Call draw on an array of gramm objects in order to have
% multiple gramms on a single figure.
% An optional logical argument can be given. Use
% gramm_object.draw(false) in order to allow the nex gramm plot to be
% superimposed on the same figure (useful for plotting with
% different grouping aesthetics). When using superimposed
% plots, the last draw() should be given without this false
% argument. Giving false as argument deactivates the automatic
% call for the redraw() function, i.e. deactivates the fancy
% axis placement.

%We set redraw() as resize callback by default
if nargin<2
    do_redraw=obj(1).layout_options.redraw;
end

%If no parent was given we use the current figure
if isempty(obj(1).parent)
    for obj_ind=1:numel(obj)
        obj(obj_ind).parent=gcf;
    end
end

%There are bugs with the hardware openGL renderer in pre 2014b version,
%on Mac OS and Win the x and y axes become invisible and on
%Windows the patch objects behave strangely.
if ~obj(1).handle_graphics
    if ismac
        warning('Mac Pre-2014b version detected, forcing to ''Painters'' renderer to show axes lines. Use set(gcf,''Renderer'',''OpenGL'') to restore transparency')
        if isprop(obj(1).parent,'Renderer')
            set(obj(1).parent,'Renderer','Painters');
        end
    else
        warning('Windows Pre-2014b version detected, forcing to software openGL which is less buggy')
        opengl software
    end
end

%Handle call of draw() on array of gramm objects by dividing the figure
%up and launching the individual draw functions of each object
if numel(obj)>1
    %Take care of the big title
    if isempty(obj(1).bigtitle)
        maxh=1;
    else
        maxh=0.94;
        tmp=axes('Position',[0.1 0.965 0.8 0.01],'Parent',obj(1).parent);
        set(tmp,'Visible','off','XLim',[-1 1],'YLim',[-1 1]);
        tmp=text(0,0,obj(1).bigtitle,...
            'FontWeight','bold',...
            'Interpreter',obj(1).text_options.interpreter,...
            'FontName',obj(1).text_options.font,...
            'FontSize',obj(1).text_options.base_size*obj(1).text_options.big_title_scaling,...
            'HorizontalAlignment','center');

        if ~isempty(obj(1).bigtitle_options)
            set(tmp,obj(1).bigtitle_options{:});
        end
    end

    nl=size(obj,1);
    nc=size(obj,2);
    for l=1:nl
        for c=1:nc
            if strcmp(obj(l,c).layout_options.position,'auto')
                obj(l,c).multi.orig=[(nl-l)*maxh/nl (c-1)/nc];
                obj(l,c).multi.size=[maxh/nl 1/nc];
            else
                obj(l,c).multi.orig=[obj(l,c).layout_options.position(2)*maxh obj(l,c).layout_options.position(1)];
                obj(l,c).multi.size=[obj(l,c).layout_options.position(4)*maxh obj(l,c).layout_options.position(3)];
            end
            obj(l,c).layout_options.redraw_gap=obj(l,c).layout_options.redraw_gap/max(nl,nc);
            obj(l,c).multi.active=true;
            draw(obj(l,c));
        end
    end

    %Set up tight redrawing for multiple plots
    if do_redraw
        redraw(obj);
        if verLessThan('matlab','8.4')
            set(gcf,'ResizeFcn',@(a,b)redraw(obj));
        else
            set(gcf,'SizeChangedFcn',@(a,b)redraw(obj));
        end
    end
    return;
end

%If x is empty the object probably is so we skip (happens for multiple graphs)
if isempty(obj.aes.x)
    return
end

%If draw() is called a second time without update() having been
%called we skip
if (~obj.updater.first_draw && ~obj.updater.updated)
    warning('Multiple draw() calls need update() calls')
    return
end

obj.aes=validate_aes(obj.aes);

%Handle multiple figures
uni_fig=unique_and_sort(obj.aes.fig,1);
if ~iscell(uni_fig)
    uni_fig=num2cell(uni_fig);
end
if numel(uni_fig)>1 %if multiple figures
    fig_pos=get(obj.parent,'Position'); %We'll create them with the same location and size as the current figure
    %Convert to a cell
    obj={obj};
    for ind_fig=2:numel(uni_fig) %Make shallow copy of original object
        obj{ind_fig}=copy(obj{1});
    end
    for ind_fig=1:numel(uni_fig) %Loop over objects
        %disp(['Plotting Fig ' num2str(ind_fig)])
        %We plot what we want by setlecting within aes on the basis of fig
        obj{ind_fig}.aes=select_aes(obj{ind_fig}.aes,multi_sel(obj{ind_fig}.aes.fig,uni_fig{ind_fig}));
        % Skip if there is no data
        if any(obj{ind_fig}.aes.subset)
            %Set title
            obj{ind_fig}.set_title([obj{ind_fig}.aes_names.fig ': ' num2str(uni_fig{ind_fig})]);
            %Create figure and do the drawing
            if ind_fig>1
                obj{ind_fig}.parent=figure('Position',fig_pos);
            end
            obj{ind_fig}.draw();
        end
    end
    return;
end

if ~any(obj.aes.subset)
    disp('Empty subset, nothing to draw');
    return;
end

%Apply subset
temp_aes=select_aes(obj.aes,obj.aes.subset);

%Convert factor data to numeric
if iscellstr(temp_aes.x)
    obj.x_factor=true;
    obj.x_ticks=unique_and_sort(temp_aes.x,obj.order_options.x);
    tempx=zeros(length(temp_aes.x),1);
    for k=1:length(obj.x_ticks)
        tempx(strcmp(temp_aes.x,obj.x_ticks{k}))=k;
    end
    temp_aes.x=tempx;
else
    obj.x_factor=false;
end

%Compute x limits
obj.var_lim.maxx=allmax(temp_aes.x);
obj.var_lim.minx=allmin(temp_aes.x);

%Depending on whether ymin of ymax are present we change what to pick for Y
%limit computation
if ~isempty(temp_aes.ymin) && ~isempty(temp_aes.ymax)
    tmp_y_for_min=temp_aes.ymin;
    tmp_y_for_max=temp_aes.ymax;
else
    tmp_y_for_min=temp_aes.y;
    tmp_y_for_max=temp_aes.y;
end

%Compute y limits
obj.var_lim.maxy=allmax(tmp_y_for_max);
obj.var_lim.miny=allmin(tmp_y_for_min);

%Compute z limits
if ~isempty(temp_aes.z)
   obj.var_lim.maxz=allmax(temp_aes.z);
   obj.var_lim.minz=allmin(temp_aes.z);
end

%Handles the multiple gramm case: we set up subtightplot so
%that it restricts the drawing (using margins) to a portion of
%the figure
mysubtightplot=@(m,n,p,gap,marg_h,marg_w)my_tightplot(m,n,p,...
    [gap(2)*obj.multi.size(1) gap(1)*obj.multi.size(2)],...
    [obj.multi.orig(1)+obj.multi.size(1)*marg_h(1) 1-obj.multi.orig(1)-obj.multi.size(1)*(1-marg_h(2))],...
    [obj.multi.orig(2)+obj.multi.size(2)*marg_w(1) 1-obj.multi.orig(2)-obj.multi.size(2)*(1-marg_w(2))],'Parent',obj.parent);

%Set subplot generation parameters and functions
if strcmp(obj.layout_options.gap,'auto')
    if obj.force_ticks  %Independent scales require more space between subplots for ticks
        mysubplot=@(nrow,ncol,row,col)mysubtightplot(nrow,ncol,(row-1)*ncol+col,[0.06 0.06],obj.layout_options.margin_height,obj.layout_options.margin_width);
    else
        mysubplot=@(nrow,ncol,row,col)mysubtightplot(nrow,ncol,(row-1)*ncol+col,[0.02 0.02],obj.layout_options.margin_height,obj.layout_options.margin_width);
    end
else
    mysubplot=@(nrow,ncol,row,col)mysubtightplot(nrow,ncol,(row-1)*ncol+col,obj.layout_options.gap,obj.layout_options.margin_height,obj.layout_options.margin_width);
end

%Subplots for wraps leave more space above axes for column
%legend
if strcmp(obj.layout_options.gap,'auto')
    mysubplot_wrap=@(ncol,col)mysubtightplot(ceil(ncol/obj.wrap_ncols),obj.wrap_ncols,col,[0.03 0.09],obj.layout_options.margin_height,obj.layout_options.margin_width);
else
    mysubplot_wrap=@(ncol,col)mysubtightplot(ceil(ncol/obj.wrap_ncols),obj.wrap_ncols,col,obj.layout_options.gap,obj.layout_options.margin_height,obj.layout_options.margin_width);
end

%Find uniques in aesthetics and sort according to options
uni_row=unique_and_sort(temp_aes.row,obj.order_options.row);
uni_column=unique_and_sort(temp_aes.column,obj.order_options.column);
uni_linestyle=unique_and_sort(temp_aes.linestyle,obj.order_options.linestyle);
uni_marker=unique_and_sort(temp_aes.marker,obj.order_options.marker);
uni_lightness=unique_and_sort(temp_aes.lightness,obj.order_options.lightness);
uni_size=unique_and_sort(temp_aes.size,obj.order_options.size);

%If the color is in a cell array of doubles, we set it as
%continuous color
if iscell(temp_aes.color) && ~iscellstr(temp_aes.color)
    obj.continuous_color_options.active = true;
else
    uni_color=unique_and_sort(temp_aes.color,obj.order_options.color);

    % If continous_color is not set and we have too many numerical values
    % for the color we switch to continuous color
    if isnan(obj.continuous_color_options.active)
        if length(uni_color)>15 && ~iscellstr(uni_color)
            obj.continuous_color_options.active = true;
        else
            obj.continuous_color_options.active = false;
        end
    elseif iscellstr(uni_color) %We can't have continuous color for cellstr input
        obj.continuous_color_options.active = false;
    end
end
if obj.continuous_color_options.active
    uni_color={1};
end

%Transform all aesthetics in cell arrays (makes it easier to handle after)
if ~iscell(uni_row)
    uni_row=num2cell(uni_row);
end
if ~iscell(uni_column)
    uni_column=num2cell(uni_column);
end
if ~iscell(uni_color)
    uni_color=num2cell(uni_color);
end
if ~iscell(uni_lightness)
    uni_lightness=num2cell(uni_lightness);
end
if ~iscell(uni_linestyle)
    uni_linestyle=num2cell(uni_linestyle);
end
if ~iscell(uni_size)
    uni_size=num2cell(uni_size);
end
if ~iscell(uni_marker)
    uni_marker=num2cell(uni_marker);
end

%Correct empty facets when wrapping
if obj.wrap_ncols>length(uni_column)
    obj.wrap_ncols=length(uni_column);
end

%The plot minimums and maximums are stored in arrays the size of the
%number of subplots
n_columns=length(uni_column);
n_rows=length(uni_row);
if obj.updater.first_draw || obj.updater.facet_updated==1
    obj.plot_lim.minx=nan(n_rows,n_columns);
    obj.plot_lim.maxx=nan(n_rows,n_columns);
    obj.plot_lim.miny=nan(n_rows,n_columns);
    obj.plot_lim.maxy=nan(n_rows,n_columns);
    obj.plot_lim.minz=nan(n_rows,n_columns);
    obj.plot_lim.maxz=nan(n_rows,n_columns);
    obj.plot_lim.minc=nan(n_rows,n_columns);
    obj.plot_lim.maxc=nan(n_rows,n_columns);
end

%Get colormap
cmap=get_colormap(length(uni_color),length(uni_lightness),obj.color_options);

%Initialize results structure (n_groups is an overestimate if
%there are redundant groups
n_groups=length(uni_row)*length(uni_column)*length(uni_marker)...
    *length(uni_size)*length(uni_linestyle)*length(uni_color)*length(unique(temp_aes.group));
aes_names_fieldnames=fieldnames(obj.aes_names);
for fn=1:length(aes_names_fieldnames)
    if ~any(strcmp(aes_names_fieldnames{fn},{'x' 'y' 'z' 'ymin' 'ymax' 'label' 'fig'})) %We ignore x y z and label
        obj.results.(aes_names_fieldnames{fn})=cell(n_groups,1);
        obj.results.(['ind_' aes_names_fieldnames{fn}])=cell(n_groups,1);
    end
end
obj.results.draw_data=cell(n_groups,1);

obj.firstrun=ones(n_rows,n_columns);
%Index in the loops
obj.result_ind=1;

%Initialize facet axes handles
if obj.updater.first_draw
    if obj.handle_graphics %Post 2014b we return objects
        obj.facet_axes_handles=gobjects(length(uni_row),length(uni_column));
    else
        obj.facet_axes_handles=zeros(length(uni_row),length(uni_column));
    end
end

%Set the dodge width across all facets
draw_data.dodge_avl_w=min(diff(unique(comb(temp_aes.x))));
if isempty(draw_data.dodge_avl_w) || draw_data.dodge_avl_w==0
    draw_data.dodge_avl_w=1;
end

obj.data_size=numel(temp_aes.x);

%% draw() looping
%Loop over rows
for ind_row=1:length(uni_row)

    sel_row=multi_sel(temp_aes.row,uni_row{ind_row});

    obj.current_row=ind_row;

    %Loop over columns
    for ind_column=1:length(uni_column)

        sel_column=sel_row & multi_sel(temp_aes.column,uni_column{ind_column});

        obj.current_column=ind_column;

        %Store limits of the subplots
        if sum(sel_column)>0
            %Here we do max/min with the current values for draw after
            %update() calls
            obj.plot_lim.maxx(obj.current_row,obj.current_column)=max(allmax(temp_aes.x(sel_column)),...
                obj.plot_lim.maxx(obj.current_row,obj.current_column));
            obj.plot_lim.minx(obj.current_row,obj.current_column)=min(allmin(temp_aes.x(sel_column)),...
                obj.plot_lim.minx(obj.current_row,obj.current_column));

            obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(allmax(tmp_y_for_max(sel_column)),...
                obj.plot_lim.maxy(obj.current_row,obj.current_column));
            obj.plot_lim.miny(obj.current_row,obj.current_column)=min(allmin(tmp_y_for_min(sel_column)),...
                obj.plot_lim.miny(obj.current_row,obj.current_column));

            if ~isempty(temp_aes.z)
                obj.plot_lim.maxz(obj.current_row,obj.current_column)=max(allmax(temp_aes.z(sel_column)),...
                    obj.plot_lim.maxz(obj.current_row,obj.current_column));
                obj.plot_lim.minz(obj.current_row,obj.current_column)=min(allmin(temp_aes.z(sel_column)),...
                    obj.plot_lim.minz(obj.current_row,obj.current_column));
            else
                obj.plot_lim.maxz(obj.current_row,obj.current_column)=0;
                obj.plot_lim.minz(obj.current_row,obj.current_column)=0;
            end
        end

        %Create subplot (much faster to do it here at the most upper level: changing subplot takes a
        %long time in Matlab)
        if obj.updater.first_draw %If we are in the normal case (first draw() call)
            if obj.wrap_ncols>0
                obj.facet_axes_handles(ind_row,ind_column)=mysubplot_wrap(length(uni_column),ind_column);
            else
                obj.facet_axes_handles(ind_row,ind_column)=mysubplot(length(uni_row),length(uni_column),ind_row,ind_column);
            end
        else %If we are in a draw() call after an update() call
            if obj.updater.facet_updated==1 %If faceting was updated from one to multiple facets
                if ind_column==1 && ind_row==1 %If we are in the first facet
                    %We don't need to create it again, it exists
                    %axes(obj.facet_axes_handles(ind_row,ind_column));
                    %Store content of first facet for copying it on the other facets !
                    first_axes_children=allchild(obj.facet_axes_handles(ind_row,ind_column));
                else
                    %We need to create the next facets because they don't exist
                    if obj.wrap_ncols>0
                        obj.facet_axes_handles(ind_row,ind_column)=mysubplot_wrap(length(uni_column),ind_column);
                    else
                        obj.facet_axes_handles(ind_row,ind_column)=mysubplot(length(uni_row),length(uni_column),ind_row,ind_column);
                    end
                    %And we copy the contents of the first facet in the new ones
                    copyobj(first_axes_children,obj.facet_axes_handles(ind_row,ind_column));
                end
            %else
                %In other cases (same facets or multiple to
                %one facet), the facets already exist
                %axes(obj.facet_axes_handles(ind_row,ind_column));
            end

            if obj.updater.facet_updated==-1 %If facets were updated from many to one facets
                if ind_column==1 && ind_row==1
                    %We store the current content of the first
                    %facet so that we can check which new
                    %things are going to be drawn in
                    first_axes_children=allchild(obj.facet_axes_handles(ind_row,ind_column));
                end
            end
        end

        %Make axes current
        axes(obj.facet_axes_handles(ind_row,ind_column));

        hold on

        %Draw polygons before plotting data, so data isn't covered up
        if obj.polygon.on
            draw_polygons(obj);
        end

        %Store all the X used for the current facet (useful for
        %correct spacing of dodged bars and boxplots when
        %missing data).
        draw_data.facet_x=temp_aes.x(sel_column);

        %Store data for advanced dodging
        [dodge_data.fallback,dodge_data.x,dodge_data.color,dodge_data.lightness,dodge_data.ind,dodge_data.n]=dodge_comp(temp_aes.x(sel_column),temp_aes.color(sel_column),temp_aes.lightness(sel_column),uni_color,uni_lightness);

        %Loop over point shapes
        for ind_marker=1:length(uni_marker)

            sel_marker=sel_column & multi_sel(temp_aes.marker,uni_marker{ind_marker});

            %Loop over line types
            for ind_linestyle=1:length(uni_linestyle)

                sel_linestyle=sel_marker & multi_sel(temp_aes.linestyle,uni_linestyle{ind_linestyle});

                %Loop over sizes
                for ind_size=1:length(uni_size)

                    sel_size=sel_linestyle & multi_sel(temp_aes.size,uni_size{ind_size});

                    %Loop over colors
                    for ind_color=1:length(uni_color)

                        if obj.continuous_color_options.active
                            sel_color=sel_size;
                        else
                            sel_color=sel_size & multi_sel(temp_aes.color,uni_color{ind_color});
                        end

                        %loop over lightness
                        for ind_lightness=1:length(uni_lightness)

                            sel_lightness=sel_color & multi_sel(temp_aes.lightness,uni_lightness{ind_lightness});

                            %We create the groups only within lightness and subplots for speed
                            uni_group=unique_and_sort(temp_aes.group(sel_lightness),0);
                            if isnumeric(uni_group)
                                uni_group=num2cell(uni_group);
                            end

                            %Select dodging parameters for current color
                            %and lightness
                            if obj.continuous_color_options.active
                                sel_dodge=true(size(dodge_data.color));
                            else
                                sel_dodge=multi_sel(dodge_data.color,uni_color{ind_color}) & multi_sel(dodge_data.lightness,uni_lightness{ind_lightness});
                            end

                            draw_data.dodge_fallback=dodge_data.fallback;
                            draw_data.dodge_x=dodge_data.x(sel_dodge);
                            draw_data.dodge_ind=dodge_data.ind(sel_dodge);
                            draw_data.dodge_n=dodge_data.n(sel_dodge);

                            %Loop over groups
                            for ind_group=1:length(uni_group)

                                sel=sel_lightness & multi_sel(temp_aes.group,uni_group{ind_group});

                                %Fill out results struct
                                obj.results.row{obj.result_ind}=uni_row{ind_row};
                                obj.results.column{obj.result_ind}=uni_column{ind_column};
                                obj.results.marker{obj.result_ind}=uni_marker{ind_marker};
                                obj.results.linestyle{obj.result_ind}=uni_linestyle{ind_linestyle};
                                obj.results.size{obj.result_ind}=uni_size{ind_size};
                                obj.results.color{obj.result_ind}=uni_color{ind_color};
                                obj.results.lightness{obj.result_ind}=uni_lightness{ind_lightness};
                                obj.results.group{obj.result_ind}=uni_group{ind_group};

                                obj.results.ind_row{obj.result_ind}=ind_row;
                                obj.results.ind_column{obj.result_ind}=ind_column;
                                obj.results.ind_marker{obj.result_ind}=ind_marker;
                                obj.results.ind_linestyle{obj.result_ind}=ind_linestyle;
                                obj.results.ind_size{obj.result_ind}=ind_size;
                                obj.results.ind_color{obj.result_ind}=ind_color;
                                obj.results.ind_lightness{obj.result_ind}=ind_lightness;
                                obj.results.ind_group{obj.result_ind}=ind_group;

                                if ~isempty(sel)

                                    %Fill up the draw_data
                                    %structure passed to the
                                    %individual geoms
                                    draw_data.x=temp_aes.x(sel);
                                    draw_data.y=temp_aes.y(sel);
                                    if ~isempty(temp_aes.ymin) && ~isempty(temp_aes.ymax)
                                        draw_data.ymin=temp_aes.ymin(sel);
                                        draw_data.ymax=temp_aes.ymax(sel);
                                    else
                                        draw_data.ymin=[];
                                        draw_data.ymax=[];
                                    end
                                    if ~isempty(temp_aes.z)
                                        draw_data.z=temp_aes.z(sel);
                                    else
                                        draw_data.z=[];
                                    end
                                    if ~isempty(temp_aes.label)
                                        draw_data.label=temp_aes.label(sel);
                                    else
                                        draw_data.label=[];
                                    end
                                    draw_data.continuous_color=temp_aes.color(sel);
                                    draw_data.color=cmap((ind_color-1)*length(uni_lightness)+ind_lightness,:);
                                    draw_data.marker=obj.point_options.markers{1+mod(ind_marker-1,length(obj.point_options.markers))};
                                    draw_data.line_style=obj.line_options.styles{1+mod(ind_linestyle-1,length(obj.line_options.styles))};
                                    if obj.line_options.use_input
                                        draw_data.line_size=obj.line_options.input_fun(uni_size{ind_size});
                                    else
                                        draw_data.line_size=obj.line_options.base_size+(ind_size-1)*obj.line_options.step_size;
                                    end
                                    if obj.point_options.use_input
                                        draw_data.point_size=obj.point_options.input_fun(uni_size{ind_size});
                                    else
                                        draw_data.point_size=obj.point_options.base_size+(ind_size-1)*obj.point_options.step_size;
                                    end
                                    draw_data.color_index=(ind_color-1)*length(uni_lightness)+ind_lightness;
                                    draw_data.n_colors=length(uni_color)*length(uni_lightness);

                                    %Loop over geoms
                                    for geom_ind=1:length(obj.geom)

                                        %Call each geom !
                                        obj.geom{geom_ind}(obj,draw_data);

                                    end
                                    obj.firstrun(obj.current_row,obj.current_column)=0;

                                    draw_data=rmfield(draw_data,{'x','y','z','continuous_color','color_index','n_colors'});
                                    obj.results.draw_data{obj.result_ind}=draw_data;

                                    %Iterate loop counter
                                    obj.result_ind=obj.result_ind+1;
                                end

                            end
                        end
                    end
                end
            end
        end

        %Set colormap of subplot if needed
        if obj.continuous_color_options.active
            colormap(obj.continuous_color_options.colormap);
        end

        %Show facet values in titles
        if obj.updater.first_draw || obj.updater.facet_updated
            if length(uni_column)>1 && obj.column_labels
                if ~isempty(obj.aes_names.column)
                    column_string=[obj.aes_names.column ': ' num2str(uni_column{ind_column})];
                else
                    column_string=num2str(uni_column{ind_column});
                end

                if ind_row==1
                    obj.facet_text_handles=[obj.facet_text_handles ...
                        text('Interpreter',obj.text_options.interpreter,'String',column_string,'Rotation',0,...
                        'Units','normalized',...
                        'Position',[0.5 1.05 2],...
                        'BackgroundColor','none',...
                        'HorizontalAlignment','Center',...
                        'VerticalAlignment','bottom',...
                        'FontWeight','bold',...
                        'FontName',obj.text_options.font,...
                        'FontSize',obj.text_options.base_size*obj.text_options.facet_scaling,...
                        'Parent',obj.facet_axes_handles(ind_row,ind_column))];
                end
            end
            if length(uni_row)>1  && obj.row_labels
                if ~isempty(obj.aes_names.row)
                    row_string=[obj.aes_names.row ': ' num2str(uni_row{ind_row})];
                else
                    row_string=num2str(uni_row{ind_row});
                end

                if ind_column==length(uni_column)
                    obj.facet_text_handles=[obj.facet_text_handles ...
                        text('Interpreter',obj.text_options.interpreter,'String',row_string,'Rotation',-90,...
                        'Units','normalized',...
                        'Position',[1.05 0.5 2],...
                        'BackgroundColor','none',...
                        'HorizontalAlignment','Center',...
                        'VerticalAlignment','bottom',...
                        'FontWeight','bold',...
                        'FontName',obj.text_options.font,...
                        'FontSize',obj.text_options.base_size*obj.text_options.facet_scaling,...
                        'Parent',obj.facet_axes_handles(ind_row,ind_column))];
                end
            end
        end

    end

end

%% draw() Title
if ~isempty(obj.title)
    if obj.updater.first_draw
        obj.title_axe_handle=axes('Position',[obj.multi.orig(2)+0.1*obj.multi.size(2)...
            obj.multi.orig(1)+0.90*obj.multi.size(1) 0.8*obj.multi.size(2) 0.05*obj.multi.size(1)],...
            'Parent',obj.parent);

        set(obj.title_axe_handle,'Visible','off','XLim',[-1 1],'YLim',[-1 1]);
        obj.title_text_handle=text(0,0,obj.title,...
            'FontWeight','bold',...
            'Interpreter',obj.text_options.interpreter,...
            'FontName',obj.text_options.font,...
            'FontSize',obj.text_options.base_size*obj.text_options.title_scaling,...
            'HorizontalAlignment','center',...
            'Parent',obj.title_axe_handle);
        if ~isempty(obj.title_options)
            set(obj.title_text_handle,obj.title_options{:});
        end
    end
else
    obj.title_axe_handle=[];
end

%% Compute continuous colormap limits

if obj.continuous_color_options.active
    if isempty(obj.continuous_color_options.CLim)
        obj.continuous_color_options.CLim = [min(min(obj.plot_lim.minc)) max(max(obj.plot_lim.maxc))];
    end
    if obj.continuous_color_options.CLim(1) == obj.continuous_color_options.CLim(2)
        if obj.continuous_color_options.CLim(1) == 0
            obj.continuous_color_options.CLim = [-1 1];
        else
            obj.continuous_color_options.CLim(1) = obj.continuous_color_options.CLim(1) - abs(obj.continuous_color_options.CLim(1))/2;
            obj.continuous_color_options.CLim(2) = obj.continuous_color_options.CLim(2) + abs(obj.continuous_color_options.CLim(2))/2;
        end
    end
end

%% draw() legends

% ensuring grouping variables are cellstrs
str_uni_color = cellfun(@num2str,uni_color,'UniformOutput',false);
str_uni_lightness = cellfun(@num2str,uni_lightness,'UniformOutput',false);
str_uni_linestyle = cellfun(@num2str,uni_linestyle,'UniformOutput',false);
str_uni_size = cellfun(@num2str,uni_size,'UniformOutput',false);
str_uni_marker = cellfun(@num2str,uni_marker,'UniformOutput',false);

%Create axes for legends
if obj.updater.first_draw
    if strcmp(obj.layout_options.legend_position,'auto')
        if strcmp(obj.layout_options.legend_width,'auto')
            obj.legend_axe_handle=axes('Position',[obj.multi.orig(2)+0.85*obj.multi.size(2)...
                obj.multi.orig(1)+0.1*obj.multi.size(1) 0.15*obj.multi.size(2) 0.8*obj.multi.size(1)],...
                'Parent',obj.parent);
        else
            obj.legend_axe_handle=axes('Position',[obj.multi.orig(2)+(1-obj.layout_options.legend_width)*obj.multi.size(2)...
                obj.multi.orig(1)+0.1*obj.multi.size(1) 0.15*obj.multi.size(2) 0.8*obj.multi.size(1)],...
                'Parent',obj.parent);
        end
    else
        obj.legend_axe_handle=axes('Position',obj.layout_options.legend_position,...
                'Parent',obj.parent);
    end
    hold on
    set(obj.legend_axe_handle,'Visible','off','NextPlot','add');
    %set(obj.legend_axe_handle,'PlotBoxAspectRatio',[1 1 1]);
else
    axes(obj.legend_axe_handle)
end

legend_y_step=1;
legend_y_additional_step=0.5;

if obj.layout_options.legend

    % Color legend
    % If the color groups are the same as marker, linestyle or size groups
    % then there will be a common legend (handled by the marker, linestyle
    % or size legend).
    if length(str_uni_color)>1 && any(strcmp(obj.color_options.legend,{'separate','separate_gray'})) % && ...
            %~(length(str_uni_color)==length(str_uni_marker) && all(strcmp(str_uni_color, str_uni_marker))) && ...
            %~(length(str_uni_color)==length(str_uni_linestyle) && all(strcmp(str_uni_color, str_uni_linestyle))) && ...
            %~(length(str_uni_color)==length(str_uni_size) && all(strcmp(str_uni_color, str_uni_size)))

        %Make a colormap with only the colors and no lightness
        color_legend_map=get_colormap(length(str_uni_color),1,obj.color_options);

        fill_legend(obj,obj.aes_names.color,...
            str_uni_color,...
            'line',...
            color_legend_map,...
            'o','-',3,0);
    end

    %Lightness legend
    if length(uni_lightness)>1 && any(strcmp(obj.color_options.legend,{'separate','separate_gray','merge'}))

        if ischar(obj.color_options.map) && strcmp(obj.color_options.map,'lch') && strcmp(obj.color_options.legend,'separate_gray')
            %With LCH we can generate a correct desaturated legend
            lightness_legend_map=pa_LCH2RGB([linspace(obj.color_options.lightness_range(1),obj.color_options.lightness_range(2),length(uni_lightness))' ...
                zeros(length(uni_lightness),1)...
                zeros(length(uni_lightness),1)]);
        else
            %Make a colormap with the first color
            lightness_legend_map=get_colormap(1,length(str_uni_lightness),obj.color_options);
        end

        fill_legend(obj,obj.aes_names.lightness,...
            str_uni_lightness,...
            'line',...
            lightness_legend_map,...
            'o','-',3,0);
    end

    if strcmp(obj.color_options.legend,'expand')
        expanded_legend_map = get_colormap(length(str_uni_color),length(str_uni_lightness),obj.color_options);
        expanded_legend_labels = cell(length(str_uni_color)*length(str_uni_lightness),1);
        for uc = 1:length(str_uni_color)
            for ul = 1:length(str_uni_lightness)
                expanded_legend_labels{ (uc-1) * length(str_uni_lightness) + ul } = [str_uni_color{uc} ' / ' str_uni_lightness{ul} ];
            end
        end

        fill_legend(obj, [obj.aes_names.color ' / ' obj.aes_names.lightness],...
            expanded_legend_labels,...
            'line',...
            expanded_legend_map,...
            'o','-',3,0);
    end

    %Continuous color legend
    if obj.continuous_color_options.active

        fill_gradient_legend(obj,obj.aes_names.color,...
            linspace(obj.continuous_color_options.CLim(1), obj.continuous_color_options.CLim(2), 3),...
            linspace(obj.continuous_color_options.CLim(1), obj.continuous_color_options.CLim(2), 5),...
            4);

    end

    %marker legend
    if length(str_uni_marker)>1

        % If marker groups are the same as color groups we combine the
        % legends by drawing each marker legend the corresponding color
        if strcmp(obj.color_options.legend,'merge') && length(str_uni_color)==length(str_uni_marker) && all(strcmp(str_uni_color, str_uni_marker))
            color_legend_map = get_colormap(length(str_uni_marker), 1, obj.color_options);
        else
            color_legend_map = [0 0 0]; %zeros(length(str_uni_marker), 3); %Otherwise in black
        end

        fill_legend(obj,obj.aes_names.marker,...
            str_uni_marker,...
            'point',...
            color_legend_map,...
            obj.point_options.markers,'-',0,obj.point_options.base_size);
    end

    %linestyle legend
    if length(str_uni_linestyle)>1

        % If linestyle groups are the same as color groups we combine the
        % legends by drawing each linestyle legend the corresponding color
        if strcmp(obj.color_options.legend,'merge') && length(str_uni_color)==length(str_uni_linestyle) && all(strcmp(str_uni_color, str_uni_linestyle))
            color_legend_map = get_colormap(length(str_uni_linestyle), 1, obj.color_options);
        else
            color_legend_map = [0 0 0]; %Otherwise in black
        end

        fill_legend(obj,obj.aes_names.linestyle,...
            str_uni_linestyle,...
            'line',...
            color_legend_map,...
            'o',obj.line_options.styles,3,0);
    end

    %Size legend
    if length(str_uni_size)>1

        % If size groups are the same as color groups we combine the
        % legends by drawing each size legend the corresponding color
        if strcmp(obj.color_options.legend,'merge') && length(str_uni_color)==length(str_uni_size) && all(strcmp(str_uni_color, str_uni_size))
            color_legend_map = get_colormap(length(str_uni_size), 1, obj.color_options);
        else
            color_legend_map = [0 0 0]; %Otherwise in black
        end

        % Handle cases where sizes are defined by the input calues instead
        % of discrete categories ('use_input' in set_line_options() and set_point_options() )
        if obj.line_options.use_input
            temp_lw=obj.line_options.input_fun([uni_size{:}]);
        else
            temp_lw=obj.line_options.base_size+(0:length(str_uni_size)-1)*obj.line_options.step_size;
        end
        if obj.point_options.use_input
            temp_ps=obj.point_options.input_fun([uni_size{:}]);
        else
            temp_ps=obj.point_options.base_size+(0:length(str_uni_size)-1)*obj.point_options.step_size;
        end

        fill_legend(obj,obj.aes_names.size,...
            str_uni_size,...
            'both',...
            color_legend_map,...
            'o','-',temp_lw,temp_ps);
    end
end

%Set size of legend axes
set(obj.legend_axe_handle,'XLim',[1 8]);
%xlim([1 8])
if obj.legend_y<0
    set(obj.legend_axe_handle,'YLim',[obj.legend_y 1]);
    %ylim([obj.legend_y 1])
end
%set(obj.legend_axe_handle,'DataAspectRatio',[1 1 1]);

%If we go from many to one facet, the new content is drawn only
%on the first facet so needs to be copied in the other one
if obj.updater.facet_updated==-1
    %Get all content from the first facet
    new_first_axes_children=allchild(obj.facet_axes_handles(1,1));
    %Check what is new there by comparing the handles with the
    %ones we stored before drawing new content
    todraw=setdiff(new_first_axes_children,first_axes_children);
    %Copy the new stuff to all other facets
    for k=2:numel(obj.facet_axes_handles)
        copyobj(todraw,obj.facet_axes_handles(k));
    end
end

%% draw() axes modifications

%Set various properties on each of the subplots
for ind_row=1:length(uni_row) %Loop over rows

    for ind_column=1:length(uni_column) %Loop over columns

        %Set current axes
        ca = obj.facet_axes_handles(ind_row,ind_column);

        set(ca,'FontName',obj.text_options.font,...
            'FontSize',obj.text_options.base_size)

        if obj.continuous_color_options.active
            %Set color limits the same way on each plot
            set(ca,'CLimMode','manual','CLim',obj.continuous_color_options.CLim);
        end

        %Ad hoc limit correction for empty facets
        obj.plot_lim.maxy(obj.plot_lim.miny==obj.plot_lim.maxy)=obj.plot_lim.maxy(obj.plot_lim.miny==obj.plot_lim.maxy)+0.01;
        obj.plot_lim.maxx(obj.plot_lim.minx==obj.plot_lim.maxx)=obj.plot_lim.maxx(obj.plot_lim.minx==obj.plot_lim.maxx)+0.01;
        obj.plot_lim.maxz(obj.plot_lim.minz==obj.plot_lim.maxz)=obj.plot_lim.maxz(obj.plot_lim.minz==obj.plot_lim.maxz)+0.01;

        if ~obj.polar.is_polar % XY Limits are only useful for non-polar plots

            %Set axes limits logic according to facet_scale and
            %wrapping. Also set up corresponding axis linking, using tip from
            % http://undocumentedmatlab.com/blog/using-linkaxes-vs-linkprop#more-5928
            if (obj.wrap_ncols>0) %in the facet_wrap case
                switch obj.facet_scale
                    case 'fixed'
                        temp_xscale='global';
                        temp_yscale='global';
                        temp_zscale='global';

                        %Both XLims and YLims are linked across
                        %all plots
                        if ind_row==1 && ind_column==1 %We only do the linking once for all plots facets
                            obj.extra.YLim_listeners=linkprop(obj.facet_axes_handles(:),'YLim');
                            obj.extra.XLim_listeners=linkprop(obj.facet_axes_handles(:),'XLim');
                            obj.extra.ZLim_listeners=linkprop(obj.facet_axes_handles(:),'ZLim');
                        end
                    case 'free_x'
                        temp_xscale='per_plot';
                        temp_yscale='global';
                        temp_zscale='global';

                        %YLims are linked across all plots
                        if ind_row==1 && ind_column==1
                            obj.extra.YLim_listeners=linkprop(obj.facet_axes_handles(:),'YLim');
                        end
                    case 'free_y'
                        temp_xscale='global';
                        temp_yscale='per_plot';
                        temp_zscale='global';

                        %XLims are linked across all plots
                        if ind_row==1 && ind_column==1
                            obj.extra.XLim_listeners=linkprop(obj.facet_axes_handles(:),'XLim');
                        end
                    case 'free'
                        temp_xscale='per_plot';
                        temp_yscale='per_plot';
                        temp_zscale='per_plot';
                    case 'independent'
                        temp_xscale='per_plot';
                        temp_yscale='per_plot';
                        temp_zscale='per_plot';
                end
            else %In the facet_grid case
                switch obj.facet_scale
                    case 'fixed'
                        temp_xscale='global';
                        temp_yscale='global';
                        temp_zscale='global';

                        %Both XLim and YLim are linked across
                        %all plots
                        if ind_row==1 && ind_column==1
                            obj.extra.YLim_listeners=linkprop(obj.facet_axes_handles(:),'YLim');
                            obj.extra.XLim_listeners=linkprop(obj.facet_axes_handles(:),'XLim');
                            obj.extra.ZLim_listeners=linkprop(obj.facet_axes_handles(:),'ZLim');
                        end
                    case 'free_x'
                        temp_xscale='per_column';
                        temp_yscale='global';
                        temp_zscale='global';

                        %YLims are linked across all plots,
                        if ind_row==1 && ind_column==1
                            obj.extra.YLim_listeners=linkprop(obj.facet_axes_handles(:),'YLim');
                        end
                        if obj.is_flipped
                            %For flipped axes, we link XLims within rows
                            if ind_column==1
                                obj.extra.XLim_listeners(ind_row)=linkprop(obj.facet_axes_handles(ind_row,:),'XLim');
                            end
                        else
                            %XLims are linked within columns for
                            %non-flipped axes
                            if ind_row==1
                                obj.extra.XLim_listeners(ind_column)=linkprop(obj.facet_axes_handles(:,ind_column),'XLim');
                            end
                        end
                    case 'free_y'
                        temp_xscale='global';
                        temp_yscale='per_row';
                        temp_zscale='global';

                        %XLims are linked across all plots,
                        if ind_row==1 && ind_column==1
                            obj.extra.XLim_listeners=linkprop(obj.facet_axes_handles(:),'XLim');
                        end
                        if obj.is_flipped
                            %For flipped axes, we link YLims within columns
                            if ind_row==1
                                obj.extra.YLim_listeners(ind_column)=linkprop(obj.facet_axes_handles(ind_column,:),'YLim');
                            end
                        else
                            %YLims are linked within rows for non-flipped
                            %axes
                            if ind_column==1
                                obj.extra.YLim_listeners(ind_row)=linkprop(obj.facet_axes_handles(ind_row,:),'YLim');
                            end
                        end
                    case 'free'
                        temp_xscale='per_column';
                        temp_yscale='per_row';

                        if obj.is_flipped
                            %Flipped case
                            %XLims are linked within rows
                            %YLims are linked within columns
                            if ind_row==1
                                obj.extra.YLim_listeners(ind_column)=linkprop(obj.facet_axes_handles(:,ind_column),'YLim');
                            end
                            if ind_column==1
                                obj.extra.XLim_listeners(ind_row)=linkprop(obj.facet_axes_handles(ind_row,:),'XLim');
                            end
                        else
                            %Non flipped case
                            %XLims are linked within columns
                            %YLims are linked within rows
                            if ind_row==1
                                obj.extra.XLim_listeners(ind_column)=linkprop(obj.facet_axes_handles(:,ind_column),'XLim');
                            end
                            if ind_column==1
                                obj.extra.YLim_listeners(ind_row)=linkprop(obj.facet_axes_handles(ind_row,:),'YLim');
                            end
                        end
                    case 'independent'
                        temp_xscale='per_plot';
                        temp_yscale='per_plot';
                        temp_zscale='per_plot';
                end
            end

            %Actually set the axes scales and presence of
            %labels depending on logic determined above
            switch temp_xscale
                case 'global'
                    temp_xlim=[min(min(obj.plot_lim.minx(:,:))) max(max(obj.plot_lim.maxx(:,:)))];
                case 'per_column'
                    if obj.is_flipped
                        temp_xlim=[min(obj.plot_lim.minx(ind_row,:)) max(obj.plot_lim.maxx(ind_row,:))];
                    else
                        temp_xlim=[min(obj.plot_lim.minx(:,ind_column)) max(obj.plot_lim.maxx(:,ind_column))];
                    end
                case 'per_plot'
                    temp_xlim=[obj.plot_lim.minx(ind_row,ind_column) obj.plot_lim.maxx(ind_row,ind_column)];
            end
            if sum(isnan(temp_xlim))==0
                set(ca,'XLim',temp_xlim+[-diff(temp_xlim)*obj.xlim_extra(1) diff(temp_xlim)*obj.xlim_extra(2)]);
            end

            switch temp_yscale
                case 'global'
                    temp_ylim=[min(min(obj.plot_lim.miny(:,:))) max(max(obj.plot_lim.maxy(:,:)))];
                case 'per_row'
                    if obj.is_flipped
                        temp_ylim=[min(obj.plot_lim.miny(:,ind_column)) max(obj.plot_lim.maxy(:,ind_column))];
                    else
                        temp_ylim=[min(obj.plot_lim.miny(ind_row,:)) max(obj.plot_lim.maxy(ind_row,:))];
                    end
                case 'per_plot'
                    temp_ylim=[obj.plot_lim.miny(ind_row,ind_column) obj.plot_lim.maxy(ind_row,ind_column)];
            end
            if sum(isnan(temp_ylim))==0
                set(ca,'YLim',temp_ylim+[-diff(temp_ylim)*obj.ylim_extra(1) diff(temp_ylim)*obj.ylim_extra(2)]);
            end

            %Is the axis at interesting locations for ticks/labels?
            on_left=(obj.wrap_ncols==-1 && ind_column==1) || ... %Or if we're in facet grid mode and are in the first column
                    (obj.wrap_ncols>0 && mod(ind_column,obj.wrap_ncols)==1) || ... %Or if we are in facet wrap mode and are in the first "column"
                    (obj.wrap_ncols==1 && ind_column==1); %Special case for single facet in facet_wrap mode

            on_bottom=(obj.wrap_ncols==-1 && ind_row==length(uni_row)) || ... %Or if we're in facet grid mode and we are in the last row
                    (obj.wrap_ncols>0 && (length(uni_column)-ind_column)<obj.wrap_ncols); %Or if we are in facet wrap mode and we are in the last facet on the "column"

            if ~isempty(temp_aes.z) %Only do the Z limit stuff if we have z data

                switch temp_zscale
                    case 'global'
                        temp_zlim=[min(min(obj.plot_lim.minz(:,:))) max(max(obj.plot_lim.maxz(:,:)))];
                    case 'per_plot'
                        temp_zlim=[obj.plot_lim.minz(ind_row,ind_column) obj.plot_lim.maxz(ind_row,ind_column)];
                end
                set(ca,'ZLim',temp_zlim+[-diff(temp_zlim)*obj.zlim_extra(1) diff(temp_zlim)*obj.zlim_extra(2)]);

                %Always have ticks if we have z data
                has_xtick=true;
                has_ytick=true;

            else %Only do optional xy ticks if we don't have z data

                %Set up logic of plot ticks presence
                has_xtick=obj.force_ticks || ... %Plot has xticks if forced
                    on_bottom ||...
                    strcmp(temp_xscale,'per_plot'); %Or if we were in a per-plot scale mode

                has_ytick=obj.force_ticks || ... %Plot has xticks if forced
                    on_left || ... %Special case for single facet in facet_wrap mode
                    strcmp(temp_yscale,'per_plot'); %Or if we were in a per-plot scale mode
            end

            %flip tick presence if coordinates are flipped
            if obj.is_flipped
                tmp=has_ytick;
                has_ytick=has_xtick;
                has_xtick=tmp;
            end

            %Remove ticks if necessary
            if ~has_xtick
                set(ca,'XTickLabel','');
            end
            if ~has_ytick
                set(ca,'YTickLabel','');
            end

            %Set appropriate x ticks if labeled
            if obj.x_factor
                temp_xlim=get(ca,'xlim');
                set(ca,'XLim',[temp_xlim(1)-0.6 temp_xlim(2)+0.6]);
                set(ca,'XTick',1:length(obj.x_ticks))
                if has_xtick
                    set(ca,'XTickLabel',obj.x_ticks)
                    if isprop(ca,'TickLabelInterpreter') %Just try it (doesn't exist in pre-2014b)
                        set(ca,'TickLabelInterpreter','none')
                    end
                end
            end

            %Add axes labels on right and botttom graphs only
            if (~obj.is_flipped && on_left) || ~isempty(temp_aes.z) || ...
               (obj.is_flipped && on_bottom)   %If coord flipped do it the other way around (y like x)
                ylabel(ca,obj.aes_names.y,...
                    'Interpreter',obj.text_options.interpreter,...
                    'FontName',obj.text_options.font,...
                    'FontSize',obj.text_options.base_size*obj.text_options.label_scaling); %,'Units','normalized','position',[-0.2 0.5 1]
            end
            if (~obj.is_flipped && on_bottom) || ~isempty(temp_aes.z) ||...
                (obj.is_flipped && on_left) %If coord flipped do it the other way around (x like y)
                xlabel(ca,obj.aes_names.x,...
                    'Interpreter',obj.text_options.interpreter,...
                    'FontName',obj.text_options.font,...
                    'FontSize',obj.text_options.base_size*obj.text_options.label_scaling)
            end
            %If we have z data
            if ~isempty(temp_aes.z)
                zlabel(ca,obj.aes_names.z,...
                    'Interpreter',obj.text_options.interpreter,...
                    'FontName',obj.text_options.font,...
                    'FontSize',obj.text_options.base_size*obj.text_options.label_scaling) %Add z label
                view(ca,3); %Reset the view so that it is a 3D view
                if ind_row==1 && ind_column==1
                    %Link the camera properties between the
                    %facets so that they all rotate together
                    obj.extra.CamPos_listeners=linkprop(obj.facet_axes_handles(:),'CameraPosition');
                    obj.extra.CamTgt_listeners=linkprop(obj.facet_axes_handles(:),'CameraTarget');
                    obj.extra.CamUp_listeners=linkprop(obj.facet_axes_handles(:),'CameraUpVector');
                end
            end

        else
            %Make polar axes
            if obj.polar.max_polar_y<0
                if strcmp(obj.facet_scale,'fixed')
                    draw_polar_axes(ca,max(max(obj.plot_lim.maxy(:,:))));
                else
                    draw_polar_axes(ca,obj.plot_lim.maxy(ind_row,ind_column));
                end
            else
                draw_polar_axes(ca,obj.polar.max_polar_y)
            end
        end

        %Set custom axes properties
        if ~isempty(obj.axe_properties)
            for ap=1:size(obj.axe_properties,1)
                if isprop(ca,obj.axe_properties{ap,1})
                    set(ca,obj.axe_properties{ap,1},obj.axe_properties{ap,2})
                else
                    warning(['Improper ''' obj.axe_properties{ap,1} ''' custom axe property'])
                end
            end
        end

        %Do the datetick
        if ~isempty(obj.datetick_params)
            for dtk=1:length(obj.datetick_params)
                %Previously used 'keepticks' by default  which could cause a truncation of the axes
                datetick(ca,obj.datetick_params{dtk}{:});
            end
        end

        %Set ablines, hlines and vlines (after axe properties in case the limits
        %are changed there
        if obj.abline.on
            %xl=get(ca,'xlim');
            for line_ind=1:length(obj.abline.intercept)
                tmp_xl=[obj.var_lim.minx obj.var_lim.maxx];
                tmp_extent=(tmp_xl(2)-tmp_xl(1))*obj.abline.extent(line_ind)/2;
                xl=[mean(tmp_xl)-tmp_extent mean(tmp_xl)+tmp_extent];
                if ~isnan(obj.abline.intercept(line_ind))
                    %abline
                    plot(xl,xl*obj.abline.slope(line_ind)+obj.abline.intercept(line_ind),obj.abline.style{line_ind},'LineWidth',obj.abline.linewidth(line_ind),'Parent',ca);
                else
                    if ~isnan(obj.abline.xintercept(line_ind))
                        %vline
                        %yl=get(ca,'ylim');
                        if obj.var_lim.miny == obj.var_lim.maxy %We are probably in a case where y wasn't provided (histogram or raster)
                            tmp_yl=[0 numel(temp_aes.x)]; %We scale y according to number of x elements
                        else
                            tmp_yl=[obj.var_lim.miny obj.var_lim.maxy];
                        end
                        tmp_extent=(tmp_yl(2)-tmp_yl(1))*obj.abline.extent(line_ind)/2;
                        yl=[mean(tmp_yl)-tmp_extent mean(tmp_yl)+tmp_extent];
                        plot([obj.abline.xintercept(line_ind) obj.abline.xintercept(line_ind)],yl,obj.abline.style{line_ind},'LineWidth',obj.abline.linewidth(line_ind),'Parent',ca);
                    else
                        if ~isnan(obj.abline.yintercept(line_ind))
                            %hline
                            plot(xl,[obj.abline.yintercept(line_ind) obj.abline.yintercept(line_ind)],obj.abline.style{line_ind},'LineWidth',obj.abline.linewidth(line_ind),'Parent',ca);
                        else
                            temp_x=linspace(xl(1),xl(2),500);
                            plot(temp_x,obj.abline.fun{line_ind}(temp_x),obj.abline.style{line_ind},'LineWidth',obj.abline.linewidth(line_ind),'Parent',ca);
                        end
                    end
                end
            end

        end

    end
end

%White background !
set(gcf,'color','w','PaperPositionMode','auto');

% Make everything tight and set the resize function so that it stays so
if do_redraw  && ~obj.multi.active %Redrawing for multiple plots is handled at the beginning of draw()
    redraw(obj);
    if verLessThan('matlab','8.4')
        set(gcf,'ResizeFcn',@(a,b)redraw(obj));
    else
        set(gcf,'SizeChangedFcn',@(a,b)redraw(obj));
    end
end

%Clean up results
result_fields=fieldnames(obj.results);
if obj.result_ind>1
    for rf=1:length(result_fields)
        %Resize
        obj.results.(result_fields{rf})=obj.results.(result_fields{rf})(1:obj.result_ind-1);
        %Transform non cell strings in arrays
        if ~iscellstr(obj.results.(result_fields{rf}))
            obj.results.(result_fields{rf})=[obj.results.(result_fields{rf}){:}]';
        end
    end
else
    obj.results=[];
end

obj.updater.first_draw=false;
end

function out=allmax(in)
%Return maximum of an array or of a cell of arrays
if iscell(in)
    out=max(cellfun(@max,in(~cellfun(@isempty,in)))); %Cellfun on non empty cells
    if isempty(out)
        out=NaN;
    end
else
    out=max(in);
end
end

function out=allmin(in)
%Return minimum of an array or of a cell of arrays
if iscell(in)
    out=min(cellfun(@min,in(~cellfun(@isempty,in))));  %Cellfun on non empty cells
    if isempty(out)
        out=NaN;
    end
else
    out=min(in);
end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function [ node_num, edge_num ] = grf_example_size ( )

%*****************************************************************************80
%
%% GRF_EXAMPLE_SIZE sizes a GRF example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Skiena,
%    Implementing Discrete Mathematics,
%    Combinatorics and Graph Theory with Mathematica,
%    Addison-Wesley, 1990.
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer EDGE_NUM, the number of edges.
%
  node_num = 28;
  edge_num = 84;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%TEST_DGT  Test DGT full window backend
%
%  This script runs a throrough test of the COMP_DGT_FAC
%  and COMP_IDGT_FAC testing them on a range of input parameters.
%
%  Use TEST_WFAC first, to verify that COMP_WFAC and COMP_IWFAC
%  are working, since this tester depends on them.

Lr=[24,16,144,108,144,24,135,35,77,20];
ar=[ 4, 4,  9,  9, 12, 6,  9, 5, 7, 1];
Mr=[ 6, 8, 16, 12, 24, 8,  9, 7,11,20];

test_failed=0;

disp('--- Used subroutines ---');

which comp_wfac
which comp_iwfac
which comp_dgt_fac
which comp_idgt_fac

for ii=1:length(Lr);

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  b=L/M;
  N=L/a;
  c=gcd(a,M);
  d=gcd(b,N);
  p=a/c;
  q=M/c;

  for W=1:3

    for R=1:3

      for rtype=1:2
	if rtype==1
	  rname='REAL ';
	  f=tester_rand(L,W);
	  g=tester_rand(L,R);
	else
	  rname='CMPLX';
	  f=tester_crand(L,W);
	  g=tester_crand(L,R);
	end;

	gf=comp_wfac(g,a,M);
	cc=comp_dgt_fac(f,gf,a,M);
	cc2=ref_dgt(f,g,a,M);

	res=norm(cc(:)-cc2(:));

	failed='';
	if res>10e-10
	  failed='FAILED';
	  test_failed=test_failed+1;
	end;

	s=sprintf('DGT  %s L:%3i W:%2i R:%2i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i %0.5g %s',rname,L,W,R,a,b,c,d,p,q,res,failed);
	disp(s)

      end;

      for rtype=1:2
	if rtype==1
	  rname='REAL ';
	  g=tester_rand(L,R);
	else
	  rname='CMPLX';
	  g=tester_crand(L,R);
	end;

	cc=tester_crand(M,N*R*W);

	gf=comp_wfac(g,a,M);
	f=comp_idgt_fac(ifft(cc)*sqrt(M),gf,L,a,M);
	f2=ref_idgt(reshape(cc,M*N*R,W),g,a,M);

	res=norm(f(:)-f2(:));

	failed='';
	if res>10e-10
	  failed='FAILED';
	  test_failed=test_failed+1;
	end;

	s=sprintf('IDGT %s L:%3i W:%2i R:%2i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i %0.5g %s',rname,L,W,R,a,b,c,d,p,q,res,failed);
	disp(s)

      end;

    end;

  end;

end;

test_failed

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ x_wave ] = wavethresh_2D( x, waveletStages, waveletFilterName, wavecoeffS, lambdaWave, L )
%
% (c) Marc Fischer, Thomas Kuestner
% ---------------------------------------------------------------------

    x_wave_helper = wavedec2(x,waveletStages,waveletFilterName);
    x_wave_helper = softthresh_real(x_wave_helper,lambdaWave/L);
    x_wave = waverec2(x_wave_helper,wavecoeffS,waveletFilterName);

end

author: thomaskuestner


text:
function hdr = makeAnalyzeHeader(datatype, fname)
%
% hdr = makeAnalyzeHeader(datatype, fname)
%
% makes an Analyze-format header struct with many default values pre-loaded.
%
%
%

if(~exist('datatype','var') | isempty(datatype))
    datatype = 4;
end
if(~exist('fname','var') | isempty(fname))
    fname = 'unknown';
end

% compute intensity range & bits per pix from datatype
%
glmin = 0;
switch(datatype)
case 1,
    bitpix = 1;
    glmax = 1;
case 2,
    bitpix = 8;
    glmax = 255;
case 4,
    bitpix = 16;
    glmax = 32767;
case 8,
    bitpix = 32;
    glmax = (2^31-1);
case 16,
    bitpix = 32;
    glmax = 1;
case 64,
    bitpix = 64;
    glmax = 1;
otherwise,
    error(['Unknown datatype ',num2str(datatype),'!']);
end

hdr.data_type = fixLen('dsr', 10);
hdr.db_name = fixLen(fname,18);
hdr.extents = 16384; % Should be 16384, the image file is created as contiguous with a minimum
hdr.session_error = 0;
hdr.regular = 'r';   % Must be `r' to indicate that all images and volumes are the same size.
hdr.hkey_un0 = 0;

% write (struct) image_dimension
%
hdr.dim = fixLen(dim,8);
hdr.vox_units = fixLen('mm',4);
hdr.cal_units = fixLen(' ',8);
hdr.unused1 = fixLen(0,1);
hdr.datatype = fixLen(datatype,1);
hdr.bitpix = fixLen(bitpix,1);
hdr.dim_un0 = fixLen(0,1);
hdr.pixdim = fixLen(pixdim,8);
hdr.vox_offset = fixLen(vox_offset,1);
hdr.funused1 = fixLen(funused1,1);
hdr.funused2 = fixLen(funused2,1);
hdr.funused3 = fixLen(funused3,1);
hdr.cal_max = fixLen(cal_max,1);
hdr.cal_min = fixLen(cal_min,1);
hdr.compressed = fixLen(compressed,1);
hdr.verified = fixLen(verified,1);
hdr.glmax = fixLen(glmax,1);
hdr.glmin = fixLen(glmin,1);

% write (struct) data_history
%---------------------------------------------------------------------------
hdr.descrip = fixLen(descrip,80),	'char');
hdr.aux_file = fixLen(aux_file,24), 'char');
% valid values for orient:
% 0 transverse unflipped, 1 coronal unflipped, 2 sagittal unflipped,
% 3 transverse flipped, 4 coronal flipped, 5 sagittal flipped
hdr.orient = fixLen(orient,1), 'char');
hdr.origin = fixLen(origin,5), 'int16');
hdr.generated = fixLen(generated,10), 'uchar');
hdr.scannum = fixLen(scannum,10), 'uchar');
hdr.patient_id = fixLen(patient_id,10), 'uchar');
hdr.exp_date = fixLen(exp_date,10), 'uchar');
hdr.exp_time = fixLen(exp_time,10), 'uchar');
hdr.hist_un0 = fixLen(hist_un0,3), 'uchar');
hdr.views = fixLen(views,1), 'int32');
hdr.vols_added = fixLen(vols_added,1), 'int32');
hdr.start_field = fixLen(start_field,1), 'int32');
hdr.field_skip = fixLen(field_skip,1), 'int32');
hdr.omax = fixLen(omax,1), 'int32');
hdr.omin = fixLen(omin,1), 'int32');
hdr.smax = fixLen(smax,1), 'int32');
hdr.smin = fixLen(smin,1), 'int32');

numOut = ftell(fid);
fclose(fid);
return;

function m = fixLen(m, len)
m = m(:)';
if(length(m) > len)
    m = m(1:len);
elseif(length(m) < len)
    if(ischar(m))
        m = [m char(zeros(1,len-length(m)))];
    else
        m = [m zeros(1,len-length(m))];
    end
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [gradISE,gradJhr,gradJrr]=GaussMixISEWeightGrads(w1,w2,sqrtw2,PDFVals12,PDFVals22)
%%GAUSSMIXISEWEIGHTGRADS Compute the gradient of the non-normalized
%               integrated squared error (ISE) between two Gaussian mixture
%               PDFs with respect to the square root of the weight terms of
%               the second Gaussian mixture. This gradient can arise when
%               optimizing over the ISE with respect to the square root of
%               the weights. The square root of the weights is typically
%               used so as to ensure that no weight ever becomes negative.
%
%INPUTS: w1 The 1Xn1 or n1X1 set of weights of the first Gaussian mixture
%           distribution. It is required that all w1>=0 and sum(w1)=1.
%        w2 The 1Xn2 or n2X1 set of weights of the second Gaussian mixture
%           distribution. It is required that all w2>=0 and sum(w2)=1.
%    sqrtw2 This is just sqrt(w2) and is the term with which the gradient
%           is being taken.
% PDFVals12 A matrix such that the value in element (i,j) is
%           N(mu1(:,i);mu2(:,j),P1(:,:,i)+P2(:,:,j)), where N indicates the
%           multivariate Gaussian PDF evaluated at the first argument with
%           the second and third arguments being the mean and covarince
%           matrix. mu1 and P1 are the means and covariance matrices of the
%           first Gaussian mixture distribution and mu2 and P2 are the same
%           for the second Gaussian mixture distribution. This parameter is
%           returned by computeGaussMixISE.
% PDFVals22 A matrix such that the value in element (i,j) is
%           N(mu2(:,i);mu2(:,j),P2(:,:,i)+P2(:,:,j)). This parameter is
%           returned by computeGaussMixISE.
%
%OUTPUTS: gradISE The 1Xn2 set of derivatives of the ISE wtih respect to
%                 the elements of q (the square roots of w2).
% gradJhr, gradJrr In Chapter 3 of [1], the ISE is expressed in terms of
%                 Jhr and Jrr terms. These are the 1Xn2 gradients of those
%                 terms.
%
%Formule for gradJhr and gradJrr are Equation 3.31 in Section 3.3.3.1 of
%[1]. They relate to the ISE via Equation 3.20. See the function
%computeGaussMixISE to compute the ISE.
%
%EXAMPLE:
%In this example with a scalar PDF, we verify that the gradient obtained
%from this function is consistent with numerical differentiation.
% w1=[0.03,0.18,0.12,0.19,0.02,0.16,0.06,0.1,0.08,0.06];
% n1=length(w1);
% mu1=[1.45,2.20,0.67,0.48,1.49,0.91,1.01,1.42,2.77,0.89];
% P1=[0.0487,0.0305,0.1171,0.0174,0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679];
% P1=reshape(P1,[1,1,n1]);
%
% %The second PDF is the first with the five least-weight components deleted.
% w2=[0.18,0.12,0.19,0.16,0.1,0.08];
% w2=w2/sum(w2);
% n2=length(w2);
% mu2=[2.20,0.67,0.48,0.91,1.42,2.77];
% P2=[0.0305,0.1171,0.0174,0.0102,0.0380,0.0115];
% P2=reshape(P2,[1,1,n2]);
%
% [ISEVal,PDFVals12,PDFVals22]=computeGaussMixISE(w1,mu1,P1,w2,mu2,P2);
% q2=sqrt(w2);
% epsVal=1e-9;
% gradISENum=zeros(1,n2);
% for k=1:n2
%     q2Cur=q2;
%     q2Cur(k)=q2Cur(k)+epsVal;
%     ISEValCur=computeGaussMixISE(w1,mu1,P1,q2Cur.^2,mu2,P2);
%     gradISENum(k)=(ISEValCur-ISEVal)/epsVal;
% end
% gradISE=GaussMixISEWeightGrads(w1,w2,PDFVals12,PDFVals22);
% RelErr=max(abs((gradISENum-gradISE)./gradISENum))
%The relative error will be about 2.1628e-6, which indicates good numeric
%agreement.
%
%REFERENCES:
%[1] J. L. Williams, "Gaussian mixture reduction for tracking multiple
%    maneuvering targets in clutter," Master's thesis, Air Force Institute
%    of Technology, Mar. 2003. [Online].
%    Available: http://www.dtic.mil/srch/doc?collection=t3&id=ADA415317
%
%May 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

q=sqrtw2;

%Both formulae are Equation 3.31 in Section 3.3.3.1 of [1].
gradJhr=2*q(:).'.*sum(bsxfun(@times,w1(:),PDFVals12),1);
gradJrr=4*q(:).'.*sum(bsxfun(@times,w2(:),PDFVals22),1);

gradISE=gradJrr-2*gradJhr;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function tests = test_ft_preproc_baselinecorrect

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_baselinecorrect

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

result = [];
result{end+1} = ft_preproc_baselinecorrect(dat, 1, 1000);
result{end+1} = ft_preproc_baselinecorrect(dat, 1, 900);
result{end+1} = ft_preproc_baselinecorrect(dat, 101, 1000);
result{end+1} = ft_preproc_baselinecorrect(dat, 101, 900);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function MimSaveAsNifti(image_to_save, path, filename, reporting)
    % MimSaveAsNifti. Writes out a PTKImage in NIFTI format
    %
    %     Syntax
    %     ------
    %
    %         MimSaveAsNifti(image_data, path, filename, data_type, orientation, reporting)
    %
    %             image_to_save   is a PTKImage (or PTKDicomImage) class containing the image
    %                             to be saved
    %             path, filename  specify the location to save the NIFTI data.
    %             reporting       an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if nargin < 4
        reporting = CoreReportingDefault();
    end

    if ~isa(image_to_save, 'PTKImage')
        reporting.Error('MimSaveAsNifti:InputMustBePTKImage', 'Requires a PTKImage as input');
    end

    image_data = image_to_save.RawImage;

    full_filename = fullfile(path, filename);

    resolution = image_to_save.VoxelSize([2, 1, 3]);

    offset = [0 0 0];
    if isa(image_data, 'PTKDicomImage')
        metadata = image_data. MetaHeader;
        if isfield(metadata, 'Offset')
            offset = metadata.Offset;
        end
    end

    image_data = permute(image_data, [2, 1, 3]);
    image_data = flip(flip(flip(image_data, 3), 2), 1);

    nii_data = make_nii(image_data, resolution, offset, [], image_to_save.Title);
    save_nii(nii_data, full_filename);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout = imqacdemo03(varargin)
% IMQACDEMO03 M-file for imqacdemo03.fig
%      IMQACDEMO03, by itself, creates a new IMQACDEMO03 or raises the existing
%      singleton*.
%
%      H = IMQACDEMO03 returns the handle to a new IMQACDEMO03 or the
%      handle to
%      the existing singleton*.
%
%      IMQACDEMO03('CALLBACK',hObject,eventData,handles,...) calls the
%      local
%      function named CALLBACK in IMQACDEMO03.M with the given input arguments.
%
%      IMQACDEMO03('Property','Value',...) creates a new IMQACDEMO03 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before imqacdemo03_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to imqacdemo03_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help imqacdemo03

% Last Modified by GUIDE v2.5 23-Jul-2004 12:45:11

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @imqacdemo03_OpeningFcn, ...
                   'gui_OutputFcn',  @imqacdemo03_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before imqacdemo03 is made visible.
function imqacdemo03_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to imqacdemo03 (see VARARGIN)

handles.output = hObject;
Monitor_Mode(hObject, eventdata, handles);
S = imread('autumn.tif');
handles.S = S;
axes(handles.axes1);
subimage(S);
S2 = imread('autumn.tif');
axes(handles.axes2);
subimage(S2);
handles.S2 = S2;
guidata(hObject, handles);

% UIWAIT makes imqacdemo03 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = imqacdemo03_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pbStart.
function pbStart_Callback(hObject, eventdata, handles)
% hObject    handle to pbStart (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'outdata')
    handles = rmfield(handles,'outdata');
end
guidata(hObject, handles);
set(handles.pbStart,'enable','off');
set(handles.pbStop,'enable','on');

set(handles.pbStop,'UserData',0);

vid = videoinput('winvideo');
handles.vid = vid;
set(handles.vid,'FramesPerTrigger',1);
set(handles.vid,'TriggerRepeat',Inf);
triggerconfig(handles.vid, 'Manual');

guidata(hObject, handles);

start(handles.vid);
trigger(handles.vid);
y = (getdata(handles.vid,1,'uint8'));
cnt_move = 0;
cnt_unmove = 0;
cntsnap = 1;
while 1
    trigger(handles.vid);
    yprev = y;
    set(handles.editTime,'string',datestr(clock));
    if get(handles.pbStop,'UserData')
        break
    else
        y = (getdata(handles.vid,1,'uint8'));
        diff = abs(y-yprev);
        abs_img = mean(diff(:));
        set(handles.editAbs,'string',abs_img);
        axes(handles.axes1);subimage(y);
        axes(handles.axes2);subimage(diff);
        out = imaqmem;
        mem_left = out.FrameMemoryLimit - out.FrameMemoryUsed;
        set(handles.editMemLeft,'string',mem_left/10^6);
        set(handles.editMemLoad,'string',out.MemoryLoad);
        if abs_img > str2num(get(handles.editThresh,'string'));
            cnt_move = cnt_move + 1;
            set(handles.editObj,'string',cnt_move);
            if cnt_move >= str2num(get(handles.editSnapCnt,'string'));
                handles.outdata.image(:,:,:,cntsnap) = y;
                handles.outdata.time{cntsnap} = datestr(clock);
                cntsnap = cntsnap + 1;
                cnt_move = 0;
                set(handles.editSnap,'string',cntsnap-1);
            end
            cnt_unmove = 0;
        else
            cnt_unmove = cnt_unmove + 1;
            cnt_move = 0;
            set(handles.editObj,'string',cnt_move);
        end

    end
end

axes(handles.axes1);
cla;
subimage(handles.S);
axes(handles.axes2);
cla;
subimage(handles.S2);
delete(handles.vid);
clear handles.vid;
imaqreset;
guidata(hObject, handles);

% --- Executes on button press in pbStop.
function pbStop_Callback(hObject, eventdata, handles)
% hObject    handle to pbStop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.pbStart,'enable','on');
set(handles.pbStop,'enable','off');

set(handles.pbStop,'UserData',1);
set(handles.pbStop,'UserData',1);
set(handles.pbStop,'UserData',1);
set(handles.pbStop,'UserData',1);
set(handles.pbStop,'UserData',1);
guidata(hObject, handles);

function editAbs_Callback(hObject, eventdata, handles)
% hObject    handle to editAbs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editAbs as text
%        str2double(get(hObject,'String')) returns contents of editAbs as a double

% --- Executes during object creation, after setting all properties.
function editAbs_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editAbs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editMemLeft_Callback(hObject, eventdata, handles)
% hObject    handle to editMemLeft (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editMemLeft as text
%        str2double(get(hObject,'String')) returns contents of editMemLeft as a double

% --- Executes during object creation, after setting all properties.
function editMemLeft_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editMemLeft (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editMemLoad_Callback(hObject, eventdata, handles)
% hObject    handle to editMemLoad (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editMemLoad as text
%        str2double(get(hObject,'String')) returns contents of editMemLoad as a double

% --- Executes during object creation, after setting all properties.
function editMemLoad_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editMemLoad (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editThresh_Callback(hObject, eventdata, handles)
% hObject    handle to editThresh (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editThresh as text
%        str2double(get(hObject,'String')) returns contents of editThresh as a double

% --- Executes during object creation, after setting all properties.
function editThresh_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editThresh (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editObj_Callback(hObject, eventdata, handles)
% hObject    handle to editObj (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editObj as text
%        str2double(get(hObject,'String')) returns contents of editObj as a double

% --- Executes during object creation, after setting all properties.
function editObj_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editObj (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editSnap_Callback(hObject, eventdata, handles)
% hObject    handle to editSnap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editSnap as text
%        str2double(get(hObject,'String')) returns contents of editSnap as a double

% --- Executes during object creation, after setting all properties.
function editSnap_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editSnap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on button press in pbPlay.
function pbPlay_Callback(hObject, eventdata, handles)
% hObject    handle to pbPlay (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if (get(hObject,'value'))

    set(hObject,'string','Play');
    recorded_data = handles.outdata;
    sz = size(recorded_data.image);
    frm_num = sz(4);
    cntstart = round(get(handles.sldFrame,'value'));

    for cnt = cntstart:frm_num
        axes(handles.axes1);
        subimage(recorded_data.image(:,:,:,cnt));
        axes(handles.axes2);
        St = edge(rgb2gray(recorded_data.image(:,:,:,cnt)));
        subimage(St);
        set(handles.sldFrame,'value',cnt);
        set(handles.editTime2,'string',recorded_data.time{cnt});
        set(handles.txtCurrent,'string',cnt);

        if ~(get(hObject,'value'))
            break
        end
        pause(0.01);
    end
else
    set(hObject,'string','Pause');
end
set(hObject,'value',0);
set(hObject,'string','Pause');

function editTime_Callback(hObject, eventdata, handles)
% hObject    handle to editTime (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editTime as text
%        str2double(get(hObject,'String')) returns contents of editTime as a double

% --- Executes during object creation, after setting all properties.
function editTime_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editTime (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
function sldFrame_Callback(hObject, eventdata, handles)
% hObject    handle to sldFrame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
%

recorded_data = handles.outdata;
cnt = round(get(handles.sldFrame,'value'));
set(handles.sldFrame,'value',cnt)

axes(handles.axes1);
subimage(recorded_data.image(:,:,:,cnt));
axes(handles.axes2);
St = edge(rgb2gray(recorded_data.image(:,:,:,cnt)));
subimage(St);
set(handles.editTime2,'string',recorded_data.time{cnt});
set(handles.txtCurrent,'string',cnt);

% --- Executes during object creation, after setting all properties.
function sldFrame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sldFrame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editTime2_Callback(hObject, eventdata, handles)
% hObject    handle to editTime2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editTime2 as text
%        str2double(get(hObject,'String')) returns contents of editTime2 as a double

% --- Executes during object creation, after setting all properties.
function editTime2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editTime2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function editSnapCnt_Callback(hObject, eventdata, handles)
% hObject    handle to editSnapCnt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editSnapCnt as text
%        str2double(get(hObject,'String')) returns contents of editSnapCnt as a double

% --- Executes during object creation, after setting all properties.
function editSnapCnt_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editSnapCnt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function Monitor_Mode(hObject, eventdata, handles)
set(handles.pbStart,'enable','on');
set(handles.pbStop,'enable','off');
set(handles.pbPlay,'enable','off');
set(handles.sldFrame,'enable','off');
set(handles.txtMin,'enable','off');
set(handles.txtMax,'enable','off');
set(handles.txtCurrent,'enable','off');

function PlayBack_Mode(hObject, eventdata, handles)
set(handles.pbStart,'enable','off');
set(handles.pbStop,'enable','off');
set(handles.pbPlay,'enable','on');
set(handles.sldFrame,'enable','on');
set(handles.txtMin,'enable','on');
set(handles.txtMax,'enable','on');
set(handles.txtCurrent,'enable','on');

% % --------------------------------------------------------------------
% function modeSelect_SelectionChangeFcn(hObject, eventdata, handles)
% % hObject    handle to modeSelect (see GCBO)
% % eventdata  reserved - to be defined in a future version of MATLAB
% % handles    structure with handles and user data (see GUIDATA)
% if get(handles.radioMonitor,'Value')==1
%     Monitor_Mode(hObject, eventdata, handles);
% else
%     PlayBack_Mode(hObject, eventdata, handles);
% end

% --- Executes on button press in radioMonitor.
function radioMonitor_Callback(hObject, eventdata, handles)
% hObject    handle to radioMonitor (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radioMonitor

selected_mode = 1;
set(hObject,'value',selected_mode);
set(handles.radioPlayBack,'value',~selected_mode);
Monitor_Mode(hObject, eventdata, handles);

% --- Executes on button press in radioPlayBack.
function radioPlayBack_Callback(hObject, eventdata, handles)
% hObject    handle to radioPlayBack (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radioPlayBack
selected_mode = 1;
set(hObject,'value',selected_mode);
set(handles.radioMonitor,'value',~selected_mode);

recorded_data = handles.outdata;
sz = size(recorded_data.image);
frm_num = sz(4);

set(handles.sldFrame,'min',1);
set(handles.sldFrame,'max',frm_num);
set(handles.sldFrame,'SliderStep',[1/frm_num 5/frm_num]);
set(handles.txtMin,'string',1);
set(handles.txtMax,'string',frm_num);

set(handles.sldFrame,'value',1);
PlayBack_Mode(hObject, eventdata, handles);

author: Sable

repo: mcbench-benchmarks


text:
function unit = ft_estimate_units(size)

% FT_ESTIMATE_UNITS tries to determine the units of a geometrical object by
% looking at its size and by relating this to the size of the human
% brain.
%
% Use as
%   unit = ft_estimate_units(size)
%
% This function will return one of the following strings
%   'm'
%   'dm'
%   'cm'
%   'mm'

% Copyright (C) 2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: ft_estimate_units.m 2885 2011-02-16 09:41:58Z roboos $

% do some magic based on the size
unit = {'m', 'dm', 'cm', 'mm'};
indx = round(log10(size)+2-0.2);

if indx>length(unit)
  indx = length(unit);
  warning('assuming that the units are "%s"', unit{indx});
end

if indx<1
  indx = 1;
  warning('assuming that the units are "%s"', unit{indx});
end

unit = unit{indx};

author: PatternRecognition

repo: OpenBMI


text:
function data = CreateEditData(fieldCells, session)
%
%  data = CreateEditData(fieldCells, session)
%
% Compare the field names in the session structure with
% those specified by fieldCells. If they match, add
% field name, label, contents, edit flat, and length to the
% output data struct array.
%
% See GetReconEdit.m for more information.
%
% DBR, 4/99

sessionFields = fieldnames(session);
iMatch = 0;

for iField=1:size(fieldCells, 1)
  field = fieldCells{iField, 1};

  if strmatch(field, sessionFields, 'exact')
    iMatch = iMatch + 1;
    data(iMatch).field = field;
    label = fieldCells{iField, 2};
    data(iMatch).label = label;
    fieldName = fieldCells{iField, 1};
    content = num2str(session.(fieldName));
    data(iMatch).content = content;
    data(iMatch).edit = fieldCells{iField, 3};
    data(iMatch).width = length(label) + length(content);
  end

end

return;


text:
function [R,current_eps]=MC_USTM_Alamouti(snrdB,T,L,epsilon,prec,filename)
%
% metaconverse bound for the Alamouti ensemble
%-------------------------------------------------------------------
%                       SET-UP PARAMETERS
%-------------------------------------------------------------------

SAVE=1;
MAT=1;

K = 2^prec; % number of monte carlo simulations (at least 100 x 1/epsilon)
rho = 10.^(snrdB/10); % SNR in linear scale
Mt=2;
Mr=2;
%-------------------------------------------------------------------
%                       MONTE CARLO SIMULATION
%-------------------------------------------------------------------
Ip = zeros(K,1); %allocate for the montecarlo runs
% Iq = zeros(K,1); %allocate for the montecarlo runs
%-------------------------------------------------------------------
%                       CONSTANTS
%-------------------------------------------------------------------
rho_tilde = T*rho/Mt;

lambda=1+rho_tilde;
lambda1=1/lambda;
lambda2=rho_tilde*lambda1;

x1=rho_tilde;
x2=rho_tilde;

D= [diag([sqrt(1+x1), sqrt(1+x2)]), zeros(Mt,T-Mt);
zeros(T-Mt,Mt), eye(T-Mt)]; % D matrix (covariance matrix of equivalent noise)

%-------------------------------------------------------------------
%                       MONTE CARLO
%-------------------------------------------------------------------

%tic

norm=sqrt(.5);

for k=1:K

        i_L = 0;

        Z = randn(T,Mr,L)*norm+1i*randn(T,Mr,L)*norm;

        for l = 1:L %Create each realization

          %COMPUTE EVERYTHING THAT HAS TO DO WITH SINGULAR VALUES
          Sigma_alt= svd(D*Z(:,:,l)).^2;
          Sigma_alt=sort(Sigma_alt,1,'descend');
          TraceZ=abs(trace(Z(:,:,l)'*Z(:,:,l)));

          Y=D*Z(:,:,l);

          Ytilde=(zeros(T,4));

          Ytilde(:,[1,3])=Y;

          Ytilde(1:2:T,[2,4])= conj(Y(2:2:T,:));

          Ytilde(2:2:T,[2,4])= -conj(Y(1:2:T,:));

          Sigma=svd(Ytilde).^2;

          Sigma=[Sigma(1),Sigma(3)];

          Sigma=Sigma*lambda2;

           if (T>4),

               M=[gammainc(Sigma(1), T-5),  gammainc(Sigma(1), T-4), exp(Sigma(2)-Sigma(1))*gammainc(Sigma(2), T-5)/(Sigma(2)/Sigma(1))^(T-4), exp(Sigma(2)-Sigma(1))*gammainc(Sigma(2), T-4)/(Sigma(2)/Sigma(1))^(T-4) ;
                   (T-2)*Sigma(1), Sigma(1)^2,(T-2)*Sigma(2), Sigma(2)^2 ;
                   T-3, Sigma(1),T-3, Sigma(2);
                  (T-4)/Sigma(1),1,(T-4)/Sigma(2),1];

              logd=log(det(M))-(T-4)*log(Sigma(1))+Sigma(1);

            else

              M=[ 1, 2*Sigma(1), 1, 0 ; ...
                1,  Sigma(1)^2, Sigma(1), 1; ...
                exp(Sigma(2)-Sigma(1)), 2*Sigma(2), 1, 0; ...
                exp(Sigma(2)-Sigma(1)), Sigma(2)^2, Sigma(2), 1];

              logd=(log(det(M)))+Sigma(1);

            end

            log_exp_sum = logd -4*log(Sigma(1)-Sigma(2));

            i = - TraceZ  +sum(Sigma_alt) - log(gamma(T)) - log_exp_sum;

            i_L = i_L + i; %add it to the total i_L

        end

        Ip(k) =  i_L; %put all computations on a pile to compute the average later
end

if (SAVE==1)
  if (MAT==1)
    save(filename,'Ip')
  else
    save(filename,'Ip','-ascii','-append')
  end
end

%---------------------------------------
%   SEARCHING THE RATE
%---------------------------------------

% load saved data (to account for append possibilities)
if (SAVE==1 && MAT==0)

    Ip=load(filename);
end

Ip=sort(Ip);

Kcurrent=length(Ip); % redefine K to account for append

current_prec=floor(log2(Kcurrent)); % actual precision

K=2^(current_prec); % round off K to avoid search errors

% first find a suitable initial point for the linear search
step=K/2;
index=step;

onevec=ones(K,1);

while(step>1),

   th=Ip(index);

   current_eps=sum(Ip<=th)/K;

   step=step/2;

   if current_eps> epsilon,

       index=index-step;

   else

       index=index+step;

   end

end

current_eps=sum(Ip<=Ip(index))/K;

if(current_eps<epsilon)
    index=index+1;
end

% now perform linear search
% Rvect=zeros(1,K-index+1);
% for ii=1:K-index+1,
% %for ii=1:K-index+1,
%
%     current_rate=Ip(ii+index-1)-log(sum(Ip<=Ip(ii+index-1))/K-epsilon);
%
%     Rvect(ii)=current_rate;
%
% end

%faster search (potentially less tight)

count=0;

R=Ip(index)-log(sum(Ip<=Ip(index))/K-epsilon);

for ii=1:K-index+1,
    current_rate=Ip(ii+index-1)-log(sum(Ip<=Ip(ii+index-1))/K-epsilon);
    if current_rate<= R
        R=current_rate;
     else
       count=count+1;
    end

    if count==20,
      break
    end
end
R=R/(L*T*log(2));

end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function value = r8_cbrt ( x )

%*****************************************************************************80
%
%% R8_CBRT computes the cube root of an R8.
%
%  Discussion:
%
%    The approximation is a generalized Chebyshev series converted
%    to polynomial form.  The approximation is nearly best in the
%    sense of relative error with 4.085 digits accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the number whose square root is desired.
%
%    Output, real VALUE, the cube root of X.
%
  persistent cbrt2
  persistent niter

  if ( isempty ( niter ) )

    cbrt2 = [ ...
      0.62996052494743658238360530363911, ...
      0.79370052598409973737585281963615, ...
      1.0, ...
      1.25992104989487316476721060727823, ...
      1.58740105196819947475170563927231 ]';

    niter = r8_aint ( 1.443 * log ( - 0.106 ...
      * log ( 0.1 * r8_mach ( 3 ) ) ) + 1.0 );

  end

  value = 0.0;

  if ( x ~= 0.0 )

    [ y, n ] = r8_upak ( abs ( x ) );
    ixpnt = floor ( n / 3 );
    irem = n - 3 * ixpnt + 3;

    value = 0.439581 + y * ( ...
            0.928549 + y * ( ...
          - 0.512653 + y * ...
            0.144586 ) );

    for iter = 1 : niter
      vsq = value * value;
      value = value + ( y - value * vsq ) / ( 3.0 * vsq );
    end

    if ( x < 0.0 )
      value = - abs ( value );
    else
      value = + abs ( value );
    end

    value = r8_pak ( cbrt2(irem) * value, ixpnt );

  end

  return
end