text:
function [data, mmap] = loadjson(fname,varargin)
%
% data=loadjson(fname,opt)
%    or
% [data, mmap]=loadjson(fname,'param1',value1,'param2',value2,...)
%
% parse a JSON (JavaScript Object Notation) file or string and return a
% matlab data structure with optional memory-map (mmap) table
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% created on 2011/09/09, including previous works from
%
%         Nedialko Krouchev: http://www.mathworks.com/matlabcentral/fileexchange/25713
%            created on 2009/11/02
%         François Glineur: http://www.mathworks.com/matlabcentral/fileexchange/23393
%            created on  2009/03/22
%         Joel Feenstra:
%         http://www.mathworks.com/matlabcentral/fileexchange/20565
%            created on 2008/07/03
%
% input:
%      fname: input file name; if fname contains "{}" or "[]", fname
%             will be interpreted as a JSON string
%      opt: (optional) a struct to store parsing options, opt can be replaced by
%           a list of ('param',value) pairs - the param string is equivallent
%           to a field in opt. opt can have the following
%           fields (first in [.|.] is the default)
%
%           SimplifyCell [1|0]: if set to 1, loadjson will call cell2mat
%                         for each element of the JSON data, and group
%                         arrays based on the cell2mat rules.
%           FastArrayParser [1|0 or integer]: if set to 1, use a
%                         speed-optimized array parser when loading an
%                         array object. The fast array parser may
%                         collapse block arrays into a single large
%                         array similar to rules defined in cell2mat; 0 to
%                         use a legacy parser; if set to a larger-than-1
%                         value, this option will specify the minimum
%                         dimension to enable the fast array parser. For
%                         example, if the input is a 3D array, setting
%                         FastArrayParser to 1 will return a 3D array;
%                         setting to 2 will return a cell array of 2D
%                         arrays; setting to 3 will return to a 2D cell
%                         array of 1D vectors; setting to 4 will return a
%                         3D cell array.
%           UseMap [0|1]: if set to 1, loadjson uses a containers.Map to
%                         store map objects; otherwise use a struct object
%           ShowProgress [0|1]: if set to 1, loadjson displays a progress bar.
%           ParseStringArray [0|1]: if set to 0, loadjson converts "string arrays"
%                         (introduced in MATLAB R2016b) to char arrays; if set to 1,
%                         loadjson skips this conversion.
%           FormatVersion [3|float]: set the JSONLab format version; since
%                         v2.0, JSONLab uses JData specification Draft 1
%                         for output format, it is incompatible with all
%                         previous releases; if old output is desired,
%                         please set FormatVersion to 1.9 or earlier.
%           Encoding ['']: json file encoding. Support all encodings of
%                         fopen() function
%           ObjectID [0|interger or list]: if set to a positive number,
%                         it returns the specified JSON object by index
%                         in a multi-JSON document; if set to a vector,
%                         it returns a list of specified objects.
%           JDataDecode [1|0]: if set to 1, call jdatadecode to decode
%                         JData structures defined in the JData
%                         Specification.
%           BuiltinJSON [0|1]: if set to 1, this function attempts to call
%                         jsondecode, if presents (MATLAB R2016b or Octave
%                         6) first. If jsondecode does not exist or failed,
%                         this function falls back to the jsonlab parser
%           MmapOnly [0|1]: if set to 1, this function only returns mmap
%           MMapInclude 'str1' or  {'str1','str2',..}: if provided, the
%                         returned mmap will be filtered by only keeping
%                         entries containing any one of the string patterns
%                         provided in a cell
%           MMapExclude 'str1' or  {'str1','str2',..}: if provided, the
%                         returned mmap will be filtered by removing
%                         entries containing any one of the string patterns
%                         provided in a cell
%
% output:
%      dat: a cell array, where {...} blocks are converted into cell arrays,
%           and [...] are converted to arrays
%      mmap: (optional) a cell array as memory-mapping table in the form of
%             {{jsonpath1,[start,length,<whitespace_pre>]},
%              {jsonpath2,[start,length,<whitespace_pre>]}, ...}
%           where jsonpath_i is a string in the JSONPath [1,2] format, and
%           "start" is an integer referring to the offset from the begining
%           of the stream, and "length" is the JSON object string length.
%           An optional 3rd integer "whitespace_pre" may appear to record
%           the preceding whitespace length in case expansion of the data
%           record is needed when using the mmap.
%
%           The format of the mmap table retruned from this function
%           follows the JSON-Mmap Specification Draft 1 [3] defined by the
%           NeuroJSON project, see https://neurojson.org/jsonmmap/draft1/
%
%           Memory-mapping table (mmap) is useful when fast reading/writing
%           specific data records inside a large JSON file without needing
%           to load/parse/overwrite the entire file.
%
%           The JSONPath keys used in mmap is largely compatible to the
%           upstream specification defined in [1], with a slight extension
%           to handle contatenated JSON files.
%
%           In the mmap jsonpath key, a '$' denotes the root object, a '.'
%           denotes a child of the preceding element; '.key' points to the
%           value segment of the child named "key" of the preceding
%           object; '[i]' denotes the (i+1)th member of the preceding
%           element, which must be an array. For example, a key
%
%           $.obj1.obj2[0].obj3
%
%           defines the memory-map of the "value" section in the below
%           hierarchy:
%             {
%                "obj1":{
%            	     "obj2":[
%                       {"obj3":value},
%                       ...
%                    ],
%                    ...
%                 }
%             }
%           Please note that "value" can be any valid JSON value, including
%           an array, an object, a string or numerical value.
%
%           To handle concatenated JSON objects (including ndjson,
%           http://ndjson.org/), such as
%
%             {"root1": {"obj1": ...}}
%             ["root2", value1, value2, {"obj2": ...}]
%             {"root3": ...}
%
%           we use '$' or '$0' for the first root-object, and '$1' refers
%           to the 2nd root object (["root2",...]) and '$2' referrs to the
%           3rd root object, and so on. Please note that this syntax is an
%           extension from the JSONPath documentation [1,2]
%
%           [1] https://goessner.net/articles/JsonPath/
%           [2] http://jsonpath.herokuapp.com/
%           [3] https://neurojson.org/jsonmmap/draft1/
%
% examples:
%      dat=loadjson('{"obj":{"string":"value","array":[1,2,3]}}')
%      dat=loadjson(['examples' filesep 'example1.json'])
%      [dat, mmap]=loadjson(['examples' filesep 'example1.json'],'SimplifyCell',0)
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

    opt=varargin2struct(varargin{:});

    if(regexp(fname,'^\s*(?:\[.*\])|(?:\{.*\})\s*$','once'))
       string=fname;
    elseif(exist(fname,'file'))
       try
           encoding = jsonopt('Encoding','',opt);
           if(isempty(encoding))
               string = fileread(fname);
           else
               fid = fopen(fname,'r','n',encoding);
               string = fread(fid,'*char')';
               fclose(fid);
           end
       catch
           try
               string = urlread(fname);
           catch
               string = urlread(['file://',fullfile(pwd,fname)]);
           end
       end
    elseif(regexpi(fname,'^\s*(http|https|ftp|file)://'))
       string = urlread(fname);
    else
       error_pos('input file does not exist');
    end

    if(jsonopt('BuiltinJSON',0,opt) && exist('jsondecode','builtin'))
        try
            newstring=regexprep(string,'[\r\n]','');
            newdata=jsondecode(newstring);
            newdata=jdatadecode(newdata,'Base64',1,'Recursive',1,varargin{:});
            data=newdata;
            return;
        catch
            warning('built-in jsondecode function failed to parse the file, fallback to loadjson');
        end
    end

    pos = 1; inputlen = length(string); inputstr = string;
    arraytokenidx=find(inputstr=='[' | inputstr==']');
    arraytoken=inputstr(arraytokenidx);

    % String delimiters and escape chars identified to improve speed:
    esc = find(inputstr=='"' | inputstr=='\' ); % comparable to: regexp(inputstr, '["\\]');
    index_esc = 1;

    opt.arraytoken_=arraytoken;
    opt.arraytokenidx_=arraytokenidx;
    opt.simplifycell=jsonopt('SimplifyCell',1,opt);
    opt.simplifycellarray=jsonopt('SimplifyCellArray',opt.simplifycell,opt);
    opt.formatversion=jsonopt('FormatVersion',3,opt);
    opt.fastarrayparser=jsonopt('FastArrayParser',1,opt);
    opt.parsestringarray=jsonopt('ParseStringArray',0,opt);
    opt.usemap=jsonopt('UseMap',0,opt);
    opt.arraydepth_=1;
    opt.mmaponly=jsonopt('MmapOnly',0,opt);

    if(jsonopt('ShowProgress',0,opt)==1)
        opt.progressbar_=waitbar(0,'loading ...');
    end

    objid=jsonopt('ObjectID',0,opt);
    maxobjid=max(objid);
    if(maxobjid==0)
        maxobjid=inf;
    end
    opt.jsonpath_='$';
    if(nargout>1 || opt.mmaponly)
        mmap={};
    end
    jsoncount=1;
    while pos <= inputlen
        [cc,pos,w1]=next_char(inputstr, pos);
        switch(cc)
            case '{'
                if(nargout>1 || opt.mmaponly)
                    mmap{end+1}={opt.jsonpath_,[pos, 0, w1]};
                    [data{jsoncount},pos,index_esc,newmmap] = parse_object(inputstr, pos, esc, index_esc,opt);
                    mmap{end}{2}(2)=pos-mmap{end}{2}(1);
                    mmap=[mmap(:);newmmap(:)];
                else
                    [data{jsoncount},pos,index_esc] = parse_object(inputstr, pos, esc, index_esc,opt);
                end
            case '['
                if(nargout>1 || opt.mmaponly)
                    mmap{end+1}={opt.jsonpath_,[pos,0,w1]};
                    [data{jsoncount},pos,index_esc,newmmap] = parse_array(inputstr, pos, esc, index_esc,opt);
                    mmap{end}{2}(2)=pos-mmap{end}{2}(1);
                    mmap=[mmap(:);newmmap(:)];
                else
                    [data{jsoncount},pos,index_esc] = parse_array(inputstr, pos, esc, index_esc,opt);
                end
            otherwise
                pos=error_pos('Outer level structure must be an object or an array',inputstr,pos);
        end
	    if(jsoncount>=maxobjid)
	        break;
	    end
        opt.jsonpath_=sprintf('$%d',jsoncount);
        jsoncount=jsoncount+1;
    end % while

    if(length(objid)>1 || min(objid)>1)
        data=data(objid(objid<=length(data)));
    end

    jsoncount=length(data);
    if(jsoncount==1 && iscell(data))
        data=data{1};
    end
    if(nargout>1 || opt.mmaponly)
        mmap=mmap';
        mmap=filterjsonmmap(mmap, jsonopt('MMapExclude',{},opt), 0);
        mmap=filterjsonmmap(mmap, jsonopt('MMapInclude',{},opt), 1);
        mmap=cellfun(@(x) {x{1},x{2}(1:(2+int8(length(x{2})>=3 && (x{2}(3)>0))))}, mmap, 'UniformOutput', false);
    end
    if(jsonopt('JDataDecode',1,varargin{:})==1)
        try
            data=jdatadecode(data,'Base64',1,'Recursive',1,varargin{:});
        catch ME
            warning(['Failed to decode embedded JData annotations, '...
                'return raw JSON data\n\njdatadecode error: %s\n%s\nCall stack:\n%s\n'], ...
                ME.identifier, ME.message, savejson('',ME.stack));
        end
    end
    if(opt.mmaponly)
        data=mmap;
    end
    if(isfield(opt,'progressbar_'))
        close(opt.progressbar_);
    end
end

%%-------------------------------------------------------------------------
%% helper functions
%%-------------------------------------------------------------------------

function [object, pos,index_esc, mmap] = parse_array(inputstr, pos, esc, index_esc, varargin) % JSON array is written in row-major order
    if(nargout>3)
        mmap={};
        origpath=varargin{1}.jsonpath_;
    end
    pos=parse_char(inputstr, pos, '[');
    object = cell(0, 1);
    arraydepth=varargin{1}.arraydepth_;
    pbar=-1;
    if(isfield(varargin{1},'progressbar_'))
        pbar=varargin{1}.progressbar_;
    end
    format=varargin{1}.formatversion;
    [cc,pos]=next_char(inputstr,pos);
    endpos=[];

    if cc ~= ']'
        try
            if((varargin{1}.fastarrayparser)>=1 && arraydepth>=varargin{1}.fastarrayparser)
                [endpos, maxlevel]=fast_match_bracket(varargin{1}.arraytoken_,varargin{1}.arraytokenidx_,pos);
                if(~isempty(endpos))
                    arraystr=['[' inputstr(pos:endpos)];
                    arraystr=sscanf_prep(arraystr);
                    if(isempty(find(arraystr=='"', 1)))
                        % handle 1D array first
                        if(maxlevel==1)
                            astr=arraystr(2:end-1);
                            astr(astr==' ')='';
                            [obj, count, errmsg, nextidx]=sscanf(astr,'%f,',[1,inf]);
                            if(nextidx>=length(astr)-1)
                                    object=obj;
                                    pos=endpos;
                                    pos=parse_char(inputstr, pos, ']');
                                    return;
                            end
                        end

                        % for N-D packed array in a nested array construct,
                        if(maxlevel>=2 && ~isempty(regexp(arraystr(2:end),'^\s*\[','once')))
                            [dims,isndarray]=nestbracket2dim(arraystr);
                            rowstart=find(arraystr(2:end)=='[',1)+1;
                            if(rowstart && isndarray)
                                [obj, nextidx]=parsendarray(arraystr,dims);
                                if(nextidx>=length(arraystr)-1)
                                    object=obj;
                                    if(format>1.9)
                                        object=permute(object,ndims(object):-1:1);
                                    end
                                    pos=endpos;
                                    pos=parse_char(inputstr, pos, ']');
                                    if(pbar>0)
                                        waitbar(pos/length(inputstr),pbar,'loading ...');
                                    end
                                    return;
                                end
                            end
                        end
                    end
                end
            end
            if(isempty(regexp(arraystr,':','once')) && isempty(regexp(arraystr,'\(','once')))
                arraystr=regexprep(arraystr,'\[','{');
                arraystr=regexprep(arraystr,'\]','}');
                if(varargin{1}.parsestringarray==0)
                    arraystr=regexprep(arraystr,'\"','''');
                end
                object=eval(arraystr);
                if(iscell(object))
                    object=cellfun(@unescapejsonstring,object,'UniformOutput',false);
                end
                pos=endpos;
            end
        catch
        end
        if(isempty(endpos) || pos~=endpos)
            w2=0;
            while 1
                varargin{1}.arraydepth_=arraydepth+1;
                if(nargout>3)
                    varargin{1}.jsonpath_=[origpath sprintf('[%d]',length(object))];
                    mmap{end+1}={varargin{1}.jsonpath_, [pos, 0, w2]};
                    [val, pos, index_esc, newmmap] = parse_value(inputstr, pos, esc, index_esc,varargin{:});
                    mmap{end}{2}(2)=pos-mmap{end}{2}(1);
                    mmap=[mmap(:);newmmap(:)];
                else
                    [val, pos,index_esc] = parse_value(inputstr, pos, esc, index_esc,varargin{:});
                end
                object{end+1} = val;
                [cc,pos]=next_char(inputstr,pos);
                if cc == ']'
                    break;
                end
                [pos, w1, w2]=parse_char(inputstr, pos, ',');
            end
        end
    end

    if(varargin{1}.simplifycell)
      if(iscell(object) && ~isempty(object) && (all(cellfun(@isnumeric, object)) || all(cellfun(@isstruct, object))))
          if(all(cellfun(@(e) isequal(size(object{1}), size(e)), object(2:end))))
              try
                  oldobj=object;
                  if(iscell(object) && length(object)>1 && ndims(object{1})>=2)
                      catdim=size(object{1});
                      catdim=ndims(object{1})-(catdim(end)==1)+1;
                      object=cat(catdim,object{:});
                      object=permute(object,ndims(object):-1:1);
                  else
                      object=cell2mat(object.').';
                  end
                  if(iscell(oldobj) && isstruct(object) && numel(object)>1 && varargin{1}.simplifycellarray==0)
                      object=oldobj;
                  end
              catch
              end
          end
      end
      if(~iscell(object) && size(object,1)>1 && ndims(object)==2)
            object=object.';
      end
    end
    pos=parse_char(inputstr, pos, ']');

    if(pbar>0)
        waitbar(pos/length(inputstr),pbar,'loading ...');
    end
end
%%-------------------------------------------------------------------------

function [pos, w1, w2]=parse_char(inputstr, pos, c)
    w1=pos;
    w2=0;
    pos=skip_whitespace(pos, inputstr);
    w1=pos-w1;
    if pos > length(inputstr) || inputstr(pos) ~= c
        pos=error_pos(sprintf('Expected %c at position %%d', c),inputstr,pos);
    else
        pos = pos + 1;
        w2=pos;
        pos=skip_whitespace(pos, inputstr);
        w2=pos-w2;
    end
end
%%-------------------------------------------------------------------------

function [c, pos, w1] = next_char(inputstr, pos)
    w1=pos;
    pos=skip_whitespace(pos, inputstr);
    w1=pos-w1;
    if pos > length(inputstr)
        c = [];
    else
        c = inputstr(pos);
    end
end

%%-------------------------------------------------------------------------
function [str, pos,index_esc, mmap] = parseStr(inputstr, pos, esc, index_esc, varargin)
    if(nargout>3)
        mmap={};
    end
    if inputstr(pos) ~= '"'
        pos=error_pos('String starting with " expected at position %d',inputstr,pos);
    else
        pos = pos + 1;
    end
    str = '';
    while pos <= length(inputstr)
        while index_esc <= length(esc) && esc(index_esc) < pos
            index_esc = index_esc + 1;
        end
        if index_esc > length(esc)
            str = [str inputstr(pos:end)];
            pos = length(inputstr) + 1;
            break;
        else
            str = [str inputstr(pos:esc(index_esc)-1)];
            pos = esc(index_esc);
        end
        nstr = length(str);
        switch inputstr(pos)
            case '"'
                pos = pos + 1;
                if(~isempty(str))
                    if(strcmp(str,'_Inf_'))
                        str=Inf;
                    elseif(strcmp(str,'-_Inf_'))
                        str=-Inf;
                    elseif(strcmp(str,'_NaN_'))
                        str=NaN;
                    end
                end
                return;
            case '\'
                if pos+1 > length(inputstr)
                    pos=error_pos('End of file reached right after escape character',inputstr,pos);
                end
                pos = pos + 1;
                switch inputstr(pos)
                    case {'"' '\' '/'}
                        str(nstr+1) = inputstr(pos);
                        pos = pos + 1;
                    case {'b' 'f' 'n' 'r' 't'}
                        str(nstr+1) = sprintf(['\' inputstr(pos)]);
                        pos = pos + 1;
                    case 'u'
                        if pos+4 > length(inputstr)
                            pos=error_pos('End of file reached in escaped unicode character',inputstr,pos);
                        end
                        str(nstr+(1:6)) = inputstr(pos-1:pos+4);
                        pos = pos + 5;
                end
            otherwise % should never happen
                str(nstr+1) = inputstr(pos);
                keyboard;
                pos = pos + 1;
        end
    end
    str=unescapejsonstring(str);
    pos=error_pos('End of file while expecting end of inputstr',inputstr,pos);
end
%%-------------------------------------------------------------------------

function [num, pos] = parse_number(inputstr, pos, varargin)
    currstr=inputstr(pos:min(pos+30,end));
    [num, one, err, delta] = sscanf(currstr, '%f', 1);
    if ~isempty(err)
        pos=error_pos('Error reading number at position %d',inputstr,pos);
    end
    pos = pos + delta-1;
end
%%-------------------------------------------------------------------------

function varargout = parse_value(inputstr, pos, esc, index_esc, varargin)
    len=length(inputstr);
    if(isfield(varargin{1},'progressbar_'))
        waitbar(pos/len,varargin{1}.progressbar_,'loading ...');
    end
    varargout{3}=index_esc;
    if(nargout>3)
            varargout{4}={};
    end
    switch(inputstr(pos))
        case '"'
            [varargout{1:nargout}] = parseStr(inputstr, pos, esc, index_esc,varargin{:});
            return;
        case '['
            [varargout{1:nargout}] = parse_array(inputstr, pos, esc, index_esc, varargin{:});
            return;
        case '{'
            [varargout{1:nargout}] = parse_object(inputstr, pos, esc, index_esc, varargin{:});
            return;
        case {'-','0','1','2','3','4','5','6','7','8','9'}
            [varargout{1:2}] = parse_number(inputstr, pos, varargin{:});
            return;
        case 't'
            if pos+3 <= len && strcmpi(inputstr(pos:pos+3), 'true')
                varargout{1} = true;
                varargout{2} = pos + 4;
                return;
            end
        case 'f'
            if pos+4 <= len && strcmpi(inputstr(pos:pos+4), 'false')
                varargout{1} = false;
                varargout{2} = pos + 5;
                return;
            end
        case 'n'
            if pos+3 <= len && strcmpi(inputstr(pos:pos+3), 'null')
                varargout{1} = [];
                varargout{2} = pos + 4;
                return;
            end
    end
    varargout{2}=error_pos('Value expected at position %d',inputstr,pos);
end

%%-------------------------------------------------------------------------
function [object, pos, index_esc, mmap] = parse_object(inputstr, pos, esc, index_esc, varargin)
    if(nargout>3)
        mmap={};
        origpath=varargin{1}.jsonpath_;
    end
    pos=parse_char(inputstr, pos, '{');
    usemap=varargin{1}.usemap;
    if(usemap)
	object = containers.Map();
    else
	object = [];
    end
    [cc,pos]=next_char(inputstr,pos);
    if cc ~= '}'
        while 1
            [str, pos, index_esc] = parseStr(inputstr, pos, esc, index_esc, varargin{:});
            if isempty(str) && ~usemap
                str='x0x0_'; % empty name is valid in JSON, decodevarname('x0x0_') restores '\0'
            end
            [pos, w1, w2]=parse_char(inputstr, pos, ':');
            if(nargout>3)
                varargin{1}.jsonpath_=[origpath,'.',str];
                mmap{end+1}={varargin{1}.jsonpath_,[pos,0,w2]};
                [val, pos,index_esc, newmmap] = parse_value(inputstr, pos, esc, index_esc, varargin{:});
                mmap{end}{2}(2)=pos-mmap{end}{2}(1);
                mmap=[mmap(:);newmmap(:)];
            else
                [val, pos,index_esc] = parse_value(inputstr, pos, esc, index_esc, varargin{:});
            end
            if(usemap)
		object(str)=val;
	    else
		object.(encodevarname(str,varargin{:}))=val;
	    end
            [cc,pos]=next_char(inputstr,pos);
            if cc == '}'
                break;
            end
            pos=parse_char(inputstr, pos, ',');
        end
    end
    pos=parse_char(inputstr, pos, '}');
end

%%-------------------------------------------------------------------------

function pos=error_pos(msg, inputstr, pos)
    poShow = max(min([pos-15 pos-1 pos pos+20],length(inputstr)),1);
    if poShow(3) == poShow(2)
        poShow(3:4) = poShow(2)+[0 -1];  % display nothing after
    end
    msg = [sprintf(msg, pos) ': ' ...
    inputstr(poShow(1):poShow(2)) '<error>' inputstr(poShow(3):poShow(4)) ];
    error( ['JSONLAB:JSON:InvalidFormat: ' msg] );
end

%%-------------------------------------------------------------------------

function newpos=skip_whitespace(pos, inputstr)
    newpos=pos;
    while newpos <= length(inputstr) && isspace(inputstr(newpos))
        newpos = newpos + 1;
    end
end

%%-------------------------------------------------------------------------
function newstr=unescapejsonstring(str)
    newstr=str;
    if(iscell(str))
        try
            newstr=cell2mat(cellfun(@(x) cell2mat(x),str(:),'un',0));
        catch
        end
    end
    if(~ischar(str))
        return;
    end
    escapechars={'\\','\"','\/','\a','\b','\f','\n','\r','\t','\v'};
    for i=1:length(escapechars)
        newstr=regexprep(newstr,regexprep(escapechars{i},'\\','\\\\'), escapechars{i});
    end
    newstr=regexprep(newstr,'\\u([0-9A-Fa-f]{4})', '${char(base2dec($1,16))}');
end

%%-------------------------------------------------------------------------
function arraystr=sscanf_prep(str)
    arraystr=str;
    if(regexp(str,'"','once'))
        arraystr=regexprep(arraystr,'"_NaN_"','NaN');
        arraystr=regexprep(arraystr,'"([-+]*)_Inf_"','$1Inf');
    end
    arraystr(arraystr==sprintf('\n'))=' ';
    arraystr(arraystr==sprintf('\r'))=' ';
end

%%-------------------------------------------------------------------------
function [obj, nextidx]=parsendarray(arraystr, dims)
    astr=arraystr;
    astr(astr=='[')=' ';
    astr(astr==']')=' ';
    astr(astr==',')=' ';
    [obj, count, errmsg, nextidx]=sscanf(astr,'%f',inf);
    if(nextidx>=length(astr)-1)
            obj=reshape(obj,dims);
            nextidx=length(arraystr)+1;
    end
end

author: fangq

repo: jsonlab

save_path: github-repos/MATLAB/fangq-jsonlab


text:
% DEMCMU35GPLVMVARGPLVMSIMPLE Run variational GPLVM with dynamics on a single
% sequence of the CMU35 data.

% VARGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

% Define constants
experimentNo = 1000;
latentDim = 9; % this is Q, the number of latent dimensions
indPoints = 100; % number of inducing points
% dynamicKern = {'matern32', 'bias', 'white'};
dynamicKern = {'rbf', 'bias', 'white'}; % kernel k_t for the GP for x(t)
initX ='ppca'; % initialize latent space with ppca

% load data
dataSetName = 'cmu35gplvm';
[Ytmp, lbls, Y, lblstest] = lvmLoadData(dataSetName);

% Set training and test sets
indicesTraining = [1:40 61:size(Y,1)];
indicesTest = setdiff(1:size(Y,1),indicesTraining);
Ytr = Y(indicesTraining,:);
Ytest = Y(indicesTest,:);

% Corresponding timestamps (artificial and equally spaced for this demo)
t = linspace(0, 2*pi, size(Y, 1)+1)';
t = t(1:end-1, 1);
timeStampsTraining = t(indicesTraining);
timeStampsTest = t(indicesTest);

Y = Ytr; clear('Ytr','Ytmp');

% Set up model
options = vargplvmOptions('dtcvar');
options.kern = {'rbfard2', 'bias', 'white'}; % Kernel k_x for the GP for f(x)
options.numActive = indPoints;
options.optimiser = 'scg';

d = size(Y, 2);
fprintf(1,'# Creating the model...\n');
model = vargplvmCreate(latentDim, d, Y, options);
model = vargplvmParamInit(model, model.m, model.X);
model.beta=1/(0.01*var(model.m(:)));
model.vardist.covars = 0.5*ones(size(model.vardist.covars)) + 0.001*randn(size(model.vardist.covars));

%-------- Add dynamics to the model -----
optionsDyn.type = 'vargpTime';
optionsDyn.t=timeStampsTraining;
optionsDyn.inverseWidth=30;
optionsDyn.initX = initX;

% Dynamic kernel:
kern = kernCreate(t, dynamicKern);
% The following is related to the expected number of
% zero-crossings.(larger inv.width numerator, rougher func)
if ~strcmp(kern.comp{1}.type,'ou')
    kern.comp{1}.inverseWidth = optionsDyn.inverseWidth./(((max(t)-min(t))).^2);
    kern.comp{1}.variance = 1;
end
optionsDyn.kern = kern;

% Fill in with default values whatever is not already set
optionsDyn = vargplvmOptionsDyn(optionsDyn);
model = vargplvmAddDynamics(model, 'vargpTime', optionsDyn, optionsDyn.t, 0, 0,optionsDyn.seq);

fprintf(1,'# Further calibration of the initial parameters...\n');
model = vargplvmInitDynamics(model,optionsDyn);
model.vardist.parallel=1;
% do not learn beta for few iterations for intitilization
model.learnBeta = 0;
display = 1;
fprintf(1,'# Intitiliazing the model (fixed beta) %d iterations...\n',100);
model = vargplvmOptimise(model, display, 100);
disp('# Saving model after optimising beta...')
modelWriteResult(model, dataSetName, experimentNo);

% Optimise the model.
model.learnBeta = 1;
iters = 1000; % Default: 1000
fprintf(1,'# Optimising the model for %d iterations...\n',iters);
model = vargplvmOptimise(model, display, iters);
% Save the results.
fprintf(1,'# Saving model after doing %d iterations\n',iters)
modelWriteResult(model, dataSetName, experimentNo);

% See the final lengthscales (see how some dimensions are switched-off).
 bar(model.kern.comp{1}.inputScales)

%% ----------------- PREDICTIONS -----------------------------
fprintf('# Only times prediction...\n');
% Prediction using the only information in the test time points
[Testmeans2 Testcovars2] = vargplvmPredictPoint(model.dynamics, timeStampsTest);
Varmu2 = vargplvmPosteriorMeanVar(model, Testmeans2, Testcovars2);

errorPredictions = mean(abs(Varmu2(:) - Ytest(:)));
fprintf(1,'# Predictions Error:%d\n', errorPredictions);

plot(Varmu2(1,:),'r'),hold on, plot(Ytest(1,:),'g')

skel = acclaimReadSkel('35.asf');
[tmpchan, skel] = acclaimLoadChannels('35_01.amc', skel);
channels{1} = demCmu35VargplvmLoadChannels(Ytest,skel);
channels{2} = demCmu35VargplvmLoadChannels(Varmu2,skel);
skelPlayData2(skel, channels,1/15,{'Ytest','Varmu2'});

author: SheffieldML

repo: vargplvm


text:
function [U,E,J] = slice_triangles(V,F,plane,varargin)
  % SLICE_TRIANGLES
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangles indices into V
  %   plane  4-long plane equation [nx ny nz -d]
  % Outputs:
  %   U  #U by 3 list of vertex positions
  %   E  #E by 2 list of edge indices into U
  %   J  #E by 1 list of birth triangle indices into F
  %

  % HACK!
  [U,E,J] = slice_tets(V,F(:,[1 2 3 3]),plane);
  % matlab fails to reutrn one face as row vecto
  if size(E,2) == 1 && size(E,1) == 3
    E = E';
  end
  if isempty(E)
    U = [];
    J = [];
    return
  end
  [U,~,I] = remove_duplicate_vertices(U,eps);
  size(E)
  E(E(:,2) == E(:,1),2) = E(E(:,2) == E(:,1),3);
  E = E(:,1:2);
  N = normals(V,F(J,:));
  N = N-sum(N.*plane(1:3),2).*plane(1:3);
  M = U(E(:,2),:)-U(E(:,1),:);
  Q = [1 plane(1:3)].*[cos(pi/4) sin(pi/4)*[1 1 1]];
  W = quatmultiply(quatmultiply(Q,[zeros(size(M,1),1) M]),Q.*[1 -1 -1 -1]);
  W = W(:,2:4);
  R = sign(sum(W.*N,2))>0;
  E(R,:) = fliplr(E(R,:));
  E = unique(E,'rows');

end

author: alecjacobson

repo: gptoolbox


text:
function [D] = floydwarshall(A, sym, w)
% % Copyright 2012 Nino Shervashidze
% Input: A - nxn adjacency matrix,
%           sym - boolean, 1 if A and w symmetric
%	    w - nxn weight matrix
% Output: D - nxn distance matrix

n = size(A,1); % number of nodes
D=zeros(n,n);

if nargin<2 % if the graph is not weighted and we have no information about sym, then
  sym=1;
  w=A;
end

if nargin<3 % if the graph is not weighted, then
  w=A;
end

D=w.*A; % if A(i,j)=1,  D(i,j)=w(i,j);
D(A+diag(repmat(Inf,n,1))==0)=Inf; % If A(i,j)~=0 and i~=j D(i,j)=Inf;
D=full(D.*(ones(n)-eye(n))); % set the diagonal to zero

%t=cputime;
if sym % then it is a bit faster
  for k=1:n
    Daux=repmat(full(D(:,k)),1,n);
    Sumdist=Daux+Daux';
    D(Sumdist<D)=Sumdist(Sumdist<D);
  end
else
  for k=1:n
    Daux1=repmat(full(D(:,k)),1,n);
    Daux2=repmat(full(D(k,:)),n,1);
    Sumdist=Daux1+Daux2;
    D(Sumdist<D)=Sumdist(Sumdist<D);
  end
end
%cputime-t
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function getVoxelHull(obj, voxelSize, varargin)
% GETVOXELHULL Compute the voxel hull of a point cloud.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% * The voxel hull is a low resolution representation of the volume occupied by
%   a point cloud. For the computation of the voxel hull the object space is
%   subdivided into a voxel structure. The voxel hull of the point cloud
%   consists of all voxels which contain at least one point of the point cloud.
%
% * Only active points are considered for the computation of the voxel hull.
%
% * The voxel hull can be used to select points, which are overlapping with
%   another point cloud. An example for this task is given in the documentation
%   of the 'select' method. For this call 'help pointCloud.select' and see the
%   examples for the 'InVoxelHull' selection strategy.
% ------------------------------------------------------------------------------
% INPUT
% 1 [voxelSize]
%   Voxel size (equal to edge length) of a single voxel.
% ------------------------------------------------------------------------------
% OUTPUT
% 1 [obj]
%   The voxel hull is attached to the following properties of the ouput object:
%   * obj.voxelHull          = n-by-3 matrix containing the x, y, z coordinates
%                              of the centers of the voxels
%   * obj.voxelHullVoxelSize = edge length of voxels (given by parameter
%                              voxelSize)
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Import a point cloud and compute its voxel hull.
%   pc = pointCloud('Lion.xyz');
%   pc.getVoxelHull(5); % voxel size is 5mm
%   % The voxel hull can now be found in pc.voxelHull and pc.voxelHullVoxelSize
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired('voxelSize', @(x) isscalar(x) && x>0);
% Undocumented
p.addParameter('Centroids', false, @islogical);
p.parse(voxelSize, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear voxelSize

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'GETVOXELHULL'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));

% Compute voxel hull -----------------------------------------------------------

% Lower left point of activated points
lim.min = min(obj.X(obj.act,:), [], 1);
% Round origin (voxel hulls have coincident voxel centers if mod(100, p.voxelSize) == 0)
lim.min = (floor(lim.min/100))*100;

% Indices of voxel cells in x, y and z direction (indices start with 0!)
idxXYZ = [floor( (obj.X(obj.act,1) - lim.min(1)) / p.voxelSize ) ...
          floor( (obj.X(obj.act,2) - lim.min(2)) / p.voxelSize ) ...
          floor( (obj.X(obj.act,3) - lim.min(3)) / p.voxelSize )];

% Remove multiple points to get unique voxels
[voxel, ~, ic] = unique(idxXYZ, 'rows'); % ic for centroids

% Transformation of indices to coordinate system and save to object
obj.voxelHull = [lim.min(1) + p.voxelSize/2 + voxel(:,1) * p.voxelSize ...
                 lim.min(2) + p.voxelSize/2 + voxel(:,2) * p.voxelSize ...
                 lim.min(3) + p.voxelSize/2 + voxel(:,3) * p.voxelSize];

% Round voxel hull (because e.g. -100+0.025+1982*0.05 ~= -0.875)
if p.voxelSize < 1
    noDigits = abs(floor(log10(p.voxelSize)-1)) + 2; % e.g. for 0.05 -> noDigits = 3+2 = 5 (+2 is arbitrary choice to be on the safe side)
    obj.voxelHull = round(obj.voxelHull, noDigits);
end

% Save voxel size to object
obj.voxelHullVoxelSize = p.voxelSize;

% Calculate centroids within each voxel ----------------------------------------

if p.Centroids

    % Add centroid matrix
    obj.voxelHull = [obj.voxelHull zeros(size(voxel,1),3)];

    % Centroid for each voxel
    % Var1
    % tic;
    % for i = 1:size(voxel,1), obj.voxelHull(i,4:6) = mean(obj.X(ic == i,:),1); end
    % toc;
    % c1 = obj.voxelHull(:,4:6);

    % Var2
    % tic;
    % cellX = cell(size(voxel,1),1);
    % for i = 1:size(voxel,1), cellX{i,1} = obj.X(ic == i,:); end
    % cellMean = cellfun(@(x) mean(x,1), cellX, 'UniformOutput', false);
    % obj.voxelHull(:,4:6) = vertcat(cellMean{:});
    % toc;
    % c2 = obj.voxelHull(:,4:6);

    % Var3
    % tic;
    % cellMean = arrayfun(@(x) mean(obj.X(ic == x,:),1), [1:size(voxel,1)]', 'UniformOutput', false);
    % obj.voxelHull(:,4:6) = vertcat(cellMean{:});
    % toc;
    % c3 = obj.voxelHull(:,4:6);

    % Var4
    % tic;
    [~, idx] = sort(ic);
    Xs = obj.X(idx,:);
    rowDist = diff(find([1; diff(ic(idx)); 1]));
    cellX = mat2cell(Xs, rowDist);
    cellMean = cellfun(@(x) mean(x,1), cellX, 'UniformOutput', false);
    obj.voxelHull(:,4:6) = vertcat(cellMean{:});
    % toc;
    % c4 = obj.voxelHull(:,4:6);

end

% End --------------------------------------------------------------------------

msg('E', {'POINTCLOUD' 'GETVOXELHULL'});

end

author: pglira


text:
function hmm = updatePCAparam (hmm,Gammasum,XXGXX,Tfactor,rangeK)

K = hmm.K; ndim = hmm.train.ndim;
if nargin < 5 || isempty(rangeK), rangeK = 1:K; end
if nargin < 4, Tfactor = 1; end
p = hmm.train.lowrank;

for k = rangeK

    % unlike Bishop's mixture of PCA, we don't have a mean vector per state here
    v = hmm.Omega.Gam_rate / hmm.Omega.Gam_shape;
    W = hmm.state(k).W.Mu_W; % posterior dist of the precision matrix
    S = XXGXX{k};% / Gammasum(k);
    %regterm = diag(hmm.state(k).beta.Gam_shape ./ hmm.state(k).beta.Gam_rate);
    SW = S * W;
    M = W'*W + v*eye(p);

    % W
    iS_W = v*eye(p) + M\W'*SW / Gammasum(k);
    S_W = inv(iS_W);
    hmm.state(k).W.Mu_W = SW * S_W / Gammasum(k);
    for n = 1:ndim
        hmm.state(k).W.iS_W(n,:,:) = iS_W;
        hmm.state(k).W.S_W(n,:,:) = S_W;
    end

    % Omega
    Wnew = hmm.state(k).W.Mu_W;
    omega_i = mean(diag(S - SW * (M \ Wnew')));
    % replace
    hmm.Omega.Gam_rate = hmm.Omega.Gam_rate - hmm.Omega.Gam_rate_state(k);
    hmm.Omega.Gam_shape = hmm.Omega.Gam_shape - hmm.Omega.Gam_shape_state(k);
    hmm.Omega.Gam_rate_state(k) = 0.5 * Tfactor * omega_i;
    hmm.Omega.Gam_shape_state(k) = 0.5  * Tfactor * Gammasum(k);
    hmm.Omega.Gam_rate = hmm.Omega.Gam_rate + hmm.Omega.Gam_rate_state(k);
    hmm.Omega.Gam_shape = hmm.Omega.Gam_shape + hmm.Omega.Gam_shape_state(k);

    % prior beta
    %hmm.state(k).beta.Gam_shape = hmm.state(k).prior.beta.Gam_shape + 0.5 * ndim;
    %hmm.state(k).beta.Gam_rate = hmm.state(k).prior.beta.Gam_rate + ...
    %    sum(hmm.state(k).W.Mu_W.^2) + ...
    %    diag(permute(sum(hmm.state(k).W.S_W,1),[2 3 1]))';

end

end

author: OHBA-analysis

repo: HMM-MAR


text:
 function [out weight] = ir_patch_avg(patches, idx, dim, varargin)
%function [out weight] = ir_patch_avg(patches, idx, dim, varargin)
%|
%| average 2D overlapping patches to form image
%| each output pixel value is average of corresponding pixel values
%| from all patches that contain that pixel, i.e, that overlapping it
%|
%| in
%|	patches	[*patch_size npatch]
%|	idx	[npatch]		from ir_im2col
%|
%| option
%|	patch_size	[2]		patch size (default: [8 8])
%|	means	[npatch]		patch means (if need subtracted)
%|						default: empty
%|	nchunk	scalar			# of patches to due in block processing
%|						default: 10000
%|
%| out
%|	out	[dim]	 image formed by averaging overlapping patches
%|
%| based on code from Sai Ravishankar
%| 2016-03-03, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(patches, 'test'), ir_patch_avg_test, return, end

arg.patch_size = [8 8];
arg.means = [];
arg.nchunk = 10000; % how many blocks of patches to do jointly
arg = vararg_pair(arg, varargin);

npatch = size(patches, 2);

if numel(arg.patch_size) ~= 2 || (prod(arg.patch_size) ~= size(patches,1))
	fail 'bad patch_size'
end

b1 = arg.patch_size(1);
b2 = arg.patch_size(2);

[rows, cols] = ind2sub(dim - [b1 b2] + 1, idx);

if ~isempty(arg.means)
	if numel(arg.means) ~= npatch, fail 'bad arg.means', end
	patches = patches + repmat(arg.means, [b1*b2 1]);
end

out = zeros(dim);
weight = zeros(dim);

for jj = 1:arg.nchunk:npatch
	jumpSize = min(jj+arg.nchunk-1, npatch);
	zz = patches(:, jj:jumpSize);
	for ii = jj:jumpSize
		col = cols(ii); row = rows(ii);
		block = reshape(zz(:, ii-jj+1), [b1 b2]);
		i1 = row:(row+b1-1);
		i2 = col:(col+b2-1);
		out(i1, i2) = out(i1, i2) + block; % +=
		weight(i1, i2) = weight(i1, i2) + 1;
	end
end

if any(weight(:) == 0)
	fail 'bug'
end
out = out ./ weight; % average

% ir_patch_avg_test
function ir_patch_avg_test

nx = 2^8 - 8*7*0; ny = 2^8;
xtrue = shepplogan(nx, ny, 1);
%xtrue = rand(nx, ny);
xtrue(end) = 5; % stress ends
xtrue(nx) = 4; % stress ends
bsize = [8 6];
stride = 3;
stride = 1;
[patches idx] = ir_im2col(xtrue, bsize, stride);
im plc 2 2
im(1, xtrue)
[xhat weight] = ir_patch_avg(patches, idx, [nx ny], 'patch_size', bsize);
im(2, xhat)
im(3, weight)
im(4, xhat - xtrue)
equivs(xhat, xtrue)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function h= testmotion
h = [];
h.plotData = @plotData;
h.setRotation = @setRotation;
h.setStepSize = @setStepSize;
h.onerot = @onerot;
h.ccrot = @ccrot;
h.cwrot = @cwrot;
h.readcell = @readcell;
% ---
rotationdirection = 1;
%1 for clockwise
%-1 for counter clockwise
%--
stepsize = 1; % in degrees

hMotor = makeMotor(stepsize,h);

    function onerot()
        hMotor.moveHand(rotationdirection);
    end

    function ccrot()
       rotationdirection = 1;
    end

    function cwrot()
       rotationdirection = -1;
    end

    function result = readcell(lednumber)
        result = hMotor.readcell(lednumber);
    end
end

function hMotor = makeMotor(setstepsize,hmotion)
hMotor = [];

stepsize = setstepsize;
hMotor.moveHand = @moveHand;
hMotor.plotData = @plotData;
hMotor.readcell = @readcell;
currentTheta = 0;
% Initiliaze the clock face
clockface = figure('name','Stepper Motor Simulator');% Figure clockface named "scotts clock"
set(clockface,'NumberTitle','off');                 %
set(clockface,'MenuBar','none');                    % Disable the menu bar
set(clockface,'color','w');                         % Set the backround color to white
set(clockface,'visible','on');                      % Set the clock to visible
set(clockface,'closerequestfcn',@closeRequestFcn);   % Specify closefunc at the bottom of the screen
set(clockface,'Resize','off');
set(clockface,'DoubleBuffer','on');
%set(clockface,'Renderer','OpenGL');
set(clockface,'UserData',hmotion);

%Draw the perimeter of the clock and position the figure
R=linspace(0,2*pi,1000);
x1=9*cos(R);
y1=9*sin(R);
plot(x1,y1,'b','linewidth',8,'color','k') %Draws a thick black circle
hold on
axis off
axis([-10 10 -10 10]) %Draws the figure with +/-10 for [Xmin Xmin Ymin Ymax]
axis equal
%== Tic Marks ==
% Define the Length of the Tic marks
TLenStart = 8.1;    % Start of the Tick mark (distance from origin)
TLenStop = 8.5;     % End of the Tick mark (distance from origin)
[STX,STY,TTX,TTY] = ticMark(TLenStart,TLenStop);
% Plot Skinny and Thick Tick marks on the clock face
plot(STX,STY, 'linewidth',1,'color','k');
%plot(TTX,TTY, 'linewidth',5,'color','k');
[HpDx,HpDy] = GetPolyData(0);

hourhand = fill(HpDx,HpDy,'k');

% Plot the numbers 1-12 on the screen
% --Declare variables to be used in the plotting the clock numbers
Clk_fSize = 12;                                  % controls the font size of the numbers of the clock
Clk_fTheta = (0:pi/4:((2*pi)-(pi/4)))';           % sets the Theta for each number position
Clk_fRad = 7;                                 % sets the Raduis for each number position
Clk_numbas = (0:45:315)';
Clk_nData = [Clk_fRad*cos(Clk_fTheta) Clk_fRad*sin(Clk_fTheta) Clk_numbas];
text(Clk_nData(:,1),Clk_nData(:,2),num2str(Clk_nData(:,3)),...
    'horizontalAlignment','center','verticalAlignment','middle','FontSize',Clk_fSize);

led1 = makeLED(5,0,0);
led2 = makeLED(0,-5,270);
led1.cursorPosition(0);
led2.cursorPosition(0);

hFace = makeFace(stepsize,hMotor,led1,led2);
firststep = hFace.getFirstStep();
currentstep = firststep;

    function result = readcell(cellnumber)
        result = 0;
        if cellnumber ==1
           result =  led1.getState();
        elseif cellnumber ==2
            result = led2.getState();
        end

    end

    function moveHand(direction)
        if direction>0
            currentstep = currentstep.getNext();
            currentstep.landedon();
        else
            currentstep = currentstep.getPrev();
            currentstep.landedon();
        end
    end

    function setRotation(rotdir)
        rotationdirection = rotdir;
    end

    function setStepSize(stepsizedeg)
        stepsize = stepsizedeg;
    end

    function onerot()
        plotData(currentTheta+stepsize*rotationdirection);
        currentTheta = currentTheta+stepsize*rotationdirection;
    end

    function plotData(angleOfArm)
        angleOfArm = angleOfArm*(pi/180);
[HpDx,HpDy] = GetPolyData(angleOfArm);
set(hourhand,'xdata',HpDx,'ydata',HpDy);

%         for idx = 1:.5:abs((angleOfArm*(180/pi))-stateTheta)
%             if currentTheta >=360
%                 currentTheta = 0;
%             elseif currentTheta<=0
%                 currentTheta = 360;
%             end
%             [HpDx,HpDy] = GetPolyData((pi/180)*(currentTheta+(rotationdirection*idx)));
%             set(hourhand,'xdata',HpDx,'ydata',HpDy);
%             led1.cursorPosition((currentTheta+(rotationdirection*idx)));
%             led2.cursorPosition((currentTheta+(rotationdirection*idx)));
%             pause(.01);
%         end
    end

    function [HpDx,HpDy] = GetPolyData(theta)

        HourHandData = [5 2/3 .2 .4];
        [HpDx,HpDy] = PolyEngine(theta,HourHandData);
    end

    function [STX,STY,TTX,TTY] = ticMark(TLenStart,TLenStop)
        %ticMark is given the distance from center to start the tick
        % marks (TLenStart) and the distance from origin to stop the
        % tick marks (TLenStop).
        %STTTheta 60 point array going clockwise skinny ticmarks
        STTheta = pi/2:-2*pi/180:-3*pi/2;
        %Calculates X Y coordinates for all 60 skinny tick marks
        STX = [TLenStart*cos(STTheta') TLenStop*cos(STTheta')]';
        STY = [TLenStart*sin(STTheta') TLenStop*sin(STTheta')]';
        %TTTheta 12 point array going around clockwise thick tic marks
        TTTheta = pi/2:-2*pi/12:-3*pi/2;
        %Calculates X Y coordinates for all 12 thick tic marks
        TTX = [TLenStart*cos(TTTheta') TLenStop*cos(TTTheta')]';
        TTY = [TLenStart*sin(TTTheta') TLenStop*sin(TTTheta')]';
    end %end ticmark function

    function [pDx,pDy] = PolyEngine(Theta,HanData)
        %PolyEngine is given the initial angle and specifications for each
        % polygon it is to generate. It then calculates the data points
        % that will makeup the polygon and passes it back in two variables
        % making up a set of X and Y coordinates that corelate to the
        % current angle. This makes it easy to plot the points as well as
        % easy to change the polygons shape. At this time it is specified
        % with 4 data points but a 5th could be added easily with only
        % changing the data in this PolyEngine function.

        %-Hand Polygon Equations
        %==================================================================
        %-Calculate the length from origin to points A and B.
        oA = sqrt(HanData(1)^2+HanData(3)^2);
        %-Calculate the length from origin to points C and D.
        oB = sqrt(HanData(2)^2+HanData(4)^2);
        %-Calculate X Y coordinates of points A B C and D.
        %-Prepare the X Y data points to be easily passed back and plotted.
        pDx = [oA*cos(Theta+atan(HanData(3)/HanData(1))), ...
            (HanData(1)+HanData(3)*5)*cos(Theta),...
            oA*cos(Theta-atan(HanData(3)/HanData(1))),...
            oB*cos(Theta+atan(HanData(4)/HanData(2))+pi),...
            oB*cos(Theta-atan(HanData(4)/HanData(2))+pi)];
        pDy = [oA*sin(Theta+atan(HanData(3)/HanData(1))) (HanData(1)+HanData(3)*5)*sin(Theta) oA*sin(Theta-atan(HanData(3)/HanData(1))) oB*sin(Theta+atan(HanData(4)/HanData(2))+pi) oB*sin(Theta-atan(HanData(4)/HanData(2))+pi)];
    end%end PolyEngine function

    function closeRequestFcn(varargin)
        closereq
    end%end closerequestfcn

end%end of motor

function hplace = makePlaceHolder(setposition,sethand)
hplace = [];
position = setposition;
hand = sethand;
hplace.setNext = @setNext;
hplace.getNext = @getNext;
hplace.setPrev = @setPrev;
hplace.getPrev = @getPrev;
hplace.addListener = @addListener;
hplace.landedon = @landedon;

theNextPlace = [];
thePrevPlace =[];
listeners = [];

    function landedon()

        hand.plotData(position);
        if ~isempty(listeners)
            listeners.turnOn()
        end
        pause(.1);
    end

    function addListener(listener)
        listeners = listener;
    end

    function setNext(nextplace)
        theNextPlace = nextplace;
    end

    function setPrev(prevplace)
        thePrevPlace = prevplace;
    end

    function result = getNext()

        result = theNextPlace;
        if ~isempty(listeners)
            listeners.turnOff()
        end
    end

    function result = getPrev()

        result = thePrevPlace;
        if ~isempty(listeners)
            listeners.turnOff()
        end
    end

end

function hface = makeFace(stepsize,hand,led1,led2)
hface = [];
hface.getFirstStep = @getFirstStep;

startstep = makePlaceHolder(0,hand);
endstep = makePlaceHolder((360-stepsize),hand);
startstep.setPrev(endstep);
startstep.addListener(led1);
endstep.setNext(startstep);
laststep = [];
beforebeforeendstep= [];
for step = (0+stepsize):stepsize:(360-(2*stepsize))
    newplace = makePlaceHolder(step,hand);
    if step==(0+stepsize)
        startstep.setPrev(endstep);
        startstep.setNext(newplace);
        newplace.setPrev(startstep);
    elseif step==(360-(2*stepsize))
        newplace.setNext(endstep);
        newplace.setPrev(laststep);
        endstep.setPrev(newplace);
        beforebeforeendstep.setNext(newplace);
    else
        laststep.setNext(newplace);
        newplace.setPrev(laststep);
    end

    if step==(360-(3*stepsize))
        beforebeforeendstep = newplace;
    end

    if step==270
       newplace.addListener(led2);
    end
    laststep = newplace;
end

    function result = getFirstStep()
        result = startstep;
    end

end

function hled = makeLED(x,y,angleOfLED)
%sensors
plot(x,y,'*k') %270 degrees
plot(x*2,y*2,'O','color','black') %the circle around the light
led = plot(x*2,y*2,'*k','color',[.75 .75 .75]); %gray it out

hled = [];
hled.cursorPosition = @cursorPosition;
hled.led = led;
hled.turnOn = @turnOn;
hled.turnOff = @turnOff;
hled.getState = @getState;
fov = 4;

if (angleOfLED==0)
    angleofled = [360:-.01:(360-fov),0:.01:(0+fov)];
else
    angleofled = (angleOfLED-fov):(angleOfLED+fov);
end

    function cursorPosition(angleOfArm)
        if any(angleOfArm==angleofled)
            turnOn();
        else
            turnOff();
        end
    end

    function result = getState()
        result = [];
        currentcolor = get(hled.led,'Color');
        if sum(currentcolor)==1
           result = 5;
        else
            result = 0;
        end
    end

    function turnOn
        set(hled.led,'Color','r');
    end

    function turnOff
        set(hled.led,'Color',[.75 .75 .75]);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function stats = ReconstructPosition(positions,spikes,phases,varargin)

%ReconstructPosition - Bayesian reconstruction of positions from spike trains.
%
% Instantaneous positions are reconstructed using a Bayesian algorithm.
% Instantaneous population firing rates can be estimated either over fixed time
% windows, or over fractions of the theta cycle (or of any other brain rhythm).
% Similarly, positions will be reconstructed either over time or phase windows.
% The model is first trained on a subset of the data, then tested on the rest.
%
% USAGE
%
%    stats = ReconstructPosition(positions,spikes,phases,<options>)
%
%    positions      linear or two-dimensional positions, in [0..1]
%    spikes         list of (t,group,cluster) triplets (obtained via e.g.
%                   GetSpikes, using full output)
%    phases         optional unwrapped phase of the LFP (see Phase)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'training'    time interval over which the model should be trained
%                   (default = first half of the position data)
%     'window'      length of the time or phase window (default = 0.020 s for
%                   time, and pi/3 for phases)
%     'type'        two letters (one for X and one for Y) indicating which
%                   coordinates are linear ('l') and which are circular ('c')
%                   - for 1D data, only one letter is used (default 'll')
%     'nBins'       firing curve or map resolution (default = [200 200])
%    =========================================================================
%
%   OUTPUT
%
%     stats.positions     real position across time or phase windows
%     stats.spikes        cell firing vector across time or phase windows
%     stats.estimations   estimated position across time or phase windows
%     stats.errors        estimation error across time or phase windows
%     stats.average       average estimation error in each phase window
%     stats.windows       time windows (possibly computed from phases)
%     stats.phases        phase windows (empty for fixed time windows)
%

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

% Defaults
wt = 0.020; % default time window
wp = pi/3; % default phase window
window = [];
nBins = 200;
training = 0.5;
type = '';
nDimensions = 1;

% Check number of parameters
if nargin < 2 || mod(length(varargin),2) ~= 0,
	builtin('error','Incorrect number of parameters (type ''help ReconstructPosition'' for details).');
end

% Optional parameter 'phases'
if nargin == 2,
	phases = [];
elseif nargin >= 3 && ischar(phases),
	varargin = {phases,varargin{:}};
	phases = [];
end

% Check parameter sizes
if ~isdmatrix(positions),
	builtin('error','Incorrect positions (type ''help ReconstructPosition'' for details).');
end
if ~isdmatrix(spikes),
	builtin('error','Incorrect spikes (type ''help ReconstructPosition'' for details).');
end
if size(positions,2) > 3,
	nDimensions = 2;
end
if ~isempty(phases) && ~isdmatrix(phases),
	builtin('error','Incorrect value for property ''phases'' (type ''help ReconstructPosition'' for details).');
end

% Check number of output parameters
if isempty(phases) && nargout == 3,
	builtin('error','Too many output parameters or missing phases (type ''help ReconstructPosition'' for details).');
end

% Parse parameters
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		builtin('error',['Parameter ' num2str(i+2) ' is not a property (type ''help ReconstructPosition'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'training',
			training = varargin{i+1};
			if ~isdvector(training,'>'),
				builtin('error','Incorrect value for property ''training'' (type ''help ReconstructPosition'' for details).');
			end
		case 'window',
			window = varargin{i+1};
			if ~isdscalar(window,'>0'),
				builtin('error','Incorrect value for property ''window'' (type ''help ReconstructPosition'' for details).');
			end
		case 'show',
			show = varargin{i+1};
			if ~isstring_FMAT(show,'on','off'),
				builtin('error','Incorrect value for property ''show'' (type ''help ReconstructPosition'' for details).');
			end
		case 'nBins',
			nBins = varargin{i+1};
			if ~isinteger(nBins),
				builtin('error','Incorrect value for property ''nBins'' (type ''help ReconstructPosition'' for details).');
			end
		case 'type',
			type = lower(varargin{i+1});
			if (nDimensions == 1 && ~isstring_FMAT(type,'cc','cl','lc','ll')) || (nDimensions == 2 && ~isstring_FMAT(type,'ccl','cll','lcl','lll','ccc','clc','lcc','llc')),
				builtin('error','Incorrect value for property ''type'' (type ''help ReconstructPosition'' for details).');
			end
		otherwise,
			builtin('error',['Unknown property ''' num2str(varargin{i}) ''' (type ''help ReconstructPosition'' for details).']);
	end
end

% Defaults
if isempty(window),
	if isempty(phases),
		window = wt;
	else
		window = wp;
		if ~isiscalar((2*pi)/window),
			builtin('error',['Incorrect phase window: not an integer fraction of 2pi (type ''help ReconstructPosition'' for details).']);
		end
	end
end
if isdscalar(training),
	training = [-Inf positions(1,1)+training*(positions(end,1)-positions(1,1))];
else
	if min([positions(1,1) spikes(1,1)]) < training(1) & min([positions(end,1) spikes(end,1)]) > training(2),
		builtin('error',['Spikes or positions occur both before and after training interval (type ''help ReconstructPosition'' for details).']);
	end
end
if isempty(type),
	if nDimensions == 2,
		type = 'lll';
	else
		type = 'll';
	end
end
nBinsX = nBins(1);
if length(nBins) > 2,
	nBinsY = nBins(2);
else
	if nDimensions == 2,
		nBinsY = nBinsX;
	else
		nBinsY = 1;
	end
end

% List units, assign them an ID (number them from 1 to N), and associate these IDs with each spike
% (IDs will be easier to manipulate than (group,cluster) pairs in subsequent computations)
[units,~,i] = unique(spikes(:,2:end),'rows');
nUnits = length(units);
index = 1:nUnits;
id = index(i)';
spikes = [spikes(:,1) id];

% Split data (training vs test)
trainingPositions = Restrict(positions,training);
trainingSpikes = Restrict(spikes,training);
testPositions = positions(~InIntervals(positions,training),:);
testSpikes = spikes(~InIntervals(spikes,training),:);

% TRAINING

% Compute occupancy probability P(x) (i.e. normalized occupancy map)
firstUnit = trainingSpikes(:,2) == 1;
s = trainingSpikes(firstUnit,1);
map = Map(trainingPositions,s,'nbins',nBins,'smooth',5,'type',type);
Px = map.time;
Px = Px ./ sum(Px(:));

% Compute average firing probability lambda for each unit (i.e. firing maps)
lambda(:,:,1) = map.z;
for i = 2:nUnits,
	nextUnit = trainingSpikes(:,2) == i;
	s = trainingSpikes(nextUnit,1);
	map = Map(trainingPositions,s,'nbins',nBins,'smooth',5,'type',type);
	lambda(:,:,i) = map.z;
end

% TEST

% Determine time windows (using unwrapped phases if necessary)
if ~isempty(phases),
	testPhases = phases(~InIntervals(phases,training),:);
	drop = testPhases(:,1) < testPositions(1,1);
	testPhases(drop,:) = [];
	startPhase = ceil(testPhases(1,2)/(2*pi))*2*pi;
	stopPhase = floor(testPhases(end,2)/(2*pi))*2*pi;
	windows = (startPhase:window:stopPhase)';
	stats.phases = windows;
	windows = Interpolate(testPhases(:,[2 1]),windows);
	windows = [windows(1:end-1,2) windows(2:end,2)];
else
	stats.phases = [];
	windows = (testPositions(1,1):window:testPositions(end,1))';
	windows = [windows(1:end-1) windows(2:end)];
end
nWindows = size(windows,1);

stats.estimations = nan(nBinsY,nBinsX,nWindows);
stats.spikes = zeros(nUnits,nWindows);
% Loop over data windows
for i = 1:nWindows,

	% Get spikes for this window
	s = Restrict(testSpikes,windows(i,:));

	if isempty(s),
		% No spikes: set uniform probability
		stats.estimations(:,:,i) = ones(nBinsY,nBinsX,1)/(nBinsX*nBinsY);
		continue;
	end

	% Population spike count vector
	stats.spikes(:,i) = Accumulate(s(:,2),1,nUnits);
	% To avoid 'for' loops, prepare for vector computation:
	% assign a spike count to each position and unit (3D array)
	n = reshape(repmat(stats.spikes(:,i),1,nBinsX*nBinsY)',nBinsY,nBinsX,nUnits);

	% For each cell i, compute P(ni|x) using a Poisson model
	dt = windows(i,2) - windows(i,1);
	Pnix = (dt*lambda).^n./factorial(n).*exp(-dt*lambda);
	% Compute P(n|x) assuming independent probabilities across units (hmm...)
	% i.e. P(n|x) = product over i of P(ni|x)
	Pnx = prod(Pnix,3);

	% Compute P(n) = sum over x of P(n|x)*P(x)
	Pn = sum(sum(Pnx.*Px));

	% Compute P(x|n) = P(n|x)*P(x)/P(n)
	Pxn = Pnx .* Px / Pn;

	% Store result
	stats.estimations(:,:,i) = Pxn;

end
stats.estimations = squeeze(stats.estimations);
stats.windows = windows;

% Estimation error

stats.errors = [];
stats.average = [];
if nDimensions == 1,
	% Bin test positions and compute distance to center
	stats.positions = Interpolate(testPositions,windows(:,1));
	stats.positions(:,2) = Bin(stats.positions(:,2),[0 1],nBinsX);
	dx = (round(nBinsX/2)-stats.positions(:,2))';
	% Shift estimated position by the real distance to center
	stats.errors = CircularShift(stats.estimations(:,1:length(dx)),dx);
	% Average over one or more cycles
	k = 2*pi/window;
	n = floor(size(stats.errors,2)/k)*k;
	stats.average = reshape(stats.errors(:,1:n),nBins,k,[]);
	stats.average = nanmean(stats.average,3);
else
	warning('Computation of estimation error not yet implemented for 2D environments');
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function open_gripper(coppelia)
    %robot_name = coppelia.robot.name;
    %gripper_name = strcat(robot_name, '_gripper');
    [res retInts retFloats retStrings retBuffer]=coppelia.sim.simxCallScriptFunction(coppelia.clientID,'R1_gripper',coppelia.sim.sim_scripttype_childscript,'open_close_gripper_function',[1],[],[],[],coppelia.sim.simx_opmode_blocking);
    %coppelia_wait(coppelia, 10)
end

repo: ARTE


text:
function triangulation_test16 ( )

%*****************************************************************************80
%
%% TEST16 tests TRIANGULATION_ORDER3_CHECK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  triangle_num = 16;
  node_num = 13;

  triangle_node(1:3,1:triangle_num) = [ ...
     3,   4,   1; ...
     3,   1,   2; ...
     3,   2,   8; ...
     2,   1,   5; ...
     8,   2,  13; ...
     8,  13,   9; ...
     3,   8,   9; ...
    13,   2,   5; ...
     9,  13,   7; ...
     7,  13,   5; ...
     6,   7,   5; ...
     9,   7,   6; ...
    10,   9,   6; ...
     6,   5,  12; ...
    11,   6,  12; ...
    10,   6,  11 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST16\n' );
  fprintf ( 1, '  For a triangulation of a set of nodes,\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_CHECK checks the triangulation.\n' );

  i4mat_transpose_print ( 3, triangle_num, triangle_node, '  Triangles:' );
%
%  Pass all tests.
%
  ierror = triangulation_order3_check ( node_num, triangle_num, triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
%
%  Fail test 1.
%
  node_num2 = 2;

  ierror = triangulation_order3_check ( node_num2, triangle_num, ...
    triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
%
%  Fail test 2.
%
  triangle_num2 = 0;

  ierror = triangulation_order3_check ( node_num, triangle_num2, ...
    triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
%
%  Fail test 3.
%
  isave = triangle_node(2,5);
  triangle_node(2,5) = 0;

  ierror = triangulation_order3_check ( node_num, triangle_num, triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
  triangle_node(2,5) = isave;
%
%  Fail test 4.
%
  isave = triangle_node(3,10);
  triangle_node(3,10) = 2 * node_num + 1;

  ierror = triangulation_order3_check ( node_num, triangle_num, triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
  triangle_node(3,10) = isave;
%
%  Fail test 5.
%
  triangle_node(3,4) = 3;
  triangle_node(3,8) = 3;
  triangle_node(3,10) = 3;
  triangle_node(3,11) = 3;
  triangle_node(2,14) = 3;

  ierror = triangulation_order3_check ( node_num, triangle_num, triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );

  triangle_node(3,4) = 5;
  triangle_node(3,8) = 5;
  triangle_node(3,10) = 5;
  triangle_node(3,11) = 5;
  triangle_node(2,14) = 5;
%
%  Fail test 6.
%
  triangle_node(1,9) = 7;
  ierror = triangulation_order3_check ( node_num, triangle_num, triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
  triangle_node(1,9) = 9;
%
%  Fail test 7.
%
  triangle_node(3,7) = 2;
  ierror = triangulation_order3_check ( node_num, triangle_num, triangle_node );

  fprintf ( 1, '  Error code = %d\n', ierror );
  triangle_node(3,7) = 9;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function val = hrfGet(prfParams,hrfParam)
% hrfGet - get the hrfParameters from each of the scans
%
%   val = hrfGet(params,val)
%
% Brief description:
%   The vistasoft scans have an implicit HRF attached to them.  This
%   function gets the HRF for the parameters attached to a particular
%   scan.
%
% Inputs
%
% Key/value pairs
%   N/A
%
% Return
%
% Example:
%
%
% 2009/03 SOD: modified from readHRFParams
if notDefined('prfParams'), error('prfParams needed'); end
if notDefined('hrfParam'),  error('hrfParam needed');  end

nScans = length(prfParams.stim);

val = [];
switch lower(hrfParam)
    case {'hrfparams'}
        val  = cell(nScans,1);
        for n = 1:nScans
            switch prfParams.stim(n).hrfType
                case 'one gamma (Boynton style)'
                    val{n} = prfParams.stim(n).hrfParams{1};
                case 'two gammas (SPM style)'
                    val{n} = prfParams.stim(n).hrfParams{2};
                case 'impulse'
                    % is this what we want??
                    val{n} = 1;
                otherwise
                    fprintf(1,'[%s]:Unknown hrf type (%s)',mfilename,prfParams.stim(n).hrfType);
            end
        end

    case {'hrftype'}
        val  = cell(nScans,1);
        for n = 1:nScans
            val{n} = prfParams.stim(n).hrfType;
        end

    otherwise
        fprintf(1,'[%s]:Unknown parameter (%s)',mfilename,hrfParam);
end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Main script
clear all, close all;

%% Global options
GO_SHOW_CHANGE = true;	% Show the difference image (DI)
GO_SHOW_MASK = true;	% Show the change map (CM)
GO_SHOW_PRETTIFIED = true;	% Show the prettified detection results
GO_SHOW_ROC_CURVE = true;	% Plot the ROC curve

GO_BAND_PRE_NORM = false;	% Perform a band-wise pre-normalization on the inputs

GO_CONFIG_ROC = {};
GO_VERBOSE = true;
GO_SAVE_RESULTS = false;
GO_OUT_FILE_PATH = './results.xls';		% Path to save the results

% PAUSE MODES:
% -1: resume next iteration when all figures closed;
% -2: press any key to continue
GO_PAUSE_MODE = 1;

PAUSE_EACH_ITER_ = GO_SHOW_CHANGE | GO_SHOW_MASK | GO_SHOW_PRETTIFIED | GO_SHOW_ROC_CURVE;

%% Opt and configure the IMPORTANT ones
%{
	Available algorithms: CVA, DPCA, ImageDiff, ImageRatio, ImageRegr, IRMAD, MAD, PCAkMeans, PCDA

	Available datasets: AirChangeDataset, BernDataset, OSCDDataset, OttawaDataset, TaizhouDataset

	Available binaryzation algorithms: FixedThre, KMeans, OTSU

	Available metrics: AUC, FMeasure, Kappa, OA, Recall, UA
%}
ALG = 'MAD'
DATASET = 'TaizhouDataset'
THRE_ALG = 'KMeans'
METRICS = {'OA', 'UA', 'Recall', 'FMeasure', 'AUC', 'Kappa'}

CONFIG_ALG = {};
CONFIG_DATASET = {'D:\\data\\CD\\Taizhou'};
CONFIG_THRE_ALG = {};
CONFIG_METRICS = {{}, {}, {}, {}, {}, {}};

% Check it
if GO_SHOW_ROC_CURVE
    [~, loc] = ismember('AUC', METRICS);
    if loc == 0
        error('AUC was not included in the desired metrics');
    end
end

%% Construct objects
alg = Algorithms.(ALG)(CONFIG_ALG{:});
dataset = Datasets.(DATASET)(CONFIG_DATASET{:});
iterDS = Datasets.CDDIterator(dataset);
threAlg = ThreAlgs.(THRE_ALG)(CONFIG_THRE_ALG{:});
nMetrics = length(METRICS);
metrics = cell(1, nMetrics);
for ii = 1:nMetrics
    metrics{ii} = Metrics.(METRICS{ii})(CONFIG_METRICS{ii}{:});
end

%% Main loop
while(iterDS.hasNext())
    % Fetch data
    [t1, t2, ref] = iterDS.nextChunk();

    if GO_BAND_PRE_NORM
        % Perform a band-wise z-score normalization before any further
        % algorithm is applied
        fcnNorm = @Utilities.normMeanStd;
        [t1, t2] = deal(fcnNorm(double(t1)), fcnNorm(double(t2)));
    end

    % Make difference image
    DI = alg.detectChange(t1, t2);
    % Segment
    CM = threAlg.segment(DI);
    % Measure
    cellfun(@(obj) obj.update(CM, ref, DI), metrics);

    if GO_VERBOSE
        for ii = 1:nMetrics
            m = metrics{ii};
            fprintf('type: %s\n', METRICS{ii});
            fprintf('\tnewest: %f\n', m.val(end));
            fprintf('\taverage: %f\n', m.avg);
        end
        fprintf('\n')
    end

    if PAUSE_EACH_ITER_
        handles = [];
        if GO_SHOW_CHANGE
            figure('Name', 'Change Map'),
            chns = size(DI, 3);
            if  chns ~= 1 && chns ~=3
                imshow(Utilities.normMinMax(Utilities.mergeAvg(DI)));
            else
                imshow(Utilities.normMinMax(DI));
            end
            handles = [handles, gcf];
        end

        if GO_SHOW_MASK
            figure('Name', 'Change Mask'),
            imshow(CM);
            handles = [handles, gcf];
        end

        if GO_SHOW_PRETTIFIED
            figure('Name', 'Prettified Change Map'),
            imshow(Utilities.pretty(DI, CM, ref));
            handles = [handles, gcf];
        end

        if GO_SHOW_ROC_CURVE
            if ~exist('aucer', 'var')
                aucer = metrics{loc};
            end
            fig = aucer.plotROC(GO_CONFIG_ROC{:});
            handles = [handles, fig];
        end

        if (iterDS.hasNext())
            if GO_PAUSE_MODE == 1
                for h = handles
                    waitfor(h);
                end
            elseif GO_PAUSE_MODE == 2
                pause
            else
                ;
            end
        end
    end
end

%% Collate and save results
results = struct('name', alg.algName, 'threAlg', threAlg.algName, 'dataset', DATASET);
for ii = 1:nMetrics
    results.(METRICS{ii}) = metrics{ii}.avg;
end

if GO_SAVE_RESULTS
    [~, ~, ext] = fileparts(GO_OUT_FILE_PATH);
    switch ext
        case '.mat'
            save(GO_OUT_FILE_PATH, 'results');
        case {'.xls', '.xlsx'}
            writetable(struct2table(results), GO_OUT_FILE_PATH);
        otherwise
            error('Unsupported type of file');
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function obj = EVD_wbmethod(obj)

global H W_mopt Nk Nt Nr  Ns Vn ifVFD n;
t1 = clock;
j = 0;
w = zeros(Nk,1);
v = zeros(Nk,1);
%init
if (ifVFD)
    W_equal = W_mopt;
else
    W_equal = exp( 1i*unifrnd(0,2*pi,Nr,Ns,Nk) );
end
V_equal = zeros(Nt,Ns,Nk);
H_equal = zeros(Ns,Ns,Nk);
m_mse = zeros(Nk,1);
for i = 1:Nk
    w(i) = trace(W_equal(:,:,i)'*W_equal(:,:,i));
end

H1 = zeros(Nt,Ns,Nk);
H2 = zeros(Nr,Ns,Nk);
trigger = 1;
m_MSE_new = 100;

%limit the iterations number by i<10
while ( trigger > 1e-4 && j<10)

    Vn1 = Vn * w;
    for i = 1: Nk
        H1(:,:,i) = H(:,:,i)'*W_equal(:,:,i);
    end
    [V_RF,V_U] = EVD_method(H1,Vn1);

    for i = 1:Nk
        V_equal(:,:,i) = V_RF * V_U(:,:,i);
        v(i) = trace(V_equal(:,:,i)'*V_equal(:,:,i));
        H2(:,:,i) = H(:,:,i)*V_equal(:,:,i);
    end
    Vn2 = Vn * v;
    [W_RF,W_B] = EVD_method(H2,Vn2);

    m_MSE_old = m_MSE_new;

    for k = 1:Nk
        W_equal(:,:,k) = W_RF * W_B(:,:,k);
        w(k) = trace(W_equal(:,:,k)'*W_equal(:,:,k));
        H_equal(:,:,k) = W_equal(:,:,k)'*H2(:,:,k);
        m_mse(k) = trace(H_equal(:,:,k) * H_equal(:,:,k)' - H_equal(:,:,k) - H_equal(:,:,k)')...
            + Vn * v(k) * w(k);
    end
    m_MSE_new = sum(m_mse)/Nk;
    trigger = m_MSE_old - m_MSE_new;
    j = j + 1;
    obj.modmse(j,n) = m_MSE_new + Ns;
end

for i = 1:Nk
    V_B(:,:,i)= V_U(:,:,i) /sqrt(v(i));
end

t2 = clock;
runtime  = etime(t2,t1);
obj.V_B = V_B;
obj.W_B = W_B;
obj.V_RF = V_RF;
obj.W_RF = W_RF;
obj.runtime = obj.runtime + runtime;
obj = get_wbmetric(obj);

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
%% Demo for semi-supervised learning
% Copyright: Andreas Damianou 2015

%% id of differerent random trials
allTrials = 1:2;
saveDir = './results/semiSupervised_master/';

%%
for trial = allTrials
    % For each trial, try different number of observed outputs
    for NtrObs = [10:10:160]
        keep('saveDir','allTrials','trial','NtrObs','res');
        % Different trials should have different seeds
        curSeed = trial + 1000;
        randn('seed', curSeed);
        rand('seed', curSeed);

        dataSetName = 'oil';
        experimentNo = curSeed;
        % This is for when we bootstrap: we try to populate the output
        % labelled space as follows: if Y_i is in the labelled set,
        % (i.e. labels_i exists), then we find this point's embeding, q(x_i) ~
        % N(mu_i, S_i), from which we sample samplesPerObserved points
        % x_samp_i,k, k=1,..., samplesPerObserved and form the pairs
        % {x_samp_i,k, labels_i}, for all k. Here we try two different
        % numbers for samplesPerObserved (a small and a big number).
        samplesPerObserved3 = 6;
        samplesPerObserved4 = 50;

        % Number of iterations for initialization of the model (helps
        % avoiding early local minima)
        initVardistIters = 500;
        % Number of optimization iterations
        itNo = [1200 1000];

        % Load data
        [Y, lbls] = lvmLoadData(dataSetName);
        labels = transformLabels(lbls)';

        % Use Nts test points randomly from the data set. Use the rest
        % NtrObs as training points.
        Nts = 700;
        perm = randperm(size(Y,1));
        ts = perm(1:Nts); t = Nts;
        trObs = perm(t+1:t+NtrObs); t = t+NtrObs;
        trMis = perm(t+1:end);

        %%%%% ATTENTION!!!! For compatibility with the rest of the code,
        %%%%% the notation here differs than from the paper:
        %%%%% Y here is Z in paper
        %%%%% labels here is Y in paper
        % SO:
        % YtrObs: Training points in the labelled set (observed and we also not know their corresponding labelsObs
        % YtrMis: Training points in the unlabelled set (observed, but we do not know their labels labelsMis (we have the below but not show them to the model))
        % Yts:    Test points (we'll infer their corresponding labels labelsTs at test time)
        YtrObs = Y(trObs,:); lblsObs = lbls(trObs, :); labelsObs = labels(trObs,:);
        YtrMis = Y(trMis,:); lblsMis = lbls(trMis, :); labelsMis = labels(trMis,:);
        Yts = Y(ts,:); lblsTs = lbls(ts, :);  labelsTs = labels(ts,:);

        % Options for bgplvm
        options = vargplvmOptions('dtcvar');
        options.kern = 'rbfardjit';
        options.numActive = min(25, size(YtrObs,1)); % number of ind. points
        options.latentDim = 8; % number of latent dimensions
        options.optimiser = 'scg2';
        options.initSNR = 100; % Initial signal to noise ratio

        %----- BGPLVM embed only from the labelled data (no semi-supervised learning here)
        % X is the embedding of only labelled data
        [X, sigma2, W, ~, modelInitVardist] = vargplvmEmbed(YtrObs, options.latentDim, options, initVardistIters, 0 );
        modelInitVardist.globalOpt.experimentNo = experimentNo;
        modelInitVardist.throwSNRError = false;
        model = vargplvmOptimiseModel(modelInitVardist, 1,1, {[], itNo}, 1);
        SNR = vargplvmShowSNR(model,false);
        [x_star_all, varx_star_all, mini] = vargplvmPredictLatent(model, Yts, [], false, 1000, false);

        %----- BGPLVM embed from the labelled and unlabelled data (semi-supervised learning)
        % X2 is a riccher embedding because it is learned from more data
        [X2, sigma22, W2, ~, modelInitVardist2] = vargplvmEmbed([YtrObs; YtrMis], options.latentDim, options, initVardistIters, 0 );
        modelInitVardist2.globalOpt.experimentNo = experimentNo;
        modelInitVardist2.throwSNRError = false;
        model2 = vargplvmOptimiseModel(modelInitVardist2, 1,1, {[], itNo}, 1);
        SNR2 = vargplvmShowSNR(model2,false);
        [x_star_all2, varx_star_all2, mini2] = vargplvmPredictLatent(model2, Yts, [], false, 1000, false);

        %----  X3 is an even richer embedding, because it is a superset of X2:
        % it starts as X2, but for every point in X2 which has a known
        % label we sample extra points and assign that label. So, X3 is a
        % larger set.
        % First copy X3 = X2 (later we'll augment it)
        X3 = model2.X(1:length(trObs),:);
        % We'll need the variances too in order to sample
        varX3 = model2.vardist.covars(1:length(trObs),:);
        % Copy the labels
        labels3 = labelsObs;
        % Xnew, labelsNew is initialized in size
        Xnew = nan(size(X3,1)*samplesPerObserved3, size(X3,2));
        labelsNew = nan(size(Xnew,1), size(labelsObs,2));
        k=1;
        % Now for each labelled point in the embedding we sample more
        % points
        for n=1:size(X3,1) % At this point X3 == X2
            for kk = 1:samplesPerObserved3
                Xnew(k,:) = X3(n,:) + randn(size(X3(n,:))).*sqrt(varX3(n,:));
                labelsNew(k,:) = labels3(n,:);
                k = k + 1;
            end
        end
        X3 = [X3; Xnew];
        labels3 = [labels3; labelsNew];
        clear 'Xnew' 'labelsNew';

        %---------- As above but obtain many more samples ----------------

        %samplesPerObserved4 = min(samplesPerObserved4, floor(10000/NtrObs));
        X4 = model2.X(1:length(trObs),:);
        varX4 = model2.vardist.covars(1:length(trObs),:);
        labels4 = labelsObs;
        Xnew = nan(size(X4,1)*samplesPerObserved4, size(X4,2));
        labelsNew = nan(size(Xnew,1), size(labelsObs,2));
        k=1;
        for n=1:size(X4,1)
            for kk = 1:samplesPerObserved4
                Xnew(k,:) = X4(n,:) + randn(size(X4(n,:))).*sqrt(varX4(n,:));
                labelsNew(k,:) = labels4(n,:);
                k = k + 1;
            end
        end
        X4 = [X4; Xnew];
        labels4 = [labels4; labelsNew];
        clear 'Xnew' 'labelsNew';

        %% %%%%%%%%%%%%%%%%%%% CLASSIFICATION STEP: %%%%%%%%%%%%%%%%
        % Now we have all 4 embeddings / models. Recap:
        % - Model 1: Embedding only from the labelled data
        % - model 2: Embedding from the labelled and unlabelled data
        % - model 3: As above but bootstrap
        % - model 4: As above but bootstrap more points
        %
        % Given these embeddings, we'll train classifiers from the embedded
        % point to the label space.
        %
        %--------- Train classifiers from X to Z
        %
        %---- Model 1
        % Training of logistic regression classifier (replace with any
        % classifier you want). Do this for each of the embeddings.
        LogRegError = util_fitLogReg(model.X, labelsObs, x_star_all, lblsTs);
        LogRegError2 = util_fitLogReg(model2.X(1:length(trObs),:), labelsObs, x_star_all2, lblsTs);
        LogRegError3 = util_fitLogReg(X3, labels3, x_star_all2, lblsTs);
        LogRegError4 = util_fitLogReg(X4, labels4, x_star_all2, lblsTs);
        % We can also classify not from the embeddings to the labels, but
        % directly from the points. However, when the dimensionality of
        % these points is too big, the following method will automatially
        % do a PCA embedding. Notice that the directly classifying from the
        % points might be better than going through an embedding, if tjhe
        % points' dimensionality is small or the data is relatively
        % simple-structured.
        LogRegErrorOutputs = util_fitLogReg(YtrObs, labelsObs, Yts, lblsTs);

        %% -------------- Save results in a structure 'res'
        tmp=['tr' num2str(trial) '_NtrObs' num2str(NtrObs)];
        res.(tmp).seed = curSeed;
        res.(tmp).LogRegError = LogRegError;
        res.(tmp).LogRegError2 = LogRegError2;
        res.(tmp).LogRegError3 = LogRegError3;
        res.(tmp).LogRegError4 = LogRegError4;
        res.(tmp).LogRegErrorOutputs = LogRegErrorOutputs;
        res.(tmp).samplesPerObserved3 = samplesPerObserved3;
        res.(tmp).samplesPerObserved4 = samplesPerObserved4;
        res.(tmp).trObs = trObs;
        res.(tmp).ts = ts;
        res.(tmp).trMis = trMis;
        res.(tmp).W = W;
        res.(tmp).W2 = W2;
        res.(tmp).SNR = SNR;
        res.(tmp).SNR2 = SNR2;
        diary(['LOG_demOilSemiSup4.txt']);
        res.(tmp)
        diary off
        save([saveDir 'demOilSemiSup4.mat'], 'res');
    end
end

%% ------------ PLOTS ------------------
close all
cd(saveDir);
ll = [];
lw = 1.7;
ms = 12;
fs = 40;
errorBars = true;
stdMult = 0.8; % Smaller error bar to be in 0

resAll = [];
doneTrials = [];
try
    load(['demOilSemiSup4.mat']);
    ff = fieldnames(res);
    for jj=1:length(ff)
        if isfield(resAll, ff{jj})
            warning('field already exists')
        end
        resAll.(ff{jj}) = res.(ff{jj});
        try
            doneTrials(end+1) = str2num(ff{jj}(3:find(ff{1}=='_')-1));
        catch ee
            error('ll')
        end
    end
catch e
    warning(['Did not find demOilSemiSup4.mat'])
end
doneTrials = unique(doneTrials);

%%% ATTENTION: Reject trials with low SNR (or repeat them with more
%%% initiVardistIters).
% If embedding results in low SNR (i.e. close to 1), then it has to be rejected. There's no
% automatic check for this now, so you can use the next array manually.
excludeTrials = [];
allTrials = setdiff(doneTrials, excludeTrials);
allNtrObs = [10:10:160];
logRegErrorMatrix = NaN(length(allTrials), length(allNtrObs));
logRegErrorSemiSupMatrix = logRegErrorMatrix;
logRegErrorSemiSupMatrix3 = logRegErrorMatrix;
logRegErrorSemiSupMatrix4 = logRegErrorMatrix;
weightKeptMatrix = logRegErrorMatrix;
weigthKeptSemiSupMatrix = logRegErrorMatrix;
SNRMatrix = logRegErrorMatrix;
SNRSemiSupMatrix = logRegErrorMatrix;

thresh = 0.0002;
%

for i=1:length(allTrials)
    for j=1:length(allNtrObs)
        trial = allTrials(i);
        NtrObs = allNtrObs(j);
        try
            logRegErrorMatrix(i, j) = resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).LogRegError;
            logRegErrorSemiSupMatrix(i, j) = resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).LogRegError2;
            logRegErrorSemiSupMatrix3(i, j) = resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).LogRegError3;
            logRegErrorSemiSupMatrix4(i, j) = resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).LogRegError4;
            weightKeptMatrix(i, j) = sum(resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).W > thresh);
            weigthKeptSemiSupMatrix(i, j) = sum(resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).W2 > thresh);
            SNRMatrix(i, j) = resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).SNR;
            SNRSemiSupMatrix(i, j) = resAll.(['tr' num2str(trial) '_NtrObs' num2str(NtrObs)]).SNR2;
        catch e
            warning(['Could not find tr' num2str(trial) '_NtrObs' num2str(NtrObs)])
        end
    end
end

inds = 1:size(logRegErrorMatrix,2);

meanlogRegErrorMatrix = nanmean(logRegErrorMatrix(:,inds));
meanlogRegErrorSemiSupMatrix = nanmean(logRegErrorSemiSupMatrix(:,inds));
meanlogRegErrorSemiSupMatrix3 = nanmean(logRegErrorSemiSupMatrix3(:,inds));
meanlogRegErrorSemiSupMatrix4 = nanmean(logRegErrorSemiSupMatrix4(:,inds));
if errorBars
    ll(end+1)=errorbar(allNtrObs(inds), meanlogRegErrorMatrix, stdMult*nanstd(logRegErrorMatrix(:,inds)), 'x-', 'LineWidth', lw, 'MarkerSize', ms); hold on;
    ll(end+1)=errorbar(allNtrObs(inds), meanlogRegErrorSemiSupMatrix, stdMult*nanstd(logRegErrorSemiSupMatrix(:,inds)), 'ro-', 'LineWidth', lw, 'MarkerSize', ms);
    ll(end+1)=errorbar(allNtrObs(inds), meanlogRegErrorSemiSupMatrix3, stdMult*nanstd(logRegErrorSemiSupMatrix3(:,inds)), 'g+-', 'LineWidth', lw, 'MarkerSize', ms);
    ll(end+1)=errorbar(allNtrObs(inds), meanlogRegErrorSemiSupMatrix4, stdMult*nanstd(logRegErrorSemiSupMatrix4(:,inds)), 'ks-', 'LineWidth', lw, 'MarkerSize', ms);
    xlabel('# Observed')
    ylabel('# Errors')
    legend('Without Z^u','With Z^u', 'With Samples', 'With More Samples');
    xlim([17,163])
    pp = ylim;
    ylim([-20 pp(2)]);
else
    plot(allNtrObs(inds), meanlogRegErrorMatrix, 'x-'); hold on; title('# Errors')
    plot(allNtrObs(inds), meanlogRegErrorSemiSupMatrix, 'ro-');
    plot(allNtrObs(inds), meanlogRegErrorSemiSupMatrix3, 'g+-');  legend('GP','Semi-supervised GP', 'SS-GP with Samples');
    plot(allNtrObs(inds), meanlogRegErrorSemiSupMatrix4, 'ks-');  legend('GP','Semi-supervised GP', 'SS-GP with More Samples');
end

% figure
% plot(allNtrObs(inds), nanmean(weightKeptMatrix(:,inds)), 'x-'); hold on; title(['Weights kept (threshold=' num2str(thresh) ')'])
% plot(allNtrObs(inds), nanmean(weigthKeptSemiSupMatrix(:,inds)), 'ro-');  legend('GP','Semi-supervised GP');
%
% figure
% plot(allNtrObs(inds), nanmean(SNRMatrix(:,inds)), 'x-'); hold on; title('SNR')
% plot(allNtrObs(inds), nanmean(SNRSemiSupMatrix(:,inds)), 'ro-');  legend('GP','Semi-supervised GP');
%


text:
%% Lucas-Kanade Optical Flow
%
% In this demo, we will:
%
% * understand the concepts of optical flow and its estimation using
%   Lucas-Kanade method.
% * use functions like |cv.calcOpticalFlowPyrLK| to track feature points
%   in a video.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d7/d8b/tutorial_py_lucas_kanade.html>
%

%% Optical Flow
%
% Optical flow is the pattern of apparent motion of image objects between two
% consecutive frames caused by the movemement of object or camera. It is 2D
% vector field where each vector is a displacement vector showing the movement
% of points from first frame to second. Consider the image below:
%
% <<https://docs.opencv.org/3.2.0/optical_flow_basic1.jpg>>
%
% It shows a ball moving in 5 consecutive frames. The arrow shows its
% displacement vector. Optical flow has many applications in areas like:
%
% * Structure from Motion
% * Video Compression
% * Video Stabilization
% * etc.
%
% Optical flow works on several assumptions:
%
% # The pixel intensities of an object do not change between consecutive frames
% # Neighbouring pixels have similar motion
%
% Consider a pixel $I(x,y,t)$ in first frame. It moves by distance $(dx,dy)$
% in next frame taken after $dt$ time. So since those pixels are the same and
% intensity does not change, we can say,
%
% $$I(x,y,t) = I(x+dx, y+dy, t+dt)$$
%
% Then take taylor series approximation of right-hand side, remove common
% terms and divide by $dt$ to get the following equation:
%
% $$f_x u + f_y v + f_t = 0$$
%
% where
%
% $$f_x = \frac{\partial f}{\partial x} \quad ; \quad
%   f_y = \frac{\partial f}{\partial y}$$
%
% $$u = \frac{dx}{dt} \quad ; \quad
%   v = \frac{dy}{dt}$$
%
% Above equation is called Optical Flow equation. In it, we can find $f_x$ and
% $f_y$, they are image gradients. Similarly $f_t$ is the gradient along time.
% But $(u,v)$ is unknown. We cannot solve this one equation with two unknown
% variables. So several methods are provided to solve this problem and one of
% them is Lucas-Kanade.
%

%% Lucas-Kanade method
%
% We have seen an assumption before, that all the neighbouring pixels will
% have similar motion. Lucas-Kanade method takes a 3x3 patch around the point.
% So all the 9 points have the same motion. We can find $(f_x, f_y, f_t)$ for
% these 9 points. So now our problem becomes solving 9 equations with two
% unknown variables which is over-determined. A better solution is obtained
% with least square fit method. Below is the final solution which is two
% equation-two unknown problem and solve to get the solution:
%
% $$
% \left[\matrix{u \cr v}\right] =
% \left[\matrix{
%   \sum_{i}{f_{x_i}}^2       & \sum_{i}{f_{x_i} f_{y_i}} \cr
%   \sum_{i}{f_{x_i} f_{y_i}} & \sum_{i}{f_{y_i}}^2
% }\right]^{-1}
% \left[\matrix{
%   - \sum_{i}{f_{x_i} f_{t_i}} \cr
%   - \sum_{i}{f_{y_i} f_{t_i}}
% }\right]
% $$
%
% (Note similarity of inverse matrix with Harris corner detector. It denotes
% that corners are better points to be tracked.)
%
% So from user point of view, idea is simple, we give some points to track,
% we receive the optical flow vectors of those points. But again there are
% some problems. Until now, we were dealing with small motions. So it fails
% when there is large motion. So again we go for pyramids. When we go up in
% the pyramid, small motions are removed and large motions becomes small
% motions. So applying Lucas-Kanade there, we get optical flow along with the
% scale.
%

%% Lucas-Kanade Optical Flow in OpenCV
%
% OpenCV provides all these in a single function, |cv.calcOpticalFlowPyrLK|.
% Here, we create a simple application which tracks some points in a video.
% To decide the points, we use |cv.goodFeaturesToTrack|. We take the first
% frame, detect some Shi-Tomasi corner points in it, then we iteratively track
% those points using Lucas-Kanade optical flow. For the function
% |cv.calcOpticalFlowPyrLK| we pass the previous frame, previous points and
% next frame. It returns next points along with some status numbers which has
% a value of 1 if next point is found, else zero. We iteratively pass these
% next points as previous points in next step. See the code below.
%
% This code doesn't check how correct are the next keypoints. So even if any
% feature point disappears in image, there is a chance that optical flow finds
% the next point which may look close to it. So actually for a robust tracking,
% corner points should be detected in particular intervals. OpenCV samples
% comes up with such a sample which finds the feature points at every 5 frames.
% It also run a backward-check of the optical flow points got to select only
% good ones.
%

%% Video
% Prepare video source
if mexopencv.require('vision')
    vid = fullfile(toolboxdir('vision'), 'visiondata', 'visiontraffic.avi');
    %cap.PosFrames = 80;  % skip first few seconds with no motion
elseif true
    vid = fullfile(mexopencv.root(), 'test', '768x576.avi');
else
    vid = fullfile(mexopencv.root(), 'test', 'sparse_optical_flow.avi');
    if exist(vid, 'file') ~= 2
        disp('Downloading video...')
        url = 'https://cdn.rawgit.com/opencv/opencv_extra/3.2.0/gpu_demos_pack/demos/sparse_optical_flow/data/sparse_optical_flow.avi';
        urlwrite(url, vid);
    end
end
if exist(vid, 'file') ~= 2, vid = 0; end
cap = cv.VideoCapture(vid);
assert(cap.isOpened(), 'Failed to initialize capturing');

%% First frame
% Grab first frame
frame = cap.read();
assert(~isempty(frame), 'Failed to read frame');
prev = cv.cvtColor(frame, 'RGB2GRAY');

%%
% Detect corners using Shi-Tomasi method (performed only once at the start)
pts0 = cv.goodFeaturesToTrack(prev, ...
    'MaxCorners',100, 'QualityLevel',0.3, 'MinDistance',7, 'BlockSize',7);
pts0 = cat(1, pts0{:});
assert(~isempty(pts0), 'No corners found');
fprintf('%d points\n', size(pts0,1));

%%
% Initialize a mask image for drawing purposes
% (on which point tracks are drawn and remembered)
mask = zeros(size(frame), class(frame));

%%
% Some random colors for plotting
N = 64;
clrs = uint8(hsv(N) * 255);  % randi([0 255], [N 3])
clrs(:,4) = 0;

%%
% Plot
hImg = imshow(frame);
title('PyrLK [Sparse]')

%% Main loop
while ishghandle(hImg)
    % Grab next frame
    frame = cap.read();
    if isempty(frame), break; end
    next = cv.cvtColor(frame, 'RGB2GRAY');

    % Calculate sparse optical flow using Lucas-Kanade method to track points
    [pts1, status] = cv.calcOpticalFlowPyrLK(prev, next, pts0, ...
        'WinSize',[15 15], 'MaxLevel',2, ...
        'Criteria',struct('type','Count+EPS', 'maxCount',10, 'epsilon',0.03));
    pts1 = cat(1, pts1{:});
    status = logical(status);

    % Keep good points (points for which the flow has been found)
    if ~any(status)
        break;
    elseif ~all(status)
        fprintf('%d points\n', nnz(status));
    end
    pts0 = pts0(status,:);
    pts1 = pts1(status,:);

    % Draw latest locations of tracked points
    clr = clrs(rem((1:size(pts1,1))-1,N)+1,:);  % cycle through colors
    frame = cv.circle(frame, pts1, 5, 'Colors',clr, 'Thickness','Filled');

    % Draw point tracks (comet-like plot)
    mask = cv.line(mask, pts1, pts0, 'Colors',clr, 'Thickness',2);
    frame = cv.addWeighted(frame, 0.5, mask, 0.5, 0.0);

    % Display result
    set(hImg, 'CData',frame);
    drawnow;

    % Next iteration: update the previous frame and previous points
    prev = next;
    pts0 = pts1;
end
cap.release();

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% spher() - return the sphering matrix for given input data
%
% Usage:
%
%        >> sphere_matrix = spher(data);
%
% Reference: T. Bell (1996) - - -
%

% S. Makeig CNL / Salk Institute, La Jolla CA 7-17-97

function sphere = spher(data)

if nargin<1 | size(data,1)<1
  help spher
  return
end

sphere = 2.0*inv(sqrtm(cov(data'))); % return the "sphering" matrix

author: PatternRecognition

repo: OpenBMI


text:
function [Fx, Fy] = contactForces(q,dq,ddq,p)
% [Fx, Fy] = contactForces(q,dq,ddq,p)
%
% This function computes the contact forces for the five-link biped
%
% INPUTS:
%   q = [5,n] = configuration
%   dq = [5,n] = rates
%   ddq = [5,n] = accelerations
%   p = parameter struct
%
% OUTPUTS:
%   Fx = [1,n] = horizontal contact force acting on robot
%   Fy = [1,n] = vertical contact force acting on robot
%

[Fx,Fy] = autoGen_contactForce(...
    q(1,:),q(2,:),q(3,:),q(4,:),q(5,:),...
    dq(1,:),dq(2,:),dq(3,:),dq(4,:),dq(5,:),...
    ddq(1,:),ddq(2,:),ddq(3,:),ddq(4,:),ddq(5,:),...
    p.m1, p.m2, p.m3, p.m4, p.m5, p.l1, p.l2, p.l3, p.l4, p.c1, p.c2, p.c3, p.c4, p.c5, p.g);

end

author: MatthewPeterKelly

repo: OptimTraj


text:
% calculate the cross product of vec1 and vec2 (scalar means in z axis)
function outputVector = cross2(inputVector1,inputVector2)
        RotMatCross = [0,-1;1,0];
        switch size(inputVector1,1)+size(inputVector2,1)
            case 2
                outputVector = zeros(1,size(inputVector1,2));
            case 3
                if size(inputVector1,1)==2
                    inputVector1Cross = -RotMatCross*inputVector1;
                    outputVector = inputVector1Cross.*inputVector2;
                else
                    inputVector2Cross = -RotMatCross*inputVector2;
                    outputVector = -inputVector2Cross.*inputVector1;
                end
            case 4
                inputVector1Cross = RotMatCross*inputVector1;
                outputVector = diag(inputVector1Cross'*inputVector2)';
            otherwise
                error('Varargin Error: Input dimension should be either 1xN or 2xN.')
        end
end

repo: paper-simulation


text:
function capproach_prt(jdtdb)

% print closest approach conditions

% Orbital Mechanics with MATLAB

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

global au drsaved iplanet1 iplanet2 pnames

% convert julian day to calendar date and time strings

[cdstr, utstr] = jd2str(jdtdb);

% create text string

textstr = horzcat('time and conditions at closest approach between ', ...
    deblank(pnames(iplanet1, 1:7)), ' and ', deblank(pnames(iplanet2, 1:7)));

fprintf('\n\n');

disp(textstr);
fprintf('===============================================================\n\n');

fprintf('calendar date                ');

disp(cdstr);

fprintf('\nTDB time                     ');

disp(utstr);

fprintf('\nTDB Julian date              %16.8f \n', jdtdb);

jdutc = tdb2utc (jdtdb);

[cdstr, utstr] = jd2str(jdutc);

fprintf('\nUTC time                     ');

disp(utstr);

fprintf('\nUTC Julian date              %16.8f \n\n', jdutc);

% display heliocentric distance (kilometers and AU)

fprintf('close approach distance      %16.6f kilometers\n', drsaved);

fprintf('                             %16.12f AU\n', drsaved / au);


text:
function r8vec_print_part ( n, a, max_print, title )

%*****************************************************************************80
%
%% R8VEC_PRINT_PART prints "part" of an R8VEC.
%
%  Discussion:
%
%    The user specifies MAX_PRINT, the maximum number of lines to print.
%
%    If N, the size of the vector, is no more than MAX_PRINT, then
%    the entire vector is printed, one entry per line.
%
%    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
%    followed by a line of periods suggesting an omission,
%    and the last entry.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines
%    to print.
%
%    Input, string TITLE, a title.
%
  if ( max_print <= 0 )
    return
  end

  if ( n <= 0 )
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  if ( n <= max_print )

    for i = 1 : n
      fprintf ( 1, '  %4d  %14g\n', i, a(i) );
    end

  elseif ( 3 <= max_print )

    for i = 1 : max_print-2
      fprintf ( 1, '  %4d  %14g\n', i, a(i) );
    end
    fprintf ( 1, '  ......  ..............\n' );
    i = n;
    fprintf ( 1, '  %4d  %14g\n', i, a(i) );

  else

    for i = 1 : max_print - 1
      fprintf ( 1, '  %4d  %14g\n', i, a(i) );
    end
    i = max_print;
    fprintf ( 1, '  %4d  %14g  ...more entries...\n', i, a(i) );

  end

  return
end


text:
function pass = test_nonlinSysDamping_C2(pref)
% Test 2x2 system (sin/cos) where damping is required
%
% Asgeir Birkisson, May 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

% Smooth domain:
d = [-pi pi];
A = chebop(@(x,u,v) [u - diff(v,2) + (1-x.^2).*u.^2; diff(u,2) + sin(v)],d);
A.lbc = @(u,v) [u-1; v+1];
A.rbc = @(u,v) [v-1/2; diff(v)];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% CHEBCOLLOC2
pref.discretization = @chebcolloc2;
[u1, u2, info] = solvebvp(A, f, pref);

% Want to check BCs as well.
bcFunLeft = chebfun(A.lbc(u1,u2));
bcFunRight = chebfun(A.rbc(u1,u2));

pass(1) = info.error < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun