text:
function planMergeGui(command, varargin)
%"planMergeGui"
%   Raises a GUI to merge two plans together, usually for the purpose of
%   fusion.
%
%JRA 2/28/04
% copyright (c) 2001-2006, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).
%
%Usage:
%   Access through the CERR Gui, OR
%   load a plan and type: planMergeGui('init');

global stateS

if nargin == 0
    command = 'INIT';
end

%Find old planMerge figures.
hFig = findobj('Tag', 'CERR_PlanMergeGui');

if isempty(hFig) && ~strcmpi(command, 'init')
    warning('planMergeGui was not initialized.  Initializing now.');
    command = 'init';
end

switch upper(command)
    case 'INIT'
        %Wipe out old plan merge figures;
        if ~isempty(hFig)
            delete(hFig);
        end

        screenSize = get(0,'ScreenSize');
        units = 'pixels';
        y = 520;
        x = 520;

        %dx for each of 2 columns.
        dx = floor((x-30)/2);

        %Create figure and UI controls.
        hFig = figure('Name','Plan Merge', 'units', 'pixels', ...
            'position',[(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], ...
            'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off',...
            'Tag', 'CERR_PlanMergeGui','CloseRequestFcn','planMergeGui(''QUIT'')');

        %Status frame.
        bartop = 40;
        statusFrame             = uicontrol(hFig, 'units',units,'Position',[10 10 x-20 30], 'style', 'frame');
        frameColor              = get(statusFrame, 'BackgroundColor');
        statusLabel             = uicontrol(hFig, 'units',units,'Position',[10 25 65 15], 'style', 'text', 'String', 'Status', 'fontweight', 'bold');
        ud.handles.statusBar    = uicontrol(hFig, 'units',units,'Position',[70 14 2*dx-165 18], 'style', 'text', 'string', 'Ready.', 'enable', 'inactive', 'backgroundcolor', frameColor, 'horizontalAlignment', 'right');

        %Load plan frame and label.
        loadbot = y-65-10;
        uicontrol(hFig, 'units',units,'Position',[10 loadbot dx*2+10 65], 'style', 'frame');
        uicontrol(hFig, 'units',units,'Position',[10 loadbot+65-15 100 15], 'style', 'text', 'String', 'Select Plan', 'fontweight', 'bold', 'backgroundcolor', frameColor);

        %Doses frame and label.
        dy = floor((loadbot-bartop-30)/2);
        uicontrol(hFig, 'units',units,'Position',[10 bartop+10 dx dy], 'style', 'frame');
        uicontrol(hFig, 'units',units,'Position',[10 bartop+10+dy-15 65 15], 'style', 'text', 'String', 'Doses', 'fontweight', 'bold', 'backgroundcolor', frameColor);

        %Scans frame and label.
        uicontrol(hFig, 'units',units,'Position',[10 bartop+20+dy dx dy], 'style', 'frame');
        uicontrol(hFig, 'units',units,'Position',[10 bartop+20+2*dy-15 65 15], 'style', 'text', 'String', 'Scans', 'fontweight', 'bold', 'backgroundcolor', frameColor);

        %Structures frame and label
        uicontrol(hFig, 'units',units,'Position',[dx+20 bartop+10 dx 2*dy+10], 'style', 'frame');
        uicontrol(hFig, 'units',units,'Position',[dx+20 bartop+20+2*dy-15 90 15], 'style', 'text', 'String', 'Structures', 'fontweight', 'bold', 'backgroundcolor', frameColor);

        ud.handles.dosesSlider   = uicontrol(hFig, 'units',units,'Position',[dx-15 bartop+10+5 20 dy-10], 'style', 'slider', 'enable', 'off');
        ud.handles.scansSlider   = uicontrol(hFig, 'units',units,'Position',[dx-15 bartop+20+5+dy 20 dy-10], 'style', 'slider', 'enable', 'off', 'tag','scanSlider','callback','scanSlider_clk(''clicked'')');
        ud.handles.structsSlider = uicontrol(hFig, 'units',units,'Position',[2*dx-5 bartop+10+5 20 2*dy], 'style', 'slider', 'enable', 'off','tag','structSlider','callback','structSlider_clk(''clicked'')');

        ud.handles.Browse       = uicontrol(hFig, 'units',units, 'Position',[20  loadbot+10 dx/2-15 25], 'style', 'pushbutton', 'string', 'Browse...', 'callback', 'planMergeGui(''BROWSE'')');
        ud.handles.filenameText = uicontrol(hFig, 'units',units, 'Position',[130 loadbot+5 375 25], 'fontsize', 8, 'style', 'text', 'string', {'Click Browse to find a plan to merge with the plan currently loaded in CERR.'}, 'backgroundcolor', frameColor);
        ud.handles.pathText     = uicontrol(hFig, 'units',units, 'Position',[130 loadbot+30 375 25], 'fontsize', 8, 'style', 'text', 'string', {''}, 'backgroundcolor', frameColor);

        %Create scan UIControls.
        uicontrol(hFig, 'units',units,'Position',[12 bartop+10+2*dy-25 12 15], 'style', 'text', 'String', '#', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[25 bartop+10+2*dy-25 80 15], 'style', 'text', 'String', 'Modality', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[120 bartop+10+2*dy-25 60 15], 'style', 'text', 'String', 'Size', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[180 bartop+10+2*dy-25 50 15], 'style', 'text', 'String', 'Merge?', 'backgroundcolor', frameColor);

        uicontrol(hFig, 'units',units,'Position',[15 bartop+10+2*dy-45 80 15], 'style', 'text', 'String', 'ALL', 'backgroundcolor', frameColor, 'fontweight', 'bold');
        uicontrol(hFig, 'units',units,'Position',[100 bartop+10+2*dy-45 60 15], 'style', 'text', 'String', '', 'backgroundcolor', frameColor);
        ud.handles.scanALLCheck = uicontrol(hFig, 'units',units,'Position',[200 bartop+10+2*dy-45 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'callback', 'planMergeGui(''ALLCHECK'')', 'tag', 'scan', 'enable', 'off');

        numRows = 6;
        for i=1:numRows
            ud.handles.scanNumber(i)  = uicontrol(hFig, 'units',units,'Position',[12 bartop+10+2*dy-45-i*20 20 15], 'style', 'text', 'String', '1', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.scanName(i)    = uicontrol(hFig, 'units',units,'Position',[30 bartop+10+2*dy-45-i*20 90 15], 'style', 'text', 'String', 'CT Scan', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.scanSize(i)    = uicontrol(hFig, 'units',units,'Position',[125 bartop+10+2*dy-45-i*20 60 15], 'style', 'text', 'String', '100', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.scanCheck(i)   = uicontrol(hFig, 'units',units,'Position',[200 bartop+10+2*dy-45-i*20 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'userdata', i, 'Tag', 'ScanCheck', 'callback', 'planMergeGui(''Check'')', 'visible', 'off');
        end
        ud.numScanRows = numRows;

        %Create dose UICONTROLS.
        uicontrol(hFig, 'units',units,'Position',[15 bartop+dy-25 80 15], 'style', 'text', 'String', 'Name', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[100 bartop+dy-25 60 15], 'style', 'text', 'String', 'Size', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[165 bartop+dy-25 60 15], 'style', 'text', 'String', 'Merge?', 'backgroundcolor', frameColor);

        uicontrol(hFig, 'units',units,'Position',[15 bartop+dy-45 80 15], 'style', 'text', 'String', 'ALL', 'backgroundcolor', frameColor, 'fontweight', 'bold');
        uicontrol(hFig, 'units',units,'Position',[100 bartop+dy-45 60 15], 'style', 'text', 'String', '', 'backgroundcolor', frameColor);
        ud.handles.doseALLCheck = uicontrol(hFig, 'units',units,'Position',[185 bartop+dy-45 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'callback', 'planMergeGui(''ALLCHECK'')', 'tag', 'dose', 'enable', 'off');

        numRows = 6;
        for i=1:numRows
            ud.handles.doseName(i)  = uicontrol(hFig, 'units',units,'Position',[15 bartop+dy-45-i*20 80 15], 'style', 'text', 'String', 'CT Scan', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.doseSize(i)  = uicontrol(hFig, 'units',units,'Position',[100 bartop+dy-45-i*20 60 15], 'style', 'text', 'String', '100', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.doseCheck(i) = uicontrol(hFig, 'units',units,'Position',[185 bartop+dy-45-i*20 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'userdata', i, 'Tag', 'DoseCheck', 'callback', 'planMergeGui(''Check'')', 'visible', 'off');
        end
        ud.numDoseRows = numRows;

        %Create struct UICONTROLS.
        uicontrol(hFig, 'units',units,'Position',[dx+10+15 loadbot-46 80 15], 'style', 'text', 'String', 'Name', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[dx+10+100 loadbot-46 60 15], 'style', 'text', 'String', 'AssocScan', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[dx+10+165 loadbot-46 60 15], 'style', 'text', 'String', 'Merge?', 'backgroundcolor', frameColor);

        uicontrol(hFig, 'units',units,'Position',[dx+10+15 loadbot-66 80 15], 'style', 'text', 'String', 'ALL', 'backgroundcolor', frameColor,'fontweight', 'bold');
        uicontrol(hFig, 'units',units,'Position',[dx+10+100 loadbot-66 60 15], 'style', 'text', 'String', '', 'backgroundcolor', frameColor);
        ud.handles.strALLCheck = uicontrol(hFig, 'units',units,'Position',[dx+10+185 loadbot-66 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'callback', 'planMergeGui(''ALLCHECK'')', 'tag', 'struct', 'enable', 'off');

        numRows = 16;
        for i=1:numRows
            ud.handles.strName(i)  = uicontrol(hFig, 'units',units,'Position',[dx+10+15 loadbot-66-i*20 80 15], 'style', 'text', 'String', 'CT Scan', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.strScan(i)  = uicontrol(hFig, 'units',units,'Position',[dx+10+100 loadbot-66-i*20 60 15], 'style', 'text', 'String', '1', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.strCheck(i) = uicontrol(hFig, 'units',units,'Position',[dx+10+185 loadbot-66-i*20 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'userdata', i, 'Tag', 'StructCheck', 'callback', 'planMergeGui(''Check'')', 'visible', 'off');
        end
        ud.numStructRows = numRows;

        ud.firstVisStruct   = 1;
        ud.firstVisScan     = 1;
        ud.firstVisDose     = 1;

        ud.handles.merge = uicontrol(hFig, 'units',units, 'Position', [2*dx-85 14 100 20], 'style', 'pushbutton', 'string', 'Merge', 'callback', 'planMergeGui(''MERGE'')', 'enable', 'off');
%         uicontrol(hFig, 'units',units, 'Position', [2*dx-100 20 dx/2-15 20], 'style', 'pushbutton', 'string', 'Cancel', 'callback', 'planMergeGui(''CANCEL'')');
%
        set(hFig, 'userdata', ud);

        stateS.handle.planMergeGui = hFig;

    case 'STATUS'
        ud = get(hFig, 'userdata');
        statusString = varargin{1};
        set(ud.handles.statusBar, 'string', statusString);
        drawnow;

    case 'CANCEL'
        close(hFig);

    case 'BROWSE'
        ud = get(hFig, 'userdata');
        [fname, pathname] = uigetfile({'*.mat;*.mat.bz2', 'CERR Plans (*.mat, *.mat.bz2)';'*.*', 'All Files (*.*)'}, 'Select CERR .mat .mat.bzip2 archive to merge.');

        if fname == 0
            return;
        else
            set(ud.handles.pathText, 'string', pathname);
            set(ud.handles.filenameText, 'string', fname, 'userdata', fullfile(pathname, fname));
        end

        planMergeGui('STATUS', 'Loading selected plan and analyzing its contents...');
        planMergeGui('LOADPLAN');

    case 'LOADPLAN'
        hFig = get(gcbo, 'parent');
        ud = get(hFig, 'userdata');
        fullName = get(ud.handles.filenameText, 'userdata');

        %Get temporary directory to extract uncompress
        optS = stateS.optS;
        if isempty(optS.tmpDecompressDir)
            tmpExtractDir = tempdir;
        elseif isfolder(optS.tmpDecompressDir)
            tmpExtractDir = optS.tmpDecompressDir;
        elseif ~isfolder(optS.tmpDecompressDir)
            error('Please specify a valid directory within CERROptions.m for optS.tmpDecompressDir')
        end

        try
            planC = loadPlanC(fullName,tmpExtractDir);
            planC = updatePlanFields(planC);
            planC = quality_assure_planC(fullName,planC);
%             indexS=planC{end};
%             for i=1:length(planC{indexS.scan}(1).scanInfo)
%                 planC{indexS.scan}(1).scanInfo(i).planDir = fullName;
%             end
        end

        if ~exist('planC','var')
            planMergeGui('status', 'UNABLE TO LOAD SPECIFIED PLAN. Try another plan.');
            return;
        end
        [planInfo, planC]= getPlanInfo(planC);
        if isempty(planInfo)
            return;
        end
        ud.planInfo = planInfo;
        setappdata(hFig, 'MERGEPLAN', planC);

        indexS = planC{end};
        nStructs = length(planC{indexS.structures});
        ud.checkedStructs = ones(nStructs, 1);
        ud.nStructs = nStructs;

        nDoses = length(planC{indexS.dose});
        ud.checkedDoses = ones(nDoses, 1);

        nScans = length(planC{indexS.scan});
        ud.checkedScans = ones(nScans, 1);

        set(ud.handles.doseALLCheck, 'enable', 'on', 'value', 1);
        set(ud.handles.scanALLCheck, 'enable', 'on', 'value', 1);
        set(ud.handles.strALLCheck, 'enable', 'on', 'value', 1);

        set(hFig, 'userdata', ud);
        planMergeGui('refresh');
        planMergeGui('STATUS', 'Plan loaded.  Select elements to merge and click MERGE button.');
        set(ud.handles.merge, 'enable', 'on')

    case 'REFRESH'
        ud = get(hFig, 'userdata');
        planInfo = ud.planInfo;

        planC = getappdata(hFig, 'MERGEPLAN');
        indexS = planC{end};

        %Refresh Struct UIelements.
        set(ud.handles.strScan,  'visible', 'off');
        set(ud.handles.strName,  'visible', 'off');
        set(ud.handles.strCheck, 'visible', 'off');
        firstVisStruct = ud.firstVisStruct;
        nStructs = length(ud.checkedStructs);

        if nStructs > 16
            hSlider = findobj('tag','structSlider');
            set(hSlider,'enable', 'on','BackgroundColor',[0 0 0]);
            set(hSlider,'min',1,'max',nStructs,'value',nStructs,'sliderstep',[1/(nStructs-1) 1/(nStructs-1)]);
            structSlider_clk('init',planC);
        else
            for i=1:length(planC{indexS.structures})
                structNum  = firstVisStruct+i-1;
                structName = planC{indexS.structures}(structNum).structureName;
                assocScan  = getStructureAssociatedScan(structNum, planC);
                set(ud.handles.strName(i), 'string', structName, 'visible', 'on');
                set(ud.handles.strScan(i), 'string', num2str(assocScan), 'visible', 'on');
                set(ud.handles.strCheck(i), 'visible', 'on', 'value', ud.checkedStructs(structNum));
            end
        end

        %Refresh scan struct UIelements.
        set(ud.handles.scanSize,  'visible', 'off');
        set(ud.handles.scanName,  'visible', 'off');
        set(ud.handles.scanCheck, 'visible', 'off');
        nScans = length(ud.checkedScans);
        firstVisScan = ud.firstVisScan;

        if nScans > 6
            hSlider = findobj('tag','scanSlider');
            set(hSlider,'enable', 'on','BackgroundColor',[0 0 0]);
            set(hSlider,'min',1,'max',nScans,'value',nScans,'sliderstep',[1/(nScans-1) 1/(nScans-1)]);
            scanSlider_clk('init', planInfo, planC);
        else
            for i=1:min(length(ud.handles.scanName), nScans)
                scanNum  = firstVisScan+i-1;
                modality = planInfo.scans(scanNum).modality;
                size     = planInfo.scans(scanNum).sizeInMB;
                set(ud.handles.scanNumber(i), 'string', num2str(scanNum), 'visible', 'on');
                set(ud.handles.scanName(i), 'string', modality, 'visible', 'on');
                set(ud.handles.scanSize(i), 'string', [num2str(round(size)) ' MB'], 'visible', 'on');
                set(ud.handles.scanCheck(i), 'visible', 'on', 'value', ud.checkedScans(scanNum));
            end
        end

        %Refresh dose struct UIelements.
        set(ud.handles.doseSize,  'visible', 'off');
        set(ud.handles.doseName,  'visible', 'off');
        set(ud.handles.doseCheck, 'visible', 'off');
        nDoses = length(ud.checkedDoses);
        firstVisDose = ud.firstVisDose;
        for i=1:min(length(ud.handles.doseName), nDoses)
            doseNum  = firstVisDose+i-1;
            doseName = planC{indexS.dose}(doseNum).fractionGroupID;
            size     = planInfo.doses(doseNum).sizeInMB;
            set(ud.handles.doseName(i), 'string', doseName, 'visible', 'on');
            set(ud.handles.doseSize(i), 'string',[num2str(round(size)) ' MB'], 'visible', 'on');
            set(ud.handles.doseCheck(i), 'visible', 'on', 'value', ud.checkedDoses(doseNum));
        end
        set(hFig, 'userdata',ud);

    case 'CHECK'
        checkNum = get(gcbo, 'userdata');
        type     = get(gcbo, 'Tag');
        value    = get(gcbo, 'value');

        ud = get(hFig, 'userdata');
        totStruct = ud.nStructs;
        switch upper(type)
            case 'DOSECHECK'
                ud.checkedDoses(checkNum+ud.firstVisDose-1) = value;
            case 'SCANCHECK'
                ud.checkedScans(checkNum+ud.firstVisScan-1) = value;
            case 'STRUCTCHECK'
                slidStat = get(findobj('tag','structSlider'),'enable');
                if strcmpi(slidStat,'off')
                    ud.checkedStructs(checkNum+ud.firstVisStruct-1) = value;
                elseif strcmpi(slidStat,'on')
                    val = get(findobj('tag','structSlider'),'value');
                    if val == totStruct
                        startPt = totStruct-val;
                    else
                        startPt = totStruct-val-1;
                    end
                    ud.checkedStructs(checkNum+startPt+ud.firstVisStruct-1) = value;
                end
        end

        set(hFig, 'userdata', ud);

    case 'ALLCHECK'
        type     = get(gcbo, 'Tag');
        value    = get(gcbo, 'value');

        ud = get(hFig, 'userdata');

        switch upper(type)
            case 'DOSE'
                if ~isempty(ud.checkedDoses)
                    ud.checkedDoses(1:end) = value;
                end
            case 'SCAN'
                if ~isempty(ud.checkedScans)
                    ud.checkedScans(1:end) = value;
                end
            case 'STRUCT'
                if ~isempty(ud.checkedStructs)
                    ud.checkedStructs(1:end) = value;
                end
        end

        set(hFig, 'userdata', ud);
        planMergeGui('REFRESH');

    case 'MERGE'
        ud = get(hFig, 'userdata');
        mergePlan = getappdata(hFig, 'MERGEPLAN');

        global planC
        fileName = fullfile(get(ud.handles.pathText,'string'),get(ud.handles.filenameText,'string'));
        planMergeGui('status', 'Merging, please wait...');
        planC = planMerge(planC, mergePlan, find(ud.checkedScans), find(ud.checkedDoses), find(ud.checkedStructs), fileName);

        clear mergePlan;
        rmappdata(hFig, 'MERGEPLAN');

        planMergeGui('CANCEL');
        stateS.planMerged = 1;

        %refresh navigation montage if it exists
        navFig = findobj('tag','navigationFigure');
        if ~isempty(navFig)
            navigationMontage('refresh')
        end

        %Refresh the Viewer
        currentAxis = stateS.currentAxis;
        for iAxis = 1:length(stateS.handle.CERRAxis)
            axView = getAxisInfo(stateS.handle.CERRAxis(iAxis),'view');
            if ismember(axView,{'transverse','sagittal','coronal'})
                stateS.currentAxis = iAxis;
                sliceCallBack('CHANGESLC','PREVSLICE')
                sliceCallBack('CHANGESLC','NEXTSLICE')
            end
        end
        stateS.currentAxis = currentAxis;

    case 'QUIT'
        stateS.handle.planMergeGui = [];
        closereq;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function a = linIt(A)
% function a = linIt(A)
%	a = A(:);
	a = A(:);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function tf = prtUtilIsMethodIncludeHidden(obj,methName)
% xxx Need Help xxx
% Internal

m = metaclass(obj);
tf = ismember(methName,cellfun(@(c)c.Name,m.Methods,'uniformoutput',false));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef prtCluster < prtAction
    % prtCluster   Base class for prt clustering objects
    %
    %   prtCluster is an abstract class and cannot be instantiated.
    %
    %   All prtCluster objects inherit all properities and methods from the
    %   prtAction object.
    %
    %   prtCluster objects have the following methods:
    %
    %   plot          - Plot the output cluster of a trained clustering
    %                   algorithm. Only valid when the clutering algorithm
    %                   has been trained with a data set of 3 dimensions or
    %                   less.
    %
    %   A prtCluster object inherits the TRAIN, RUN, CROSSVALIDATE and KFOLDS
    %   functions from the prtAction object.
    %
    %   See also prtClusterGmm, prtClusterKmeans

    properties (Abstract)
        nClusters  % The number of clusters
    end
    properties (Dependent)
        internalDecider  % Optional prtDecider object for making decisions
    end
    properties (SetAccess = private, GetAccess = private, Hidden=true)
        internalDeciderDepHelper = [];
    end
    properties (Dependent)
        includesDecision % Flag indicating if result includes a decision
    end
    properties (SetAccess=protected, Hidden = true)
        yieldsMaryOutput = nan; % Determined in trainProcessing()
    end
    properties (Hidden = true)
        plotOptions = prtClass.initializePlotOptions();
    end
    properties (SetAccess = protected)
        isSupervised = false;  % False
        isCrossValidateValid = true; % True
	end

	methods (Hidden = true)
        function featureNameModificationFunction = getFeatureNameModificationFunction(self)
            if ~self.includesDecision
				featureNameModificationFunction = prtUtilFeatureNameModificationFunctionHandleCreator('%s Membership in cluster #index#', self.nameAbbreviation);
			else
				featureNameModificationFunction = prtUtilFeatureNameModificationFunctionHandleCreator('Class Label');
            end
        end
	end

    methods
        function self = prtCluster()
            % As an action subclass we must set the properties to reflect
            % our dataset requirements
            self.classTrain = 'prtDataSetClass'; % Can't this be prtDataSetStandard?
            self.classRun = 'prtDataSetStandard';
            self.classRunRetained = true;
        end

        function val = get.internalDecider(self)
            val = self.internalDeciderDepHelper;
        end
        function self = set.internalDecider(self,val)
            if ~isempty(val) && ~isa(val,'prtDecision')
                error('prtClass:internalDecider','internalDecider must be an empty vector ([]) of type prtDecision, but input is a %s',class(val));
            end

            if ~isempty(val) && ~val.isTrained && self.isTrained
                % We are adding a non-trained decision to a trained
                % classifier if we have verboseStorage = true we can train
                % the decider. If we don't we have to untrain and warn

                if ~isempty(self.dataSet)
                    % We have the dataset
                    self.internalDeciderDepHelper = val;
                    self = postTrainProcessing(self, self.dataSet);
                    self.yieldsMaryOutput = false;
                else
                    %Make self.isTrained false:
                    self.isTrained = false;
                    %warn user?
                    warning('prt:internalDecider:isTrained','Setting the internalDecider property of a prtClass object with an un-trained internalDecider sets the prtClass object''s isTrained to false');
                end
            else
                self.internalDeciderDepHelper = val;
            end

        end

        function has = get.includesDecision(self)
            has = ~isempty(self.internalDecider);
        end

        function varargout = plot(self)
            % PLOT  Plot the output of the prtCluster object
            %
            %   OBJ.plot() plots the output confidence of a prtClass
            %   object. This function only operates when the dimensionality
            %   of dataset is 3 or less. When verboseStorage is set to
            %   'true', the training data points are also displayed on the
            %   plot.
            %
            %   See also: prtClass\plotDecision

            assert(self.isTrained,'Clusterer must be trained before it can be plotted.');
            assert(self.dataSetSummary.nFeatures < 4, 'nFeatures in the training dataset must be less than or equal to 3');

            if self.yieldsMaryOutput
                self = trainAutoDecision(self);
            end

            handleStructure = plotBinaryClusterConfidence(self);

            if ~isempty(self.dataSet) && ~isempty(self.dataSet.name)
                title(sprintf('%s (%s)',self.name,self.dataSet.name));
            else
                title(self.name);
            end

            varargout = {};
            if nargout > 0
                varargout = {handleStructure};
            end
        end
    end

    methods (Access = protected, Hidden = true)

        function self = postTrainProcessing(self,ds)
             if ~isempty(self.internalDecider)
                tempself = self;
                tempself.internalDecider = [];
                yOut = tempself.run(ds);
                self.internalDecider = self.internalDecider.train(yOut);
                self.internalDecider.classList = 1:self.nClusters;
            end
        end

        function self = preTrainProcessing(self, ds)
            % Overload preTrainProcessing() so that we can determine mary
            % output status

            self.yieldsMaryOutput = ~self.includesDecision;

            self = preTrainProcessing@prtAction(self, ds);
        end

        function outputDs = postRunProcessing(self, inputDs, outputDs)
            % Overload postRunProcessing (from prtAction) so that we can
            % enforce the internaldecider

            if ~isempty(self.internalDecider)
                outputDs = self.internalDecider.run(outputDs);
            end

            outputDs = postRunProcessing@prtAction(self, inputDs, outputDs);
        end

        % Plotting functions
        function [outputDs, linGrid, gridSize] = runClustererOnGrid(self, upperBounds, lowerBounds)

            if nargin < 3 || isempty(lowerBounds)
                lowerBounds = self.dataSetSummary.lowerBounds;
            end

            if nargin < 2 || isempty(upperBounds)
                upperBounds = self.dataSetSummary.upperBounds;
            end

            [linGrid, gridSize] = prtPlotUtilGenerateGrid(upperBounds, lowerBounds, self.plotOptions.nSamplesPerDim);

            outputDs = run(self,prtDataSetClass(linGrid));
        end

        function handleStructure = plotBinaryClusterConfidence(self)

            [outputDs, linGrid, gridSize] = runClustererOnGrid(self);

            %internalDeciders* output the right colors:
            imageHandle = prtPlotUtilPlotGriddedEvaledClassifier(outputDs.getObservations(), linGrid, gridSize, prtPlotUtilLightenColors(self.plotOptions.colorsFunction(self.nClusters)),self);

            if ~isempty(self.dataSet)
                hold on;
                handles = plot(self.dataSet);
                hold off;
                handleStructure.Axes = struct('imageHandle',{imageHandle},'handles',{handles});
            else
                handleStructure.Axes = struct('imageHandle',{imageHandle},'handles',{[]});
            end
        end

        function self = trainAutoDecision(self)
            if ~isempty(self.dataSet)
                % warning('prt:prtCluster:plot:autoDecision','prtCluster.plot() requires a prtCluster with an internal decider. A prtDecisionMap has been trained and set as the internalDecider to enable plotting.');

                % It is much easier (and safer) to just train the
                % prtDecisionMap
                self.internalDecider = prtDecisionMap;
                self = postTrainProcessing(self, self.dataSet);
                self.yieldsMaryOutput = false;
            else
                % warning('prt:prtCluster:plot:autoDecision','prtCluster.plot() requires a prtCluster with an internal decider. A prtDecisionMap has been trained and set as the internalDecider to enable plotting.');
                % error('prt:prtCluster:plot','prtCluster.plot() requires a prtCluster with an internal decider. A prtDecisionMap cannot be trained and set as the internalDecider to enable plotting because this cluster object does not have verboseStorege turned on and therefore the dataSet used to train the clusterer is unknow. To enable plotting, set an internalDecider and retrain the clusterer.');

                % Since we don't have the dataset we just mock train
                % This is slightly dangerous in the case when
                % prtDecisionMap changes to include more parameters. or
                % dataSetSummary changes. (It is currently only a struct)
                tempInternalDecider = prtDecisionMap;
                tempInternalDecider.isTrained = true;
                tempInternalDecider.classList = 1:self.nClusters;
                tempInternalDecider.dataSetSummary = self.dataSetSummary;
                tempInternalDecider.dataSetSummary.nFeatures = self.nClusters;

                self.internalDecider = tempInternalDecider;
            end
        end
    end

    methods (Static, Hidden = true)
        function plotOptions = initializePlotOptions()
            if ~isdeployed
                plotOptions = prtOptionsGet('prtOptionsClusterPlot');
            else
                plotOptions = prtOptions.prtOptionsClusterPlot;
            end
        end
    end
end


text:
function p = prior_sqinvunif(varargin)
%PRIOR_SQINVUNIF  Uniform prior structure for the square inverse of the parameter
%
%  Description
%    P = PRIOR_SQINVUNIF creates uniform prior structure for the
%    square inverse of the parameter.
%
%  See also
%    PRIOR_*

% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010,2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_SQINVUNIFORM';
  ip.addOptional('p', [], @isstruct);
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'SqInv-Uniform';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'SqInv-Uniform')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  if init
    % set functions
    p.fh.pak = @prior_sqinvunif_pak;
    p.fh.unpak = @prior_sqinvunif_unpak;
    p.fh.lp = @prior_sqinvunif_lp;
    p.fh.lpg = @prior_sqinvunif_lpg;
    p.fh.recappend = @prior_sqinvunif_recappend;
  end

end

function [w, s] = prior_sqinvunif_pak(p, w)
  w=[];
  s={};
end

function [p, w] = prior_sqinvunif_unpak(p, w)
  w = w;
  p = p;
end

function lp = prior_sqinvunif_lp(x, p)
  lJ = -log(x)*3 + log(2);  % log(-2/x^3) log(|J|) of transformation
  lp = sum(lJ);
end

function lpg = prior_sqinvunif_lpg(x, p)
  lJg = -3./x;              % gradient of log(|J|) of transformation
  lpg = lJg;
end

function rec = prior_sqinvunif_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
N = 25;
x = 0:N-1;
H = zeros(N,1);

for i = 1:N
    if i == (N+1)/2
        H(i) = 1/6;
    else
        H(i) = sin(pi*(i-(N+1)/2)/6)/(pi*(i-(N+1)/2));
    end;
end;
H = H.*blackman(N);
newN = 1000;
H = [H',zeros(1,newN-length(H))];
Y = fft(H);
Y = fftshift(Y);
W = calculateDiscreteFrequencyIndex(length(Y));
subplot(2,1,1);
%plot(W,abs(Y));
plot(W,10*log10(abs(Y)));
subplot(2,1,2);
plot(W,angle(Y));

author: VipaiLab

repo: Signals-and-Systems-course


text:
function checkMeshAdjacentFaces(vertices, edges, faces)
%CHECKMESHADJACENTFACES Check if adjacent faces of a mesh have similar orientation.
%
%   checkMeshAdjacentFaces(VERTICES, EDGES, FACES)
%   The functions returns no output, but if two faces share a common edge
%   with the same direction (meaning that adjacent faces have normals in
%   opposite direction), a warning is displayed.
%
%   Example
%   [v e f] = createCube();
%   checkMeshAdjacentFaces(v, e, f);
%   % no output -> all faces have normal outwards of the cube
%
%   v = [0 0 0; 10 0 0; 0 10 0; 10 10 0];
%   e = [1 2;1 3;2 3;2 4;3 4];
%   f = [1 2 3; 2 3 4];
%   checkMeshAdjacentFaces(v, e, f);
%      Warning: Faces 1 and 2 run through the edge 3 (2-3) in the same direction
%
%   See also
%     meshes3d, trimeshMeanBreadth
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2010-10-06, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

% the message pattern that is displayed when an inconsistency is encountered
pattern = 'Faces %d and %d run through the edge %d (%d-%d) in the same direction';

% If edges are not specified, compute them
if nargin == 2
    faces = edges;
    edges = meshEdges(vertices, faces);
end

% compute edges to faces map
edgeFaces = meshEdgeFaces(vertices, edges, faces);
Ne = size(edgeFaces, 1);

for i = 1:Ne
    % indices of extreimty vertices
    v1 = edges(i, 1);
    v2 = edges(i, 2);

    % index of adjacent faces
    indF1 = edgeFaces(i, 1);
    indF2 = edgeFaces(i, 2);

    % if one of the faces has index 0, then the edge is at the boundary
    if indF1 == 0 || indF2 == 0
        continue;
    end
    % vertices of adjacent faces
    face1 = meshFace(faces, indF1);
    face2 = meshFace(faces, indF2);

    % position of vertices in face vertex array
    ind11 = find(face1 == v1);
    ind12 = find(face1 == v2);
    ind21 = find(face2 == v1);
    ind22 = find(face2 == v2);

    % check if edge is traveled forward or backard
    direct1 = (ind12 == ind11+1) | (ind12 == 1 & ind11 == length(face1));
    direct2 = (ind22 == ind21+1) | (ind22 == 1 & ind21 == length(face2));

    % adjacent faces should travel the edge in opposite direction
    if direct1 == direct2
        warning(pattern, indF1, indF2, i, v1, v2); %#ok<WNTAG>
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function V = hogDraw( H, w, fhog )
% Create visualization of hog descriptor.
%
% USAGE
%  V = hogDraw( H, [w], [fhog] )
%
% INPUTS
%  H          - [m n oBin*4] computed hog features
%  w          - [15] width for each glyph
%  fhog       - [0] if true draw features returned by fhog
%
% OUTPUTS
%  V          - [m*w n*w] visualization of hog features
%
% EXAMPLE
%
% See also hog, fhog
%
% Piotr's Image&Video Toolbox      Version 3.23
% Copyright 2013 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% fold normalizations
if(nargin<3 || isempty(fhog)), fhog=0; end
m=size(H,3); if(fhog), m=(m-5)/3; H=H(:,:,1:m*3); m=3; else m=4; end
s=size(H); s(3)=s(3)/m; w0=H; H=zeros(s);
for o=0:m-1, H=H+w0(:,:,(1:s(3))+o*s(3)); end;

% construct a "glyph" for each orientaion
if(nargin<2 || isempty(w)), w=15; end
bar=zeros(w,w); bar(:,round(.45*w):round(.55*w))=1;
bars=zeros([size(bar) s(3)]);
for o=1:s(3), bars(:,:,o)=imrotate(bar,-(o-1)*180/s(3),'crop'); end

% make pictures of positive weights by adding up weighted glyphs
H(H<0)=0; V=zeros(w*s(1:2));
for r=1:s(1), rs=(1:w)+(r-1)*w;
  for c=1:s(2), cs=(1:w)+(c-1)*w;
    for o=1:s(3), V(rs,cs)=V(rs,cs)+bars(:,:,o)*H(r,c,o); end
  end
end

author: Cloud-CV

repo: object-proposals


text:
function [net,stats] = cnn_train_dicnn_dag(net, imdb, getBatch, varargin)
%CNN_DICNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
addpath(fullfile(vl_rootnn, 'examples'));

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.epochSize = inf;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = [] ;  % Empty array means use the default SGD solver
[opts, varargin] = vl_argparse(opts, varargin) ;
if ~isempty(opts.solver)
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver; expected a function handle with two outputs.') ;
  % Call without input arguments, to get default options
  opts.solverOpts = opts.solver() ;
end

opts.momentum = 0.9 ;
opts.saveSolverState = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isscalar(opts.train) && isnumeric(opts.train) && isnan(opts.train)
  opts.train = [] ;
end
if isscalar(opts.val) && isnumeric(opts.val) && isnan(opts.val)
  opts.val = [] ;
end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.train = params.train(1:min(opts.epochSize, numel(opts.train)));
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end

  if ~isempty(opts.postEpochFn)
    if nargout(opts.postEpochFn) == 0
      opts.postEpochFn(net, params, state) ;
    else
      lr = opts.postEpochFn(net, params, state) ;
      if ~isempty(lr), opts.learningRate = lr; end
      if opts.learningRate == 0, break; end
    end
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
  state.solverState(:) = {0} ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gpuArray(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gpuArray, s, 'UniformOutput', false) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    inputs = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (s ~= 1) ;
      net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    state = accumulateGradients(net, state, params, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gather(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gather, s, 'UniformOutput', false) ;
    end
  end
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

den = params.numSubBatches * max(numGpus,1) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/den/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/den, parDer, ...
                           thisDecay, net.params(p).value) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{p} = vl_taccum(...
            params.momentum, state.solverState{p}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{p} - parDer ;
          else
            delta = state.solverState{p} ;
          end

          % Update parameters.
          net.params(p).value = vl_taccum(...
            1,  net.params(p).value, thisLR, delta) ;

        else
          % call solver function to update weights
          [net.params(p).value, state.solverState{p}] = ...
            params.solver(net.params(p).value, state.solverState{p}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end
    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
  if net.layers(sel(i)).block.ignoreAverage, continue; end;
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end

author: hbilen

repo: dynamic-image-nets


text:
function [P]=tesSmoothPosNeg(TES,V,IND_V,cPar)

%% CONTROL PARAMETERS

if isfield(cPar,'LambdaSmooth')
    LambdaSmooth=cPar.LambdaSmooth;
else
    LambdaSmooth=0.5; %DEFAULT
end

if isfield(cPar,'n')
    nMax=cPar.n;
else
    nMax=1; %DEFAULT
end

if isfield(cPar,'RigidConstraints')
    indRigid=cPar.RigidConstraints;
else
    indRigid=[]; %DEFAULT
end

if isfield(cPar,'Tolerance')
    SSQD_Tol=cPar.Tolerance;
else
    SSQD_Tol=[]; %DEFAULT
end

%%

if isempty(IND_V)
    [~,IND_V]=patchIND(TES,V);
end
logicValid=IND_V>0;

nDims=size(V,2); %Number of dimensions

VP=NaN(size(IND_V,1),size(IND_V,2),nDims);

P=V;
PP=V;
Q=V;
if ~isempty(SSQD_Tol)
    SSQD_old=[];
    SSQD_ratio=0;
end

for qIter=1:nMax;

    %% SIMPLE LAPLACIAN SMOOTHENING

    %Loop for all dimensions
    for qDim=1:1:nDims
        Xp=VP(:,:,qDim);
        Xp(logicValid)=P(IND_V(logicValid),qDim);
        Xp=gnanmean(Xp,2);
        PP(:,qDim)=Xp;
    end

    %Switch sign every iteration to partialy avoid shrinkage
    if iseven(qIter)
        wFac=1;
    else
        wFac=1;
    end

    P=P+wFac.*LambdaSmooth.*(PP-P);

    %%

    %Put back constrained points
    if ~isempty(indRigid)
       P(indRigid,:)=V(indRigid,:);
    end

    if ~isempty(SSQD_Tol)
        %Compute sum of squared differences with respect to previous iteration
        SSQD_new=gnansum((P(:)-Q(:)).^2);
        if ~isempty(SSQD_old)
            SSQD_ratio=SSQD_new./SSQD_old;
        end

        %Store current metrics
        Q=P;
        SSQD_old=SSQD_new;
        if abs(1-SSQD_ratio)<=SSQD_Tol
           break %STOP SMOOTHING LOOP IF TOLERANCE IS REACHED
        end
    end

end

%%
% _*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:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% dslam_open is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with dslam_open. If not, see <http://www.gnu.org/licenses/>.

function [num_frames, connected_ate, decentr_state] = getConnectedAte(...
    decentr_state, group)
%% Poses to positions ground truth.
p_gt_C = [];
for i = group
    p_gt_C = [p_gt_C; cell2mat(cellfun(@(x) x(1:3, 4)', ...
        decentr_state{i}.gt_T_W_C, 'UniformOutput', false))];
end
%% Poses to positions group estimate.
p_gW_C = [];
for i = group
    T_gW_C = cellfun(@(x) decentr_state{i}.Sim_W_O * x, ...
        decentr_state{i}.Sim_O_C, 'UniformOutput', false);
    p_gW_C = [p_gW_C; cell2mat(cellfun(@(x) x(1:3, 4)', ...
        T_gW_C, 'UniformOutput', false))];
end

%% Optimize: min error between positions, func of T_gt_O_ate
[decentr_state{min(group)}.T_gt_O_ate, connected_ate] = ...
    alignTrajs(p_gt_C', p_gW_C', decentr_state{min(group)}.T_gt_O_ate);
num_frames = size(p_gW_C, 1);

end

author: uzh-rpg


text:
function randlc_test01 ( )

%*****************************************************************************80
%
%% RANDLC_TEST01 tests RANDLC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2010
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDLC_TEST01\n' );
  fprintf ( 1, '  RANDLC computes pseudorandom values\n' );
  fprintf ( 1, '  in the interval [0,1].\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         I          RANDLC\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 10
    [ x, seed ] = randlc ( seed );
    fprintf ( 1, '  %8d  %14f\n', i, x );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Nonnegative matrix factorization
% Argyris Zymnis, Joelle Skaf, Stephen Boyd
%
% We are given a matrix A in R^{m*n}
% and are interested in solving the problem:
%
% minimize    ||A - Y*X||_F
% subject to  Y >= 0, X >= 0
%
% where Y in R{m*k} and X in R{k*n}.
% This script generates a random matrix A and obtains an
% *approximate* solution to the above problem by first generating
% a random initial guess for Y and the alternatively minimizing
% over X and Y for a fixed number of iterations.

% Generate data matrix A
rstate = rand('state');
m = 10; n = 10; k = 5;
A = rand(m,k)*rand(k,n);

% Initialize Y randomly
Y = rand(m,k);

% Perform alternating minimization
MAX_ITERS = 30;
residual = zeros(1,MAX_ITERS);
for iter = 1:MAX_ITERS
    cvx_begin quiet
        if mod(iter,2) == 1
            variable X(k,n)
            X >= 0;
        else
            variable Y(m,k)
            Y >= 0;
        end
        minimize(norm(A - Y*X,'fro'));
    cvx_end
    fprintf(1,'Iteration %d, residual norm %g\n',iter,cvx_optval);
    residual(iter) = cvx_optval;
end

% Plot residuals
plot(residual);
xlabel('Iteration Number');
ylabel('Residual Norm');

% Display results
disp( 'Original matrix:' );
disp( A );
disp( 'Left factor Y:' );
disp( Y );
disp( 'Right factor X:' );
disp( X );
disp( 'Residual A - Y * X:' );
disp( A - Y * X );
fprintf( 'Residual after %d iterations: %g\n', iter, cvx_optval );

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [sys,x0,str,ts]=ADRC_2(t,x,u,flag,h,TD,ESO,NLSEF,b0)

switch flag
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes(h);
    case 2
        sys=mdlUpdate(x,u,h,TD,ESO,b0);
    case 3
        sys=mdlOutputs(x,NLSEF,b0);
    case 4,
        sys=mdlGetTimeOfNextVarHit(t,h);
    case {1,9}
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes(h)
    sizes=simsizes;
    sizes.NumContStates=0;
    sizes.NumDiscStates=5;
    sizes.NumOutputs=2;
    sizes.NumInputs=3;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0;0;0;0];
    str=[];
    ts=[h 0];
function sys=mdlUpdate(x,u,h,TD,ESO,b0)
    e1=x(1)-u(1);
    fh=fhan(e1,x(2),TD(1),TD(2));
    sys(1)=x(1)+h*x(2);
    sys(2)=x(2)+h*fh;
    e2=x(3)-u(2);
    sys(3)=x(3)+h*(x(4)-ESO(1)*e2);
    sys(4)=x(4)+h*(x(5)-ESO(2)*fal(e2,0.5,ESO(4))+b0*u(3));
    sys(5)=x(5)+h*(-ESO(3)*fal(e2,0.25,ESO(4)));
%     fe=fal(e,0.5,Delta);
%     fe1=fal(e,0.25,Delta);
%     sys(1,1)=x(1)+h2*(x(2)-BB(1)*e);
%     sys(2,1)=x(2)+h2*(x(3)-BB(2)*fe+u(1));
%     sys(3,1)=x(3)+h2*(-BB(3)*fe1);
function sys=mdlOutputs(x,NLSEF,b0)
    e3=x(1)-x(3);
    e4=x(2)-x(4);
    sys(1)=-fhan(e3,NLSEF(1)*e4,NLSEF(2),NLSEF(3))-x(5)/b0;
    sys(2)=x(1);
function sys=mdlGetTimeOfNextVarHit(t,h)
    sys=t+h;

function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end

function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function an=ea_bounding(phi,theta,psi)
% bounding of euler angles
% 0<=phi<2*pi
% 0<=theta<pi
% 0<=psi<2*pi

ifchange=false; % if anles come out from boundaries

if (0<=phi)&&(phi<2*pi)

else
    ifchange=true;
end

if (0<=theta)&&(theta<pi)

else
    ifchange=true;
end

if (0<=psi)&&(psi<2*pi)

else
    ifchange=true;
end

%theta=mod(theta+pi,2*pi)-pi; % -pi<=theta<pi
theta=-(mod(theta*(-1)+pi,2*pi)-pi); % -pi<theta<=pi

if theta<0
    theta=-theta;
    phi=pi+phi;
    psi=pi+psi;
end

phi=mod(phi,2*pi);
psi=mod(psi,2*pi);

an=[ifchange phi theta psi];

author: Sable

repo: mcbench-benchmarks


text:
function hRD = plugin_rd_rxtx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.adrv9371x.common.plugin_rd('ZCU102', 'Rx & Tx');
AnalogDevices.adrv9371x.zcu102.add_io(hRD, 'Rx & Tx');

author: analogdevicesinc


text:
function [ x, y, z, w ] = ld0230 ( )

%*****************************************************************************80
%
%% LD0230 computes the 230 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(230,1);
  y = zeros(230,1);
  z = zeros(230,1);
  w = zeros(230,1);
  a = 0.0;
  b = 0.0;
  v = -0.5522639919727325E-01;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.4450274607445226E-02;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.4492044687397611;
  v = 0.4496841067921404E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2520419490210201;
  v = 0.5049153450478750E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6981906658447242;
  v = 0.3976408018051883E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6587405243460960;
  v = 0.4401400650381014E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4038544050097660E-01;
  v = 0.1724544350544401E-01;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5823842309715585;
  v = 0.4231083095357343E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.3545877390518688;
  v = 0.5198069864064399E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.2272181808998187;
  b = 0.4864661535886647;
  v = 0.4695720972568883E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end


text:
function cc_level_compose_animate ( )

%*****************************************************************************80
%
%% CC_LEVEL_COMPOSE_ANIMATE displays the grids that compose one level.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVEL_COMPOSE_ANIMATE:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display the nested 2D Clenshaw-Curtis grids\n' );
  fprintf ( 1, '  that compose one level.\n' );

  dim_num = 2;

  while ( 1 )
%
%  Get user input.
%
    fprintf ( 1, '\n' );
    level = input ( 'Enter the LEVEL or RETURN to exit;' );

    if ( isempty ( level ) )
      break
    end
%
%  Generate the entire set of points.
%
    [ grid_num, point_num ] = cc_levels_minmax_size ( dim_num, ...
      level, level );

    [ grid_level, grid_order, grid_points ] = cc_levels_minmax ( dim_num, ...
      level, level,  grid_num, point_num );

    [ dim_num, grid_points_num ] = size ( grid_points );
%
%  Display the full set of points as filled blue circles.
%
    clf
    axes_handle = axes;
    handle_new = scatter ( grid_points(1,:), grid_points(2,:), 'b', 'filled' );
    axis square
    grid on
    set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'xticklabel', [] );
    set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'yticklabel', [] );
    axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
    s = sprintf ( 'Entire grid for LEVEL %d', level );
    title ( s );
    fprintf ( 1, 'Press return\n', level );
    pause

    hold off
%
%  Display the full set of points in gray.
%  Display each contributing grid, one at a time, in red.
%
    for gridd = 1 : grid_num

      order_1d(1:2) = grid_order(1:2,gridd)
      order_nd = prod ( order_1d(1:2) );
      grid_points_new = cc_grid ( dim_num, order_1d, order_nd );

      clf
%
%  We have to name the axes in order to control the grid.
%
      axes_handle = axes;
%
%  Plot the OLD points.
%
      handle_old = scatter ( grid_points(1,:), grid_points(2,:), 'bo' );

      hold on

      handle_new = scatter ( grid_points_new(1,:), grid_points_new(2,:), ...
        'r', 'filled' );
%
%  Force the plotting region to be square, not rectangular.
%
      axis square
%
%  Request grid lines.
%
      grid on
%
%  Specify the location of the grid lines, and suppress labeling.
%
      set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
      set ( axes_handle, 'xticklabel', [] );
      set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
      set ( axes_handle, 'yticklabel', [] );
%
%  Make the plotting region slightly bigger than the data.
%
      axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
%
%  Title
%
      s = sprintf ( '+ %d * %d CC grid', grid_level(1,gridd), grid_level(2,gridd) );
      title ( s );

      fprintf ( 1, '+ %d*%d CC grid, Press return\n', ...
        grid_level(1,gridd), grid_level(2,gridd) );
      pause

      hold off

    end
%
%  Display the full set of points as filled blue circles.
%
    clf
    axes_handle = axes;
    handle_new = scatter ( grid_points(1,:), grid_points(2,:), 'b', 'filled' );
    axis square
    grid on
    set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'xticklabel', [] );
    set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'yticklabel', [] );
    axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
    s = sprintf ( 'Entire grid for LEVEL %d', level );
    title ( s );
    fprintf ( 1, 'Press return\n', level );
    pause

    hold off

    fprintf ( 1, 'Press return to CLEAR the grid!\n' );
    pause
    clf

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVEL_COMPOSE_ANIMATE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

sigmoid = sigmoid(z);
g = sigmoid .* (1 .- sigmoid);

% =============================================================

end

repo: MachineLearning


text:
function polynomial_compress_test ( )

%*****************************************************************************80
%
%% POLYNOMIAL_COMPRESS_TEST tests POLYNOMIAL_COMPRESS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  o = 10;

  c = [ 7.0, - 5.0, 5.0, 9.0, 11.0, 3.0, 6.0, 0.0, - 13.0, 1.0E-20 ];
  e = [ 1, 2, 2, 4, 5, 5, 5, 12, 33, 35 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYNOMIAL_PRINT_TEST\n' );
  fprintf ( 1, '  POLYNOMIAL_PRINT prints a polynomial.\n' );

  fprintf ( 1, '\n' );
  title = '  Uncompressed P(X) =';
  polynomial_print ( m, o, c, e, title );

  [ o2, c2, e2 ] = polynomial_compress ( o, c, e );

  fprintf ( 1, '\n' );
  title = '  Compressed P(X) =';
  polynomial_print ( m, o2, c2, e2, title );

  return
end


text:
function [y, m ] = month_borrow_common ( y, m )

%*****************************************************************************80
%
%% MONTH_BORROW_COMMON borrows a year of months on the Common calendar.
%
%  Discussion:
%
%    If the month index is legal, nothing is done.  If the month index
%    is too small, then one or more years are "cashed in" to bring the
%    month index up to a legal value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the YM date.
%    On input, M might be negative.  On output, Y should have decreased by
%    one, and M gone up by the number of months in the year that we
%    "cashed in".  The routine knows there was no year 0.
%
  while ( m <= 0 )

    months = year_length_months_common ( y );

    m = m + months;
    y = y - 1;

    if ( y == 0 )
      y = - 1;
    end

  end

  return
end


text:
classdef NegativeBinomialD
%%NEGATIVEBINOMIALD Functions to handle the negative binomial
%           distribution. The Polya distribution and the Pascal
%           distribution are special cases of this distribution. The
%           negative binomial distribution is the number of sucesses in
%           independent trials before a specified number of failures
%           occurs.
%Implemented methods are: mean, var, PMF, CDF, rand
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(r,p)
%%MEAN Obtain the mean of the negative binomial distribution.
%
%INPUTS: r This is the integer number of failed trials until sampling stops
%          (defining the distribution).
%        p This is the probability of success for each trial.
%
%OUTPUTS: val The mean of the distribution.
%
%The negative binomial distribution is discussed in Chapter 4 of [1].
%
%EXAMPLE:
%We verify the computed mean by comparing it to the sample mean.
% p=0.8;
% r=4;
% meanVal=NegativeBinomialD.mean(r,p)
% numSamp=1e5;
% sampMeanVal=mean(NegativeBinomialD.rand([1,numSamp],r,p))
%One will find both values are about 16.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=p*r/(1-p);
end

function val=var(r,p)
%%VAR Obtain the variance of the negative binomial distribution.
%
%INPUTS: r This is the integer number of failed trials until sampling stops
%          (defining the distribution).
%        p This is the probability of success for each trial.
%
%OUTPUTS: val The variance of the distribution.
%
%The negative binomial distribution is discussed in Chapter 4 of [1].
%
%EXAMPLE:
%We verify the computed variance by comparing it to the sample variance.
% p=0.8;
% r=4;
% varVal=NegativeBinomialD.var(r,p)
% numSamp=1e5;
% sampVarVal=var(NegativeBinomialD.rand([1,numSamp],r,p))
%One will find both values are about 80.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=p*r/(1-p)^2;
end

function val=PMF(k,r,p)
%%PMF Evaluate the probability mass function (PMF) of the negative binomial
%     distribution at given points.
%
%INPUTS: x The point(s) at which the negative binomial PMF is to be
%          evaluated. x is an integer. For nonzero PMF values, x>=r.
%        r This is the integer number of failed trials until sampling stops
%          (defining the distribution).
%        p This is the probability of success for each trial.
%
%OUTPUTS: val The value(s) of the negative binomial PMF.
%
%The negative binomial distribution is discussed in Chapter 4 of [1].

%EXAMPLE:
%Here, we validate the PMF by generating random samples and comparing the
%PMF plot with a histogram of the random samples.
% p=0.8;
% r=4;
% numSamples=1e5;
% figure(1)
% clf
% histogram(NegativeBinomialD.rand([numSamples,1],r,p),'BinWidth',1,'Normalization','pdf')
% hold on
% x=0:80;
% vals=NegativeBinomialD.PMF(x,r,p);
% stem(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(k));
    numEls=numel(val);

    constVal=(1-p)^r;
    for curEl=1:numEls
        if(k(curEl)>=0)
            val(curEl)=binomial(k(curEl)+r-1,k(curEl))*constVal*p^k(curEl);
        end
    end
end

function val=CDF(k,r,p)
%%CDF Evaluate the cumulative distribution function (CDF) of the negative
%     binomial distribution at given points.
%
%INPUTS: x The point(s) at which the negative binomial CDF is to be
%          evaluated. x is an integer. For nonzero CDF values, x>=r.
%        r This is the integer number of failed trials until sampling stops
%          (defining the distribution).
%        p This is the probability of success for each trial.
%
%OUTPUTS: val The value(s) of the negative binomial CDF.
%
%The negative binomial distribution is discussed in Chapter 4 of [1]. The
%CDF can be expressed in terms of a regularized incomplete beta function.
%
%EXAMPLE:
%We validate the CDF value by comparing it to a value computed from random
%samples.
% p=0.8;
% r=5;
% x=16;
% numSamples=1e5;
% prob=NegativeBinomialD.CDF(x,r,p)
% probSamp=mean(NegativeBinomialD.rand([numSamples,1],r,p)<=x)
%One will find the values ot both be about 0.414.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=1-betainc(p,k+1,r);
end

function vals=rand(N,r,p)
%%RAND Generate negative binomial random variables.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%        r This is the integer number of failed trials until sampling
%          stops (defining the distribution); r>0.
%        p This is the probability of success for each trial.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated negative binomial random variables.
%
%The negative binomial distribution arises from a Poisson distribution
%whose mean is a random variable having a central gamma distribution. This
%relationship is briefly mentioned in [1] and can be easily derived. Thus,
%this function generates gamma random variables and then generates poisson
%random variables.
%
%REFERENCES:
%[1] M. H. Quenouille, "A relation between the logarithmic, Poisson, and
%    negative binomial distributions," Biometrics, vol. 5, no. 2, pp.
%    162-164, Jun. 1949.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N, N];
    else
        dims=N;
    end

    lambda=GammaD.rand(dims,r,p/(1-p));

    numEls=numel(lambda);

    vals=zeros(size(lambda));
    for curEl=1:numEls
        vals(curEl)=PoissonD.rand(1,lambda(curEl));
    end
end
end
end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [gx,dgdx,dgdp] = g_classif(x,P,u,in)

gx = sss([in.X',ones(size(in.X,2),1)]*P);
dgdx = [];
dgdp = diag(gx.*(1-gx))*[in.X',ones(size(in.X,2),1)];
dgdp = dgdp';

function sx = sss(x)
sx = 1./(1+exp(-x));
sx(sx < 1e-8) = 1e-8;
sx(sx > 1-1e-8) = 1-1e-8;

author: MBB-team

repo: VBA-toolbox


text:
function [Timg, tThreshFDR, n_signif,index_signif, pvals] = dtiTTestImage(Tvalues, DISTR, df, mask)
%
% Basic interpretations of the dtiTTest results to produce images and other
% values related to false discovery rate
%
%  [Timg, tThreshFDR, n_signif,index_signif, pvals] = ...
%         dtiTTestImage(Tvalues, DISTR, df, mask)
%

% Create an image of the results
Timg = dtiIndToImg(Tvalues, mask, NaN);

% When you view the montage, this sets how many slices we use.
% showSlices = [20:60];
% showSlices = [10:50];
% All the brain containing slices
% showSlices = [25:50];  % Small number for debugging
% showSlices = [15:62];
% figure; imagesc(makeMontage(Timg,showSlices)); axis image; colormap cool; colorbar;
% set(gcf,'Name','FA test'); title(sprintf('tthresh, no FDR (p<10^-^4) = %0.1f',tThresh));

% Figure out the statistical significance now.
tThresh = tinv(1-10^-4, df(1));
tMax = tinv(1-10^-12, df(1));
Timg(abs(Timg)>tMax) = tMax;
tMax = max(Timg(:));

% Perform an FDR analysis for the FA test
%
fdrVal = 0.05;   % This is the p-value we are using.
fdrType = 'general';
Tvalues(isnan(Tvalues)) = 0;
pvals = 1-tcdf(Tvalues, df(1));
[n_signif,index_signif] = fdr(pvals,fdrVal,fdrType,'mean');

% Convert back to an fThreshold.  Needs more comments.  We think that this
% function returns the t-value needed to achieve a significance fdrVal, say
% 0.05 or 0.01.  It is possible that tThreshFDR/tMax is the fThreshold,
% though the comment in the printf doesn't say that.  It just puts the
% value in a parenthesis.
if n_signif > 0
    tThreshFDR = tinv(1-max(pvals(index_signif)), df(1));
    fprintf('t-threshold for FDR (%s) when p < %0.3f: %0.2f (%0.3f).\n',...
        fdrType,fdrVal,tThreshFDR,tThreshFDR/tMax);
else
    tThreshFDR = [];
    fprintf('Nothing returned as significant');
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function UFpages (figures, list)
%UFPAGES create web page for each matrix in UF Sparse Matrix Collection
% Usage: UFpages (figures, list)
%
% figures: 1 if figures are to be created, 0 otherwise
% list: list of matrix id's to process.  All matrices in the collection are
%   processed if not present.
%
% Example:
%
%   UFpages             % create all the pages
%   UFpages (0)         % create all the pages, but not the figures
%   UFpages (1,1:10)    % create pages for just matrices 1 to 10
%   UFpages (0,1:10)    % ditto, but do not create the figures
%
% See also UFpage, UFget.

% Copyright 2006-2007, Timothy A. Davis

if (nargin < 1)
    figures = 1 ;
end
index = UFget ;
if (nargin < 2)
    list = 1:length (index.nrows) ;
end
for i = list
    UFpage (i, index, figures) ;
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Test different Hawkes process simulation methods.
% Given synthetic data generated via different simulation methods, we
% estimate the parameters via MLE, and plot the relative estimation errors
% w.r.t. the number of training sequences.
%
% Provider:
% Hongteng Xu @ Georgia Tech
% June 14, 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

options.N = 50; % the number of sequences
options.Nmax = 100; % the maximum number of events per sequence
options.Tmax = 100; % the maximum size of time window
options.tstep = 0.2;% the step length for computing sup intensity
options.M = 50; % the number of steps
options.GenerationNum = 5; % the number of generations
D = 4; % the dimension of Hawkes processes
nTest = 5;
nSeg = 5;
nNum = options.N/nSeg;

disp('Fast simulation of Hawkes processes with exponential kernel')
para1.mu = rand(D,1)/D;
para1.A = rand(D, D);
para1.A = 0.25 * para1.A./max(abs(eig(para1.A)));
para1.A = reshape(para1.A, [D, 1, D]);
para1.w = 1;
Seqs1 = SimulationFast_Thinning_ExpHP(para1, options);

disp('Thinning-based simulation of Hawkes processes with exponential kernel')
para2 = para1;
para2.kernel = 'exp';
para2.landmark = 0;
Seqs2 = Simulation_Thinning_HP(para2, options);

disp('Approximate simulation of Hawkes processes via branching process')
para3 = para2;
Seqs3 = Simulation_Branch_HP(para3, options);

disp('Learning Hawkes processes from synthetic data')
alg.LowRank = 0;
alg.Sparse = 0;
alg.GroupSparse = 0;
alg.outer = 4;
alg.rho = 0.1;
alg.inner = 4;
alg.thres = 1e-5;
alg.Tmax = [];
alg.storeErr = 0;
alg.storeLL = 0;

disp('Evaluation of quality of synthetic data.')

Err = zeros(3,nSeg,nTest);
for n = 1:nTest
    for i = 1:nSeg
        % initialize
        model.A = rand(D,1,D)./(D^2);
        model.mu = rand(D,1)./D;
        model.kernel = 'exp';
        model.w = 1;
        model.landmark = 0;

        model1 = model;
        model2 = model;
        model3 = model;
        model1 = Learning_MLE_Basis( Seqs1(1:i*nNum), model1, alg );
        model2 = Learning_MLE_Basis( Seqs2(1:i*nNum), model2, alg );
        model3 = Learning_MLE_Basis( Seqs3(1:i*nNum), model3, alg );

        Err(1,i,n) = norm([model1.mu; model1.A(:)] - [para1.mu; para1.A(:)])/...
            norm([para1.mu; para1.A(:)]);
        Err(2,i,n) = norm([model2.mu; model2.A(:)] - [para2.mu; para2.A(:)])/...
            norm([para2.mu; para2.A(:)]);
        Err(3,i,n) = norm([model3.mu; model3.A(:)] - [para3.mu; para3.A(:)])/...
            norm([para3.mu; para3.A(:)]);
    end
end

Error = mean(Err,3);
Std = std(Err, 0, 3);
figure
hold on
for i = 1:3
    errorbar(nNum:nNum:options.N, Error(i,:), Std(i,:), 'o-');
end
hold off
axis tight
xlabel('The number of training sequences');
ylabel('Relative estimation error')
legend('FastThinning', 'Thinning', 'Branching')
title('Learning results based on different simulation methods')

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
%GETOBJSIZE Get object size
%
%    OBJSIZE = GETOBJSIZE(A,N)
%    [OBJSIZE1,OBJSIZE2,OBJSIZE3] = GETOBJSIZE(A)
%
% Note that the object size can be a scalar, or, in case of feature images,
% an image size vector.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% Optimal doping profile optimization
% Boyd, Kim, Vandenberghe, and Hassibi, "A tutorial on geometric programming"
% Joshi, Boyd, and Dutton, "Optimal doping profiles via geometric programming"
% Written for CVX by Almir Mutapcic 02/08/06
% (a figure is generated)
%
% Determines the optimal doping profile that minimizes base transit
% time in a (homojunction) bipolar junction transistor.
% This problem can be posed as a GP:
%
%   minimize   tau_B
%       s.t.   Nmin <= v <= Nmax
%              y_(i+1) + v_i^const1 <= y_i
%              w_(i+1) + v_i^const2 <= w_i, etc...
%
% where variables are v_i, y_i, and w_i.

% discretization size
M = 50;
% M = 1000; % takes a few minutes to process constraints

% problem constants
g1 = 0.42;
g2 = 0.69;
Nmax = 5*10^18;
Nmin = 5*10^16;
Nref = 10^17;
Dn0 = 20.72;
ni0 = 1.4*(10^10);
WB = 10^(-5);
C =  WB^2/((M^2)*(Nref^g1)*Dn0);

% exponent powers
pwi = g2 -1;
pwj = 1+g1-g2;

% optimization variables
cvx_begin gp
  variables v(M) y(M) w(M)

  % objective function is the base transmit time
  tau_B = C*w(1);

  minimize( tau_B )
  subject to
    % problem constraints
    v >= Nmin;
    v <= Nmax;

    for i = 1:M-1
      if( mod(i,100) == 0 ), fprintf(1,'progress counter: %d\n',i), end;
      y(i+1) + v(i)^pwj <= y(i);
      w(i+1) + y(i)*v(i)^pwi <= w(i);
    end

    y(M) == v(M)^pwj;
    w(M) == y(M)*v(M)^pwi;
cvx_end

% plot the basic optimal doping profile
figure, clf
nbw = 0:1/M:1-1/M;
semilogy(nbw,v,'LineWidth',2);
axis([0 1 1e16 1e19]);
xlabel('base');
ylabel('doping');
text(0,Nmin,'Nmin ', 'HorizontalAlignment','right');
text(0,Nmax,'Nmax ', 'HorizontalAlignment','right');
disp('Optimal doping profile is plotted.')


text:
function indx = multigrid_index0 ( dim_num, order_1d, order_nd )

%*****************************************************************************80
%
%% MULTIGRID_INDEX0 returns an indexed multidimensional grid.
%
%  Discussion:
%
%    For dimension DIM, the second index of INDX may vary from
%    0 to ORDER_1D(DIM)-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the points.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the
%    rule in each dimension.
%
%    Input, integer ORDER_ND, the product of the entries of ORDER_1D.
%
%    Output, integer INDX(DIM_NUM,ORDER_ND), the indices of the points in
%    the grid.  The second dimension of this array is equal to the
%    product of the entries of ORDER_1D.
%
  p = 0;
  indx = zeros(dim_num,order_nd);

  a = [];
  more = 0;

  while ( 1 )

    [ a, more ] = vec_colex_next2 ( dim_num, order_1d, a, more );

    if ( ~more )
      break
    end

    p = p + 1;

    indx(1:dim_num,p) = a(1:dim_num)';

  end

  return
end


text:
function [xUpdate,LUpdate,PUpdate,innov,Pzz,W]=separatedCovUpdateWithPred(z,R,zPred,PzPred,otherInfo,LPred,c)
%%SEPARATEDCOVUPDATEWITHPRED Given the output of the measurement prediction
%           step from separatedCovMeasPred and a measurement, complete the
%           measurement update step of the separated covariance filter.
%           Separating the measurement prediction step from the rest of the
%           update step can make the creation of multiple measurement
%           association hypotheses from a single target prediction more
%           efficient. The full measurement update function is
%           separatedCovUpdate.
%
%INPUTS: z The zDimX1 vector measurement.
%        R The zDimXzDim measurement covariance matrix in the native
%          coordinate system of the measurement.
%    zPred The zDimXnumComp measurement predictions from the filter.
%   PzPred The zDimXzDimXnumComp covariance matrices associated with zPred.
% otherInfo The intermediate results returned in the otherInfo output of
%          the separatedCovMeasPred function.
%    LPred The xDimXzDimXnumComp predicted delay vectors (defined before
%          Equation 12 in [1]). The use of multiple columns represents a
%          choice in how the algorithm was generalized to multiple
%          dimensions.
%        c The confidence region under consideration by the filter. 0<c<1.
%          If this parameter is omitted or an empty matrix is passed, the
%          default value of c=0.99 is used.
%
%OUTPUTS: xUpdate The xDimXnumComp updated state vectors.
%         LUpdate The xDimXzDimXnumComp updated delay vectors.
%         PUpdate The xDimXxDimXnumComp covariance matrices of the state
%                 estimates. This is a combination of LUpdate and PUpdate.
%      innov, Pzz The zDimXnumComp innovations and a zDimXzDimXnumComp
%                 set of matrices Pzz that are akin to an innovation
%                 covariance matrices are returned in case one wishes to
%                 analyze the consistency of the estimator or use those
%                 values in gating or likelihood evaluation.
%               W The xDimXzDimXnumComp gains used in the update.
%
%See the comments to the function separatedCovMeasPred for an example of
%usage of this function. See the comments to separatedCovUpdate for more
%information on the algorithm.
%
%REFERENCES:
%[1] G. J. Portmann, J. R. Moore, and W. G. Bath, "Separated covariance
%    filtering," in Proceedings of the IEEE International Radar Conference,
%    Arlington, VA, 7-10 May 1990, pp. 456-460.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<7||isempty(c))
    c=0.99;
end

H=otherInfo.H;
TPred=otherInfo.TPred;
xPred=otherInfo.xPred;
Pxz=otherInfo.Pxz;

xDim=size(xPred,1);
numComp=size(xPred,2);
zDim=size(z,1);

xUpdate=zeros(xDim,numComp);
LUpdate=zeros(xDim,zDim,numComp);
PUpdate=zeros(xDim,xDim,numComp);
innov=zeros(zDim,numComp);
Pzz=zeros(zDim,zDim,numComp);
W=zeros(xDim,zDim,numComp);
for k=1:numComp
    Pzz(:,:,k)=PzPred(:,:,k)+c^2*R;
    W(:,:,k)=Pxz(:,:,k)/Pzz(:,:,k);%The gain

    innov(:,k)=z-zPred(:,k);%The innovation
    xUpdate(:,k)=xPred(:,k)+W(:,:,k)*innov(:,k);

    xDim=size(xPred,1);
    diff=eye(xDim,xDim)-W(:,:,k)*H;

    LUpdate(:,:,k)=diff*LPred(:,:,k);
    TUpdate=diff*TPred(:,:,k);

    PUpdate(:,:,k)=(TUpdate-LUpdate(:,:,k)*LUpdate(:,:,k)')/c^2;
end
end

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


text:
function [GHz] = Hz2GHz(Hz)
% Convert frequency from hertz to gigahertz.
% Chad A. Greene 2012
GHz = Hz*1e-9;


text:
%% background model around vehicle
function [bgr, stack]  = mdl(stack, pts, mat, st)                          % background model

%% fill stack: shift previous + insert new
stack         = stck(stack, pts, mat, st);                                 % transform previous mats and put new mat in stack
%% background modeling
bgr.mat       = zeros(st.vx.ix, st.vx.iy);
bgr.con       = zeros(st.vx.ix, st.vx.iy);                                 % confidence
for i         = 1 : size(stack.mat, 1)
   for j      = 1 : size(stack.mat, 2)
      if sum(stack.ind(i, j, :) ~= 0)                                      % if that cell has any valid data on it
      ind.val = stack.ind(i, j, :) ~= 0;                                   % find indexes of those valid cells
      if sum(ind.val) > st.stc.in                                          % if it is more than st.stc.in, just take first st.stc.in measurements
      ind.ind = find(ind.val == 1);                                        % keep first st.stc.in indexes
      ind.val(ind.ind(st.stc.in) + 1 : end) = 0;                           % to average over first st.stc.in filled stack's cells
      end
      bgr.mat(i, j) = mean(stack.mat(i, j, ind.val ~= 0));                 % average over all
      bgr.con(i, j) = sum(double(ind.val));                                % confidence (number of measurements)
      end
   end
end
bgr.mat       = (bgr.con >= 5) .* bgr.mat;                                 % add to option! cell must have been observed for minimum number of 5
%% matrix to vector (bgr.pts, bgr.ptn)
ivx           = 1 : size(bgr.mat, 1);                                      % vxl pts at the origin
ivy           = 1 : size(bgr.mat, 2);                                      % index vector
[iy, ix]      = meshgrid(ivy, ivx);                                        % notice! ivy ivx
bgr.pts       = [(ix(:) - 1) * st.vx.x + st.vm.xb, ...                     % map vector at origin
                 (iy(:) - 1) * st.vx.y + st.vm.yr, bgr.mat(:)];
ptm           = pts.rtn * bgr.pts(:, 1:3)';                                % vxl pts in the world coordinates
bgr.ptn       = [ptm(1,:)' + pts.trn(1), ptm(2,:)' + pts.trn(2), ...
                 bgr.pts(:, 3)];                                           % map vector at real coord

end

%% fill stack: shift previous + insert new
function stack  = stck(stack, pts, mat, st)            % fill stack

%% stack: shift
for i     = st.stc.sz - 1 : -1 : 1       % stack process from bottom to top
bef.mat   = stack.mat(:, :, i);          % load previous mat
bef.pts   = stack.pts(:, :, i);
bef.ptn   = stack.ptn(:, :, i);
aftm      = trns(bef, pts, st);          % transformation: mat [frm - 1], pts [frm]
bef.mat   = stack.ind(:, :, i);          % put index instead of mat
afti      = trns(bef, pts, st);          % do the same for index
stack.mat(:, :, i + 1) = aftm.mat;       % shift stack
stack.ind(:, :, i + 1) = afti.mat;
stack.pts(:, :, i + 1) = aftm.pts;
stack.ptn(:, :, i + 1) = aftm.ptn;
end
%% stack: insert
stack.mat(:, :, 1)     = mat.mat;        % insert new observation
stack.ind(:, :, 1)     = mat.ind;
stack.pts(:, :, 1)     = mat.pts;
stack.ptn(:, :, 1)     = mat.ptn;

end

%% transform map frm - 1 (on the map frm)
function mat  = trns(mat, pts, st)                                         % mat (frm - 1), pts [trans] (frm)

%% move the map [frm] (mat.ptn) to the origin by transformation of [frm + 1] (pts.trn, pts.rtn)
nts           = zeros(size(mat.pts));
nts(:, 1:3)   = mat.ptn - repmat(pts.trn', size(mat.pts, 1), 1);           % trajectory at frm + 1: translation
nts(:, 1:3)   = (pts.rtn \ nts(:, 1:3)')';                                 %                        rotation
nts(:, 3)     = mat.mat(:);                                                % put map values on corresponding locations:[x y values]
%%  move the map [x, y] to the matrix coordinate [i, j]
nts(:, 1:2)   = nts(:, 1:2) ./ repmat([st.vx.x, st.vx.y], ...              % compensate voxel size
                 size(nts, 1), 1);
nts(:, 1:2)   = nts(:, 1:2) + repmat([(-st.vm.xb / st.vx.x + 1), ...       % bias (start index from 1)
                (-st.vm.yr / st.vx.y + 1)], size(nts, 1), 1);
nts(:, 1:2)   = round(nts(:, 1:2));                                        % interpolate: use round!
%% keep valid values [valid indexes]
vind          = ((nts(:,1) >= 1) & (nts(:,1) <= (st.vm.xf - st.vm.xb) / st.vx.x) & ...
                 (nts(:,2) >= 1) & (nts(:,2) <= (st.vm.yl - st.vm.yr) / st.vx.y));
nts(~vind, :) = [];                                                        % filter points inside grid [remove points outsdie grid]
%% convert points (vectors) to matrix with valid index
mat.mat       = zeros(size(mat.mat));                                      % matrix
for i         = 1 : size(nts, 1)
if nts(i, 3) ~= 0                                                          % for object cells
mat.mat(nts(i, 1), nts(i, 2)) = nts(i, 3);                                 % vector to matrix
end
end
%% matrix to vector (mat.pts, mat.ptn)
ivx           = 1 : size(mat.mat, 1);                                      % vxl pts at the origin
ivy           = 1 : size(mat.mat, 2);                                      % index vector
[iy, ix]      = meshgrid(ivy, ivx);                                        % notice! ivy ivx
mat.pts       = [(ix(:) - 1) * st.vx.x + st.vm.xb, ...                     % map vector at origin
                 (iy(:) - 1) * st.vx.y + st.vm.yr, mat.mat(:)];
ptm           = pts.rtn * mat.pts(:, 1:3)';                                % vxl pts in the world coordinates
mat.ptn       = [ptm(1,:)' + pts.trn(1), ptm(2,:)' + pts.trn(2), ...
                 mat.pts(:, 3)];                                           % map vector at real coord

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function y=roundTies2EvenOdd(x,mode)
%%ROUNDTIES2EVENODD Round the values in x to the nearest integer. Ties will
%                   either round to even or odd integers.
%
%INPUTS:   x A scalar, vector, matrix of values to round to the nearest
%            integer.
%       mode A parameter specifying how ties (when the value is an integer
%           +0.5) are handled. Possible values are
%           0 (The default if omitted or an empty matrix is passed) Round
%             ties to even integers.
%           1 Round ties to odd integers.
%
%OUTPUTS: y x rounded with the selected handling of ties.
%
%EXAMPLE:
%As an example of rounding to even versus odd integers, consider:
% x=[-1.5,-2.5,-3.5,-1.1,-1.9,-2.1,-2.9,0,1.5,2.5,3.5,1.1,1.9,2.1,2.9];
% y1=roundTies2EvenOdd(x,0)
% y2=roundTies2EvenOdd(x,1)
%One will find that
%y1=[-2,-2,-4,-1,-2,-2,-3,0,2,2,4,1,2,2,3]
%y2=[-1,-3,-3,-1,-2,-2,-3,0,1,3,3,1,2,2,3]
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<2||isempty(mode))
        mode=0;
    end

    switch(mode)
        case 0%Round with ties going to even
            y=round(x);
            sel1=abs(x-fix(x))==0.5 & mod(fix(x),2)==0;
            y(sel1)=fix(x(sel1));
        case 1%Round with ties going to odd.
            y=round(x);
            sel1=abs(x-fix(x))==0.5 & mod(fix(x),2)==1;
            y(sel1)=fix(x(sel1));
        otherwise
            error('Wrong Mode specified')
    end
end

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