text:
function [radius,subsettrunc,outstate] = ncorr_gui_setstrainradius(reference,current,rois_ref,rois_cur,plots_u_ref,plots_v_ref,plots_u_cur,plots_v_cur,spacing,pixtounits,units,pos_parent,params_init)
% This is a GUI for setting the strain radius.
%
% Inputs -----------------------------------------------------------------%
%   reference - ncorr_class_img; used for displaying the background image.
%   current - ncorr_class_img; used for displaying the background image.
%   rois_ref - ncorr_class_roi; contains ROIs WRT the reference image. Note
%   that rois_ref has already been reduced.
%   rois_cur - ncorr_class_roi; contains ROIs WRT the current image. Note
%   that rois_ref has already been reduced.
%   plots_u_ref - cell; contains displacement data WRT the reference image.
%   Uses real units. Note that displacement plots are reduced by default.
%   plots_v_ref - cell; contains displacement data WRT the reference image.
%   Uses real units. Note that displacement plots are reduced bym default.
%   plots_u_cur - cell; contains displacement data WRT the current image.
%   Uses real units. Note that displacement plots are reduced by default.
%   plots_v_cur - cell; contains displacement data WRT the current image.
%   Uses real units. Note that displacement plots are reduced by default.
%   spacing - double; spacing parameter
%   pixtounits - double; units/pixels conversion factor
%   units - string; this a string containing the units
%   pos_parent - integer array; this is the position of the parent figure
%   which determines where to position this figure
%   params_init - integer; contains [radius subsettrunc] if they have been set before.
%
% Outputs ----------------------------------------------------------------%
%   radius - integer; strain radius
%   subsettrunc - logical; true if subset truncation is set
%   outstate - integer; returns either out.cancelled, out.failed, or
%   out.success.
%
% Note that some regions for the input ROIs can be empty. This is handled
% properly.

    % Data ---------------------------------------------------------------%
    % Initialize outputs
    outstate = out.cancelled;
    radius = [];
    subsettrunc = [];
    % Get GUI handles
    handles_gui = init_gui();
    % Run c-tor
    feval(ncorr_util_wrapcallbacktrycatch(@constructor,handles_gui.figure));

    % Callbacks and functions --------------------------------------------%
    function constructor()
        % Set zoom and pan
        handle_zoom = zoom(handles_gui.figure);
        handle_pan = pan(handles_gui.figure);

        % Allow rotations - For some reason when rotations are enabled, the
        % impoint is no longer draggable, so only allow rotations when the
        % cursor is over the least squares plane fit axes.
        handle_rotate3d = rotate3d(handles_gui.axes_leastsquares);
        set(handles_gui.figure,'WindowButtonMotionFcn',ncorr_util_wrapcallbacktrycatch(@callback_moveplot,handles_gui.figure));

        % Set buffers
        % Min radius and spacing
        min_radius = 1;
        max_radius = 100;

        % Create radius function this will take the input from the slider
        % and return a radius
        get_radius = @(x)round(x*50+min_radius);

        % Initialize strain radius and slider
        if (isempty(params_init))
            radius_prelim = 15;
            radius_slider_prelim = min(1,max(0,fzero(@(x) get_radius(x)-radius_prelim,0)));
            subsettrunc_prelim = false;
        else
            radius_prelim = params_init(1);
            radius_slider_prelim = min(1,max(0,fzero(@(x) get_radius(x)-radius_prelim,0)));
            subsettrunc_prelim = params_init(2);
        end

        % Set data
        setappdata(handles_gui.figure,'min_radius',min_radius);
        setappdata(handles_gui.figure,'max_radius',max_radius);
        setappdata(handles_gui.figure,'radius_prelim',radius_prelim);
        setappdata(handles_gui.figure,'radius_slider_prelim',radius_slider_prelim);
        setappdata(handles_gui.figure,'subsettrunc_prelim',subsettrunc_prelim);
        setappdata(handles_gui.figure,'get_radius',get_radius);
        setappdata(handles_gui.figure,'num_cur',length(current)-1);
        setappdata(handles_gui.figure,'handle_points',[]);
        setappdata(handles_gui.figure,'handle_preview',[]);
        setappdata(handles_gui.figure,'handle_rotate3d',handle_rotate3d);
        setappdata(handles_gui.figure,'val_popupmenu_disp',1);
        setappdata(handles_gui.figure,'val_popupmenu_lore',1);
        setappdata(handles_gui.figure,'pos_points',[]);
        setappdata(handles_gui.figure,'num_region',0);
        setappdata(handles_gui.figure,'handle_zoom',handle_zoom);
        setappdata(handles_gui.figure,'handle_pan',handle_pan);

        % Update
        update_axes('set');
        update_sidemenu();

        % Set Visible
        set(handles_gui.figure,'Visible','on');
    end

    function callback_popupmenu_lore(hObject,eventdata) %#ok<INUSD>
        % Get data
        val_popupmenu_lore = get(handles_gui.popupmenu_lore,'value');

        % Clear position of points
        pos_points = [];

        % Set data
        setappdata(handles_gui.figure,'pos_points',pos_points);
        setappdata(handles_gui.figure,'val_popupmenu_lore',val_popupmenu_lore);

        % Update
        update_axes('set');
        update_sidemenu();
    end

    function callback_popupmenu_disp(hObject,eventdata) %#ok<INUSD>
        % Get data
        val_popupmenu_disp = get(handles_gui.popupmenu_disp,'value');

        % Set data
        setappdata(handles_gui.figure,'val_popupmenu_disp',val_popupmenu_disp);

        % Update
        update_axes('set');
        update_sidemenu();
    end

    function callback_slider_radius(hObject,eventdata) %#ok<INUSD>
        % Get data
        get_radius = getappdata(handles_gui.figure,'get_radius');

        % Get slider value
        radius_slider_prelim = get(handles_gui.slider_radius,'value');

        % Get radius_prelim
        radius_prelim = get_radius(radius_slider_prelim);

        % Set data
        setappdata(handles_gui.figure,'radius_prelim',radius_prelim);
        setappdata(handles_gui.figure,'radius_slider_prelim',radius_slider_prelim);

        % Update
        update_axes('update');
        update_sidemenu();
    end

    function callback_checkbox_subsettrunc(hObject,eventdata) %#ok<INUSD>
        % Get value
        subsettrunc_prelim = get(handles_gui.checkbox_subsettrunc,'Value');

        % Set data
        setappdata(handles_gui.figure,'subsettrunc_prelim',subsettrunc_prelim);

        % Update
        update_axes('update');
        update_sidemenu();
    end

    function callback_edit_radius(hObject,eventdata) %#ok<INUSD>
        min_radius = getappdata(handles_gui.figure,'min_radius');
        max_radius = getappdata(handles_gui.figure,'max_radius');
        radius_prelim = getappdata(handles_gui.figure,'radius_prelim');
        radius_slider_prelim = getappdata(handles_gui.figure,'radius_slider_prelim');
        get_radius = getappdata(handles_gui.figure,'get_radius');

        % Get radius
        radius_buffer = str2double(get(handles_gui.edit_radius,'string'));
        if (ncorr_util_isintbb(radius_buffer,min_radius,max_radius,'Radius') == out.success)
            % Set radius
            radius_prelim = radius_buffer;
            radius_slider_prelim = min(1,max(0,fzero(@(x)get_radius(x)-radius_prelim,0)));
        end

        % Set data
        setappdata(handles_gui.figure,'radius_prelim',radius_prelim);
        setappdata(handles_gui.figure,'radius_slider_prelim',radius_slider_prelim);

        % Update
        update_axes('update');
        update_sidemenu();
    end

    function callback_edit_imgnum(hObject,eventdata) %#ok<INUSD>
        % Get Data
        num_cur = getappdata(handles_gui.figure,'num_cur');
        pos_points = getappdata(handles_gui.figure,'pos_points');

        % Get img num
        num_cur_prelim = str2double(get(handles_gui.edit_imgnum,'string'));
        if (ncorr_util_isintbb(num_cur_prelim,1,length(current),'Current image number') == out.success)
            % Clear position of points
            pos_points = [];

            % num_cur_prelim is zero based indexed
            num_cur = num_cur_prelim-1;
        end

        % Set Data
        setappdata(handles_gui.figure,'num_cur',num_cur);
        setappdata(handles_gui.figure,'pos_points',pos_points);

        % Update
        update_axes('set');
    end

    function callback_button_zoom(hObject,eventdata) %#ok<INUSD>
        % Get data
        handle_zoom = getappdata(handles_gui.figure,'handle_zoom');

        if (strcmp(get(handle_zoom,'Enable'),'on'))
            % Zoom is already enabled; disable it
            set(handle_zoom,'Enable','off');
        else
            % Zoom not enabled; enable it
            set(handle_zoom,'Enable','on');
        end

        % Set data
        setappdata(handles_gui.figure,'handle_zoom',handle_zoom);

        % Update
        update_sidemenu();
    end

    function callback_button_pan(hObject,eventdata) %#ok<INUSD>
        % Get data
        handle_pan = getappdata(handles_gui.figure,'handle_pan');

        if (strcmp(get(handle_pan,'Enable'),'on'))
            % Pan is already enabled; disable it
            set(handle_pan,'Enable','off');
        else
            % Pan not enabled; enable it
            set(handle_pan,'Enable','on');
        end

        % Set data
        setappdata(handles_gui.figure,'handle_pan',handle_pan);

        % Update
        update_sidemenu();
    end

    function callback_button_finish(hObject,eventdata) %#ok<INUSD>
        % Get data
        radius_prelim = getappdata(handles_gui.figure,'radius_prelim');
        subsettrunc_prelim = getappdata(handles_gui.figure,'subsettrunc_prelim');

        if (subsettrunc_prelim)
            text_st = 'Enabled';
        else
            text_st = 'Disabled';
        end
        contbutton = questdlg({['Strain radius is set to: ' num2str(radius_prelim) '.'], ...
                               ['Subset Truncation is: ' text_st '.'], ...
                               'Is this correct?'}, ...
                               'Continue Operation','Yes','No','No');
        if (strcmp(contbutton,'Yes'))
            % Set Outputs
            radius = radius_prelim;
            subsettrunc = subsettrunc_prelim;
            outstate = out.success;

            % Exit
            close(handles_gui.figure);
        end
    end

    function callback_button_cancel(hObject,eventdata) %#ok<INUSD>
        close(handles_gui.figure);
    end

    function callback_button_left(hObject,eventdata) %#ok<INUSD>
        % Get Data
        num_cur = getappdata(handles_gui.figure,'num_cur');
        pos_points = getappdata(handles_gui.figure,'pos_points');

        % Check for overshoot
        if (num_cur > 0)
            % Clear position of points
            pos_points = [];

            num_cur = num_cur-1;
        end

        % Set Data
        setappdata(handles_gui.figure,'num_cur',num_cur);
        setappdata(handles_gui.figure,'pos_points',pos_points);

        % Update
        update_axes('set');
        update_sidemenu();
    end

    function callback_button_right(hObject,eventdata) %#ok<INUSD>
        % Get Data
        num_cur = getappdata(handles_gui.figure,'num_cur');
        pos_points = getappdata(handles_gui.figure,'pos_points');

        % Check for overshoot
        if (num_cur < length(current)-1)
            % Clear position of points
            pos_points = [];

            num_cur = num_cur+1;
        end

        % Set Data
        setappdata(handles_gui.figure,'num_cur',num_cur);
        setappdata(handles_gui.figure,'pos_points',pos_points);

        % Update
        update_axes('set');
        update_sidemenu();
    end

    function callback_impoint(pos,num_region)
        % Get data
        pos_points = getappdata(handles_gui.figure,'pos_points');

        % Store position point
        pos_points(num_region+1,:) = round(pos)-1;

        % Set Data
        setappdata(handles_gui.figure,'pos_points',pos_points);
        setappdata(handles_gui.figure,'num_region',num_region);

        % Update
        update_axes('update');
    end

    function callback_moveplot(hObject,eventdata) %#ok<INUSD>
        % Get data
        handle_rotate3d = getappdata(handles_gui.figure,'handle_rotate3d');
        handle_zoom = getappdata(handles_gui.figure,'handle_zoom');
        handle_pan = getappdata(handles_gui.figure,'handle_pan');

        % Enable rotation if zoom and pan are disabled
        pos_mouse = get(handles_gui.figure,'CurrentPoint');
        pos_group = get(handles_gui.group_preview,'Position');
        pos_axes = get(handles_gui.axes_leastsquares,'Position');
        inleastsquares = (pos_mouse(1) > pos_group(1)+pos_axes(1) && pos_mouse(1) < pos_group(1)+pos_axes(1)+pos_axes(3) && pos_mouse(2) > pos_group(2)+pos_axes(2) && pos_mouse(2) < pos_group(2)+pos_axes(2)+pos_axes(4));
        if (inleastsquares)
            set(handle_rotate3d,'Enable','on');
        end

        % If pan and zoom are both disabled and cursor is outside of least
        % square axes, then disable rotation
        if (strcmp(get(handle_zoom,'Enable'),'off') && strcmp(get(handle_pan,'Enable'),'off') && ~inleastsquares)
            set(handle_rotate3d,'Enable','off');
        end

        % Update
        update_sidemenu();
    end

    function update_sidemenu()
        % Get data
        radius_prelim = getappdata(handles_gui.figure,'radius_prelim');
        radius_slider_prelim = getappdata(handles_gui.figure,'radius_slider_prelim');
        val_popupmenu_disp = getappdata(handles_gui.figure,'val_popupmenu_disp');
        val_popupmenu_lore = getappdata(handles_gui.figure,'val_popupmenu_lore');
        handle_zoom = getappdata(handles_gui.figure,'handle_zoom');
        handle_pan = getappdata(handles_gui.figure,'handle_pan');
        subsettrunc_prelim = getappdata(handles_gui.figure,'subsettrunc_prelim');

        % Update slider
        set(handles_gui.slider_radius,'value',radius_slider_prelim);

        % Update edit field
        set(handles_gui.edit_radius,'String',num2str(radius_prelim))

        % Update popupmenus
        set(handles_gui.popupmenu_disp,'value',val_popupmenu_disp);
        set(handles_gui.popupmenu_lore,'value',val_popupmenu_lore);

        if (strcmp(get(handle_pan,'Enable'),'on'))
            set(handles_gui.button_pan,'FontWeight','bold');
        else
            set(handles_gui.button_pan,'FontWeight','normal');
        end

        if (strcmp(get(handle_zoom,'Enable'),'on'))
            set(handles_gui.button_zoom,'FontWeight','bold');
        else
            set(handles_gui.button_zoom,'FontWeight','normal');
        end

        set(handles_gui.checkbox_subsettrunc,'Value',subsettrunc_prelim);
    end

    function update_axes(action)
        % Get data
        radius_prelim = getappdata(handles_gui.figure,'radius_prelim');
        num_cur = getappdata(handles_gui.figure,'num_cur');
        handle_points = getappdata(handles_gui.figure,'handle_points');
        handle_preview = getappdata(handles_gui.figure,'handle_preview');
        val_popupmenu_disp = getappdata(handles_gui.figure,'val_popupmenu_disp');
        val_popupmenu_lore = getappdata(handles_gui.figure,'val_popupmenu_lore');
        pos_points = getappdata(handles_gui.figure,'pos_points');
        num_region = getappdata(handles_gui.figure,'num_region');
        subsettrunc_prelim = getappdata(handles_gui.figure,'subsettrunc_prelim');

        % Get img, data plot, and ROI
        if (val_popupmenu_disp == 1 && val_popupmenu_lore == 1)
            img_bg = reference;
            plot_data = plots_u_ref{num_cur+1};
            roi_data = rois_ref(num_cur+1);
        elseif (val_popupmenu_disp == 2 && val_popupmenu_lore == 1)
            img_bg = reference;
            plot_data = plots_v_ref{num_cur+1};
            roi_data = rois_ref(num_cur+1);
        elseif (val_popupmenu_disp == 1 && val_popupmenu_lore == 2)
            img_bg = current(num_cur+1);
            plot_data = plots_u_cur{num_cur+1};
            roi_data = rois_cur(num_cur+1);
        else
            img_bg = current(num_cur+1);
            plot_data = plots_v_cur{num_cur+1};
            roi_data = rois_cur(num_cur+1);
        end

        if (strcmp(action,'set'))
            % Get reduced img
            img_reduced = img_bg.reduce(spacing);

            % Initialize displacements
            imshow(img_reduced.get_img(),[img_reduced.min_gs img_reduced.max_gs],'Parent',handles_gui.axes_preview);
            hold(handles_gui.axes_preview,'on');

            % Get bounds
            cmin = prctile(plot_data(roi_data.mask),5)-1e-5;
            cmax = prctile(plot_data(roi_data.mask),95)+1e-5;
            if (abs(cmin-cmax)<1e-4)
                cmax = cmax + 1e-5; % add a small increment to insure they are not the same number
                cmin = cmin - 1e-5;
            end

            handle_preview = imshow(plot_data,'Parent',handles_gui.axes_preview);
            set(handles_gui.axes_preview,'CLim',[cmin cmax]);
            hold(handles_gui.axes_preview,'off');

            % Set colorbar
            handle_colorbar = colorbar('peer',handles_gui.axes_preview);
            set(handle_colorbar,'UIContextMenu','');
            set(get(handle_colorbar,'child'),'YData',[cmin cmax]);
            set(handle_colorbar,'YLim',[cmin cmax]);
            set(handle_colorbar,'Units','Pixels');

            % Set Impoints - only set impoints in regions which are full - use
            % a cell array here since it is not gauranteed that regions are
            % full.
            if (isempty(pos_points))
                setnewpos = true;
            else
                setnewpos = false;
            end

            handle_points = cell(1,length(roi_data.region));
            for i = 0:length(roi_data.region)-1
                if (roi_data.region(i+1).totalpoints > 0)
                    % Form point
                    handle_points{i+1} = impoint(handles_gui.axes_preview,round(size(roi_data.mask,2)/2),round(size(roi_data.mask,1)/2));
                    setColor(handle_points{i+1},'g');
                    set(handle_points{i+1},'UIContextMenu','');

                    % Set constraint function
                    api_point = iptgetapi(handle_points{i+1});
                    constrainfcn = ncorr_util_formregionconstraint(roi_data.region(i+1));
                    api_point.setPositionConstraintFcn(constrainfcn);

                    % Initialize point location - use 1 based indexing for
                    % input to constrainfcn
                    x_center = round(roi_data.region(i+1).leftbound+(roi_data.region(i+1).rightbound-roi_data.region(i+1).leftbound)/2)+1;
                    y_center = round(roi_data.region(i+1).upperbound+(roi_data.region(i+1).lowerbound-roi_data.region(i+1).upperbound)/2)+1;

                    % See if points already exist
                    if (setnewpos)
                        % constrainfcn returns points based on 1 based
                        % indexing, so convert it to zero based.
                        pos_points(i+1,1:2) = constrainfcn([x_center y_center])-1;
                    end
                    % Set position based on 1 based indexing
                    setPosition(handle_points{i+1},pos_points(i+1,1:2)+1);

                    % Assign Position callback
                    addNewPositionCallback(handle_points{i+1},ncorr_util_wrapcallbacktrycatch(@(pos)callback_impoint(pos,i),handles_gui.figure));

                    % Set num_region
                    if (setnewpos)
                        num_region = i;
                    end
                end
            end

            % Set colormap
            ncorr_util_colormap(handles_gui.figure);

            % Update buttons
            set(handles_gui.edit_imgnum,'String',num2str(num_cur+1));
            if (length(current) == 1)
                set(handles_gui.button_left,'Enable','off');
                set(handles_gui.button_right,'Enable','off');
                set(handles_gui.edit_imgnum,'Enable','off');
            elseif (num_cur == 0)
                set(handles_gui.button_left,'Enable','off');
                set(handles_gui.button_right,'Enable','on');
                set(handles_gui.edit_imgnum,'Enable','on');
            elseif (num_cur == length(current)-1)
                set(handles_gui.button_left,'Enable','on');
                set(handles_gui.button_right,'Enable','off');
                set(handles_gui.edit_imgnum,'Enable','on');
            else
                set(handles_gui.button_left,'Enable','on');
                set(handles_gui.button_right,'Enable','on');
                set(handles_gui.edit_imgnum,'Enable','on');
            end
        end

        if (strcmp(action,'set') || strcmp(action,'update'))
            % Get cirroi
            cirroi = roi_data.get_cirroi(pos_points(num_region+1,1),pos_points(num_region+1,2),radius_prelim,subsettrunc_prelim);

            % Initialize matrix and vector to zeros
            mat_LS = zeros(3);
            z_vec_LS = zeros(3,1);

            % Initialize alphamask
            alpha_preview = 0.5*roi_data.mask;

            % Preallocate arrays
            plotpoints_x = zeros(cirroi.region.totalpoints,1);
            plotpoints_y = zeros(cirroi.region.totalpoints,1);
            plotpoints_z = zeros(cirroi.region.totalpoints,1);

            % Initialize counter
            counter = 0;
            for i = 0:size(cirroi.region.noderange,1)-1
                x = i+(cirroi.x-cirroi.radius);
                for j = 0:2:cirroi.region.noderange(i+1)-1
                    y_vec = cirroi.region.nodelist(i+1,j+1):cirroi.region.nodelist(i+1,j+2);

                    % Form cirroi.x,cirroi.y,z vectors
                    plotpoints_x(counter+1:counter+length(y_vec)) = x*ones(size(y_vec))*(pixtounits*(spacing+1));
                    plotpoints_y(counter+1:counter+length(y_vec)) = y_vec*(pixtounits*(spacing+1));
                    plotpoints_z(counter+1:counter+length(y_vec)) = plot_data(y_vec+1,x+1)';

                    % Do least squares analysis to find plane
                    x_LS = (x*ones(size(y_vec))-cirroi.x)*(pixtounits*(spacing+1));
                    y_LS = (y_vec-cirroi.y)*(pixtounits*(spacing+1));
                    z_LS = plot_data(y_vec+1,x+1)';

                    % Do the matrix first
                    mat_LS(1,1) = mat_LS(1,1)+sum(x_LS.^2);
                    mat_LS(2,1) = mat_LS(2,1)+sum(x_LS.*y_LS);
                    mat_LS(3,1) = mat_LS(3,1)+sum(x_LS);
                    mat_LS(1,2) = mat_LS(1,2)+sum(x_LS.*y_LS);
                    mat_LS(2,2) = mat_LS(2,2)+sum(y_LS.^2);
                    mat_LS(3,2) = mat_LS(3,2)+sum(y_LS);
                    mat_LS(1,3) = mat_LS(1,3)+sum(x_LS);
                    mat_LS(2,3) = mat_LS(2,3)+sum(y_LS);
                    mat_LS(3,3) = mat_LS(3,3)+sum(ones(size(y_vec)));

                    % Do vec next
                    z_vec_LS(1) = z_vec_LS(1)+sum(x_LS.*z_LS);
                    z_vec_LS(2) = z_vec_LS(2)+sum(y_LS.*z_LS);
                    z_vec_LS(3) = z_vec_LS(3)+sum(z_LS);

                    % Update Preview Mask
                    alpha_preview(y_vec+1,x+1) = 1;

                    % Update counter
                    counter = counter + length(y_vec);
                end
            end

            % Set Alpha Mask
            set(handle_preview,'AlphaData',alpha_preview);

            if (abs(det(mat_LS)) > 10^-10) % Make sure matrix is not singular
                % Find dp/dx dp/cirroi.y
                du_vec = linsolve(mat_LS,z_vec_LS);

                % Plot points first
                plot3(plotpoints_x,plotpoints_y,plotpoints_z,'o','Parent',handles_gui.axes_leastsquares);
                zlabel(handles_gui.axes_leastsquares,['Displacement - ' units]);
                hold(handles_gui.axes_leastsquares,'on');

                % Plot the plane next
                [plane_x,plane_y] = meshgrid((-cirroi.radius:cirroi.radius)*(pixtounits*(spacing+1)),(-cirroi.radius:cirroi.radius)*(pixtounits*(spacing+1)));
                plane_z = du_vec(3)+du_vec(2)*plane_y+du_vec(1)*plane_x;
                mesh(plane_x+cirroi.x*(pixtounits*(spacing+1)),plane_y+cirroi.y*(pixtounits*(spacing+1)),plane_z,'Parent',handles_gui.axes_leastsquares);

                hold(handles_gui.axes_leastsquares,'off');
            else
                cla(handles_gui.axes_leastsquares);
            end

            % Set Name
            set(handles_gui.text_ref,'String',['Reference Name: ' reference.name(1:min(end,22))]);
            set(handles_gui.text_cur,'String',['Current Name: ' current(num_cur+1).name(1:min(end,22))]);
        end

        % Set data
        setappdata(handles_gui.figure,'handle_points',handle_points);
        setappdata(handles_gui.figure,'handle_preview',handle_preview);
        setappdata(handles_gui.figure,'pos_points',pos_points);
        setappdata(handles_gui.figure,'num_region',num_region);
    end

    function handles_gui = init_gui()
    % GUI controls -------------------------------------------------------%
        % Figure
        handles_gui.figure = figure( ...
            'Tag', 'figure', ...
            'Units', 'characters', ...
            'Position', ncorr_util_figpos(pos_parent,[32.5 167]), ...
            'Name', 'Set Strain Parameters', ...
            'MenuBar', 'none', ...
            'NumberTitle', 'off', ...
            'Color', get(0,'DefaultUicontrolBackgroundColor'), ...
            'handlevisibility','off', ...
            'DockControls','off', ...
            'Resize','off', ...
            'WindowStyle','modal', ...
            'Visible','off', ...
            'IntegerHandle','off', ...
            'Interruptible','off');

        % Panels
        handles_gui.group_strainoptions = uipanel( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_strainoptions', ...
            'Units', 'characters', ...
            'Position', [2 25.5 35 6.4], ...
            'Title', 'Strain Options', ...
            'Interruptible','off');

        handles_gui.group_viewoptions = uipanel( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_viewoptions', ...
            'Units', 'characters', ...
            'Position', [2 17.8 35 7.2], ...
            'Title', 'View Options', ...
            'Interruptible','off');

        handles_gui.group_subsettrunc = uipanel( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_subsettrunc', ...
            'Units', 'characters', ...
            'Position', [2 12.8 35 4.3], ...
            'Title', 'Discontinuous Analysis', ...
            'Interruptible','off');

        handles_gui.group_zoompan = uibuttongroup( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_zoompan', ...
            'Units', 'characters', ...
            'Position', [2 8.0 35 4.5], ...
            'Title', 'Zoom/Pan', ...
            'Interruptible','off');

        handles_gui.group_menu = uipanel( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_menu', ...
            'Units', 'characters', ...
            'Position', [2 1 35 6.7], ...
            'Title', 'Menu', ...
            'Interruptible','off');

        handles_gui.group_preview = uipanel( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_preview', ...
            'Units', 'characters', ...
            'Position', [39 0.75 126 31.2], ...
            'Title', 'Preview', ...
            'Interruptible','off');

        % Axes
        handles_gui.axes_preview = axes( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'axes_preview', ...
            'Units', 'characters', ...
            'Visible','off', ...
            'Position', [2 4 57 23.5], ...
            'Interruptible','off');

        handles_gui.axes_leastsquares = axes( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'axes_leastsquares', ...
            'Units', 'characters', ...
            'Visible', 'off', ...
            'Position', [74.8 6 38 19.5], ...
            'Interruptible','off');

        % Drop-down Menu
        handles_gui.popupmenu_lore = uicontrol( ...
            'Parent', handles_gui.group_viewoptions, ...
            'Tag', 'popupmenu_lore', ...
            'Style', 'popupmenu', ...
            'Units', 'characters', ...
            'Position', [2.4 3.6 29.1 1.3], ...
            'BackgroundColor', [1 1 1], ...
            'String', {'Lagrangian','Eulerian'}, ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_popupmenu_lore,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.popupmenu_disp = uicontrol( ...
            'Parent', handles_gui.group_viewoptions, ...
            'Tag', 'popupmenu_disp', ...
            'Style', 'popupmenu', ...
            'Units', 'characters', ...
            'Position', [2.4 1.5 29.1 1.3], ...
            'BackgroundColor', [1 1 1], ...
            'String', {'U-Displacement','V-Displacement'}, ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_popupmenu_disp,handles_gui.figure), ...
            'Interruptible','off');

        % Static Texts
        handles_gui.text_radius = uicontrol( ...
            'Parent', handles_gui.group_strainoptions, ...
            'Tag', 'text_radius', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.5 2.6 16 1.3], ...
            'String', 'Strain Radius: ', ...
            'HorizontalAlignment', 'left', ...
            'Interruptible','off');

        handles_gui.text_ref = uicontrol( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'text_ref', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [3.0 1.8 56.9 1.3], ...
            'String', 'Reference Name: ', ...
            'HorizontalAlignment', 'left', ...
            'Interruptible','off');

        handles_gui.text_cur = uicontrol( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'text_cur', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [3.0 0.6 56.9 1.3], ...
            'String', 'Current Name: ', ...
            'HorizontalAlignment', 'left', ...
            'Interruptible','off');

        % Sliders
        handles_gui.slider_radius = uicontrol( ...
            'Parent', handles_gui.group_strainoptions, ...
            'Tag', 'slider_radius', ...
            'Style', 'slider', ...
            'Units', 'characters', ...
            'Position', [2.5 1.2 28.9 1.3], ...
            'BackgroundColor', [0.9 0.9 0.9], ...
            'String', {'Slider'}, ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_slider_radius,handles_gui.figure), ...
            'Interruptible','off');

        % Checkbox
        handles_gui.checkbox_subsettrunc = uicontrol( ...
            'Parent', handles_gui.group_subsettrunc, ...
            'Tag', 'checkbox_subsettrunc', ...
            'Style', 'checkbox', ...
            'Units', 'characters', ...
            'Position', [2.3 1.0 26.7 1.3], ...
            'String', 'Subset Truncation', ...
            'Value', 0, ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_checkbox_subsettrunc,handles_gui.figure), ...
            'Interruptible','off');

        % Edit texts
        handles_gui.edit_radius = uicontrol( ...
            'Parent', handles_gui.group_strainoptions, ...
            'Tag', 'edit_radius', ...
            'Style', 'edit', ...
            'Units', 'characters', ...
            'Position', [20.2 2.8 11.1 1.3], ...
            'String', '', ...
            'HorizontalAlignment', 'left', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_edit_radius,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.edit_imgnum = uicontrol( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'edit_imgnum', ...
            'Style', 'edit', ...
            'Units', 'characters', ...
            'Position', [108.2 1.3 7 1.6], ...
            'String', '', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_edit_imgnum,handles_gui.figure), ...
            'Enable', 'off', ...
            'Interruptible','off');

        % Pushbuttons
        handles_gui.button_zoom = uicontrol( ...
            'Parent', handles_gui.group_zoompan, ...
            'Tag', 'button_zoom', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [2.1 1 14.4 1.7], ...
            'String', 'Zoom', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_zoom,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_pan = uicontrol( ...
            'Parent', handles_gui.group_zoompan, ...
            'Tag', 'button_pan', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [17.4 1 14.4 1.7], ...
            'String', 'Pan', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_pan,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_finish = uicontrol( ...
            'Parent', handles_gui.group_menu, ...
            'Tag', 'button_finish', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [2.1 3 29.7 1.7], ...
            'String', 'Finish', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_finish,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_cancel = uicontrol( ...
            'Parent', handles_gui.group_menu, ...
            'Tag', 'button_cancel', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [2.1 1 29.7 1.7], ...
            'String', 'Cancel', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_cancel,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_left = uicontrol( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'button_left', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [101.2 1.2 6 1.8], ...
            'String', '<', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_left,handles_gui.figure), ...
            'Interruptible','off', ...
            'Enable', 'off');

        handles_gui.button_right = uicontrol( ...
            'Parent', handles_gui.group_preview, ...
            'Tag', 'button_right', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [116.2 1.2 6 1.8], ...
            'String', '>', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_right,handles_gui.figure), ...
            'Interruptible','off', ...
            'Enable', 'off');
    end

    % Pause until figure is closed ---------------------------------------%
    waitfor(handles_gui.figure);
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [X] = DCA(networks, dim, rsp, maxiter)
	Q = [];
	for i = 1 : length(networks)
		fileID = char(strcat('../network/', networks(i), '.txt'));
		net = load(fileID);
		tQ = diffusionRWR(net, maxiter, rsp);
		Q = [Q, tQ];
	end

	nnode = size(Q, 1);
	alpha = 1 / nnode;
	Q = log(Q + alpha) - log(alpha);

	Q = Q * Q';
	[U, S] = svds(Q, dim);
	X = U * sqrt(sqrt(S));
end

author: luoyunan

repo: DTINet

save_path: github-repos/MATLAB/luoyunan-DTINet


text:
function range = p13_range ( )

%*****************************************************************************80
%
%% P13_RANGE returns an interval bounding the root for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = - 10000000000000.0;
  range(2) =   10000000000000.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fx = p02_fun ( n, x )

%*****************************************************************************80
%
%% P02_FUN evaluates the integrand for problem 2.
%
%  Discussion:
%
%    The integrand is discontinuous at X = 0.3.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    if ( x < 0.3 )
%      f(x) = 0
%    else
%      f(x) = 1
%
%  Antiderivative:
%
%    if ( x < 0.3 )
%      g(x) = 0
%    else
%      g(x) = X - 0.3
%
%  Exact Integral:
%
%    0.7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner,
%    Comparison of Numerical Quadrature Formulas,
%    in Mathematical Software, edited by John R Rice,
%    Academic Press, 1971.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = ( 0.3 <= x );

  return
end


text:
function plot_data_against_probabilities(sv,type,specs)

if nargin<3

    specs=[3,3];

    if nargin<2

        type='state';

    end

end

% myLimits=[-sqrt(eps),1+sqrt(eps)];

[f,the_regimes,endog,data,tex]=load_filters(sv);

switch type

    case 'state'

        probs=f.smoothed_state_probabilities;

    case 'regime'

        probs=f.smoothed_regime_probabilities;

    otherwise

    error('second input must be either "state" or "regime"')

end

r0=specs(1);

c0=specs(2);

all_regimes=fieldnames(probs);

all_regimes(strcmp(all_regimes,'const_1'))=[];

for iv=1:numel(endog)

    vn=endog{iv};

    utils.plot.multiple(@plotfuncr,all_regimes,...
        [vn,' against smoothed ',type,' probabilities'],r0,c0);

    if ~isempty(tex)

        [~,h]=sup_label(tex.(vn),'t');

        set(h,'fontsize',12)

    end

end

    function [mytex,leg]=plotfuncr(vname)

        dd=data.(vn)(probs.(vname).date_numbers);

        [AX,H1,H2]=plotyy(probs.(vname),dd,'linewidth',2);

        axis(AX,'tight')

        if isempty(tex)

            mytex=vname;

        else

            mytex=tex.(vname);

        end

        if strcmp(type,'regime')

            vname(1:numel('regime'))=[];

            vname=strrep(vname,'_','');

            mytex=sprintf('%s(%s)',mytex,the_regimes{str2double(vname)});

        end

        leg='';

    end

end

author: jmaih


text:
function d = p08_fd ( p )

%*****************************************************************************80
%
%% P08_FD is a signed distance function for problem 08.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Parameters:
%
%    Input, real P, one or more points.
%
%    Output, real D, the signed distance of each point to the boundary of the region.
%
  g1 = dcircle ( p, 0, 0, 1 );
  g2 = drectangle ( protate ( p, pi/12 ), -1, 1, 0, 1 );
  g3 = drectangle ( protate ( p, -pi/12 ), -1, 1, -1, 0 );
  g4 = drectangle ( protate ( pshift ( p, -0.9, 0 ), -pi/4 ), 0, 0.2, 0, 0.2 );
  g5 = dcircle ( p, 0.6, 0, 0.1 );

  d = ddiff ( ddiff ( ddiff ( ddiff ( g1, g2 ), g3 ), g4 ), g5 );

  return
end


text:
function[h]=nocontours(h)
%NOCONTOURS  Removes contours from a CONTOURF plot.
%
%   NOCONTOURS(H) where H is an array of handles to lines, output by
%   CONTOURF, sets the linestyle for all the lines to 'none'.  This
%   removes the contours from the CONTOURF plot.
%
%   NOCONTOURS with no input arguments sets the linestyle of all patch
%   objects (for Matlab versions < 2014b) or all contour objects (for
%   versions 2014b and later) in the current axes to 'none'.
%
%   H=NOCONTOURS also returns the handles to the objects.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2009--2015 J.M. Lilly --- type 'help jlab_license' for details

if nargin==0
   h=contourhandles(gca);
end

set(h,'linestyle','none');

if nargout ==0
  clear h
end

function[h]=contourhandles(axh)
%CONTOURHANDLES  Finds all contour handles from a given set of axes.
%
%   H=CONTOURHANDLES returns handles to all contoures associated with the
%   current axes.
%
%   H=CONTOURHANDLES(AX) returns handles to all contoures associated with
%   the set of axes whose handle is AX.
%
%   See also LINEHANDLES, PATCHHANDLES, AXESHANDLES.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014 J.M. Lilly --- type 'help jlab_license' for details

if nargin==0
   axh=gca;
end

h=get(axh,'children');
bool=false(size(h));
for j=1:length(h)
    if strcmpi(get(h(j),'type'),'contour')
       bool(j)=true;
    end
end
h=h(bool);

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function rule_num = triangle_nco_rule_num ( )

%*****************************************************************************80
%
%% TRIANGLE_NCO_RULE_NUM returns the number of NCO rules available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Output, integer RULE_NUM, the number of rules available.
%
  rule_num = 9;

  return
end


text:
% LLE ALGORITHM (using K nearest neighbors)
%
% [Y] = lle(X,K,dmax)
%
% X = data as D x N matrix (D = dimensionality, N = #points)
% K = number of neighbors
% dmax = max embedding dimensionality
% Y = embedding as dmax x N matrix

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

function [Y] = lle(X,K,d)

[D,N] = size(X);
fprintf(1,'LLE running on %d points in %d dimensions\n',N,D);

% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS
fprintf(1,'-->Finding %d nearest neighbours.\n',K);

X2 = sum(X.^2,1);
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

[sorted,index] = sort(distance);
neighborhood = index(2:(1+K),:);

% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.\n');

if(K>D)
  fprintf(1,'   [note: K>D; regularization will be used]\n');
  tol=1e-3; % regularlizer in case constrained fits are ill conditioned
else
  tol=0;
end

W = zeros(K,N);
for ii=1:N
   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % shift ith pt to origin
   C = z'*z;                                        % local covariance
   C = C + eye(K,K)*tol*trace(C);                   % regularlization (K>D)
   W(:,ii) = C\ones(K,1);                           % solve Cw=1
   W(:,ii) = W(:,ii)/sum(W(:,ii));                  % enforce sum(w)=1
end;

% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W)
fprintf(1,'-->Computing embedding.\n');

% M=eye(N,N); % use a sparse matrix with storage for 4KN nonzero elements
M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N);
for ii=1:N
   w = W(:,ii);
   jj = neighborhood(:,ii);
   M(ii,jj) = M(ii,jj) - w';
   M(jj,ii) = M(jj,ii) - w;
   M(jj,jj) = M(jj,jj) + w*w';
end;

% CALCULATION OF EMBEDDING
options.disp = 0; options.isreal = 1; options.issym = 1;
[Y,eigenvals] = eigs(M,d+1,0,options);
Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0

fprintf(1,'Done.\n');

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

% other possible regularizers for K>D
%   C = C + tol*diag(diag(C));                       % regularlization
%   C = C + eye(K,K)*tol*trace(C)*K;                 % regularlization

author: jindongwang

repo: activityrecognition


text:
function px = gmmEval(x, gmm)
% Computes the probability of the data given a Guassian mixture model
% x(ndata, ndim) - the data to be evaluated
% gmm.mu{K}, gmm.sigma{K}, gmm.priors{K} - the mixture parameters

K = length(gmm.priors);
[ndata, ndim] = size(x);

px = zeros(ndata, 1);

if K ==0
    return;
end

for k = 1:K
    invSigma = inv(gmm.sigma{k});
    for i = 1:ndata
        px(i) = px(i) + 1/(sqrt(2*pi)*sqrt(det(gmm.sigma{k}))) * ...
            exp(-1/2 * (x(i, :)-gmm.mu{k}) * invSigma * (x(i, :)-gmm.mu{k})') ...
            * gmm.priors(k);
    end
end

author: zouchuhang

repo: LayoutNet


text:
function bool = is_dim_custom(varargin)
%IS_DIM_CUSTOM true for a custom grid
%
%   Usage: bool = is_grid_custom(x1,x2,...)
%
%   Input parameters:
%       x1,x2,... - axis / m; single value or [xmin,xmax] or nD-array
%
%   Output parameters:
%       bool      - array of logical indicating whether each input is an
%                   nD-array
%
%   IS_DIM_CUSTOM(x1,x2,..) checks if we have a custom grid by checking if any
%   of the given x,y,z values is a nD-array.
%
%   See also: xyz_grid, xyz_axes_selection, plot_sound_field

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking input parameters =======================================
isargnumeric(varargin{:});

%% ===== Computation =====================================================
bool = cellfun(@(x) numel(x)>2, varargin);

author: sfstoolbox

repo: sfs-matlab


text:
function [AnDartest] = AnDartest(x,alpha)
%ANDARTEST Anderson-Darling test for assessing normality of a sample data.
% The Anderson-Darling test (Anderson and Darling, 1952) is used to test if
% a sample of data comes from a specific distribution. It is a modification
% of the Kolmogorov-Smirnov (K-S) test and gives more weight to the tails
% than the K-S test. The K-S test is distribution free in the sense that the
% critical values do not depend on the specific distribution being tested.
% The Anderson-Darling test makes use of the specific distribution in calculating
% critical values. This has the advantage of allowing a more sensitive test
% and the disadvantage that critical values must be calculated for each
% distribution.
% The Anderson-Darling test is only available for a few specific distributions.
% The test is calculated as:
%
%        AD2 = integral{[F_o(x)-F_t(x)]^2/[F_t(x)(1-F_t(x)0]}dF_t(x)
%
%        AD2a = AD2*a
%
% Note that for a given distribution, the Anderson-Darling statistic may be
% multiplied by a constant, a (which usually depends on the sample size, n).
% These constants are given in the various papers by Stephens (1974, 1977a,
% 1977b, 1979, 1986). This is what should be compared against the critical
% values. Also, be aware that different constants (and therefore critical
% values) have been published. You just need to be aware of what constant
% was used for a given set of critical values (the needed constant is typically
% given with the critical values).
% The critical values for the Anderson-Darling test are dependent on the
% specific distribution that is being tested. Tabulated values and formulas
% have been published for a few specific distributions (normal, lognormal,
% exponential, Weibull, logistic, extreme value type 1). The test is a one-sided
% test and the hypothesis that the distribution is of a specific form is
% rejected if the test statistic, AD2a, is greater than the critical value.
% Here we develop the m-file for detecting departure from normality. It is
% one of the most powerful statistics for test this.
%
% Syntax: function AnDartest(X)
%
%     Input:
%          x - data vector
%      alpha - significance level (default = 0.05)
%
%     Output:
%            - Complete Anderson-Darling normality test
%
% Example: From the Table 1 base data of secondary coating diameter, SCD
%    (micrometers) for an optical fiber manufacturing and testing [Application
%    and acceptance sampling in testing of optical fiber (Bhaumik and Bhargava,
%    2005)], taked from Internet:
%    http://www.sterliteoptical.com/pdf/ProdWithContents/Application%20of%
%    20acceptance%20sampling%20in%20testing%20of%20optical%20fiber.pdf
%    The measurements of the sample are analyzed for the normality using the
%    Anderson-Darling test with a significance of 0.05.
%
% Data vector is:
%  x=[245.3;245.6;245.0;244.7;244.6;244.6;245.9;245.2;245.5;246.1;246.5;246.7;
%  246.0;245.8;245.3;245.3;246.0;245.5;246.7;245.3;245.8;245.8;245.2;245.7;
%  244.8;246.7;246.5;245.5;246.0;244.8;244.8;244.4;244.8;245.6;245.5;244.0;
%  245.2;244.8;245.4;246.1;245.9;245.6;245.2];
%
% Calling on Matlab the function:
%            AnDartest(x)
%
% Answer is:
%
% Sample size: 43
% Anderson-Darling statistic: 0.2858
% Anderson-Darling adjusted statistic: 0.2912
% Probability associated to the Anderson-Darling statistic = 0.6088
% With a given significance = 0.050
% The sampled population is normally distributed.
% Thus, this sample have been drawn from a normal population with a mean & variance = 245.4814    0.4139
%
% Created by A. Trujillo-Ortiz, R. Hernandez-Walls, K. Barba-Rojo and
%             A. Castro-Perez
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%
% Copyright. April 20, 2007.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls, K. Barba-Rojo and A. Castro-Perez. (2007).
%   AnDartest:Anderson-Darling test for assessing normality of a sample data.
%   A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/
%   fileexchange/loadFile.do?objectId=14807
%
% References:
%  Anderson, T. W. and Darling, D. A. (1952), Asymptotic theory of certain
%      'goodness-of-fit' criteria based on stochastic processes. Annals of
%      Mathematical Statistics, 23:193-212.
%  Stephens, M. A. (1974), EDF Statistics for goodness of fit and some
%      comparisons. Journal of the American Statistical Association,
%      69:730-737.
%  Stephens, M. A. (1976), Asymptotic Results for goodness-of-fit statistics
%      with unknown parameters. Annals of Statistics, 4:357-369.
%  Stephens, M. A. (1977a), Goodness of fit for the extreme value distribution.
%      Biometrika, 64:583-588.
%  Stephens, M. A. (1977b), Goodness of fit with special reference to tests
%      for exponentiality. Technical Report No. 262, Department of Statistics,
%      Stanford University, Stanford, CA.
%  Stephens, M. A. (1979), Tests of fit for the logistic distribution based
%      on the empirical distribution function. Biometrika, 66:591-595.
%  Stephens, M. A. (1986), Tests based on EDF statistics. In: D'Agostino,
%      R.B. and Stephens, M.A., eds.: Goodness-of-Fit Techniques. Marcel
%      Dekker, New York.
%

switch nargin
    case{2}
        if isempty(x) == false && isempty(alpha) == false
            if (alpha <= 0 || alpha >= 1)
                fprintf('Warning: Significance level error; must be 0 < alpha < 1 \n');
                return;
            end
        end
    case{1}
        alpha = 0.05;
    otherwise
        error('Requires at least one input argument.');
end

n = length(x);

if n < 7,
    disp('Sample size must be greater than 7.');
    return,
else
    x = x(:);
    x = sort(x);
    fx = normcdf(x,mean(x),std(x));
    i = 1:n;

    S = sum((((2*i)-1)/n)*(log(fx)+log(1-fx(n+1-i))));
    AD2 = -n-S;

    AD2a = AD2*(1 + 0.75/n + 2.25/n^2); %correction factor for small sample sizes: case normal

    if (AD2a >= 0.00 && AD2a < 0.200);
        P = 1 - exp(-13.436 + 101.14*AD2a - 223.73*AD2a^2);
    elseif (AD2a >= 0.200 && AD2a < 0.340);
        P = 1 - exp(-8.318 + 42.796*AD2a - 59.938*AD2a^2);
    elseif (AD2a >= 0.340 && AD2a < 0.600);
        P = exp(0.9177 - 4.279*AD2a - 1.38*AD2a^2);
    else (AD2a >= 0.600 && AD2a <= 13);
        P = exp(1.2937 - 5.709*AD2a + 0.0186*AD2a^2);
    end
end

disp(' ')
fprintf('Sample size: %i\n', n);
fprintf('Anderson-Darling statistic: %3.4f\n', AD2);
fprintf('Anderson-Darling adjusted statistic: %3.4f\n', AD2a);
fprintf('Probability associated to the Anderson-Darling statistic = %3.4f\n', P);
fprintf('With a given significance = %3.3f\n', alpha);
if P >= alpha;
   disp('The sampled population is normally distributed.');
   s = [mean(x),var(x)];
   fprintf('Thus, this sample have been drawn from a normal population with a mean & variance = %6.4f  %8.4f\n',s);
else
   disp('The sampled population is not normally distributed.');
end

return,

author: Sable

repo: mcbench-benchmarks


text:
% Convert adjacency matrix (nxn) to edge list (mx3)
%
% INPUTS: adjacency matrix: nxn
% OUTPUTS: edge list: mx3
%
% GB: last updated, Sep 24, 2012

function el=adj2edgeL(adj)

n=length(adj); % number of nodes
edges=find(adj>0); % indices of all edges

el=[];
for e=1:length(edges)
  [i,j]=ind2sub([n,n],edges(e)); % node indices of edge e
  el=[el; i j adj(i,j)];
end

author: aeolianine

repo: octave-networks-toolbox


text:
%% ricianpdf
% Below is a demonstration of the features of the |ricianpdf| function

%%
clear; close all; clc;

%% Syntax
% |P=ricianpdf(A,M,S);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function Z = replace(X,Y,W,expand)
%REPLACE Substitutes variables
%
%Z = REPLACE(X,Y,W)  Replaces any occurence of the SDPVAR object Y
%                    in the SDPVAR object X with the expression W
%
% Example
%  x = sdpvar(1,1);
%  t = sdpvar(1,1);
%  Y = [1+t;1+x+t];
%  Y = replace(Y,x,2) generates Y=[1+t;3+t]

if nargin<4
    expand = 1;
end

if ~isa(X,'sdpvar')
    Z = X;
    return
end
if ~isa(Y,'sdpvar')
    error('Second arguments must be an sdpvar object')
end

if ~is(Y,'linear')
    error('Second arguments must be linear')
end

if prod(size(W)) == 1
    W = repmat(W,size(Y));
end

if ~isequal(size(Y),size(W))
    if isequal(fliplr(size(Y)),size(W))
        W = W';
    else
        error('Both arguments must have same size')
    end
end

if max(size(Y))>1
    if ~isequal(Y.basis(:,2:end),speye(prod(Y.dim)))
        [Y,keptY] = unique(reshape(Y,[],1));
        W = extsubsref(W,keptY);
    end
end

if isa(W,'sdpvar')
    % This is tricky...
    Z = variable_replace(X,Y,W);
    return
end

if ~isnumeric(W)
    error('Third arguments must be numeric')
end

% Replace with NaN   destroys everything, assume it should be cleared
W(isnan(W)) = 0;

y_lmi_variables = Y.lmi_variables;
b = W(:)-Y.basis(:,1);
A = Y.basis(:,2:end);
feas_var = A\b;
if norm(A*feas_var-b)>sqrt(eps)
    error('Inconsistent assignment')
end

x_lmi_variables = X.lmi_variables;
n = X.dim(1);
m = X.dim(2);

[monomtable,variabletype] = yalmip('monomtable');
if all(variabletype(x_lmi_variables)==0) % is(X,'linear')
    Z = X.basis(:,1);
    %v = [];
    v1 = [];
    v2 = [];
    i1 = [];
    i2 = [];
    for i = 1:length(x_lmi_variables)
        j = find(x_lmi_variables(i) == y_lmi_variables);
        if isempty(j)
          %  v = [v ;recover(x_lmi_variables(i))];
            v1 = [v1 x_lmi_variables(i)];
            i1 = [i1 i];
            %Z = Z + recover(x_lmi_variables(i))*X.basis(:,i+1);
        else
          %  v = [v ;feas_var(j)];
            v2 = [v2 j];
            i2 = [i2 i];
            %Z = Z + feas_var(j)*X.basis(:,i+1);
        end
    end
    v = sparse(i1,ones(length(i1),1),recover(v1),length(x_lmi_variables),1);
    v = v + sparse(i2,ones(length(i2),1),feas_var(v2),length(x_lmi_variables),1);
    Z = Z + X.basis(:,2:end)*v;
else
    base = getbase(Y);base = base(:,2:end);
    [i,j,k] = find(base);
    replaced_vars = getvariables(Y);
    replaced_vars = replaced_vars(i);
    %for i = 1:length(Y)
    %    replaced_vars(i) = getvariables(extsubsref(Y,i));
    %end
    % used_variables = getvariables(X);
    used_variables = x_lmi_variables;
    %  monomtable = yalmip('monomtable');
    local_monom = monomtable(used_variables,replaced_vars);
    W = W(:)';
    gain = zeros(length(used_variables),1);
    for i = 1:length(used_variables)
        % F**N 6.5 0^sparse(0) and 0^0 differ
        gain(i) = prod(W.^full(local_monom(i,:)));
    end

    local_monoms_left = monomtable(used_variables,:);
    local_monoms_left(:,replaced_vars) = 0;
    used_left = find(sum(local_monoms_left,1));
    base = recovermonoms(local_monoms_left(:,used_left),recover(used_left));
    base = base.*gain(:);
    Z = X.basis(:,1);
    Z = Z + X.basis(:,2:end)*base;
end

if expand
    Xvariables = getvariables(Z);
    extvar = yalmip('extvariables');
    Xext = find(ismember(Xvariables,extvar));
    if ~isempty(Xext)
        %We must dig down in extended operators to see if they use the replaced
        %set of variables
        for i = 1:length(Xext)
            extstruct = yalmip('extstruct',Xvariables(Xext(i)));
            anychanged = 0;
            for j = 1:length(extstruct.arg)
                if isa(extstruct.arg{j},'sdpvar')
                    XinY = find(ismembc(getvariables(extstruct.arg{j}),y_lmi_variables));
                    if ~isempty(XinY)
                        anychanged = 1;
                        extstruct.arg{j} = replace(extstruct.arg{j},Y,W);
                    else
                    end
                end
            end
            if anychanged
                if isequal(extstruct.fcn,'pwa_yalmip') | isequal(extstruct.fcn,'pwq_yalmip')
                    % Change data in MPT structure to allow plotting of PWA
                    % and PWQ functions (and improve numerics)
                    % FIXME : Generalize code from sdpvar/plot
                    % extstruct = derive_new_mpt_data(extstruct);
                end
                Zi = yalmip('define',extstruct.fcn,extstruct.arg{:});
                Xvariables(Xext(i)) = getvariables(Zi);
            end
        end
        % And now recover this sucker
        Z = struct(Z);
        Z.lmi_variables = Xvariables;
        % Fucked up order (lmi_variables should be sorted)
        if any(diff(Z.lmi_variables)<0)
            [i,j]=sort(Z.lmi_variables);
            Z.basis = [Z.basis(:,1) Z.basis(:,j+1)];
            Z.lmi_variables = Z.lmi_variables(j);
        end
        Z = sdpvar(Z.dim(1),Z.dim(2),[],Z.lmi_variables,Z.basis);
    end
end

if isa(Z,'sdpvar')
    Z.dim(1) = n;
    Z.dim(2) = m;
    Z.typeflag = X.typeflag;
    Z.extra = X.extra;
    % Reset info about conic terms
    Z.conicinfo = [0 0];
else
    Z = reshape(full(Z),n,m);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function STATS=mwwtest(x1,x2)
% Mann-Whitney-Wilcoxon non parametric test for two unpaired groups.
% This file execute the non parametric Mann-Whitney-Wilcoxon test to evaluate the
% difference between unpaired samples. If the number of combinations is less than
% 20000, the algorithm calculates the exact ranks distribution; else it
% uses a normal distribution approximation. The result is not different from
% RANKSUM MatLab function, but there are more output informations.
% There is an alternative formulation of this test that yields a statistic
% commonly denoted by U. Also the U statistic is computed.
%
% Syntax: 	STATS=MWWTEST(X1,X2)
%
%     Inputs:
%           X1 and X2 - data vectors.
%     Outputs:
%           - T and U values and p-value when exact ranks distribution is used.
%           - T and U values, mean, standard deviation, Z value, and p-value when
%           normal distribution is used.
%        If STATS nargout was specified the results will be stored in the STATS
%        struct.
%
%      Example:
%
%         X1=[181 183 170 173 174 179 172 175 178 176 158 179 180 172 177];
%
%         X2=[168 165 163 175 176 166 163 174 175 173 179 180 176 167 176];
%
%           Calling on Matlab the function: mwwtest(X1,X2)
%
%           Answer is:
%
% MANN-WHITNEY-WILCOXON TEST
%
%                        Group_1    Group_2
%                        _______    _______
%
%     Numerosity          15         15
%     Sum_of_Rank_W      270        195
%     Mean_Rank           18         13
%     Test_variable_U     75        150
%
% Sample size is large enough to use the normal distribution approximation
%
%     Mean       SD        Z       p_value_one_tail    p_value_two_tails
%     _____    ______    ______    ________________    _________________
%
%     112.5    24.047    1.5386    0.061947            0.12389
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2009). MWWTEST: Mann-Whitney-Wilcoxon non parametric test for two unpaired samples.
% http://www.mathworks.com/matlabcentral/fileexchange/25830
%Input Error handling
p = inputParser;
addRequired(p,'x1',@(x) validateattributes(x,{'numeric'},{'row','real','finite','nonnan','nonempty'}));
addRequired(p,'x2',@(x) validateattributes(x,{'numeric'},{'row','real','finite','nonnan','nonempty'}));
parse(p,x1,x2);
%set the basic parameter
n1=length(x1); n2=length(x2); NP=n1*n2; N=n1+n2; N1=N+1; k=min([n1 n2]);
[A,B]=tiedrank([x1(:); x2(:)]); %compute the ranks and the ties
R1=A(1:n1); R2=A(n1+1:end);
T1=sum(R1); T2=sum(R2);
U1=NP+(n1*(n1+1))/2-T1; U2=NP-U1;
disp('MANN-WHITNEY-WILCOXON TEST')
disp(' ')
disp(table([n1;T1;T1/n1;U1],[n2;T2;T2/n2;U2],...
    'VariableNames',{'Group_1' 'Group_2'},...
    'RowNames',{'Numerosity' 'Sum_of_Rank_W' 'Mean_Rank' 'Test_variable_U'}))
if nargout
    STATS.n=[n1 n2];
    STATS.W=[T1 T2];
    STATS.mr=[T1/n1 T2/n2];
    STATS.U=[U1 U2];
end
if round(exp(gammaln(N1)-gammaln(k+1)-gammaln(N1-k))) > 20000
    mU=NP/2;
    if B==0
        sU=realsqrt(NP*N1/12);
    else
        sU=realsqrt((NP/(N^2-N))*((N^3-N-2*B)/12));
    end
    Z1=(abs(U1-mU)-0.5)/sU;
    p=1-normcdf(Z1); %p-value
    disp('Sample size is large enough to use the normal distribution approximation')
    disp(' ')
    disp(table(mU,sU,Z1,p,2*p,'VariableNames',{'Mean' 'SD' 'Z' 'p_value_one_tail' 'p_value_two_tails'}))
    if nargout
        STATS.method='Normal approximation';
        STATS.mU=mU;
        STATS.sU=sU;
        STATS.Z=Z1;
        STATS.p=[p 2*p];
    end
else
    disp('Sample size is small enough to use the exact Mann-Whitney-Wilcoxon distribution')
    disp(' ')
    if n1<=n2
        w=T1;
    else
        w=T2;
    end
    pdf=sum(nchoosek(A,k),2);
    P = [sum(pdf<=w) sum(pdf>=w)]./length(pdf);
    p = min(P);
    disp(table(w,p,2*p,'VariableNames',{'W' 'p_value_one_tail' 'p_value_two_tails'}))
    if nargout
        STATS.method='Exact distribution';
        STATS.T=w;
        STATS.p=[p 2*p];
    end
end
disp(' ')

author: kavli-ntnu


text:
function [output1, output2] = functionTemplate(model, arg1, optArg1, optArg2, optArg3)
% (Description of the function) Template for a main utility function for
% the COBRA toolbox. Allow passing COBRA-toolbox-wide parameters (e.g., feasTol, optTol, printLevel)
% and solver-specific parameter structure (e.g., struct('lpmethod', 1) for Cplex, struct('Presolve', 1) for Gurobi)
% to the function and other functions called by the function.
%
% USAGE:
%     (positional argument inputs)
%     [output1, output2] = functionTemplate(model, arg1, optArg1, optArg2, optArg3, solverParams)
%     (positional + name/value argument inputs)
%     [output1, output2] = functionTemplate(model, arg1, 'name', value, ..., solverParams)
%     (parameter structure input)
%     [output1, output2] = functionTemplate(model, ..., paramStruct)
%
% INPUT(S):
% (must be inputted as positional arguments)
%     model:               COBRA model
%     arg1:                argument 1
%
% OPTIONAL INPUT(S):
% (can be inputted as positional arguments, or name/value arguments, or a single parameter structure)
%     optArg1:            optional argument 1
%     optArg2:            optional argument 2
%     optArg3:            optional argument 3
%     ...
%
% OUTPUT(S):
%     output1:            output argument 1
%     output2:            output argument 2
%
% EXAMPLE:
%
% NOTE:
%
% Author(s):
%

% optional arguments that can be inputted as positional arguments,
% name/value arguments, or parameter structure, in the same order as
% specified in USAGE
optArgin =      {'optArg1', 'optArg2', 'optArg3'};
% defaultValues for the optional arguments above
defaultValues = {1,         'test',    {'test'}};
% validators for the optional arguments above, i.e., the inputs when passed
% to the validators must return true
validator = {...
    @(x) isscalar(x) & isnumeric(x) & x >= 0, ... % optArg1 (must be non-negative numbers)
    @ischar, ...                                  % optArg2 (must be character)
    @(x) true, ...                                % optArg3 (accept anything)
    };

% types of optimization problems to be solved in this function and the
% functions called by this function
problemTypes = {'LP', 'QP'};

% provide the optional argument name for the solver-specific parameter
% structure (to be passed to solvers like Cplex, Gurobi in e.g., solveCobraLP)
% This is intended for backward compatibility only and NOT RECOMMENDED TO USE
% because as a convention among cobra functions, solver-specific parameter
% structure can be inputted without keyword
keyForSolverParams = '';

% treat empty inputs ([]) for positional arguments as using default values or not.
% For example, if emptyForDefault = true, functionTemplate(model, [], 'test')
% takes arg1 as the default value 1. But if emptyForDefault = false, this
% returns error since '[]' does not pass the validator @(x) isscalar(x) & isnumeric(x) & x >= 0
emptyForDefault = true;

% run the input parser to get (1) function arguments, (2) cobra parameters
% and (3) solver-specific parameters
[funParams, cobraParams, solverVarargin] = parseCobraVarargin(varargin, ...
    optArgin, defaultValues, validator, problemTypes, keyForSolverParams, emptyForDefault);

% get all the optional function arguments
[optArg1, optArg2, optArg3] = deal(funParams{:});

%%%%% Do what you want to do with the function. For example, construct
%%%%% optimization problems and solve

LPproblem = struct();
[LPproblem.A, LPproblem.b, LPproblem.c, LPproblem.lb, LPproblem.ub, ...
    LPproblem.osense, LPproblem.csense] = deal(speye(2), [1; 1], [1; 0], ...
    [0; 0], [1; 1], -1, 'LL');

% pass the solver-specific parameters to optimization problem
% make sure to use {:}
solLP = solveCobraLP(LPproblem, solverVarargin.LP{:});
LPproblem.F = speye(2);
solQP = solveCobraQP(LPproblem, solverVarargin.QP{:});

% use the cobra parameters (e.g., manually check feasibility)
if ~(checkSolFeas(solLP) <= cobraParams.LP.feasTol)
    [output1, output2] = deal(0);
else
    [output1, output2] = deal(1);
end

help functionTemplate

end

author: opencobra

repo: cobratoolbox


text:
function [AF,BF]=uwfbtbounds(wt,varargin)
%UWFBTBOUNDS Frame bounds of Undecimated WFBT
%   Usage: fcond=uwfbtbounds(wt,L);
%          [A,B]=uwfbtbounds(wt,L);
%          [...]=uwfbtbounds(wt);
%
%   `uwfbtbounds(wt,L)` calculates the ratio $B/A$ of the frame bounds
%   of the undecimated filterbank specified by *wt* for a system of length
%   *L*. The ratio is a measure of the stability of the system.
%
%   `uwfbtbounds({w,J,'dwt'},L)` calculates the ratio $B/A$ of the frame
%   bounds of the undecimated DWT (|UFWT|) filterbank specified by *w* and
%   *J* for a system of length *L*.
%
%   `uwfbtbounds(wt)` does the same thing, but *L* is the length of the
%   longest filter in the identical filterbank.
%
%   `[A,B]=uwfbtbounds(...)` returns the lower and upper frame bounds
%   explicitly.
%
%   See |wfbt| for explanation of parameter *wt* and |fwt| for explanation
%   of parameters *w* and *J*.
%
%   The function supports the following flags:
%
%   `'sqrt'`(default),`'noscale'`,`'scale'`
%       The filters in the filterbank tree are scaled to reflect the
%       behavior of |uwfbt| and |iuwfbt| with the same flags.
%
%   See also: uwfbt, filterbankbounds

warning('UWFBTBOUNDS is deprecated. Please use WFBTBOUNDS with a appropriate flag.');
complainif_notenoughargs(nargin,1,'UWFBTBOUNDS');

definput.keyvals.L = [];
definput.import = {'uwfbtcommon'};
[flags,~,L]=ltfatarghelper({'L'},definput,varargin);

if nargout<2
   AF = wfbtbounds(wt,L,flags.scaling);
elseif nargout == 2
   [AF,BF] = wfbtbounds(wt,L,flags.scaling);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
I=imnoise(I, 'salt & pepper', 0.05);
J=medfilt2(I, [3, 3]);
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function inside = p05_inside ( m, n, point )

%*****************************************************************************80
%
%% P05_INSIDE reports if a point is inside the region in problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, logical INSIDE(N), is TRUE if the point is in the region.
%
  center1 = [  0.0, 0.0 ];
  center2 = [ -0.4, 0.0 ];
  r1 =   1.00;
  r2 =   0.55;

  inside(1:n) =                                                   ...
       center1(2) <=   point(2,1:n)                               ...
       &                                                          ...
                     ( point(1,1:n) - center1(1) ).^2             ...
                   + ( point(2,1:n) - center1(2) ).^2  <= r1 * r1 ...
       &                                                          ...
       r2 * r2 <=    ( point(1,1:n) - center2(1) ).^2             ...
                   + ( point(2,1:n) - center2(2) ).^2;

  return
end


text:
function [sol,times,ocp] = bouncingball

  before_contact = ocl.Stage([], @before_contact_vars, @before_contact_ode, 'N', 3, 'd', 2);
  after_contact = ocl.Stage(1, @after_contact_vars, @after_contact_ode, ...
                            @after_contact_cost, 'N', 5, 'd', 2);

  before_contact.setInitialStateBounds('s', 1);
  before_contact.setInitialStateBounds('v', 0);
  before_contact.setEndStateBounds('s', 0);

  after_contact.setEndStateBounds('s', 1);

  ocp = ocl.MultiStageProblem({before_contact, after_contact}, {@stage_transition});

  [sol,times] = ocp.solve(ocp.getInitialGuess());

  % stage 1
  figure;
  subplot(1,2,1)
  hold on; grid on;
  ocl.plot(times{1}.states, sol{1}.states.s)
  ocl.plot(times{1}.states, sol{1}.states.v)
  legend({'s','v'})
  xlabel('time [s]');
  ylim([-5 3])
  yticks(-5:3)
  title('stage 1')

  % stage 2
  subplot(1,2,2)
  hold on; grid on;
  ocl.plot(times{2}.states, sol{2}.states.s)
  ocl.plot(times{2}.states, sol{2}.states.v)
  ocl.stairs(times{2}.controls, sol{2}.controls.F)
  legend({'s','v','F'})
  xlabel('time [s]');
  ylim([-5 3])
  yticks(-5:3)
  title('stage 2')

end

function before_contact_vars(sh)
  sh.addState('s');
  sh.addState('v');
end

function before_contact_ode(sh,x,~,~,~)
  sh.setODE('s', x.v);
  sh.setODE('v', -10);
end

function after_contact_vars(sh)
  sh.addState('s');
  sh.addState('v');
  sh.addControl('F');
end

function after_contact_ode(sh,x,~,u,~)
  sh.setODE('s', x.v);
  sh.setODE('v', -10 + 10*u.F);
end

function after_contact_cost(ch,~,~,u,~)
  ch.add( u.F^2 );
end

function stage_transition(ch, x0, xF)
  % x0 current stage
  % xF previous stage
  ch.add(x0.s, '==', xF.s);
  ch.add(x0.v, '==', -xF.v/2);
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function soln = particleSwarmOptimization(problem)
% soln = particleSwarmOptimization(problem)
%
% This function computes the solution to the optimization problem using
% particle swarm optimization
%
% INPUTS:
%   problem = struct that describes the optimization problem
%       .options.maxIter = number of iterations for the search
%       .options.populationCount = number of particles to use in search
%       .options.w = search parameter - damping
%       .options.pl = search paramter - local search
%       .options.pg = search parameter - global search
%       .xLow = [nDim, 1] = column vector for lower bound on initial search space
%       .xUpp = [nDim, 1] = column vector for upper bound on initial search space
%       .objFun = function handle:
%           f = objFun(x)
%               x = [nDim, nPop] = vectorized objective function
%               f = [1,nPop] = cost for each column in x
%
% OUTPUTS:
%   soln = output struct with solution and history log
%       .x = best point found
%       .f = objective function at best point found
%       .log = history of the search, used for plotting and analysis
%       .problem = copy of the problem struct
%
% NOTES:
%   This is just a quick test function that I wrote to study how the
%   optimization parameters effect the behavior of the optimization, as
%   well as to figure out how it works with a noisy objective function. For
%   more rigerous optimization the code should be improved, most notably by
%   adding a convergence criteria and exit flag.
%
%

maxIter = problem.options.maxIter;
nPop = problem.options.populationCount;
nDim = length(problem.xLow);

% Memory allocation
Gx = zeros(nDim,1,maxIter);  %Global best solution point
Gf = zeros(1,1,maxIter);  %Global best solution value
Lx = zeros(nDim,nPop,maxIter);  %Local best solution
Lf = zeros(1,nPop,maxIter);  %Local best value
X = zeros(nDim,nPop,maxIter);  %Particle location
F = zeros(1,nPop,maxIter);  %objFun(X);
V = zeros(nDim,nPop,maxIter);  %Particle velocity

% Linear transform for vectorized random numbers
a = problem.xLow*ones(1,nPop);
b = (problem.xUpp-problem.xLow)*ones(1,nPop);

% Initialization:
X(:,:,1) = a + b.*rand(nDim,nPop);
V(:,:,1) = -abs(b) + 2*abs(b).*rand(nDim,nPop);
F(1,:,1) = problem.objFun(X(:,:,1));
Lx(:,:,1) = X(:,:,1);
Lf(1,:,1) = F(1,:,1);
[Gf(1,1,1),idx] = min(F(1,:,1));
Gx(:,1,1) = Lx(:,idx,1);

% Run the optimization here:
pl = problem.options.pl;
pg = problem.options.pg;
w = problem.options.w;
for i=2:maxIter
    rl = rand(nDim,nPop);  %Local random search dist
    rg = rand(nDim,nPop);  %global random search dist

    %%% Update all points
    V(:,:,i) = w*V(:,:,i-1) + ...    % Inertial terms
        pl.*rl.*(Lx(:,:,i-1)-X(:,:,i-1)) + ...  % Local search
        pg.*rg.*(Gx(:,1,i-1)*ones(1,nPop)-X(:,:,i-1));   %Global search
    X(:,:,i) = X(:,:,i-1) + V(:,:,i);  %Position update
    F(1,:,i) = problem.objFun(X(:,:,i));  %Objective function evaluation

    %%% Update local best:
    update = F(1,:,i) < Lf(1,:,i-1);  %Points that have new best values
    Lf(1,update,i) = F(1,update,i);
    Lf(1,~update,i) = Lf(1,~update,i-1);
    Lx(:,update,i) = X(:,update,i);
    Lx(:,~update,i) = Lx(:,~update,i);

    %%% Update global best:
    [Gf(1,1,i),idx] = min(Lf(1,:,i));
    Gx(:,1,i) = Lx(:,idx,i);

end

% Store the final solution
soln.x = reshape(Gx(:,1,end),nDim,1);
soln.f = Gf(end);

% Store all of the search history for analysis
soln.log.Gx = Gx;
soln.log.Gf = Gf;
soln.log.Lx = Lx;
soln.log.Lf = Lf;
soln.log.X = X;
soln.log.F = F;
soln.log.V = V;
soln.problem = problem;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
 function I_sr = superresolve_bvsr(slidingWindow, magFactor)

    vid_l_bic = slidingWindow.frames;

    % Reshape input images from 3-D array to cell array.
    if ~iscell(vid_l_bic)
        if ndims(vid_l_bic) == 3
            for j = 1:size(vid_l_bic,3)
                I{j}(:,:,1) = vid_l_bic(:,:,j);
                I{j}(:,:,2) = vid_l_bic(:,:,j);
                I{j}(:,:,3) = vid_l_bic(:,:,j);
            end
            vid_l_bic = I;
        end
    end

    % Set parameters
    param.LINUX = 1;
    param.MOTION_ESTIMATION = 1;
    param.OPTICAL_FLOW = 1;
    param.BLUR_ESTIMATION = 1;
    param.NOISE_ESTIMATION = 1;
    param.SHOW_IMAGE = 0;
    param.SAVE_RESULT = 1;
    param.MAKE_VIDEO = 0;

    opt(1).nref = round((length(vid_l_bic) - 1) / 2); % num of reference frames (one-side)
    opt(1).res = magFactor; % for scale factor
    opt(1).maxit = 5; % maximum iteration number
    opt(1).eps = 0.0005;         % iteration stop criteria
    opt(1).eps_out = 0.0001;    % iteration stop criteria
    opt(1).eps_blur = 0.0001;    % iteration stop criteria
    opt(1).eta = 20; % for derivative of image
    opt(1).xi = 0.7; % for derivative of kernel
    opt(1).alpha = 1;  % for noise
    opt(1).beta = 0.1; % for noise
    opt(1).degrad = 0; % 0:imresize, 1:LPF+downsampling
    opt(1).hsigma = 0.4*magFactor; % for degradation (blur kernel)
    opt(1).hsize = round(6*magFactor*0.4); % for degradation (blur kernel)
    if mod(opt(1).hsize, 2) == 0
        opt(1).hsize = opt(1).hsize + 1;
    end
    opt(1).noisesd = 0; % for degradation (noise st.dev.) 0,0.01,0.03,0.05
    if param.BLUR_ESTIMATION % Initialization of estimated blur kernel (h_2d)
        opt(1).hmode_init = 1; % 1:gaussian, 2:uniform
        opt(1).hsigma_init = 0.4*magFactor;
    end

    for l = 1:length(opt)

        opt(l).nFrames = length(vid_l_bic);
        opt(l).maxnbf = min(opt(l).nFrames, opt(l).nref); % maximum frame numbers for reference

        [opt(l).M,opt(l).N] = size(vid_l_bic{1}(:,:,1)); % height and width of high-res image
        opt(l).M = opt(1).res * opt(l).M;
        opt(l).N = opt(1).res * opt(l).N;
        [opt(l).m,opt(l).n] = size(vid_l_bic{1}(:,:,1)); % height and widht of low-res image
        M = opt(l).M;
        N = opt(l).N;

        h_2d_sim = fspecial('gaussian', [opt(l).hsize opt(l).hsize], opt(l).hsigma); % blur kernel
        if param.BLUR_ESTIMATION
            [h_1d,h_2d_init] = create_h1(opt(l));
            h_2d = h_2d_init;
        else
            h_2d = h_2d_sim;
        end

        J = extract_y(vid_l_bic);
        I_init = create_initial_I(J,opt(l).res, h_2d);

        % Initialize the variables
        W0 = zeros(opt(l).m,opt(l).n); % for weight matrix for high-res img
        Ws = zeros(M,N); % for weight matrix for derivative
        Wk = zeros(M,N); % for weight matrix for kernel
        Wi = []; % for weight matrix for high-res img (neighboring frames)
        I_sr = I_init; % initialization for super-resolved image

        % Loop for each frame
        for j = round((opt(l).nFrames + 1) / 2)

            n_back = min(opt(l).maxnbf, j-1);
            n_for = min(opt(l).maxnbf, opt(l).nFrames-j);
            th = ones(n_for+n_back+1, 1); % for noise estimation

            % Coordinate descent algorithm
            I = I_init{j}; % current frame (high-res)
            J0 = J{j}; % current frame (low-res)
            for i = -n_back:n_for
                FI{j+i} = I_init{j+i};
            end

            %  Outer iteration
            for k = 1:opt(l).maxit % Loop for each sweep of the algorithm

                I_old_out = I;

                % (1) Estimate motion
                % IRLS
                % motion estimation with optical flow algorithm
                if param.MOTION_ESTIMATION && param.OPTICAL_FLOW
                    for i = -n_back:n_for
                        if i == 0
                            u{j+i} = zeros(size(I));
                            v{j+i} = zeros(size(I));
                            ut{j+i} = zeros(size(I));
                            vt{j+i} = zeros(size(I));
                        else
                            %[u{j+i},v{j+i}] = opticalFlow_CLG_TV(I,I_sr{j+i}); % I->Ii(J_bic)
                            oflParams = [0.025, ...  % Regularization weight
                                0.5, ...    % Downsample ratio
                                20, ...     % Width of image pyramid coarsest level
                                7, ...      % Number of outer fixed point iterations
                                1, ...      % Number of inner fixed point iterations
                                20];        % Number of SOR iterations
                            [u{j+i},v{j+i}] = Coarse2FineTwoFrames(I, I_sr{j+i}, oflParams);
                            %[ut{j+i},vt{j+i}] = Coarse2FineTwoFrames(I_sr{j+i}, I, oflParams);
                            ut{j+i} = -u{j+i};
                            vt{j+i} = -v{j+i};
                            FI{j+i} = warped_img(I,u{j+i},v{j+i});
                        end
                    end
                end

                % (2) Estimate noise
                Nq = opt(l).m*opt(l).n;
                if k == 1
                    for i = -n_back:n_for
                        th(j+i) = max(1,max(n_back,n_for)) / (abs(i)+1);
                    end
                else
                    for i = -n_back:n_for
                        if i == 0
                            KI = cconv2d(h_2d,I);
                            SKI = down_sample(KI,opt(l).res);
                            x_tmp = sum(sum(abs(J{j+i}-SKI))) / Nq;
                            th(j+i) = (opt(l).alpha+Nq-1)/(opt(l).beta+Nq*x_tmp);
                        else
                            KFI = cconv2d(h_2d,FI{j+i});
                            SKFI = down_sample(KFI,opt(l).res);
                            x_tmp = sum(sum(abs(J{j+i}-SKFI)))/Nq;
                            th(j+i) = (opt(l).alpha+Nq-1)/(opt(l).beta+Nq*x_tmp);
                        end
                    end
                end

                % (3) Estimate high-res img
                % IRLS
                for m = 1:opt(l).maxit
                    I_old_in = I;
                    % compute W0,Ws,Wi
                    W0 = compute_W0(I,J0,h_2d,opt(l).res);
                    Ws = compute_Ws(I);
                    if param.MOTION_ESTIMATION
                        for i = -n_back:n_for
                            if i ~= 0
                                FI{j+i} = warped_img(I,u{j+i},v{j+i});
                                Wi{j+i} = compute_Wi(FI{j+i},J{j+i},h_2d,opt(l).res);
                            end
                        end
                    end

                    % estimat I
                    AI = zeros(M,N);
                    b = zeros(M,N);
                    if param.MOTION_ESTIMATION
                        AI = compute_Ax_h(I,W0,Ws,th,h_2d,opt(l),j,n_back,n_for,FI,Wi,ut,vt,param);
                        b = compute_b_h(J,W0,th,h_2d,opt(l),j,n_back,n_for,Wi,ut,vt,param);
                        I = conj_gradient_himg(AI,I,b,W0,Ws,th,h_2d,opt(l),j,n_back,n_for,FI,Wi,ut,vt,param);
                    else
                        AI = compute_Ax(I,W0,Ws,th(j),h_2d,opt(l));
                        b = compute_b1(J0,W0,th(j),h_2d,opt(l));
                        I = conj_gradient_himg_0(AI,I,b,W0,Ws,th(j),h_2d,opt(l),param);
                    end

                    % stop criteria
                    diff_in = norm(I-I_old_in)/norm(I_old_in);
                    if diff_in < opt(l).eps
                        break;
                    end
                end

                % (4) Estimate kernel
                % IRLS
                if param.BLUR_ESTIMATION
                    K = otf2psf(psf2otf(h_2d,size(I))); % 2d kernel (Kx x Ky)
                    for m = 1:opt(l).maxit-2
                        K_old = K;
                        % compute W0,Ws,Wi
                        Wk = compute_Wk(K,I,J0,opt(l).res);
                        Ws = compute_Ws(K);
                        % estimat Kx
                        AK = compute_Ax_k(K,Wk,Ws,th(j),I,opt(l));
                        b = compute_b1_k(J0,Wk,th(j),I,opt(l));
                        K = conj_gradient_kernel(AK,K,b,Wk,Ws,th(j),I,opt(l),param);
                        % stop criteria
                        diff_K = norm(K-K_old)/norm(K_old);
                        if diff_K < opt(l).eps_blur
                            break;
                        end
                    end
                    half = floor(opt(l).hsize/2);
                    h_2d = K(M/2-half+1:M/2+half+1,N/2-half+1:N/2+half+1);
                    h_2d = h_2d / sum(sum(h_2d));
                end

                % (5) Check convergence
                diff_out = norm(I-I_old_out)/norm(I_old_out);
                if diff_out < opt(l).eps_out
                    break;
                end

            end

            I_sr{j} = I;

        end

        % Correct offset
        if opt(l).degrad == 0
            I_sr = shift_adjust(I_sr,opt(l).res,-1);
        end
        I_sr = I_sr{round((opt(l).nFrames + 1) / 2)}(:,:,1);

    end

    clear W0 Ws Wk Wi;

author: thomas-koehler

repo: SupER


text:
function [tfr,t,f] = tfrwv(x,t,N,trace);
%TFRWV	Wigner-Ville time-frequency distribution.
%	[TFR,T,F]=TFRWV(X,T,N,TRACE) computes the Wigner-Ville distribution
%	of a discrete-time signal X,
%	or the cross Wigner-Ville representation between two signals.
%
%	X     : signal if auto-WV, or [X1,X2] if cross-WV.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR   : time-frequency representation. When called without
%	        output arguments, TFRWV runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); tfrwv(sig);
%
%	See also all the time-frequency representations listed in
%	the file CONTENTS (TFR*)

%	F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least one parameter required');
end;
[xrow,xcol] = size(x);

if (nargin == 1),
 t=1:xrow; N=xrow ; trace=0;
elseif (nargin == 2),
 N=xrow ; trace=0;
elseif (nargin == 3),
 trace = 0;
end;

if (N<0),
 error('N must be greater than zero');
end;

[trow,tcol] = size(t);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

tfr= zeros (N,tcol);
if trace, disp('Wigner-Ville distribution'); end;
for icol=1:tcol,
 ti= t(icol); taumax=min([ti-1,xrow-ti,round(N/2)-1]);
 tau=-taumax:taumax; indices= rem(N+tau,N)+1;
 tfr(indices,icol) = x(ti+tau,1) .* conj(x(ti-tau,xcol));
 tau=round(N/2);
 if (ti<=xrow-tau)&(ti>=tau+1),
  tfr(tau+1,icol) = 0.5 * (x(ti+tau,1) * conj(x(ti-tau,xcol))  + ...
                           x(ti-tau,1) * conj(x(ti+tau,xcol))) ;
 end;
 if trace, disprog(icol,tcol,10); end;
end;
tfr= fft(tfr);
if (xcol==1), tfr=real(tfr); end ;

if (nargout==0),
 tfrqview(tfr,x,t,'tfrwv');
elseif (nargout==3),
 f=(0.5*(0:N-1)/N)';
end;

author: HeLiangHIT


text:
function [tag, val, typ, dep, chk, cj] = harvest(item, cj, dflag, rflag)

% function [tag, val, typ, dep, chk, cj] = harvest(item, cj, dflag, rflag)
% Generic harvest function, suitable for all const/entry items.
% The configuration tree cj is passed unmodified. If rflag is true and a
% dependency can be resolved, the resolved value will be returned,
% otherwise the cfg_dep object will be returned in val and dep.
% Input arguments:
% item  - item to be harvested
% cj    - configuration tree (passed unmodified)
% dflag - if true, harvest defaults tree, otherwise filled tree
% rflag - if true, resolve dependencies in leaf nodes
% Output arguments:
% tag - tag of harvested item
% val - harvested value
% typ - class of harvested item (currently unused)
% dep - list of unresolved dependencies
% chk - meaningful if ~dflag and all dependencies are resolved. Then it
%       returns success status of this items .check function and its
%       childrens check functions. A job is ready to run if all
%       dependencies are resolved and chk status is true.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: harvest.m 5678 2013-10-11 14:58:04Z volkmar $

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

typ = class(item);
tag = item.tag;
val = '<UNDEFINED>';
dep = []; % Placeholder for dependencies. Will be classified during
          % dependency evaluation
chk = ~dflag && rflag;

if ~isempty(item.val)
    if isa(item.val{1},'cfg_dep')
        if dflag % do not harvest references if defaults are requested
            val = '<UNDEFINED>';
        else
            if rflag
                [val, sts] = resolve_deps(item, cj);
            else
                sts = false;
            end;
            if ~sts % deps not resolved
                dep = item.val{1}; % dep is a array of cfg_dep objects
                for k = 1:numel(dep) % we may have multiple dependencies
                    dep(k).tname = item.name; % set target name
                end;
                val = dep; % return deps also in val for saving
            end;
        end;
    else
        val = item.val{1};
    end;
end;
chk = chk && isempty(dep);
if chk
    chk = docheck(item, val);
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% div = GRAPH_DIV(S, A, ki, kj)
%
%Divergence of signal residing on edges of a graph. The output resides on
%nodes. The result should be such that
%   graph_div(graph_grad(s)) = L * s,
%where L is the graph Laplacian associated with the adjacency matrix A.
%
%OUTPUT:
%   div: divergence of columns of S defined on vertices
%
%INPUTS:
%   S: signal defined on edges e
%   A: adjacency matrix (weighted or not)
%   ki, kj: indices of edges by graph_adj2vec
%
%
% TODO: does not handle normalized Laplacian case!!
%
%
%note that the norm of the operator is
%
% ||D||^2_2 = ||L||_2, where G is the divergence operator and L is the
% Laplacian used.
%
%
%see also: graph_grad, graph_adj2vec, sgwt_laplacian
%
%code author: Vassilis Kalofolias
%date: Aug 2013

function div = gsp_div_old(G,s)

if ~strcmp(G.lap_type,'combinatorial')
    error('Not implemented yet. However ask Nathanael it is very easy');
end

if size(s,1) ~= G.Ne
    error('Signal size not equal to number of edges');
end

k = size(s,2);

% What happens with loops? (diagonal elements of G, A)
%
% I think both div and grad is zero on the diagonal elements. But then is
% it consistent with the Laplacian? (it has to be div(grad(s)) = 2*L*s) for
% a signal s defined on the vertices.

div = zeros(G.N, k);

for ii = 1 : k
    tmp = sparse(G.v_in, G.v_out, s(:,ii), G.N,G.N);
    tmp = tril(tmp,1) - tril(tmp,1)';

    % the divergence
    div(:,ii) = sum(tmp.* sqrt(G.W), 2);
end

end

repo: gspbox


text:
function y= unique_unsort(x, n)
%UNIQUE - Like MATLAB's unique, but without sorting
%
%Synopsis:
%  Y= union_unsort(X, N)
%
%Description:
%  Returns the first N unique elements of X.
%
%Arguments:
%  X - can be a vector or a cell array of strings
%
%% TODO: should be written more generally for matrices

if nargin<2,
  n= inf;
end

y= x([]);  %% empty matrix or empty cell
ii= 0;
while ii<length(x) & length(y)<n,
  ii= ii + 1;
  if ~ismember(x(ii), y,'legacy'),
    y= [y, x(ii)];
  end
end

author: bbci


text:
function toms743_test03 ( nbits, xmin, xmax, n )

%*****************************************************************************80
%
%% TOMS743_TEST03 tests WAPR(X) when X is the argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Andrew Barry, S. J. Barry,
%    Patricia Culligan-Hensley.
%    This MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer NBITS, the number of bits in the mantissa.
%
%    Input, real XMIN, XMAX, the range.
%
%    Input, integer N, the number of equally spaced values
%    in the range at which arguments are to be chosen.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS743_TEST03\n' );
  fprintf ( 1, '  Input X is the argument.\n' );

  ifmt = 1;
  l = 0;

  if ( xmax < xmin )
    temp = xmin;
    xmin = xmax;
    xmax = temp;
  end

  dx = ( xmax - xmin ) / n;
  xmin = xmin - dx;

  if ( xmax <= 0.0 )
    iw = 1;
    fprintf ( 1, '  Both branches of the W function will be checked.\n' );
  else
    iw = 0;
    fprintf ( 1, '  Wp has been selected (maximum x is > 0)\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Results for Wp(x):\n' );

  if ( ifmt == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '   Offset x    W(x) (WAPR)' );
    fprintf ( 1, '     W(x) (BISECT)  Digits Correct\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '     x     W(x) (WAPR)' );
    fprintf ( 1, '     W(x) (BISECT)  Digits Correct\n' );
  end
  fprintf ( 1, '\n' );

  for i = 1 : n + 1

    x = xmin + i * dx;
    [ w, nerror ] = wapr ( x, 0, l );

    if ( nerror == 1 )

      fprintf ( 1, '  The value of X = %g is out of range.\n', x );

    else

      [ we, ner ] = bisect ( x, 0, l );

      if ( ner == 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, ' BISECT did not converge for x = %g\n', x );
        fprintf ( 1, '  Try reducing NBITS.\n' );
      end

      if ( w == we )
        nd = floor ( log10 ( 2.0 ^ nbits ) + 0.5 );
      else
        nd = floor ( log10 ( abs ( we / ( w - we ))) + 0.5 );
      end

      fprintf ( 1, '%17.8g      %17.8g      %17.8g      %3d\n', x, w, we, nd );

    end

  end

  if ( iw == 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Results for Wm(x):\n' );

    if ( ifmt == 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '   Offset x    W(x) (WAPR)' );
      fprintf ( 1, '     W(x) (BISECT)  Digits Correct\n' );
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, '     x     W(x) (WAPR)' );
      fprintf ( 1, '     W(x) (BISECT)  Digits Correct\n' );
    end

    for i = 1 : n + 1

      x = xmin + i * dx;
      [ w, nerror ] = wapr ( x, 1, l );

      if ( nerror == 1 )

        fprintf ( 1, '  The value of X = %g is out of range.\n', x );

      else

        [ we, ner ] = bisect ( x, 1, l );

        if ( ner == 1 )
          fprintf ( 1, '\n' );
          fprintf ( 1, ' BISECT did not converge for x = %g\n', x );
          fprintf ( 1, '  Try reducing NBITS.\n' );
        end

        if ( w == we )
          nd = log10 ( 2.0 ^ nbits ) + 0.5;
        else
          nd = log10 ( abs ( we / ( w - we ))) + 0.5;
        end

        fprintf ( 1, '%17.8g      %17.8g      %17.8g      %3d\n', x, w, we, nd );

      end

    end

  end

  return
end

author: johannesgerer