text:
classdef SnirfClass < matlab.mixin.Copyable

    properties
        formatVersion
        metaDataTags
        data
        stim
        probe
        aux
    end

    properties (Access = private)
        fid
        gid
        location
        nirsdatanum
        nirs_tb
        stim0
        filename
        fileformat
        dataStorageScheme
    end

    methods

        % -------------------------------------------------------
        function obj = SnirfClass(varargin)
            %
            % Syntax:
            %   obj = SnirfClass()
            %   obj = SnirfClass(filename);
            %   obj = SnirfClass(filename, nirsdatanum);
            %   obj = SnirfClass(filename, nirsdatanum, dataStorageScheme);
            %   obj = SnirfClass(filename, dataStorageScheme);
            %   obj = SnirfClass(dotnirs);
            %   obj = SnirfClass(dotnirs, numdatabllocks);
            %   obj = SnirfClass(data, stim);
            %   obj = SnirfClass(data, stim, probe);
            %   obj = SnirfClass(data, stim, probe, aux);
            %   obj = SnirfClass(d, t, SD, aux, s);
            %   obj = SnirfClass(d, t, SD, aux, s, CondNames);
            %
            %   Also for debugging/simulation of time bases
            %
            %   obj = SnirfClass(dotnirs, tfactors);
            %
            % Example 1:
            %
            %   % Save .nirs file in SNIRF format
            %   snirf1 = SnirfClass(load('neuro_run01.nirs','-mat'));
            %   snirf1.Save('neuro_run01.snirf');
            %   snirf1.Info()
            %
            %   % Check that the file was saved correctly
            %   snirf2 = SnirfClass();
            %   snirf2.Load('neuro_run01.snirf');
            %   snirf2.Info()
            %
            % Example 2:
            %
            %   Nirs2Snirf('Simple_Probe1.nirs');
            %   obj = SnirfClass('Simple_Probe1.snirf');
            %
            %   Here's some of the output:
            %
            %   obj(1).data ====>
            %
            %       DataClass with properties:
            %
            %           data: [1200x8 double]
            %           time: [1200x1 double]
            %           measurementList: [1x8 MeasListClass]
            %

            % Initialize properties from SNIRF spec
            obj.Initialize()

            % Set class properties NOT part of the SNIRF format
            obj.fileformat = 'hdf5';
            obj.location = '/nirs';
            obj.nirsdatanum = 1;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Between 1 and 4 arguments covers the following syntax variants
            %
            % obj = SnirfClass(filename);
            % obj = SnirfClass(filename, nirsdatanum);
            % obj = SnirfClass(filename, nirsdatanum, dataStorageScheme);
            % obj = SnirfClass(filename, dataStorageScheme);
            % obj = SnirfClass(dotnirs);
            % obj = SnirfClass(dotnirs, numdatabllocks);
            % obj = SnirfClass(data, stim);
            % obj = SnirfClass(data, stim, probe);
            % obj = SnirfClass(data, stim, probe, aux);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if nargin>0 && nargin<5

                % obj = SnirfClass(filename);
                if isa(varargin{1}, 'SnirfClass')
                    obj.Copy(varargin{1});
                    return;
                end

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % obj = SnirfClass(filename, nirsdatanum);
                % obj = SnirfClass(filename, nirsdatanum, dataStorageScheme);
                % obj = SnirfClass(filename, dataStorageScheme);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                if ischar(varargin{1})
                    obj.filename = varargin{1};
                    if nargin>1
                        % obj = SnirfClass(filename, nirsdatanum);
                        if isnumeric(varargin{2})
                            obj.nirsdatanum = varargin{2};

                            % obj = SnirfClass(filename, nirsdatanum, dataStorageScheme);
                            if nargin>2
                                obj.dataStorageScheme = varargin{3};
                            end

                            % obj = SnirfClass(filename, dataStorageScheme);
                        elseif ischar(varargin{2})
                            obj.dataStorageScheme = varargin{2};

                        end
                    end

                    % Load Snirf file here ONLY if data storage scheme is 'memory'
                    if strcmpi(obj.dataStorageScheme, 'memory')
                        obj.Load(varargin{1});
                    end

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % obj = SnirfClass(dotnirs);
                    % obj = SnirfClass(dotnirs, numdatabllocks);
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                elseif isstruct(varargin{1}) || isa(varargin{1}, 'NirsClass')

                    % obj = SnirfClass(dotnirs);
                    tfactors = 1;    % Debug simulation parameter

                    % obj = SnirfClass(dotnirs, numdatabllocks);
                    if nargin==2
                        tfactors = varargin{2};
                    end
                    dotnirs = varargin{1};
                    obj.GenSimulatedTimeBases(dotnirs, tfactors);
                    for ii=1:length(tfactors)
                        obj.data(ii) = DataClass(obj.nirs_tb(ii).d, obj.nirs_tb(ii).t(:), obj.nirs_tb(ii).SD.MeasList);
                    end

                    for ii=1:size(dotnirs.s,2)
                        if isfield(dotnirs, 'CondNames')
                            obj.stim(ii) = StimClass(dotnirs.s(:,ii), dotnirs.t(:), dotnirs.CondNames{ii});
                        else
                            obj.stim(ii) = StimClass(dotnirs.s(:,ii), dotnirs.t(:), num2str(ii));
                        end
                    end
                    obj.probe      = ProbeClass(dotnirs.SD);
                    for ii=1:size(dotnirs.aux,2)
                        obj.aux(ii) = AuxClass(dotnirs.aux(:,ii), dotnirs.t(:), sprintf('aux%d',ii));
                    end

                    % Add metadatatags
                    obj.metaDataTags   = MetaDataTagsClass();

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % obj = SnirfClass(data, stim);
                    % obj = SnirfClass(data, stim, probe);
                    % obj = SnirfClass(data, stim, probe, aux);
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                elseif isa(varargin{1}, 'DataClass')

                    % obj = SnirfClass(data, stim);
                    data = varargin{1};
                    obj.SetData(data);
                    stim = varargin{2};
                    obj.SetStim(stim);

                    % obj = SnirfClass(data, stim, probe);
                    if nargin>2
                        probe = varargin{3};
                        obj.SetSd(probe);
                    end

                    % obj = SnirfClass(data, stim, probe, aux);
                    if nargin>3
                        aux = varargin{4};
                        obj.SetAux(aux);
                    end

                end

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Between 5 and 6 arguments covers the following syntax variants
                %
                % obj = SnirfClass(d, t, SD, aux, s);
                % obj = SnirfClass(d, t, SD, aux, s, CondNames);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            elseif nargin>4

                % obj = SnirfClass(d, t, SD, aux, s);
                d         = varargin{1};
                t         = varargin{2}(:);
                SD        = varargin{3};
                aux       = varargin{4};
                s         = varargin{5};
                CondNames = {};

                % obj = SnirfClass(d, t, SD, aux, s, CondNames);
                if nargin>5
                    CondNames = varargin{6};
                end

                obj.data(1) = DataClass(d, t(:), SD.MeasList);
                for ii=1:size(s,2)
                    if nargin==5
                        condition = num2str(ii);
                    else
                        condition = CondNames{ii};
                    end
                    obj.stim(ii) = StimClass(s(:,ii), t(:), condition);
                end
                obj.probe      = ProbeClass(SD);
                for ii=1:size(aux,2)
                    obj.aux(ii) = AuxClass(aux, t(:), sprintf('aux%d',ii));
                end

                % Add metadatatags
                obj.metaDataTags   = MetaDataTagsClass();

            end

        end

        % -------------------------------------------------------
        function Initialize(obj)
            obj.formatVersion      = '1.0';
            obj.metaDataTags      = MetaDataTagsClass().empty();
            obj.data              = DataClass().empty();
            obj.stim              = StimClass().empty();
            obj.probe             = ProbeClass().empty();
            obj.aux               = AuxClass().empty();

            obj.stim0             = StimClass().empty();
            obj.dataStorageScheme = 'memory';
        end

        % -------------------------------------------------------
        function err = Copy(obj, obj2)
            err=0;
            if ~isa(obj2, 'SnirfClass')
                err=1;
                return;
            end
            obj.formatVersion = obj2.formatVersion;
            obj.metaDataTags  = CopyHandles(obj2.metaDataTags);
            obj.data          = CopyHandles(obj2.data);
            obj.stim          = CopyHandles(obj2.stim);
            obj.probe         = CopyHandles(obj2.probe);
            obj.aux           = CopyHandles(obj2.aux);

            try
                obj.stim0     = CopyHandles(obj2.stim0);
            catch
            end

            if ~isempty(obj2.filename)
                obj.filename = obj2.filename;
            end
        end

        % -------------------------------------------------------
        function objnew = CopyMutable(obj, options)
            if nargin==1
                options = '';
            end

            % Load mutable data from Snirf file here ONLY if data storage scheme is 'files'
            if strcmpi(obj.dataStorageScheme, 'files')
                obj.LoadStim(obj.filename);
            end

            % Generate new instance of SnirfClass
            objnew = SnirfClass();

            objnew.filename = obj.filename;

            % Copy mutable properties to new object instance;
            objnew.stim = CopyHandles(obj.stim);

            if strcmp(options, 'extended')
                t = obj.GetTimeCombined();
                objnew.data = DataClass([],t,[]);
            end
        end

        % -------------------------------------------------------
        function SortStims(obj)
            if isempty(obj.stim)
                return;
            end
            temp = CopyHandles(obj.stim);
            delete(obj.stim);
            names = cell(length(temp),1);
            for ii=1:length(temp)
                names{ii} = temp(ii).name;
            end
            [~,idx] = sort(names);
            obj.stim = temp(idx).copy;
        end

        % -------------------------------------------------------
        function err = SetLocation(obj)
            err = 0;
            gid1 = HDF5_GroupOpen(obj.fid, sprintf('%s%d', obj.location, obj.nirsdatanum));
            gid2 = HDF5_GroupOpen(obj.fid, obj.location);

            if gid1.double > 0
                obj.location = sprintf('%s%d', obj.location, obj.nirsdatanum);
                return;
            elseif gid2.double > 0
                return;
            end
            err = -1;
        end

        % -------------------------------------------------------
        function err = LoadFormatVersion(obj)
            err = 0;
            formatVersionFile = HDF5_DatasetLoad(obj.gid, 'formatVersion'); %#ok<*PROPLC>
            formatVersionFile = str2double(formatVersionFile);
            formatVersionCurr = str2double(obj.formatVersion);
            if formatVersionFile < formatVersionCurr
                fprintf('Warning: Current SNIRF version is %0.1f. Cannot load older version (%0.1f) file. Backward compatibility not yet implemented ...\n', formatVersionCurr, formatVersionFile)
                err = -2;
                return
            end
        end

        % -------------------------------------------------------
        function err = LoadMetaDataTags(obj, fileobj)
            err = 0;
            obj.metaDataTags = MetaDataTagsClass();
            if obj.metaDataTags.LoadHdf5(fileobj, [obj.location, '/metaDataTags']) < 0
                err = -1;
            end
        end

        % -------------------------------------------------------
        function err = LoadData(obj, fileobj)
            err = 0;
            ii=1;
            while 1
                if ii > length(obj.data)
                    obj.data(ii) = DataClass;
                end
                if obj.data(ii).LoadHdf5(fileobj, [obj.location, '/data', num2str(ii)]) < 0
                    obj.data(ii).delete();
                    obj.data(ii) = [];
                    if ii==1
                        err = -1;
                    end
                    break;
                end
                ii=ii+1;
            end
        end

        % -------------------------------------------------------
        function err = LoadStim(obj, fileobj)
            err = 0;

            ii=1;
            while 1
                if ii > length(obj.stim)
                    obj.stim(ii) = StimClass;
                end
                if obj.stim(ii).LoadHdf5(fileobj, [obj.location, '/stim', num2str(ii)]) < 0
                    obj.stim(ii).delete();
                    obj.stim(ii) = [];
                    if ii==1
                        err = -1;
                    end
                    break;
                end
                ii=ii+1;
            end
            obj.SortStims();

            % Load original, unedited stims, if they exist
            ii=1;
            while 1
                if ii > length(obj.stim0)
                    obj.stim0(ii) = StimClass;
                end
                if obj.stim0(ii).LoadHdf5(fileobj, [obj.location, '/stim0', num2str(ii)]) < 0
                    obj.stim0(ii).delete();
                    obj.stim0(ii) = [];
                    break;
                end
                ii=ii+1;
            end

        end

        % -------------------------------------------------------
        function err = LoadProbe(obj, fileobj, ~)
            obj.probe = ProbeClass();
            err = obj.probe.LoadHdf5(fileobj, [obj.location, '/probe']);
        end

        % -------------------------------------------------------
        function err = LoadAux(obj, fileobj)
            err = 0;
            ii=1;
            while 1
                if ii > length(obj.aux)
                    obj.aux(ii) = AuxClass;
                end
                if obj.aux(ii).LoadHdf5(fileobj, [obj.location, '/aux', num2str(ii)]) < 0
                    obj.aux(ii).delete();
                    obj.aux(ii) = [];
                    if ii==1
                        err = -1;
                    end
                    break;
                end
                ii=ii+1;
            end
        end

        % -------------------------------------------------------
        function err = LoadHdf5(obj, fileobj, ~)
            err = 0;

            % Arg 1
            if ~exist('fileobj','var') || ~exist(fileobj,'file')
                fileobj = '';
            end

            % Error checking
            if ~isempty(fileobj) && ischar(fileobj)
                obj.filename = fileobj;
            elseif isempty(fileobj)
                fileobj = obj.filename;
            end
            if isempty(fileobj)
                err = -1;
                return;
            end

            % Don't reload if not empty
            if ~obj.IsEmpty()
                return;
            end

            %%%%%%%%%%%% Ready to load from file

            try

                % Open group
                [obj.gid, obj.fid] = HDF5_GroupOpen(fileobj, '/');

                if obj.SetLocation() < 0
                    err = -1;
                    return
                end

                %%%% Load formatVersion
                if obj.LoadFormatVersion() < 0
                    err = -2;
                end

                %%%% Load metaDataTags
                if obj.LoadMetaDataTags(obj.fid) < 0
                    err = -3;
                end

                %%%% Load data
                if obj.LoadData(obj.fid) < 0
                    err = -4;
                end

                %%%% Load stim
                if obj.LoadStim(obj.fid)
                    err = -5;
                end

                %%%% Load probe
                if obj.LoadProbe(obj.fid)
                    err = -6;
                end

                %%%% Load aux. This is an optional field, therefore error must
                %%%% be less then -1 (-1 means aux is not in SNIRF file) to be
                %%%% error for whole SNIRF file
                if obj.LoadAux(obj.fid)<-1
                    err = -7;
                end

                % Close group
                HDF5_GroupClose(fileobj, obj.gid, obj.fid);

            catch

                err = -1;

            end

            if obj.fid>0
                H5F.close(obj.fid);
            end

        end

        % -------------------------------------------------------
        function err = Load(obj, fileobj)
            err = LoadHdf5(obj, fileobj);
        end

        % -------------------------------------------------------
        function SaveMetaDataTags(obj, fileobj)
            obj.metaDataTags.SaveHdf5(fileobj, [obj.location, '/metaDataTags']);
        end

        % -------------------------------------------------------
        function SaveData(obj, fileobj)
            for ii=1:length(obj.data)
                obj.data(ii).SaveHdf5(fileobj, [obj.location, '/data', num2str(ii)]);
            end
        end

        % -------------------------------------------------------
        function SaveStim(obj, fileobj)
            for ii=1:length(obj.stim)
                obj.stim(ii).SaveHdf5(fileobj, [obj.location, '/stim', num2str(ii)]);
            end
            if isempty(obj.stim0)
                obj.stim0 = obj.stim.copy();
                for ii=1:length(obj.stim0)
                    obj.stim0(ii).SaveHdf5(fileobj, [obj.location, '/stim0', num2str(ii)]);
                end
            end
        end

        % -------------------------------------------------------
        function SaveProbe(obj, fileobj)
            obj.probe.SaveHdf5(fileobj, [obj.location, '/probe']);
        end

        % -------------------------------------------------------
        function SaveAux(obj, fileobj)
            for ii=1:length(obj.aux)
                obj.aux(ii).SaveHdf5(fileobj, [obj.location, '/aux', num2str(ii)]);
            end
        end

        % -------------------------------------------------------
        function SaveHdf5(obj, fileobj, ~)
            % Arg 1
            if ~exist('fileobj','var') || isempty(fileobj)
                error('Unable to save file. No file name given.')
            end

            % Args
            if exist(fileobj, 'file')
                delete(fileobj);
            end
            obj.fid = H5F.create(fileobj, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
            H5F.close(obj.fid);

            %%%%% Save this object's properties

            % Save formatVersion
            if isempty(obj.formatVersion)
                obj.formatVersion = '1.1';
            end
            hdf5write_safe(fileobj, '/formatVersion', obj.formatVersion);

            % Save metaDataTags
            obj.SaveMetaDataTags(fileobj);

            % Save data
            obj.SaveData(fileobj);

            % Save stim
            obj.SaveStim(fileobj);

            % Save sd
            obj.SaveProbe(fileobj);

            % Save aux
            obj.SaveAux(fileobj);
        end

        % -------------------------------------------------------
        function Save(obj, fileobj)
            SaveHdf5(obj, fileobj);
        end

        % -------------------------------------------------------
        function [stimFromFile, changes] = UpdateStim(obj, fileobj)
            flags = zeros(length(obj.stim), 1);

            % Load stim from file and update it
            snirfFile = SnirfClass();
            snirfFile.LoadStim(fileobj);

            % Update stims from file with edited stims
            for ii = 1:length(obj.stim)
                for jj = 1:length(snirfFile.stim)
                    if strcmp(obj.stim(ii).GetName(), snirfFile.stim(jj).GetName())
                        if obj.stim(ii) ~= snirfFile.stim(jj)
                            snirfFile.stim(jj).Copy(obj.stim(ii));
                        end
                        flags(ii) = 1;
                        break;
                    end
                end
                if ~flags(ii)
                    % We have new stimulus condition added
                    if ~obj.stim(ii).IsEmpty()
                        snirfFile.stim(jj+1) = StimClass(obj.stim(ii));
                    end
                end
            end

            % If stims were edited then update snirf file with new stims
            changes = sum(flags);
            if changes
                snirfFile.SaveStim(fileobj);
            end
            stimFromFile = snirfFile.stim;
        end

        % -------------------------------------------------------
        function CopyStim(obj, obj2)
            for ii = 1:length(obj2.stim)
                if ii > length(obj.stim)
                    obj.stim(ii) = StimClass(obj2.stim(ii));
                else
                    obj.stim(ii).Copy(obj2.stim(ii));
                end
            end
        end

        % -------------------------------------------------------
        function changes = StimChangesMade(obj)

            flags = zeros(length(obj.stim), 1);

            % Load stims from file
            snirf = SnirfClass();
            snirf.LoadStim(obj.filename);
            stimFromFile = snirf.stim;

            % Update stims from file with edited stims
            for ii = 1:length(obj.stim)
                for jj = 1:length(stimFromFile)
                    if strcmp(obj.stim(ii).GetName(), stimFromFile(jj).GetName())
                        if obj.stim(ii) ~= stimFromFile(jj)
                            flags(ii) = 1;
                        else
                            flags(ii) = -1;
                        end
                        break;
                    end
                end
                if flags(ii)==0
                    % We have new stimulus condition added
                    if ~obj.stim(ii).IsEmpty()
                        flags(ii) = 1;
                    end
                end
            end
            flags(flags ~= 1) = 0;
            changes = sum(flags)>0;
        end

        % -------------------------------------------------------
        function b = DataModified(obj)
            b = obj.StimChangesMade();
        end

        % -------------------------------------------------------
        function err = SaveMutable(obj, fileobj)
            if isempty(obj)
                return
            end

            % Arg 1
            if ~exist('fileobj','var') || ~exist(fileobj,'file')
                fileobj = '';
            end

            % Error checking
            if ~isempty(fileobj) && ischar(fileobj)
                obj.filename = fileobj;
            elseif isempty(fileobj)
                fileobj = obj.filename;
            end
            if isempty(fileobj)
                err = -1;
                return;
            end

            % Update original stims and save back to file
            obj.UpdateStim(fileobj);
        end

        % -------------------------------------------------------
        function B = eq(obj, obj2)
            B = false;
            if ~strcmp(obj.formatVersion, obj2.formatVersion)
                return;
            end
            if length(obj.data) ~= length(obj2.data)
                return;
            end
            for ii = 1:length(obj.data)
                if ~(obj.data(ii) == obj2.data(ii))
                    return;
                end
            end
            if length(obj.stim)~=length(obj2.stim)
                return;
            end
            for ii = 1:length(obj.stim)
                flag = false;
                for jj = 1:length(obj2.stim)
                    if obj.stim(ii) == obj2.stim(jj)
                        flag = true;
                        break;
                    end
                end
                if flag==false
                    return;
                end
            end
            if ~(obj.probe == obj2.probe)
                return;
            end
            if length(obj.aux) ~= length(obj2.aux)
                return;
            end
            for ii=1:length(obj.aux)
                if ~(obj.aux(ii) == obj2.aux(ii))
                    return;
                end
            end
            if length(obj.metaDataTags) ~= length(obj2.metaDataTags)
                return;
            end
            for ii = 1:length(obj.metaDataTags)
                if ~(obj.metaDataTags(ii) == obj2.metaDataTags(ii))
                    return;
                end
            end
            B = true;
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Basic methods to Set/Get native variable
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ---------------------------------------------------------
        function val = GetFormatVersion(obj)
            val = obj.formatVersion;
        end

        % ---------------------------------------------------------
        function val = GetFormatVersionString(obj)
            val = sprintf('SNIRF v%s', obj.formatVersion);
        end

        % ---------------------------------------------------------
        function SetData(obj, val)
            obj.data = CopyHandles(val);
        end

        % ---------------------------------------------------------
        function val = GetData(obj)
            val = obj.data;
        end

        % ---------------------------------------------------------
        function SetStim(obj, val)
            obj.stim = CopyHandles(val);
        end

        % ---------------------------------------------------------
        function val = GetStim(obj)
            val = obj.stim;
        end

        % ---------------------------------------------------------
        function SetSd(obj, val)
            obj.probe = CopyHandles(val);
        end

        % ---------------------------------------------------------
        function val = GetSd(obj)
            val = obj.probe;
        end

        % ---------------------------------------------------------
        function SetAux(obj, val)
            obj.aux = CopyHandles(val);
        end

        % ---------------------------------------------------------
        function val = GetAux(obj)
            val = obj.aux;
        end

        % ---------------------------------------------------------
        function SetMetaDataTags(obj, val)
            obj.metaDataTags = val;
        end

        % ---------------------------------------------------------
        function val = GetMetaDataTags(obj)
            val = obj.metaDataTags;
            if isempty(obj)
                return;
            end
            if isempty(obj.metaDataTags)
                return;
            end
            val = obj.metaDataTags.Get();
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods that must be implemented as a child class of AcqDataClass
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ---------------------------------------------------------
        function t = GetTime(obj, iBlk)
            t = [];
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk=1;
            end
            if iBlk>length(obj.data)
                return;
            end
            t = obj.data(iBlk).GetTime();
        end

        % ---------------------------------------------------------
        function datamat = GetDataTimeSeries(obj, options, iBlk)
            datamat = [];
            if ~exist('options','var')
                options = '';
            end
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk=1;
            end
            if iBlk>length(obj.data)
                return;
            end
            datamat = obj.data(iBlk).GetDataTimeSeries(options);
        end

        % ---------------------------------------------------------
        function datamat = GetAuxDataMatrix(obj)
            datamat = [];
            if isempty(obj.aux)
                return;
            end
            for ii=1:length(obj.aux)
                datamat(:,ii) = obj.aux(ii).GetDataTimeSeries();
            end
        end

        % ---------------------------------------------------------
        function names = GetAuxNames(obj)
            names = {};
            for ii=1:length(obj.aux)
                names{ii} = obj.aux(ii).GetName();
            end
        end

        % ---------------------------------------------------------
        function aux = GetAuxiliary(obj)
            aux = struct('names',{{}}, 'data', obj.aux);
            aux.names = obj.GetAuxNames();
            aux.data = obj.GetAuxDataMatrix();
        end

        % ---------------------------------------------------------
        function ml = GetMeasList(obj, iBlk)
            ml = [];
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk=1;
            end
            if iBlk>length(obj.data)
                return;
            end
            ml = obj.data(iBlk).GetMeasList();
        end

        % ---------------------------------------------------------
        function wls = GetWls(obj)
            wls = obj.probe.GetWls();
        end

        % ---------------------------------------------------------
        function SetStims_MatInput(obj, s, t)
            if nargin<2
                return
            end
            if isempty(t)
                return;
            end
            for ii=1:size(s,2)
                tidxs = find(s(:,ii)~=0);
                for jj=1:length(tidxs)
                    if ~obj.stim(ii).Exists(t(tidxs(jj)))
                        obj.stim(ii).AddStims(t(tidxs(jj)));
                    else
                        obj.stim(ii).EditValue(t(tidxs(jj)), s(tidxs(jj),ii));
                    end
                end
            end
        end

        % ---------------------------------------------------------
        function s = GetStims(obj, t)
            s = zeros(length(t), length(obj.stim));
            for ii=1:length(obj.stim)
                [ts, v] = obj.stim(ii).GetStim();
                [~, k] = nearest_point(t, ts);
                if isempty(k)
                    continue;
                end
                s(k,ii) = v;
            end
        end

        % ----------------------------------------------------------------------------------
        function SetConditions(obj, CondNames)
            if nargin==1
                return;
            end
            CondNamesLocal = unique({obj.stim.name});
            stimnew = StimClass().empty;
            for ii=1:length(CondNames)
                k = find(strcmp(CondNamesLocal, CondNames{ii}));
                if ~isempty(k)
                    stimnew(ii) = StimClass(obj.stim(k));
                else
                    stimnew(ii) = StimClass(CondNames{ii});
                end
            end
            obj.stim = stimnew.copy;
        end

        % ---------------------------------------------------------
        function CondNames = GetConditions(obj)
            CondNames = cell(1,length(obj.stim));
            for ii=1:length(obj.stim)
                CondNames{ii} = obj.stim(ii).GetName();
            end
        end

        % ---------------------------------------------------------
        function SD = GetSDG(obj)
            SD = [];
            if isempty(obj)
                return;
            end
            if isempty(obj.probe)
                return;
            end
            SD.Lambda = obj.probe.GetWls();
            SD.SrcPos = obj.probe.GetSrcPos();
            SD.DetPos = obj.probe.GetDetPos();
        end

        % ---------------------------------------------------------
        function srcpos = GetSrcPos(obj)
            srcpos = obj.probe.GetSrcPos();
        end

        % ---------------------------------------------------------
        function detpos = GetDetPos(obj)
            detpos = obj.probe.GetDetPos();
        end

        % ----------------------------------------------------------------------------------
        function n = GetDataBlocksNum(obj)
            n = length(obj.data);
        end

        % ----------------------------------------------------------------------------------
        function [iDataBlks, ich] = GetDataBlocksIdxs(obj, ich0)
            iDataBlks=[];
            ich={};
            if nargin==1
                ich0=[];
            end
            if isempty(ich0)
                iDataBlks=1:length(obj.data);
                return;
            end

            % Get channel matrix for whole probe
            mlAll = [];
            nDataBlks = length(obj.data);
            for iBlk = 1:nDataBlks
                mlAll = [mlAll; obj.GetMeasList(iBlk)];
            end

            iSrc = mlAll(ich0,1);
            iDet = mlAll(ich0,2);

            % Now search block by block for the selecdted channels
            ich = cell(nDataBlks,1);
            for iBlk=1:nDataBlks
                ml = obj.GetMeasList(iBlk);
                for ii=1:length(ich0)
                    k = find(ml(:,1)==iSrc(ii) & ml(:,2)==iDet(ii));
                    if ~isempty(k)
                        iDataBlks = [iDataBlks; iBlk];
                        ich{iBlk} = [ich{iBlk}, k(1)];
                    end
                end
            end

            % Important: make sure iDataBlks is row vector (: + transpose does that) .
            % For some reason a for-loop traversing through empty column vector doesn't work properly
            iDataBlks = sort(unique(iDataBlks(:)'));

        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Pubic interface for .nirs processing stream
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function d = Get_d(obj, iBlk)
            d = [];
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk = 1;
            end
            if iBlk>length(obj.data)
                return;
            end
            d = obj.data(iBlk).GetDataTimeSeries();
        end

        % ----------------------------------------------------------------------------------
        function t = Get_t(obj, iBlk)
            t = [];
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk = 1;
            end
            if iBlk>length(obj.data)
                return;
            end
            t = obj.data(iBlk).GetTime();
        end

        % ----------------------------------------------------------------------------------
        function SD = Get_SD(obj, iBlk)
            SD = [];
            if isempty(obj.probe)
                return;
            end
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk = 1;
            end
            if iBlk>length(obj.data)
                return;
            end
            SD.Lambda   = obj.probe.GetWls();
            SD.SrcPos   = obj.probe.GetSrcPos();
            SD.DetPos   = obj.probe.GetDetPos();
            SD.MeasList = obj.data(iBlk).GetMeasList();
            SD.MeasListAct = ones(size(SD.MeasList,1),1);
        end

        % ----------------------------------------------------------------------------------
        function s = Get_s(obj, iBlk)
            s = [];
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk = 1;
            end
            if iBlk>length(obj.data)
                return;
            end
            t = obj.data(iBlk).GetTime();
            s = zeros(length(t), length(obj.stim));
            for ii=1:length(obj.stim)
                [ts, v] = obj.stim(ii).GetStim();
                [~, k] = nearest_point(t, ts);
                if isempty(k)
                    continue;
                end
                s(k,ii) = v;
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % All other public methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function AddStims(obj, tPts, condition)
            % Try to find existing condition to which to add stims.
            for ii=1:length(obj.stim)
                if strcmp(condition, obj.stim(ii).GetName())
                    obj.stim(ii).AddStims(tPts);
                    return;
                end
            end

            % Otherwise we have a new condition to which to add the stims.
            obj.stim(end+1) = StimClass(tPts, condition);
            obj.SortStims();
        end

        % ----------------------------------------------------------------------------------
        function DeleteStims(obj, tPts, condition)
            % Find all stims for any conditions which match the time points.
            for ii=1:length(obj.stim)
                obj.stim(ii).DeleteStims(tPts);
            end
        end

        % ----------------------------------------------------------------------------------
        function ToggleStims(obj, tPts, condition)
            % Find all stims for any conditions which match the time points.
            for ii=1:length(obj.stim)
                obj.stim(ii).ToggleStims(tPts);
            end
        end

        % ----------------------------------------------------------------------------------
        function MoveStims(obj, tPts, condition)
            if ~exist('tPts','var') || isempty(tPts)
                return;
            end
            if ~exist('condition','var') || isempty(condition)
                return;
            end

            % Find the destination condition to move stims (among the time pts in tPts)
            % to
            j = [];
            for ii=1:length(obj.stim)
                if strcmp(condition, obj.stim(ii).GetName())
                    j=ii;
                    break;
                end
            end

            % If no destination condition found among existing conditions,
            % then create a new condition to move stims to
            if isempty(j)
                j = length(obj.stim)+1;

                % Otherwise we have a new condition to which to add the stims.
                obj.stim(j) = StimClass([], condition);
                obj.SortStims();

                % Recalculate j after sort
                for ii=1:length(obj.stim)
                    if strcmp(condition, obj.stim(ii).GetName())
                        j=ii;
                        break;
                    end
                end
            end

            % Find all stims for any conditions which match the time points.
            for ii=1:length(tPts)
                for kk=1:length(obj.stim)
                    d = obj.stim(kk).GetData();
                    if isempty(d)
                        continue;
                    end
                    k = find(d(:,1)==tPts(ii));
                    if ~isempty(k)
                        if kk==j
                            continue;
                        end

                        % If stim at time point tPts(ii) exists in stim
                        % condition kk, then move stim from obj.stim(kk) to
                        % obj.stim(j)
                        obj.stim(j).AddStims(tPts(ii), d(k(1),2), d(k(1),3));

                        % After moving stim from obj.stim(kk) to
                        % obj.stim(j), delete it from obj.stim(kk)
                        d(k(1),:)=[];
                        obj.stim(kk).SetData(d);

                        % Move on to next time point
                        break;
                    end
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function SetStimTpts(obj, icond, tpts)
            obj.stim(icond).SetTpts(tpts);
        end

        % ----------------------------------------------------------------------------------
        function tpts = GetStimTpts(obj, icond)
            if icond>length(obj.stim)
                tpts = [];
                return;
            end
            tpts = obj.stim(icond).GetTpts();
        end

        % ----------------------------------------------------------------------------------
        function SetStimDuration(obj, icond, duration)
            obj.stim(icond).SetDuration(duration);
        end

        % ----------------------------------------------------------------------------------
        function duration = GetStimDuration(obj, icond)
            if icond>length(obj.stim)
                duration = [];
                return;
            end
            duration = obj.stim(icond).GetDuration();
        end

        % ----------------------------------------------------------------------------------
        function SetStimValues(obj, icond, vals)
            obj.stim(icond).SetValues(vals);
        end

        % ----------------------------------------------------------------------------------
        function vals = GetStimValues(obj, icond)
            if icond>length(obj.stim)
                vals = [];
                return;
            end
            vals = obj.stim(icond).GetValues();
        end

        % ----------------------------------------------------------------------------------
        function RenameCondition(obj, oldname, newname)
            if ~exist('oldname','var') || ~ischar(oldname)
                return;
            end
            if ~exist('newname','var')  || ~ischar(newname)
                return;
            end
            k=[];
            for ii=1:length(obj.stim)
                if strcmp(obj.stim(ii).GetName(), oldname)
                    k = ii;
                    break;
                end
            end
            if isempty(k)
                return;
            end
            obj.stim(k).SetName(newname);
            obj.SortStims();
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return;
            end
            if isempty(obj.data) || (isa(obj.data, 'DataClass') && obj.data(1).IsEmpty())
                return;
            end
            if isempty(obj.data) || (isa(obj.probe, 'ProbeClass') && obj.probe.IsEmpty())
                return;
            end
            b = false;
        end

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            nbytes = nbytes + sizeof(obj.formatVersion);
            for ii=1:length(obj.metaDataTags)
                nbytes = nbytes + obj.metaDataTags(ii).MemoryRequired();
            end
            for ii=1:length(obj.data)
                nbytes = nbytes + obj.data(ii).MemoryRequired();
            end
            for ii=1:length(obj.stim)
                nbytes = nbytes + obj.stim(ii).MemoryRequired();
            end
            if ~isempty(obj.probe)
                nbytes = nbytes + obj.probe.MemoryRequired();
            end
            for ii=1:length(obj.aux)
                nbytes = nbytes + obj.aux(ii).MemoryRequired();
            end
        end

        % ----------------------------------------------------------------------------------
        function Info(obj)

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Read formatVersion
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            fprintf('\n');
            fprintf('    FormatVersion:\n');
            fv = obj.GetFormatVersion();
            fprintf('        Format version: %s\n', fv);
            fprintf('\n');

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Load meta data tags from file and extract the tag names and values for display
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            fprintf('    MetaDataTags:\n');
            tags = obj.GetMetaDataTags();
            for ii=1:length(tags)
                fprintf('        Tag #%d: {''%s'', ''%s''}\n', ii, tags(ii).key, tags(ii).value);
            end
            fprintf('\n');

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Load data from file and extract .nirs-style d and ml matrices for display
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            fprintf('    Data (.nirs-style display):\n');
            for ii=1:length(obj.data)

                % Display data matrix dimensions and data type
                d = obj.data(ii).GetDataTimeSeries();
                pretty_print_struct(d, 8, 1);

                % Display meas list dimensions and data type
                ml = obj.data(ii).GetMeasList();
                pretty_print_struct(ml, 8, 1);

            end
            fprintf('\n');

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Load probe and extract .nirs-style SD structure
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            fprintf('    Probe (.nirs-style display):\n');
            SD = obj.GetSDG();
            pretty_print_struct(SD, 8, 1);
            fprintf('\n');

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Load stim from file and extract it for display
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            fprintf('    Stim (.snirf-style display):\n');
            for ii=1:length(obj.stim)
                fprintf('        stim(%d): {name = ''%s'', data = [', ii, obj.stim(ii).name);
                for jj=1:size(obj.stim(ii).data,1)
                    if jj==size(obj.stim(ii).data,1)
                        fprintf('%0.1f', obj.stim(ii).data(jj,1));
                    else
                        fprintf('%0.1f, ', obj.stim(ii).data(jj,1));
                    end
                end
                fprintf(']}\n');
            end
            fprintf('\n');

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Load aux from file and extract nirs-style data for display
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            fprintf('    Aux (.nirs-style display):\n');
            auxl = obj.GetAuxiliary();
            pretty_print_struct(auxl, 8, 1);
            fprintf('\n');

        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Private methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = private)

        function GenSimulatedTimeBases(obj, nirs, tfactors)
            obj.nirs_tb = struct('SD', struct('Lambda',nirs.SD.Lambda, 'MeasList',[], 'SrcPOos',[], 'DetPOos',[]), ...
                't',[], ...
                'd',[] ...
                );

            if length(tfactors)==1 && tfactors==1
                obj.nirs_tb = nirs;
                return;
            end

            % a) Subdivide data and measurement list among time bases
            % b) Resample time and data
            % c) Throw out source and detectors that don't belong in each time base
            nCh = size(nirs.SD.MeasList,1)/length(nirs.SD.Lambda);
            nTimeBases = length(tfactors);
            baseSize = round(nCh/nTimeBases);
            irows = 1:baseSize:nCh;

            % Assign channels for time bases
            obj.nirs_tb = repmat(obj.nirs_tb, nTimeBases,1);
            for iWl=1:length(nirs.SD.Lambda)
                iBase = 1;
                for ii=irows+(iWl-1)*nCh
                    istart = ii;
                    if ii+baseSize-1 <= iWl*nCh
                        iend = ii+baseSize-1;
                    else
                        iend = iWl*nCh;
                    end
                    nChAdd = iend-istart+1;
                    obj.nirs_tb(iBase).d(:,end+1:end+nChAdd) = nirs.d(:,istart:iend);
                    obj.nirs_tb(iBase).SD.MeasList(end+1:end+nChAdd,:) = nirs.SD.MeasList(istart:iend,:);
                    if iBase<nTimeBases
                        iBase = iBase+1;
                    end
                end
            end

            % Resample data time and throw out optodes that don't belong in each time base
            for iBase=1:length(obj.nirs_tb)
                % Resample time
                [n,d] = rat(tfactors(iBase));
                obj.nirs_tb(iBase).t = resample(nirs.t, n, d);

                % Resample data
                obj.nirs_tb(iBase).d = resample(obj.nirs_tb(iBase).d, n, d);

                % Throw out source and detectors that don't belong in each time base
                iSrc = unique(obj.nirs_tb(iBase).SD.MeasList(:,1));
                obj.nirs_tb(iBase).SD.SrcPos = nirs.SD.SrcPos(iSrc);
                iDet = unique(obj.nirs_tb(iBase).SD.MeasList(:,2));
                obj.nirs_tb(iBase).SD.DetPos = nirs.SD.DetPos(iDet);
            end

        end

    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x, y, h] = draw_graph(adj, labels, node_t, x, y, varargin)
% DRAW_LAYOUT		Draws a layout for a graph
%
%  [X, Y, H] = DRAW_LAYOUT(ADJ, <LABELS, ISBOX, X, Y>)
%
% Inputs :
%	ADJ : Adjacency matrix (source, sink)
%       LABELS : Cell array containing labels <Default : '1':'N'>
%       ISBOX : 1 if node is a box, 0 if oval <Default : zeros>
%       X, Y, : Coordinates of nodes on the unit square <Default : calls make_layout>
%
% Outputs :
%	X, Y : Coordinates of nodes on the unit square
%       H    : Object handles
%
% Usage Example : [x, y] = draw_layout([0 1;0 0], {'Hidden','Visible'}, [1 0]');
%
% h(i,1) is the text handle - color
% h(i,2) is the circle handle - facecolor
%
% See also MAKE_LAYOUT

% Change History :
% Date		Time		Prog	Note
% 13-Apr-2000	 9:06 PM	ATC	Created under MATLAB 5.3.1.29215a (R11.1)
%
% ATC = Ali Taylan Cemgil,
% SNN - University of Nijmegen, Department of Medical Physics and Biophysics
% e-mail : cemgil@mbfys.kun.nl

adj = double(adj);
N = size(adj,1);
if nargin<2,
  labels = cellstr(int2str((1:N)'));
end

if nargin<3,
  node_t = zeros(N,1);
else
  node_t = node_t(:);
end;

axis([0 1 0 1]);
set(gca,'XTick',[],'YTick',[],'box','on');
% axis('square');
%colormap(flipud(gray));

if nargin<4,
  [x y] = make_layout(adj);
end;

idx1 = find(node_t==0); h1 = []; wd1=[];
if ~isempty(idx1)
  [h1 wd1] = textoval(x(idx1), y(idx1), labels(idx1), varargin{:});
end;

idx2 = find(node_t~=0); h2 = []; wd2 = [];
if ~isempty(idx2)
  [h2 wd2] = textbox(x(idx2), y(idx2), labels(idx2), varargin{:});
end;

wd = zeros(size(wd1,1)+size(wd2,1),2);
if ~isempty(idx1), wd(idx1, :) = wd1;  end;
if ~isempty(idx2), wd(idx2, :) = wd2; end;

% coloring
color.box = 'black';
color.text = color.box;
color.edge = [1 1 1]*3/4;
%color.edge = 'green';
if ~isempty(idx1)
  set(h1(:,1),'Color',color.text)
  set(h1(:,2),'EdgeColor',color.box)
end
if ~isempty(idx2)
  set(h2(:,1),'Color',color.text)
  set(h2(:,2),'EdgeColor',color.box)
end

% bug: this code assumes [x y] is the center of each box and oval, which
% isn't exactly true.
h_edge = [];
for i=1:N,
  j = find(adj(i,:)==1);
  for k=j,
    if x(k)-x(i)==0,
	sign = 1;
	if y(i)>y(k), alpha = -pi/2; else alpha = pi/2; end;
    else
	alpha = atan((y(k)-y(i))/(x(k)-x(i)));
	if x(i)<x(k), sign = 1; else sign = -1; end;
    end;
    dy1 = sign.*wd(i,2).*sin(alpha);   dx1 = sign.*wd(i,1).*cos(alpha);
    dy2 = sign.*wd(k,2).*sin(alpha);   dx2 = sign.*wd(k,1).*cos(alpha);
    if adj(k,i)==0, % if directed edge
      h = arrow([x(i)+dx1 y(i)+dy1],[x(k)-dx2 y(k)-dy2],'BaseAngle',30);
      set(h, 'FaceColor', color.edge)
      set(h, 'EdgeColor', color.edge)
    else
      h = line([x(i)+dx1 x(k)-dx2],[y(i)+dy1 y(k)-dy2]);
      set(h, 'Color', color.edge)
      adj(k,i)=-1; % Prevent drawing lines twice
    end;
    h_edge = [h_edge h];
  end;
end;

if nargout>2,
  h = zeros(length(wd),2);
  if ~isempty(idx1),
    h(idx1,:) = h1;
  end;
  if ~isempty(idx2),
    h(idx2,:) = h2;
  end;
end;

%%%%%

function [t, wd] = textoval(x, y, str, varargin)
% TEXTOVAL		Draws an oval around text objects
%
%  [T, WIDTH] = TEXTOVAL(X, Y, STR)
%  [..] = TEXTOVAL(STR)  % Interactive
%
% Inputs :
%    X, Y : Coordinates
%    TXT  : Strings
%
% Outputs :
%    T : Object Handles
%    WIDTH : x and y Width of ovals
%
% Usage Example : [t] = textoval('Visit to Asia?');
%
%
% Note     :
% See also TEXTBOX

% Uses :

% Change History :
% Date		Time		Prog	Note
% 15-Jun-1998	10:36 AM	ATC	Created under MATLAB 5.1.0.421
% 12-Mar-2004   10:00 AM        minka   Changed placement/sizing.
%
% ATC = Ali Taylan Cemgil,
% SNN - University of Nijmegen, Department of Medical Physics and Biophysics
% e-mail : cemgil@mbfys.kun.nl

temp = [];
textProperties = {'BackgroundColor','Color','FontAngle','FontName','FontSize','FontUnits','FontWeight','Rotation'};
varargin = argfilter(varargin,textProperties);

if nargin == 1
  str = x;
end
if ~isa(str,'cell') str=cellstr(str); end;
N = length(str);
wd = zeros(N,2);
for i=1:N,
  if nargin == 1
    [x, y] = ginput(1);
  end
  tx = text(x(i),y(i),str{i},'HorizontalAlignment','center',varargin{:});
  % minka
  [ptc wx wy] = draw_oval(tx);
  wd(i,:) = [wx wy];
  % draw_oval will paint over the text, so need to redraw it
  delete(tx);
  tx = text(x(i),y(i),str{i},'HorizontalAlignment','center',varargin{:});
  temp = [temp;  tx ptc];
end
if nargout>0, t = temp; end;

%%%%%%%%%

function [ptc, wx, wy] = draw_oval(tx, x, y)
% Draws an oval box around a tex object
sz = get(tx,'Extent');
% minka
wy = 2/3*sz(4);
wx = 2/3*sz(3);
x = sz(1)+sz(3)/2;
y = sz(2)+sz(4)/2;
ptc = ellipse(x, y, wx, wy);
set(ptc, 'FaceColor','w');

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

function [p] = ellipse(x, y, rx, ry, c)
% ELLIPSE		Draws Ellipse shaped patch objects
%
%  [<P>] = ELLIPSE(X, Y, Rx, Ry, C)
%
% Inputs :
%    X : N x 1 vector of x coordinates
%    Y : N x 1 vector of y coordinates
%    Rx, Ry : Radii
%    C : Color index
%
%
% Outputs :
%    P = Handles of Ellipse shaped path objects
%
% Usage Example : [] = ellipse();
%
%
% Note     :
% See also

% Uses :

% Change History :
% Date		Time		Prog	Note
% 27-May-1998	 9:55 AM	ATC	Created under MATLAB 5.1.0.421

% ATC = Ali Taylan Cemgil,
% SNN - University of Nijmegen, Department of Medical Physics and Biophysics
% e-mail : cemgil@mbfys.kun.nl

if (nargin < 2) error('Usage Example : e = ellipse([0 1],[0 -1],[1 0.5],[2 0.5]); '); end;
if (nargin < 3) rx = 0.1; end;
if (nargin < 4) ry = rx; end;
if (nargin < 5) c = 1; end;

if length(c)==1, c = ones(size(x)).*c; end;
if length(rx)==1, rx = ones(size(x)).*rx; end;
if length(ry)==1, ry = ones(size(x)).*ry; end;

n = length(x);
p = zeros(size(x));
t = 0:pi/30:2*pi;
for i=1:n,
	px = rx(i)*cos(t)+x(i);
	py = ry(i)*sin(t)+y(i);
	p(i) = patch(px,py,c(i));
end;

if nargout>0, pp = p; end;

%%%%%

function [t, wd] = textbox(x,y,str,varargin)
% TEXTBOX	Draws A Box around the text
%
%  [T, WIDTH] = TEXTBOX(X, Y, STR)
%  [..] = TEXTBOX(STR)
%
% Inputs :
%    X, Y : Coordinates
%    TXT  : Strings
%
% Outputs :
%    T : Object Handles
%    WIDTH : x and y Width of boxes
%%
% Usage Example : t = textbox({'Ali','Veli','49','50'});
%
%
% Note     :
% See also TEXTOVAL

% Uses :

% Change History :
% Date		Time		Prog	Note
% 09-Jun-1998	11:43 AM	ATC	Created under MATLAB 5.1.0.421
% 12-Mar-2004   10:00 AM        minka   Changed placement/sizing.
%
% ATC = Ali Taylan Cemgil,
% SNN - University of Nijmegen, Department of Medical Physics and Biophysics
% e-mail : cemgil@mbfys.kun.nl

temp = [];
textProperties = {'BackgroundColor','Color','FontAngle','FontName','FontSize','FontUnits','FontWeight','Rotation'};
varargin = argfilter(varargin,textProperties);

if nargin == 1
  str = x;
end
if ~isa(str,'cell') str=cellstr(str); end;
N = length(str);
wd = zeros(N,2);
for i=1:N,
  if nargin == 1
    [x, y] = ginput(1);
  end
  tx = text(x(i),y(i),str{i},'HorizontalAlignment','center',varargin{:});
  % minka
  [ptc wx wy] = draw_box(tx);
  wd(i,:) = [wx wy];
  % draw_box will paint over the text, so need to redraw it
  delete(tx);
  tx = text(x(i),y(i),str{i},'HorizontalAlignment','center',varargin{:});
  temp = [temp; tx ptc];
end;

if nargout>0, t = temp; end;

function [ptc, wx, wy] = draw_box(tx)
% Draws a box around a text object
sz = get(tx,'Extent');
% minka
wy = 1/2*sz(4);
wx = 1/2*sz(3);
x = sz(1)+sz(3)/2;
y = sz(2)+sz(4)/2;
ptc = patch([x-wx x+wx x+wx x-wx], [y+wy y+wy y-wy y-wy],'w');
set(ptc, 'FaceColor','w');

function args = argfilter(args,keep)
%ARGFILTER  Remove unwanted arguments.
% ARGFILTER(ARGS,KEEP), where ARGS = {'arg1',value1,'arg2',value2,...},
% returns a new argument list where only the arguments named in KEEP are
% retained.  KEEP is a character array or cell array of strings.

% Written by Tom Minka

if ischar(keep)
  keep = cellstr(keep);
end
i = 1;
while i < length(args)
  if ~ismember(args{i},keep)
    args = args(setdiff(1:length(args),[i i+1]));
  else
    i = i + 2;
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ y, m, d, f ] = jed_to_ymdf_soor_san ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_SOOR_SAN converts a JED to a Soor San YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%
  jed_epoch = epoch_to_jed_soor_san ( );

  j = floor ( jed - jed_epoch );
  f = ( jed - jed_epoch ) - j;

  d = 1 + j;
  m = 1;
  y = 1;
%
%  Account for the number of completed 4 year cycles of 1461 days.
%
  n = floor ( ( d - 1 ) / 1461 );
  y = y + 4 * n;
  d = d - n * 1461;
%
%  Account for the number of completed 365 day years.
%
  n = floor ( ( d - 1 ) / 365 );
  y = y + n;
  d = d - n * 365;
%
%  Account for the number of completed 30 day months.
%
  n = floor ( ( d - 1 ) / 30 );
  m = m + n;
  d = d - n * 30;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef SAMSO < ALGORITHM
% <single> <real/integer> <large/none> <expensive>
% Multiswarm-assisted expensive optimization

%------------------------------- Reference --------------------------------
% F. Li, X. Cai, L. Gao, and W. Shen, A surrogate-assisted multiswarm
% optimization algorithm for high-dimensional computationally expensive
% problems, IEEE Transactions on Cybernetics, 2021, 51(3): 1390-1402.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            assert(~isempty(ver('Optim')),'The execution of SAMPSO requires the Optimization Toolbox.');

            %% Parameter setting
            [Wnc,Pr] = Algorithm.ParameterSet(1,0.5);
            eta      = min(sqrt(0.001^2*Problem.D),5e-4*min(Problem.upper-Problem.lower));

            %% Initialze DB
            if Problem.D > 50
                N = 80;
                K = 2*Problem.D;
            else
                N = 40;
                K = N;
            end
            PopDec = UniformPoint(K,Problem.D,'Latin');
            Population = Problem.Evaluation(repmat(Problem.upper-Problem.lower,K,1).*PopDec+repmat(Problem.lower,K,1));

            %% Initialize swarm
            % Determine position
            [~,idx]  = sort(Population.objs,'ascend');
            Select   = idx(1:N);
            Position = [Population(Select).decs,Population(Select).objs];
            % velocity
            Vmax     = 0.5*(Problem.upper-Problem.lower);
            Vmin     = -0.5*Vmax;
            Velocity = rand(N,Problem.D).*(repmat(Vmax-Vmin,N,1)) + repmat(Vmin,N,1);
            % Pbest and Gbest
            Pbest    = Position;
            Gbest    = Position(1,:);
            maxFES   = Problem.maxFE - K;

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % Build RBF surrogate model
                [model,~] = rbf_build(Population.decs,Population.objs);

                % Find the minimum of the surrogate
                srgtMin = FindOpt(model,Population,Problem.upper,Problem.lower);

                % Calculate distance
                dist  = pdist2(Population.decs,srgtMin);
                dxRBF = min(dist);
                if dxRBF > eta
                    optSrgt    = Problem.Evaluation(srgtMin);
                    Population = [Population,optSrgt];
                    if optSrgt.objs < Gbest(:,end)
                        [model,~] = rbf_build(Population.decs,Population.objs);
                        Gbest     = [optSrgt.decs,optSrgt.objs];
                    end
                end
                currFES = Problem.FE - K;
                [Population,Position,Velocity,Gbest,Pbest] = UpdatePosition(Problem,Population,Position,Velocity,Pbest,Gbest,currFES,maxFES,Wnc,Pr,model,eta);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function outputImage = shear(this, shear_dr)
% Geometric shear of MrImage
%
%   Y = MrImage()
%   Y.shear(inputs)
%
% This is a method of class MrImage.
%
% NOTE: This is a method of MrImage rather than MrImageGeometry, because
%       the latter is composed on the fly from affineTransformation and
%       dimInfo to integrate both information and sustain consistency.
%       Thus, the effect will only be visible using the world space plot
%       options, i.e. Y.plot('plotType', 'spmi').
% NOTE: The rotation is *added* to the rotation already defined in the
%       affine transformation matrix, equivalent to the rotation in SPM
%       Display.
% IN
%   shear_dr            [1,3] shear components, i.e., [sx, sy, sz]
%
% OUT
%
% EXAMPLE
%   shear
%
%   See also MrImage

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-12-09
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

outputImage = this.copyobj();

[outputImage.dimInfo, outputImage.affineTransformation] = ...
    this.geometry.perform_world_space_operation('shear', shear_dr, this.dimInfo);

end

author: translationalneuromodeling

repo: tapas


text:
function compiler = mexcompiler
% mexcompiler returns the name of the compiler used by mex.

% Written by Tom Minka

mexopts = fullfile(prefdir,'mexopts.bat');
if ~exist(mexopts,'file')
  compiler = '';
  return
end
fid = fopen(mexopts);
while 1
  txt = fgetl(fid);
  if ~ischar(txt), break, end
  if strmatch('set COMPILER=',txt)
    compiler = txt(14:end);
    break
  end
end
fclose(fid);

author: Cloud-CV

repo: object-proposals


text:
function determ = conex3_determinant ( n )

%*****************************************************************************80
%
%% CONEX3_DETERMINANT returns the determinant of the CONEX3 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = - 1.0;

  return
end


text:
%This Matlab script can be used to reproduce Figure 5.13 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Load SE simulation data, generated using the code from Section 4
load section5_Mvarying_K10_20;

%Number of UEs per BS
k_index = 2; %Selecting K = 10 from the loaded SE results
K = Krange(k_index);

%Fractions of data samples used for UL and DL
ULfraction = 1/3;
DLfraction = 2/3;

%Compute joint UL/DL sum SE using the fractions of UL/DL data
sumSE_MMMSE = ULfraction*sumSE_MMMSE_UL(:,k_index) + DLfraction*sumSE_MMMSE_DL(:,k_index);
sumSE_SMMSE = ULfraction*sumSE_SMMSE_UL(:,k_index) + DLfraction*sumSE_SMMSE_DL(:,k_index);
sumSE_RZF = ULfraction*sumSE_RZF_UL(:,k_index) + DLfraction*sumSE_RZF_DL(:,k_index);
sumSE_ZF = ULfraction*sumSE_ZF_UL(:,k_index) + DLfraction*sumSE_ZF_DL(:,k_index);
sumSE_MR = ULfraction*sumSE_MR_UL(:,k_index) + DLfraction*sumSE_MR_DL(:,k_index);

%Number of BSs
L = 16;

%Communication bandwidth
B = 20e6;

%PA efficiency UEs and BSs
mu_UE = 0.4;
mu_BS = 0.5;

%Define the pilot reuse factor
f = 1;

%Select length of coherence block
tau_c = 200;

%Compute length of pilot sequences
tau_p = f*K;

%Transmit power per UE in W
p = 0.1;

%Compute total effective transmit power
ETP_total =  K*p*(tau_p/mu_UE + (tau_c-tau_p)*(ULfraction/mu_UE + DLfraction/mu_BS))/tau_c;

%% Go through the two value sets of the CP model
for valueset = 1:2

    %Compute the total CP with different schemes
    [P_MR,P_RZF,P_MMMSE,P_ZF,P_SMMSE] = functionCPcomputation(Mrange,K,L,B,tau_c,tau_p,valueset,sumSE_MR,sumSE_RZF,sumSE_MMMSE,sumSE_ZF,sumSE_SMMSE);

    %Compute EE with M-MMSE
    EE_MMMSE = (B*sumSE_MMMSE)./(ETP_total + P_MMMSE);

    %Compute EE with S-MMSE
    EE_SMMSE = (B*sumSE_SMMSE)./(ETP_total + P_SMMSE);

    %Compute EE with RZF
    EE_RZF = (B*sumSE_RZF)./(ETP_total + P_RZF);

    %Compute EE with ZF
    EE_ZF = (B*sumSE_ZF)./(ETP_total + P_ZF);

    %Compute EE with MR
    EE_MR = (B*sumSE_MR)./(ETP_total + P_MR);

    %Plot simulation results
    figure;

    plot(B*sumSE_MMMSE/10^6,EE_MMMSE/10^6, 'rd-','LineWidth',1);hold on;
    plot(B*sumSE_SMMSE/10^6,EE_SMMSE/10^6,'b:','LineWidth',1);
    plot(B*sumSE_RZF/10^6,EE_RZF/10^6,'k-.','LineWidth',1);
    plot(B*sumSE_ZF/10^6,EE_ZF/10^6,'r--','LineWidth',1);
    plot(B*sumSE_MR/10^6,EE_MR/10^6,'bs-','LineWidth',1);

    xlabel('Throughput [Mbit/s/cell]');
    ylabel('EE [Mbit/Joule/cell]');
    legend('M-MMSE','S-MMSE','RZF','ZF','MR','Location','NorthWest');

end

author: emilbjornson

repo: massivemimobook


text:
function tapas_hgf_whichworld_plotTraj(r)
% Plots trajectories estimated by fitModel for the hgf_whichworld perceptual model
% Usage:  est = tapas_fitModel(responses, inputs); tapas_hgf_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Check whether we have a configuration structure
if ~isfield(r,'c_prc')
    error('tapas:hgf:ConfigRequired', 'Configuration required: before calling tapas_hgf_whichworld_plotTraj, tapas_hgf_whichworld_config has to be called.');
end

% Number of worlds
nw = r.c_prc.nw;

% Define colors
colors = [1 0 0; 0.67 0 1; 0 0.67 1; 0.67 1 0];

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name','HGF trajectories');

% Number of trials
t = size(r.u,1);

% Optional plotting of standard deviations (true or false)
plotsd2 = true;
plotsd3 = true;

% Subplots
subplot(3,1,1);

if plotsd3 == true
    upper3prior = r.p_prc.mu3_0 +sqrt(r.p_prc.sa3_0);
    lower3prior = r.p_prc.mu3_0 -sqrt(r.p_prc.sa3_0);
    upper3 = [upper3prior; r.traj.mu(:,3)+sqrt(r.traj.sa(:,3))];
    lower3 = [lower3prior; r.traj.mu(:,3)-sqrt(r.traj.sa(:,3))];

    plot(0, upper3prior, 'ob', 'LineWidth', 1);
    hold all;
    plot(0, lower3prior, 'ob', 'LineWidth', 1);
    fill([0:t, fliplr(0:t)], [(upper3)', fliplr((lower3)')], ...
         'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
plot(0:t, [r.p_prc.mu3_0; r.traj.mu(:,3)], 'b', 'LineWidth', 2);
hold all;
plot(0, r.p_prc.mu3_0, 'ob', 'LineWidth', 2); % prior
xlim([0 t]);
title('Posterior expectation \mu_3 of log-volatility of tendency x_3', 'FontWeight', 'bold');
xlabel('Trial number');
ylabel('\mu_3');

subplot(3,1,2);
if plotsd2 == true
    for j=1:nw
    upper2prior = r.p_prc.mu2_0(j) +sqrt(r.p_prc.sa2_0(j));
    lower2prior = r.p_prc.mu2_0(j) -sqrt(r.p_prc.sa2_0(j));
    upper2 = [upper2prior; r.traj.mu(:,2,j)+sqrt(r.traj.sa(:,2,j))];
    lower2 = [lower2prior; r.traj.mu(:,2,j)-sqrt(r.traj.sa(:,2,j))];

    plot(0, upper2prior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
    hold all;
    plot(0, lower2prior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
    fill([0:t, fliplr(0:t)], [(upper2)', fliplr((lower2)')], ...
         colors(j,:), 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
end
for j=1:nw
    plot(0:t, [r.p_prc.mu2_0(j); r.traj.mu(:,2,j)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu2_0(j), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
xlim([0 t]);
title('Posterior expectations \mu_2 of tendencies x_2', 'FontWeight', 'bold');
xlabel({'Trial number', ' '}); % A hack to get the relative subplot sizes right
ylabel('\mu_2');
hold off;

subplot(3,1,3);
for j=1:nw
    plot(0:t, [tapas_sgm(r.p_prc.mu2_0(j), 1); tapas_sgm(r.traj.mu(:,2,j), 1)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, tapas_sgm(r.p_prc.mu2_0(j), 1), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
plot(1:t, r.u(:,1), '.', 'Color', [0 0 0]); % inputs
if ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    y = r.y(:,1);
    if ~isempty(find(strcmp(fieldnames(r),'irr')))
        y(r.irr) = NaN; % weed out irregular responses
        plot(r.irr,  1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
    end
    if ~any(y>1) && ~any(y<0)
        plot(1:t, y, '.', 'Color', [1 0.7 0]); % responses
    else
        for j=1:nw
            plot(find(y==j), 1.08*ones([1 length(find(y==j))]), '.', 'Color', colors(j,:)); % responses
        end
    end
    title(['Response y, input u (black), and posterior probability of worlds s(\mu_2) for \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ', \vartheta=', num2str(r.p_prc.th)], ...
          'FontWeight', 'bold');
    ylabel('y, u, s(\mu_2)');
    axis([0 t -0.1 1.15]);
else
    title(['Input u (black) and posterior probability of worlds s(\mu_2) for m_3=', ...
           num2str(r.p_prc.m), ', \phi_3=', num2str(r.p_prc.phi), ', \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ', \vartheta=', num2str(r.p_prc.th)], ...
      'FontWeight', 'bold');
    ylabel('u, s(\mu_2)');
    axis([0 t -0.1 1.1]);
end
plot(1:t, 0.5, 'k');
xlabel('Trial number');
hold off;


text:
classdef TestLinearPolar
    %TestLinearPolar

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
            dst = cv.linearPolar(img, [100 100], 500);
            validateattributes(dst, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.linearPolar();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function nuclides
% Solution for a chain of radionuclides
%    using MATLAB expm
%   $Ekkehard Holzbecher  $Date: 2006/28/12 $
%--------------------------------------------------------------------------

T = 10;                % maximum time
lambda = [1; 0.1; 0.5];% decay rates
c0 = [1; 0; 0];        % initial concentrations
q = [0.1; 0; 0];       % source rates
N = 60;                % discretization of time

t = linspace (0,T,N);
B = -diag(lambda);
for i = 2:size(lambda,1)
    B(i,i-1) = lambda(i-1);
end
c = c0;

for i = 2:N
    E = expm(B*t(i));
    c = [c E*c0-(eye(size(B,1))-E)*inv(B)*q];
end
plot (t,c');
legend ('mother','daughter 1','daughter 2');
text (T/2,0.8,'Steady state:'); text (T/2,0.7,num2str(-(inv(B)*q)'));
xlabel ('time');

author: Sable

repo: mcbench-benchmarks


text:
% GP_COV_DAMPEDCOS - Damped cosine covariance function.
%
% c(r) = exp(-r/lengthscale) * cos( 2*pi*r/wavelength )
%
% COVFUNC = GP_COV_DAMPEDCOS(D2)
%
% [K, DK] = COVFUNC(THETA)
%
% [N_THETA, N1, N2] = COVFUNC()

function covfunc = gp_cov_dampedcos(D, varargin)

options = struct('wavelength', [],...
                 'lengthscale', [] );
[options, errmsg] = argparse(options, varargin{:});
error(errmsg);

n_theta = isempty(options.wavelength) + isempty(options.lengthscale);

covfunc = @get_covariance;

  function varargout = get_covariance(theta)

  varargout = cell(nargout,1);

  if nargin == 0
    % Return the number of parameters and the size of the covariance matrix
    if nargout >= 1
      varargout{1} = n_theta; % number of parameters
      if nargout >= 2
        varargout{2} = size(D,1); % number of rows
        if nargout == 3
          varargout{3} = size(D,2); % number of columns
        else
          error('Too many outputs');
        end
      end
    end
    return
  end

  if numel(theta) ~= n_theta
    error('Wrong number of parameters (%d), should be %d', numel(theta), n_theta);
  end

  % Parse parameters
  if isempty(options.wavelength)
      wl = theta(1);
      if isempty(options.lengthscale)
          ls = theta(2);
      else
          ls = options.lengthscale;
      end
  else
      wl = options.wavelength;
      if isempty(options.lengthscale)
          ls = theta(1);
      else
          ls = options.lengthscale;
      end
  end

  % Covariance matrix
  K = exp(-D/ls) .* cos( 2*pi*D/wl );
  varargout{1} = K;

  % Gradient for hyperparameters
  if nargout >= 2
      dK = cell(n_theta,1);
      n = 1;
      if isempty(options.wavelength)
          dK{n} = exp(-D/ls) .* sin( 2*pi*D/wl ) .* ( 2*pi*D*wl^(-2) );
          n = n + 1;
      end
      if isempty(options.lengthscale)
          dK{n} = exp(-D/ls) .* ( D*ls^(-2) ) .* cos( 2*pi*D/wl );
      end
      varargout{2} = dK;
  end

  end % function get_covariance

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [y,deriv] = logdet_chol(w)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
% y = log(det(W)), where W is positive definite and W = reshape(w,...)

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)logdet_chol(w);
    return;
end

if isa(w,'function_handle')
    outer = logdet_chol([]);
    y = compose_mv(outer,w,[]);
    return;
end

dim = sqrt(length(w));
W = reshape(w,dim,dim);

if nargout>1
    %[inv_map,bi_inv_map,logdet,iW] = invchol2(W);
    [inv_map,bi_inv_map,logdet,iW] = invchol_or_lu(W);
    y = logdet;
    deriv = @(dy) deriv_this(dy,bi_inv_map,iW);
else
    %[inv_map,bi_inv_map,logdet] = invchol2(W);
    [inv_map,bi_inv_map,logdet] = invchol_or_lu(W);
    y = logdet;
end

function [g,hess,linear] = deriv_this(dy,bi_inv_map,iW)
G = iW.';
grad = G(:);
g = dy*grad;
linear = false;
hess = @(d) hess_this(grad,bi_inv_map,dy,d);

function [h,Jd] = hess_this(grad,bi_inv_map,dy,d)
dim = sqrt(length(d));
D = reshape(d,dim,dim);
H = - dy*bi_inv_map(D).';
h = H(:);
if nargout>1
    Jd = grad.'*d(:);
end

function test_this()
m = 3;
n = 10;

w = [];
A = UtU(w,n,m);

f = logdet_chol(A);
w = randn(m*n,1);

test_MV2DF(f,w,true);

author: nesl


text:
%% Transient diffusion equation evaluated using FVTool
%
% adapted from A.A. Eftekhari by M.H.V. Werts (2020)
%
% GNU Octave 4.2.2 was used for development
%
% This script is intended to be run from the command line interface.
% It is called in this manner from within the accompanying Jupyter
% Python notebook.
%
% It should be inside the 'FVTool' directory tree as downloaded/cloned
% from Github, under
%   './Examples/External/Diffusion1DSpherical_Analytic-vs-FVTool-vs-Fipy'
%
% Script calculates diffusion in a 1D spherical geometry for an 'infinite'
% medium, with the initial condition that all mass at $t = 0$ is
% homogeneously confined inside a sphere of radius $a$. This is sometimes
% called a 'spherical initial condition'.
%
% see J. Crank (1975) "The Mathematics of Diffusion", 2nd Ed.,
%      Clarendon Press (Oxford), pages 29-30
%      Equation 3.8, Figure 3.1
%
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc),
% $D$ is the diffusion coefficient, and $\alpha$ is a constant (1, here).
%
%
clc; clear;
more off;
run('../../../FVToolStartUp.m')

%% Define the domain and create a mesh structure
% Here we work in a 1D spherical coordinate system (r coordinate)
L = 10.0;  % domain length
Nx = 2000; % number of cells
m = createMeshSpherical1D(Nx, L);

%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure

%% define the transfer coeffs
D = createCellVariable(m, 1.0);
alfa = createCellVariable(m, 1.0);

%% define initial condition
c_init = 0;
c_old = createCellVariable(m, c_init, BC); % initial values
r = c_old.domain.cellcenters.x;
c_old.value(r<1.0) = 1.0;

%% calculate volumes of FV cellslices
%  We use this for demonstrating mass conservation
cellA = m.facecenters.x(1:end-1);
cellB = m.facecenters.x(2:end);
cellvol = 4/3 .* pi .* (cellB.^3 - cellA.^3);
cellsum = sum(cellvol)

c = c_old; % assign the old value of the cells to the current values

t = 0.0; % master time
deltat = 0.0625/20; % time step

% output total mass in the system
m_tot = sum(c.value(2:end-1) .* cellvol);
t,m_tot

%% loop for "time-stepping" the solution
% It outputs the spatial profile C(r) after
% 20, 80 and 320 time-steps
% This corresponds to t=0.0625, t=0.25 and t=1, respectively.
ti = 0
for s=[20,60,240]
  for n=1:s
      [M_trans, RHS_trans] = transientTerm(c, deltat, alfa);
      Dave = harmonicMean(D);
      Mdiff = diffusionTerm(Dave);
      [Mbc, RHSbc] = boundaryCondition(BC);
      M = M_trans-Mdiff+Mbc;
      RHS = RHS_trans+RHSbc;
      c = solvePDE(m,M, RHS);
      t += deltat;
      c_old = c;
  endfor
  m_tot = sum(c.value(2:end-1) .* cellvol);
  n,t,m_tot
  % The following writes the result to a file
  %  adapted from visualizeCells with domain.dimension = 1.8
  x = [c.domain.facecenters.x(1); c.domain.cellcenters.x; c.domain.facecenters.x(end)];
  cval = [0.5*(c.value(1)+c.value(2)); c.value(2:end-1); 0.5*(c.value(end-1)+c.value(end))];
  ti += s;
  filename = ["diffusion1Dspherical_FVTool_tstep",num2str(ti),".mat"]
  save('-6',filename,'x','cval');
endfor

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function signedcos = acos2d(alpha, beta)
%ACOS2D      4-quadrant arccosine function, in degrees.
%
%   ACOS2D(alpha, beta) computes the four-quadrant arccosine of the amgle
%   [alpha]. For arguments |alpha| > 1, the result is NaN. The resulting
%   angle is not uniquely determined by alpha, nor by the lengths or
%   order of the sides of the triangle (as in ATAN2), so an additional
%   argument [beta] is required. If [beta] < 180�, the small angle
%   (0� <= alpha <= 180�) is returned. If [beta] > 180�, the large angle
%   (180� < alpha < 360�) is returned.
%
%   See also acos2.

% Rody P.S. Oldenhuis
% Delft University of Technology
% oldenhuis@gmail.com
%
% Crated     : 23/Feb/2009
% Last edited: 30/Nov/2012

    % compute 4-quadrant cosine by calling ACOS2 directly
    signedcos = acos2(alpha, beta)*180/pi;

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% DOF
%
% Files
%   dof3BDM1  - dof structure for BDM1 element in 3-D.
%   dof3edge  - dof structure for the lowest order edge elements.
%   dof3face  - dof sturctrue for the lowest order Raviart-Thomas element in 3-D.
%   dof3P2    - dof structure for P2 element in 3-D.
%   dof3RT0   - dof sturctrue for the lowest order face element in 3-D.
%   dofBDM1   - dof structure for BDM1 element
%   dofedge   - dof structure for edges.
%   dofP2     - dof structure for P2 element.
%   dofRT0    - dof structure for the lowest order Ravairt-Thomas element.
%   sortelem3 - sort elem in ascend ordering

repo: ifem


text:
function H=calcSpherConvHessian(zSpher,systemType,useHalfRange,lTx,lRx,M)
%%CALCSPHERCONVHESSIAN Calculate the Hessian (a matrix of second partial
%            derivatives) for a monostatic or bistatic range and a
%            spherical direction measurement in 3D, ignoring atmospheric
%            effects, with respect to Cartesian position. This type of
%            Jacobian is useful when performing tracking using Cartesian-
%            converted measurements where the clutter density is specified
%            in the measurement coordinate system, not the converted
%            measurement coordinate system.
%
%INPUTS: zSpher A 3X1 point in range and azimuth and angle in the format
%           [range;azimuth;angle], where the angles are given in radians.
% systemType An optional parameter specifying the axes from which the
%           angles are measured in radians. Possible values are
%           0 (The default if omitted) Azimuth is measured
%             counterclockwise from the x-axis in the x-y plane. Elevation
%             is measured up from the x-y plane (towards the z-axis). This
%             is consistent with common spherical coordinate systems for
%             specifying longitude (azimuth) and geocentric latitude
%             (elevation).
%           1 Azimuth is measured counterclockwise from the z-axis in the
%             z-x plane. Elevation is measured up from the z-x plane
%             (towards the y-axis). This is consistent with some spherical
%             coordinate systems that use the z axis as the boresight
%             direction of the radar.
%           2 This is the same as 0 except instead of being given
%             elevation, one desires the angle away from the z-axis, which
%             is (pi/2-elevation).
%           3 This is the same as 0 except azimuth is measured clockwise
%             from the y-axis in the x-y plane instead of counterclockwise
%             from the x-axis. This coordinate system often arises when
%             given "bearings" in a local East-North-Up coordinate system,
%             where the bearing directions are measured East of North.
% useHalfRange An optional boolean value specifying whether the bistatic
%           (round-trip) range value has been divided by two. This normally
%           comes up when operating in monostatic mode (the most common
%           type of spherical coordinate system), so that the range
%           reported is a one-way range (or just half a bistatic range).
%           The default if this parameter is not provided is false if lTx
%           is provided and true if it is are omitted (monostatic). If no
%           range values are provided, an empty matrix can be passed.
%       lTx The 3X1 [x;y;z] location vector of the transmitter in global
%           Cartesian coordinates. If this parameter is omitted or an
%           empty matrix is passed, then the transmitter is assumed to be
%           at the origin.
%       lRx The 3X1 [x;y;z] location vector of the receiver in Cartesian
%           coordinates. If this parameter is omitted or an empty matrix
%           is passed, then the receiver is assumed to be at the origin.
%         M A 3X3 rotation matrix to go from the alignment of the global
%           coordinate system to that at the receiver. The z-axis of the
%           local coordinate system of the receiver is the pointing
%           direction of the receiver. If omitted or an empty matrix is
%           passed, then it is assumed that the local coordinate system is
%           aligned with the global and M=eye(3) --the identity matrix is
%           used.
%
%OUTPUTS: H The 3X3X3 Hessian matrix, where H(:,:,1) is the Hessian for the
%           range component, H(:,:,2) is the Hessian for the azimuth
%           component, and H(:,:,3) is the Hessian for the elevation
%           component. The ordering of the derivatives in each matrix is:
%                  [d^2/(dxdx), d^2/(dxdy), d^2/(dxdz);
%                   d^2/(dydx), d^2/(dydy), d^2/(dydz);
%                   d^2/(dzdx), d^2/(dzdy), d^2/(dzdz)];
%          note that each matrix is symmetric (i.e.
%                   d^2/(dydx)=d^2/(dxdy) ).
%
%This function just calls spher2Cart and then calcSpherHessian.
%
%June 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(M))
	M=eye(3,3);
end

if(nargin<5)
   lRx=[];
end

if(nargin<4)
    lTx=[];
end

if(isempty(lTx)&&(nargin<3||isempty(useHalfRange)))
    useHalfRange=true;
elseif(~isempty(lTx)&&(nargin<3||isempty(useHalfRange)))
    useHalfRange=false;
end

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

x=spher2Cart(zSpher,systemType,useHalfRange,lTx,lRx,M);
H=calcSpherHessian(x,systemType,useHalfRange,lTx,lRx,M);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function adapted = adapted_ranks(ranks, different, mode)
% adapted_ranks Performs rank adaptation on a set of ranks
%
% Performs different types of rank adaptation on a set of ranks and a matrix of determined actual differences of ranked entities.
%
% The function looks for sets of entities that are considered equal according to the `different` matrix and adapts their ranks
% in one of the following ways:
% - none: no adaptation, keep the original ranks.
% - mean: a mean rank of a set is assigned to all trackers in the set.
% - median: a median rank of a set is assigned to all trackers in the set.
% - best: a minimum rank in a set is assigned to all trackers in the set.
%
% Input:
% - ranks (double vector): A vector of `N` ranks.
% - different (boolean matrix): A `N x N` matrix that denotes which entities should actually be equal and which not.
% - mode (string): Type of rank adaptation.
%
% Output
% - adapted (double vector): Adapted ranks

adapted = zeros(1, length(ranks)) ;

switch mode
	case 'none'
        adapted = ranks;

	case 'mean'

		for tracker = 1:length(ranks)
			adapted(tracker) = mean(ranks(find(~different(tracker,:)))) ; %#ok<FNDSB>
		end

	case 'median'

		for tracker = 1:length(ranks)
			adapted(tracker) = median(ranks(find(~different(tracker,:)))) ; %#ok<FNDSB>
		end

	case 'best'

        [sorted_ranks, idx] = sort(ranks);
        [~, unsort_idx] = sort(idx);

        sorted_matrix = different(idx,idx);

        for tracker = 2:length(ranks)
            if sorted_matrix(tracker,tracker-1)
                sorted_ranks(tracker) = sorted_ranks(tracker-1);
            end;
        end;
        adapted = sorted_ranks(unsort_idx);

	%	for tracker = 1:length(ranks)
	%		adapted(tracker) = min(ranks(find(~different(tracker,:)))) ; %#ok<FNDSB>
	%	end

end;

author: votchallenge

repo: toolkit-legacy


text:
function marker = get_markers

marker = { '-*k', '-or', '-*b', '-ob', '-og', '-+m', '-*r', '-+g', '', '', '-*k', '-or', '-+b' };

end

author: andrewssobral

repo: lrslibrary


%% NISwGSP-parrington
% %{
imfolder = 'images\NISwGSP-parrington';
im_n = 18;
imfile = cell(im_n,1);
for ii = 1:im_n
    imfile{ii} = sprintf('%s\\prtn%02d.jpg', imfolder, ii-1);
end

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = zeros(im_n,2);
ei = 0;
for ii = 1:im_n-1
    ei = ei + 1;
    edge_list(ei,:) = [ii,ii+1];
end
edge_list(im_n,:) = [im_n,1];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 0, 'equi', 0.02, imfolder );
%}


text:
% File: RecognizeActions.m
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function [accuracy, predicted_labels] = RecognizeActions(datasetTrain, datasetTest, G, maxIter)

% INPUTS
% datasetTrain: dataset for training models, see PA for details
% datasetTest: dataset for testing models, see PA for details
% G: graph parameterization as explained in PA decription
% maxIter: max number of iterations to run for EM

% OUTPUTS
% accuracy: recognition accuracy, defined as (#correctly classified examples / #total examples)
% predicted_labels: N x 1 vector with the predicted labels for each of the instances in datasetTest, with N being the number of unknown test instances

% Train a model for each action
% Note that all actions share the same graph parameterization and number of max iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[P1 logli1 ClassProb1 PairProb1] = EM_HMM(datasetTrain(1).actionData, datasetTrain(1).poseData, G, datasetTrain(1).InitialClassProb, datasetTrain(1).InitialPairProb, maxIter);

[P2 logli2 ClassProb2 PairProb2] = EM_HMM(datasetTrain(2).actionData, datasetTrain(2).poseData, G, datasetTrain(2).InitialClassProb, datasetTrain(2).InitialPairProb, maxIter);

[P3 logli3 ClassProb3 PairProb3] = EM_HMM(datasetTrain(3).actionData, datasetTrain(3).poseData, G, datasetTrain(3).InitialClassProb, datasetTrain(3).InitialPairProb, maxIter);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Classify each of the instances in datasetTest
% Compute and return the predicted labels and accuracy
% Accuracy is defined as (#correctly classified examples / #total examples)
% Note that all actions share the same graph parameterization

accuracy = 0;
predicted_labels = [];
P = [P1,P2,P3];
N = size(datasetTest.poseData,1);
I = size(datasetTest.labels,1);
K = 3;
loglikelihood = zeros(I,3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for iclass = 1:3
  logEmissionProb = zeros(N,K);

  for i = 1:N
	  data = reshape(datasetTest.poseData(i,:,:),10,3);
	  for j = 1:10
		  if G(j,1) == 1
			  parent = data(G(j,2),:);
			  for k = 1:K
				  theta = P(iclass).clg(j).theta(k,:);
				  mu_y = sum(theta(1:4).*[1,parent]);
				  mu_x = sum(theta(5:8).*[1,parent]);
				  mu_a = sum(theta(9:12).*[1,parent]);
				  logEmissionProb(i,k) += lognormpdf(data(j,1),mu_y,P(iclass).clg(j).sigma_y(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,2),mu_x,P(iclass).clg(j).sigma_x(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,3),mu_a,P(iclass).clg(j).sigma_angle(k));
			  end
		  else
			  for k = 1:K
				  logEmissionProb(i,k) += lognormpdf(data(j,1),P(iclass).clg(j).mu_y(k),P(iclass).clg(j).sigma_y(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,2),P(iclass).clg(j).mu_x(k),P(iclass).clg(j).sigma_x(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,3),P(iclass).clg(j).mu_angle(k),P(iclass).clg(j).sigma_angle(k));
			  end
		  end
	  end
  end

  for i = 1:I
	  % construct all the three types of factors for each action and do inference to fill ClassProb and PairProb
	  m = length(datasetTest.actionData(i).marg_ind); % 1 to m represents S variables
	  F = repmat(struct('var',[],'card',[],'var',[]),1,2*m); % 1 is P(S1), 2 to m is P(Si/Si-1) and m+1 to 2m P(S/O)
	  F(1).var = 1; F(1).card = [K];F(1).val = log(P(iclass).c);
	  temp = log(reshape(P(iclass).transMatrix',1,9));
	  for j = 2:m
		  F(j).var = [j j-1]; F(j).card = [K K]; F(j).val = temp;
	  end
	  for j = 1:m
		  F(j+m).var = [j];F(j+m).card = [K];
		  F(j+m).val = logEmissionProb(datasetTest.actionData(i).marg_ind(j),:);
	  end
	  [M, PCalibrated] = ComputeExactMarginalsHMM(F);
	  loglikelihood(i,iclass) += logsumexp(PCalibrated.cliqueList(1).val);
  end
end
[temp, predicted_labels] = max(loglikelihood,[],2);
accuracy = sum(predicted_labels==datasetTest.labels)/I;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: anhncs

repo: Probabilistic-Graphical-Models


text:
%SAMMONM Sammon mapping
%
%   W = SAMMONM(A,K,MAX)
%   W = A*SAMMONM([],K,MAX)
%   W = A*SAMMONM(K,MAX)
%   D = B*W
%
% INPUT
%   A    Dataset, used for training the mapping
%   B    Dataset, same dimensionality as A, to be mapped
%   K    Target dimension of mapping (default 2)
%   MAX  Maximum number of iterations, default 100
%
% OUTPUT
%   W    Trained mapping
%   D    K-dimensional dataset
%
% DESCRIPTION
% This is a simplified interface to the more complex MDS routine for high
% dimensional data visualisation. The output is a non-linear projection of
% the original vector space to a K-dimensional target space.
%
% The main differences with MDS are that SAMMONM operates on feature based
% datasets, while MDS expects dissimilarity matrices; MDS maps new objects
% by a second optimization procedures minimizing the stress for the test
% objects, while SAMMONM uses a linear mapping between dissimilarities and
% the target space. See also PREX_MDS for examples. A different procedure
% for the same purpose is TSNEM.
%
% EXAMPLE
% prdatasets;            % make sure prdatasets is in the path
% a = satellite;         % 36D dataset, 6 classes, 6435 objects
% [x,y] = gendat(a,0.5); % split in train and test set
% w = x*sammonm;         % compute mapping
% figure; scattern(x*w); % show trainset mapped to 2D: somewhat overtrained
% figure; scattern((x+randn(size(x))*1e-5)*w): % some noise helps
% figure; scattern(y*w); % show test set mapped to 2D
%
% REFERENCES
% 1. JW Sammon Jr A nonlinear mapping for data structure analysis,
%   IEEE Transactions on Computers C-18, pp. 401-409,1969.
% 2. E. Pekalska, D. de Ridder, R.P.W. Duin, and M.A. Kraaijveld, A new
%    method of generalizing Sammon mapping with application to algorithm
%    speed-up, ASCI99, Proc. 5th Annual ASCI Conf., 1999, 221-228. [pdf]
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, PCAM, MDS, TSNEM, PREX_MDS, SCATTERD, SCATTERN

% Copyright: E. Pekalska, R.P.W. Duin, r.p.w.duin@37steps.com

function out = sammonm(varargin)

  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],2,100,0.01);
  if mapping_task(argin,'definition')
    out = define_mapping(argin,'untrained');
  elseif mapping_task(argin,'training')
    [a,k,max_iter,pow] = deal(argin{:});
    d = sqrt(distm(+a));
    opt.maxiter = max_iter;
    w = mds(d,a*pcam(a,k),opt);
    y = getdata(w,1);
    % refer to distances with a low power to get some smoothness
    v = prpinv(d.^pow)*y;
    out = trained_mapping(a,{v,a,pow},2);
  elseif mapping_task(argin,'trained execution')
    [b,w] = deal(argin{1:2});
    [v,a,pow] = getdata(w);
    out = setdata(b,(sqrt(distm(b,a)).^pow)*v);
  else
    error('Illegal call');
  end

return

author: marianux

repo: ecg-kit

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