text:
function [streams,fileheader] = load_xdf(filename,varargin)
% Import an XDF file.
% [Streams,FileHeader] = load_xdf(Filename, Options...)
%
% This is a MATLAB importer for multi-stream XDF (Extensible Data Format) recordings. All
% information covered by the XDF 1.0 specification is imported, plus any additional meta-data
% associated with streams or with the container file itself.
%
% See http://code.google.com/p/xdf/ for more information on XDF.
%
% The function supports several further features, such as compressed XDF archives, robust
% time synchronization, support for breaks in the data, as well as some other defects.
%
% In:
%   Filename : name of the file to import (*.xdf or *.xdfz)
%
%   Options... : A list of optional name-value arguments for special use cases. The allowed names
%                are listed in the following:
%
%                Parameters that control various processing features:
%
%                'Verbose' : Whether to print verbose diagnostics. (default: false)
%
%                'HandleClockSynchronization' : Whether to enable clock synchronization based on
%                                               ClockOffset chunks. (default: true)
%
%                'HandleJitterRemoval' : Whether to perform jitter removal for regularly sampled
%                                        streams. (default: true)
%
%                'OnChunk' : Function that is called for each chunk of data as it
%                            is being retrieved from the file; the function is allowed to modify the
%                            data (for example, sub-sample it). The four input arguments are 1) the
%                            matrix of [#channels x #samples] values (either numeric or 2d cell
%                            array of strings), 2) the vector of unprocessed local time stamps (one
%                            per sample), 3) the info struct for the stream (same as the .info field
%                            in the final output, buth without the .effective_srate sub-field), and
%                            4) the scalar stream number (1-based integers). The three return values
%                            are 1) the (optionally modified) data, 2) the (optionally modified)
%                            time stamps, and 3) the (optionally modified) header (default: []).
%
%                Parameters for advanced failure recovery in clock synchronization:
%
%                'HandleClockResets' : Whether the importer should check for potential resets of the
%                                      clock of a stream (e.g. computer restart during recording, or
%                                      hot-swap). Only useful if the recording system supports
%                                      recording under such circumstances. (default: true)
%
%                'ClockResetThresholdStds' : A clock reset must be accompanied by a ClockOffset
%                                            chunk being delayed by at least this many standard
%                                            deviations from the distribution. (default: 5)
%
%                'ClockResetThresholdSeconds' : A clock reset must be accompanied by a ClockOffset
%                                               chunk being delayed by at least this many seconds.
%                                               (default: 5)
%
%                'ClockResetThresholdOffsetStds' : A clock reset must be accompanied by a
%                                                  ClockOffset difference that lies at least this many
%                                                  standard deviations from the distribution. (default: 10)
%
%                'ClockResetThresholdOffsetSeconds' : A clock reset must be accompanied by a
%                                                     ClockOffset difference that is at least this
%                                                     many seconds away from the median. (default: 1)
%
%                'ClockResetMaxJitter' : Maximum tolerable jitter (in seconds of error) for clock
%                                        reset handling. (default: 5)
%
%                Parameters for jitter removal in the presence of data breaks:
%
%                'JitterBreakThresholdSeconds' : An interruption in a regularly-sampled stream of at least this
%                                                many seconds will be considered as a potential break (if also
%                                                the BreakThresholdSamples is crossed) and multiple segments
%                                                will be returned. Default: 1
%
%                'JitterBreakThresholdSamples' : An interruption in a regularly-sampled stream of at least this
%                                                many samples will be considered as a potential break (if also
%                                                the BreakThresholdSeconds is crossed) and multiple segments
%                                                will be returned. Default: 500
%
% Out:
%   Streams : cell array of structs, one for each stream; the structs have the following content:
%             .time_series field: contains the stream's time series [#Channels x #Samples]
%                                 this matrix is of the type declared in .info.channel_format
%             .time_stamps field: contains the time stamps for each sample (synced across streams)
%
%             .info field: contains the meta-data of the stream (all values are strings)
%               .name: name of the stream
%               .type: content-type of the stream ('EEG','Events', ...)
%               .channel_format: value format ('int8','int16','int32','int64','float32','double64','string')
%               .nominal_srate: nominal sampling rate of the stream (as declared by the device);
%                               zero for streams with irregular sampling rate
%               .effective_srate: effective (measured) sampling rate of the stream, if regular
%                                 (otherwise omitted)
%               .desc: struct with any domain-specific meta-data declared for the stream; see
%                      www.xdf.org for the declared specifications
%
%             .segments field: struct array containing segment ranges for regularly sampled
%                              time series with breaks (not present if the stream is irregular)
%               .index_range: 1st and last index of the segment within the .time_series/.time_stamps
%                             arrays
%               .t_begin: time of the 1st sample in the segment, in seconds
%               .t_end: time of the last sample in the segment, in seconds
%               .duration: duration of the segment, in seconds
%               .num_samples: number of samples in the segment
%               .effective_srate: effective (i.e. measured) sampling rate within the segment
%
%   FileHeader : struct with file header contents in the .info field
%
% Examples:
%   % load the streams contained in a given XDF file
%   streams = load_xdf('C:\Recordings\myrecording.xdf')
%
% License:
%     This file is covered by the BSD license.
%
%     Copyright (c) 2012, Christian Kothe
%     Portions Copyright (c) 2010, Wouter Falkena
%     All rights reserved.
%
%     Redistribution and use in source and binary forms, with or without
%     modification, are permitted provided that the following conditions are
%     met:
%
%         * Redistributions of source code must retain the above copyright
%           notice, this list of conditions and the following disclaimer.
%         * Redistributions in binary form must reproduce the above copyright
%           notice, this list of conditions and the following disclaimer in
%           the documentation and/or other materials provided with the distribution
%
%     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%     AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%     IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%     ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%     LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%     CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%     SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%     INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%     CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%     ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%     POSSIBILITY OF SUCH DAMAGE.
%
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-04-22
%
%                                Contains portions of xml2struct Copyright (c) 2010, Wouter Falkena,
%                                ASTI, TUDelft, 21-08-2010
%
%                                version 1.14

% check inputs
opts = cell2struct(varargin(2:2:end),varargin(1:2:end),2);
if ~isfield(opts,'OnChunk')
    opts.OnChunk = []; end
if ~isfield(opts,'Verbose')
    opts.Verbose = false; end
if ~isfield(opts,'HandleClockSynchronization')
    opts.HandleClockSynchronization = true; end
if ~isfield(opts,'HandleClockResets')
    opts.HandleClockResets = true; end
if ~isfield(opts,'HandleJitterRemoval')
    opts.HandleJitterRemoval = true; end
if ~isfield(opts,'JitterBreakThresholdSeconds')
    opts.JitterBreakThresholdSeconds = 1; end
if ~isfield(opts,'JitterBreakThresholdSamples')
    opts.JitterBreakThresholdSamples = 500; end
if ~isfield(opts,'ClockResetThresholdSeconds')
    opts.ClockResetThresholdSeconds = 5; end
if ~isfield(opts,'ClockResetThresholdStds')
    opts.ClockResetThresholdStds = 5; end
if ~isfield(opts,'ClockResetThresholdOffsetSeconds')
    opts.ClockResetThresholdOffsetSeconds = 1; end
if ~isfield(opts,'ClockResetThresholdOffsetStds')
    opts.ClockResetThresholdOffsetStds = 10; end
if ~isfield(opts,'WinsorThreshold')
    opts.WinsorThreshold = 0.0001; end
if ~isfield(opts,'ClockResetMaxJitter')
    opts.ClockResetMaxJitter = 5; end
if ~exist(filename,'file')
    error(['The file "' filename '" does not exist.']); end

if opts.Verbose
    disp(['Importing XDF file ' filename '...']); end

% uncompress if necessary (note: "bonus" feature, not part of the XDF 1.0 spec)
[p,n,x] = fileparts(filename);
if strcmp(x,'.xdfz')
    % idea for this type of approach by Michael Kleder
    import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
    src = java.io.FileInputStream(filename);
    flt = java.util.zip.InflaterInputStream(src);
    filename = [p filesep n '_temp_uncompressed' x];
    dst = java.io.FileOutputStream(filename);
    copier = InterruptibleStreamCopier.getInterruptibleStreamCopier;
    copier.copyStream(flt,dst);
    dst.close();
    src.close();
end

streams = {};                                   % cell array of returned streams (in the order of appearance in the file)
idmap = sparse(2^31-1,1);                       % remaps stream id's onto indices in streams
temp = struct();                                % struct array of temporary per-stream information
fileheader = struct();                          % the file header
f = fopen(filename,'r','ieee-le.l64');          % file handle
closer = onCleanup(@()close_file(f,filename));  % object that closes the file when the function exits

filesize = getfield(dir(filename),'bytes');

% there is a fast C mex file for the inner loop, but it's
% not necessarily available for every platform
have_mex = exist('load_xdf_innerloop','file');
if ~have_mex
    disp('NOTE: apparently you are missing a compiled binary version of the inner loop code. Using the slow MATLAB code instead.'); end

% ======================================================
% === parse the file ([SomeText] refers to XDF Spec) ===
% ======================================================

% read [MagicCode]
if ~strcmp(fread(f,4,'*char')','XDF:')
    error(['This is not a valid XDF file (' filename ').']); end

% for each chunk...
while 1
    % read [NumLengthBytes], [Length]
    len = double(read_varlen_int(f));
    if ~len
        if ftell(f) < filesize-1024
            fprintf('  got zero-length chunk, scanning forward to next boundary chunk.\n');
            scan_forward(f);
            continue;
        else
            if opts.Verbose
                fprintf('  reached end of file.\n'); end
            break;
        end
    end
    % read [Tag]
    tag = fread(f,1,'uint16');
    if opts.Verbose
        fprintf('  read tag: %i at %d bytes, length=%d\n',tag,ftell(f),len); end
    switch tag
        case 3 % read [Samples] chunk
            try
                % read [StreamId]
                id = idmap(fread(f,1,'uint32'));
                if have_mex
                    % read the chunk data at once
                    data = fread(f,len-6,'*uint8');
                    % run the mex kernel
                    [values,timestamps] = load_xdf_innerloop(data, temp(id).chns, temp(id).readfmt, temp(id).sampling_interval, temp(id).last_timestamp);
                    temp(id).last_timestamp = timestamps(end);
                else % fallback MATLAB implementation
                    % read [NumSampleBytes], [NumSamples]
                    num = read_varlen_int(f);
                    % allocate space
                    timestamps = zeros(1,num);
                    if strcmp(temp(id).readfmt,'*string')
                        values = cell(temp(id).chns,num);
                    else
                        values = zeros(temp(id).chns,num);
                    end
                    % for each sample...
                    for s=1:num
                        % read or deduce time stamp
                        if fread(f,1,'*uint8')
                            timestamps(s) = fread(f,1,'double');
                        else
                            timestamps(s) = temp(id).last_timestamp + temp(id).sampling_interval;
                        end
                        % read the values
                        if strcmp(temp(id).readfmt,'*string')
                            for v = 1:size(values,1)
                                values{v,s} = fread(f,double(read_varlen_int(f)),'*char')'; end
                        else
                            values(:,s) = fread(f,size(values,1),temp(id).readfmt);
                        end
                        temp(id).last_timestamp = timestamps(s);
                    end
                end
                % optionally send through the OnChunk function
                if ~isempty(opts.OnChunk)
                    [values,timestamps,streams{id}] = opts.OnChunk(values,timestamps,streams{id},id); end %#ok<*AGROW>
                % append to the time series...
                temp(id).time_series{end+1} = values;
                temp(id).time_stamps{end+1} = timestamps;
            catch e
                % an error occurred (perhaps a chopped-off file): emit a warning
                % and scan forward to the next recognized chunk
                fprintf('  got error "%s" (%s), scanning forward to next boundary chunk.\n',e.identifier,e.message);
                scan_forward(f);
            end
        case 2 % read [StreamHeader] chunk
            % read [StreamId]
            streamid = fread(f,1,'uint32');
            id = length(streams)+1;
            idmap(streamid) = id; %#ok<SPRIX>
            % read [Content]
            header = parse_xml_struct(fread(f,len-6,'*char')');
            streams{id} = header;
            if opts.Verbose
                fprintf(['  found stream ' header.info.name '\n']); end
            % generate a few temporary fields
            temp(id).chns = str2num(header.info.channel_count); %#ok<*ST2NM>
            temp(id).srate = str2num(header.info.nominal_srate);
            temp(id).last_timestamp = 0;
            temp(id).time_series = {};
            temp(id).time_stamps = {};
            temp(id).clock_times = [];
            temp(id).clock_values = [];
            if temp(id).srate > 0
                temp(id).sampling_interval = 1/temp(id).srate;
            else
                temp(id).sampling_interval = 0;
            end
            % fread parsing format for data values
            temp(id).readfmt = ['*' header.info.channel_format];
            if strcmp(temp(id).readfmt,'*double64') && ~have_mex
                temp(id).readfmt = '*double'; end % for fread()
        case 6 % read [StreamFooter] chunk
            % read [StreamId]
            id = idmap(fread(f,1,'uint32'));
            % read [Content]
            footer = parse_xml_struct(fread(f,len-6,'*char')');
            streams{id} = hlp_superimposedata(footer,streams{id});
        case 1 % read [FileHeader] chunk
            fileheader = parse_xml_struct(fread(f,len-2,'*char')');
        case 4 % read [ClockOffset] chunk
            % read [StreamId]
            id = idmap(fread(f,1,'uint32'));
            % read [CollectionTime]
            temp(id).clock_times(end+1) = fread(f,1,'double');
            % read [OffsetValue]
            temp(id).clock_values(end+1) = fread(f,1,'double');
        case 5 % read [Boundary] chunk
            fread(f,len-2,'*uint8');
        otherwise
            % skip other chunk types
            fread(f,len-2,'*uint8');
    end
end

% concatenate the signal across chunks
for k=1:length(temp)
    try
        temp(k).time_series = [temp(k).time_series{:}];
        temp(k).time_stamps = [temp(k).time_stamps{:}];
    catch e
        disp(['Could not concatenate time series for stream ' streams{k}.info.name '; skipping.']);
        disp(['Reason: ' e.message]);
        temp(k).time_series = [];
        temp(k).time_stamps = [];
    end
end

% ===================================================================
% === perform (fault-tolerant) clock synchronization if requested ===
% ===================================================================

if opts.HandleClockSynchronization
    if opts.Verbose
        disp('  performing clock synchronization...'); end
    for k=1:length(temp)
        if ~isempty(temp(k).time_stamps)
            try
                clock_times = temp(k).clock_times;
                clock_values = temp(k).clock_values;
                if isempty(clock_times)
                    error('No clock offset values present.'); end
            catch
                disp(['No clock offsets were available for stream "' streams{k}.info.name '"']);
                continue;
            end

            % detect clock resets (e.g., computer restarts during recording) if requested
            % this is only for cases where "everything goes wrong" during recording
            % note that this is a fancy feature that is not needed for normal XDF compliance
            if opts.HandleClockResets
                % first detect potential breaks in the synchronization data; this is only necessary when the
                % importer should be able to deal with recordings where the computer that served a stream
                % was restarted or hot-swapped during an ongoing recording, or the clock was reset otherwise
                time_diff = diff(clock_times);
                value_diff = abs(diff(clock_values));
                % points where a glitch in the timing of successive clock measurements happened
                time_glitch = (time_diff < 0 | (((time_diff - median(time_diff)) ./ mad(time_diff,1)) > opts.ClockResetThresholdStds & ...
                    ((time_diff - median(time_diff)) > opts.ClockResetThresholdSeconds)));
                % points where a glitch in successive clock value estimates happened
                value_glitch = (value_diff - median(value_diff)) ./ mad(value_diff,1) > opts.ClockResetThresholdOffsetStds & ...
                    (value_diff - median(value_diff)) > opts.ClockResetThresholdOffsetSeconds;
                % points where both a time glitch and a value glitch co-occur are treated as resets
                resets_at = time_glitch & value_glitch;
                % determine the [begin,end] index ranges between resets
                if any(resets_at)
                    tmp = find(resets_at)';
                    tmp = [tmp tmp+1]';
                    tmp = [1 tmp(:)' length(resets_at)];
                    ranges = num2cell(reshape(tmp,2,[])',2);
                    if opts.Verbose
                        disp(['  found ' num2str(nnz(resets_at)) ' clock resets in stream ' streams{k}.info.name '.']); end
                else
                    ranges = {[1,length(clock_times)]};
                end
            else
                % otherwise we just assume that there are no clock resets
                ranges = {[1,length(clock_times)]};
            end

            % calculate clock offset mappings for each data range
            mappings = {};
            for r=1:length(ranges)
                idx = ranges{r};
                if idx(1) ~= idx(2)
                    % to accomodate the Winsorizing threshold (in seconds) we rescale the data (robust_fit sets it to 1 unit)
                    mappings{r} = robust_fit([ones(idx(2)-idx(1)+1,1) clock_times(idx(1):idx(2))']/opts.WinsorThreshold, clock_values(idx(1):idx(2))'/opts.WinsorThreshold);
                else
                    mappings{r} = [clock_values(idx(1)) 0]; % just one measurement
                end
            end

            if length(ranges) == 1
                % apply the correction to all time stamps
                temp(k).time_stamps = temp(k).time_stamps + (mappings{1}(1) + mappings{1}(2)*temp(k).time_stamps);
            else
                % if there are data segments measured with different clocks we need to
                % determine, for any time stamp lying between two segments, to which of the segments it belongs
                clock_segments = zeros(size(temp(k).time_stamps));  % the segment index to which each stamp belongs
                begin_of_segment = 1;                               % first index into time stamps that belongs to the current segment
                end_of_segment = NaN; %#ok<NASGU>                   % last index into time stamps that belongs to the current segment
                for r=1:length(ranges)-1
                    cur_end_time = clock_times(ranges{r}(2));       % time at which the current segment ends
                    next_begin_time = clock_times(ranges{r+1}(1));  % time at which the next segment begins
                    % get the data that is not yet processed
                    remaining_indices = begin_of_segment:length(temp(k).time_stamps);
                    if isempty(remaining_indices)
                        break; end
                    remaining_data = temp(k).time_stamps(remaining_indices);
                    if next_begin_time > cur_end_time
                        % clock jumps forward: the end of the segment is where the data time stamps
                        % lie closer to the next segment than the current in time
                        end_of_segment = remaining_indices(min(find([abs(remaining_data-cur_end_time) > abs(remaining_data-next_begin_time),true],1)-1,length(remaining_indices)));
                    else
                        % clock jumps backward: the end of the segment is where the data time stamps
                        % jump back by more than the max conceivable jitter (as any negative delta is jitter)
                        end_of_segment = remaining_indices(min(find([diff(remaining_data) < -opts.ClockResetMaxJitter,true],1),length(remaining_indices)));
                    end
                    % assign the segment of data points to the current range
                    % go to next segment
                    clock_segments(begin_of_segment:end_of_segment) = r;
                    begin_of_segment = end_of_segment+1;
                end
                % assign all remaining time stamps to the last segment
                clock_segments(begin_of_segment:end) = length(ranges);
                % apply corrections on a per-segment basis
                for r=1:length(ranges)
                    temp(k).time_stamps(clock_segments==r) = temp(k).time_stamps(clock_segments==r) + (mappings{r}(1) + mappings{r}(2)*temp(k).time_stamps(clock_segments==r)); end
            end
        end
    end
end

% ===========================================
% === perform jitter removal if requested ===
% ===========================================

if opts.HandleJitterRemoval
    % jitter removal is a bonus feature that yields linearly increasing timestamps from data
    % where samples had been time stamped with some jitter (e.g., due to operating system
    % delays)
    if opts.Verbose
        disp('  performing jitter removal...'); end
    for k=1:length(temp)
        if ~isempty(temp(k).time_stamps) && temp(k).srate
            % identify breaks in the data
            diffs = diff(temp(k).time_stamps);
            breaks_at = diffs > max(opts.JitterBreakThresholdSeconds,opts.JitterBreakThresholdSamples*temp(k).sampling_interval);
            if any(breaks_at)
                % turn the break mask into a cell array of [begin,end] index ranges
                tmp = find(breaks_at)';
                tmp = [tmp tmp+1]';
                tmp = [1 tmp(:)' length(breaks_at)];
                ranges = num2cell(reshape(tmp,2,[])',2);
                if opts.Verbose
                    disp(['  found ' num2str(nnz(breaks_at)) ' data breaks in stream ' streams{k}.info.name '.']); end
            else
                ranges = {[1,length(temp(k).time_stamps)]};
            end

            % process each segment separately
            segments = repmat(struct(),1,length(ranges));
            for r=1:length(ranges)
                range = ranges{r};
                segments(r).num_samples = range(2)-range(1)+1;
                segments(r).index_range = range;
                if segments(r).num_samples > 0
                    indices = segments(r).index_range(1):segments(r).index_range(2);
                    % regress out the jitter
                    mapping = temp(k).time_stamps(indices) / [ones(1,length(indices)); indices];
                    temp(k).time_stamps(indices) = mapping(1) + mapping(2) * indices;
                end
                % calculate some other meta-data about the segments
                segments(r).t_begin = temp(k).time_stamps(range(1));
                segments(r).t_end = temp(k).time_stamps(range(2));
                segments(r).duration = segments(r).t_end - segments(r).t_begin;
                segments(r).effective_srate = segments(r).num_samples / segments(r).duration;
            end

            % calculate the weighted mean sampling rate over all segments
            temp(k).effective_rate = sum(bsxfun(@times,[segments.effective_srate],[segments.num_samples]/sum([segments.num_samples])));

            % transfer the information into the output structs
            streams{k}.info.effective_srate = temp(k).effective_rate;
            streams{k}.segments = segments;
        end
    end
else
    % calculate effective sampling rate
    for k=1:length(temp)
        if ~isempty(temp(k).time_stamps)
            temp(k).effective_srate = length(temp(k).time_stamps) / (temp(k).time_stamps(end) - temp(k).time_stamps(1));
        else
            temp(k).effective_srate = NaN;
        end
    end
end

% copy the information into the output
for k=1:length(temp)
    streams{k}.time_series = temp(k).time_series;
    streams{k}.time_stamps = temp(k).time_stamps;
end
end

% ========================
% === helper functions ===
% ========================

% scan forward through the file until we find the boundary signature
function scan_forward(f)
    block_len = 2^20;
    signature = uint8([67 165 70 220 203 245 65 15 179 14 213 70 115 131 203 228]);
    while 1
        curpos = ftell(f);
        block = fread(f,block_len,'*uint8');
        matchpos = strfind(block',signature);
        if ~isempty(matchpos)
            fseek(f,curpos+matchpos+15,'bof');
            fprintf('  scan forward found a boundary chunk.\n');
            break;
        end
        if length(block) < block_len
            fprintf('  scan forward reached end of file with no match.\n');
            break;
        end
    end
end

% read a variable-length integer
function num = read_varlen_int(f)
try
    switch fread(f,1,'*uint8')
        case 1
            num = fread(f,1,'*uint8');
        case 4
            num = fread(f,1,'*uint32');
        case 8
            num = fread(f,1,'*uint64');
        otherwise
            error('Invalid variable-length integer encountered.');
    end
catch %#ok<*CTCH>
    num = 0;
end
end

% close the file and delete temporary data
function close_file(f,filename)
fclose(f);
if strfind(filename,'_temp_uncompressed.xdf')
    delete(filename); end
end

% parse a simplified (attribute-free) subset of XML into a MATLAB struct
function result = parse_xml_struct(str)
import org.xml.sax.InputSource
import javax.xml.parsers.*
import java.io.*
tmp = InputSource();
tmp.setCharacterStream(StringReader(str));
result = parseChildNodes(xmlread(tmp));

% this is part of xml2struct (slightly simplified)
    function [children,ptext] = parseChildNodes(theNode)
        % Recurse over node children.
        children = struct;
        ptext = [];
        if theNode.hasChildNodes
            childNodes = theNode.getChildNodes;
            numChildNodes = childNodes.getLength;
            for count = 1:numChildNodes
                theChild = childNodes.item(count-1);
                [text,name,childs] = getNodeData(theChild);
                if (~strcmp(name,'#text') && ~strcmp(name,'#comment'))
                    if (isfield(children,name))
                        if (~iscell(children.(name)))
                            children.(name) = {children.(name)}; end
                        index = length(children.(name))+1;
                        children.(name){index} = childs;
                        if(~isempty(text))
                            children.(name){index} = text; end
                    else
                        children.(name) = childs;
                        if(~isempty(text))
                            children.(name) = text; end
                    end
                elseif (strcmp(name,'#text'))
                    if (~isempty(regexprep(text,'[\s]*','')))
                        if (isempty(ptext))
                            ptext = text;
                        else
                            ptext = [ptext text];
                        end
                    end
                end
            end
        end
    end

% this is part of xml2struct (slightly simplified)
    function [text,name,childs] = getNodeData(theNode)
        % Create structure of node info.
        name = char(theNode.getNodeName);
        if ~isvarname(name)
            name = regexprep(name,'[-]','_dash_');
            name = regexprep(name,'[:]','_colon_');
            name = regexprep(name,'[.]','_dot_');
        end
        [childs,text] = parseChildNodes(theNode);
        if (isempty(fieldnames(childs)))
            try
                text = char(theNode.getData);
            catch
            end
        end
    end
end

function x = robust_fit(A,y,rho,iters)
% Perform a robust linear regression using the Huber loss function.
% x = robust_fit(A,y,rho,iters)
%
% Input:
%   A : design matrix
%   y : target variable
%   rho : augmented Lagrangian variable (default: 1)
%   iters : number of iterations to perform (default: 1000)
%
% Output:
%   x : solution for x
%
% Notes:
%   solves the following problem via ADMM for x:
%     minimize 1/2*sum(huber(A*x - y))
%
% Based on the ADMM Matlab codes also found at:
%   http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-03-04

if ~exist('rho','var')
    rho = 1; end
if ~exist('iters','var')
    iters = 1000; end
Aty = A'*y;
L = sparse(chol(A'*A,'lower')); U = L';
z = zeros(size(y)); u = z;
for k = 1:iters
    x = U \ (L \ (Aty + A'*(z - u)));
    d = A*x - y + u;
    z = rho/(1+rho)*d + 1/(1+rho)*max(0,(1-(1+1/rho)./abs(d))).*d;
    u = d - z;
end
end

function res = hlp_superimposedata(varargin)
% Merge multiple partially populated data structures into one fully populated one.
% Result = hlp_superimposedata(Data1, Data2, Data3, ...)
%
% The function is applicable when you have cell arrays or structs/struct arrays with non-overlapping
% patterns of non-empty entries, where all entries should be merged into a single data structure
% which retains their original positions. If entries exist in multiple data structures at the same
% location, entries of later items will be ignored (i.e. earlier data structures take precedence).
%
% In:
%   DataK : a data structure that should be super-imposed with the others to form a single data
%           structure
%
% Out:
%   Result : the resulting data structure
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-08-19

% first, compactify the data by removing the empty items
compact = varargin(~cellfun('isempty',varargin));
% start with the last data structure, then merge the remaining data structures into it (in reverse
% order as this avoids having to grow arrays incrementally in typical cases)
res = compact{end};
for k=length(compact)-1:-1:1
    res = merge(res,compact{k}); end
end

function A = merge(A,B)
% merge data structures A and B
if iscell(A) && iscell(B)
    % make sure that both have the same number of dimensions
    if ndims(A) > ndims(B)
        B = grow_cell(B,size(A));
    elseif ndims(A) < ndims(B)
        A = grow_cell(A,size(B));
    end
    % make sure that both have the same size
    if all(size(B)==size(A))
        % we're fine
    elseif all(size(B)>=size(A))
        % A is a minor of B: grow A
        A = grow_cell(A,size(B));
    elseif all(size(A)>=size(B))
        % B is a minor of A: grow B
        B = grow_cell(B,size(A));
    else
        % A and B have mixed sizes... grow both as necessary
        M = max(size(A),size(B));
        A = grow_cell(A,M);
        B = grow_cell(B,M);
    end
    % find all non-empty elements in B
    idx = find(~cellfun(@(x)isequal(x,[]),B));
    if ~isempty(idx)
        % check if any of these is occupied in A
        clean = cellfun('isempty',A(idx));
        if ~all(clean)
            % merge all conflicting items recursively
            conflicts = idx(~clean);
            for k=conflicts(:)'
                A{k} = merge(A{k},B{k}); end
            % and transfer the rest
            if any(clean)
                A(idx(clean)) = B(idx(clean)); end
        else
            % transfer all to A
            A(idx) = B(idx);
        end
    end
elseif isstruct(A) && isstruct(B)
    % first make sure that both have the same fields
    fnA = fieldnames(A);
    fnB = fieldnames(B);
    if isequal(fnA,fnB)
        % we're fine
    elseif isequal(sort(fnA),sort(fnB))
        % order doesn't match -- impose A's order on B
        B = orderfields(B,fnA);
    elseif isempty(setdiff(fnA,fnB))
        % B has a superset of A's fields: add the remaining fields to A, and order them according to B
        remaining = setdiff(fnB,fnA);
        for fn = remaining'
            A(1).(fn{1}) = []; end
        A = orderfields(A,fnB);
    elseif isempty(setdiff(fnB,fnA))
        % A has a superset of B's fields: add the remaining fields to B, and order them according to A
        remaining = setdiff(fnA,fnB);
        for fn = remaining'
            B(1).(fn{1}) = []; end
        B = orderfields(B,fnA);
    else
        % A and B have incommensurable fields; add B's fields to A's fields, add A's fields to B's
        % and order according to A's fields
        remainingB = setdiff(fnB,fnA);
        for fn = remainingB'
            A(1).(fn{1}) = []; end
        remainingA = setdiff(fnA,fnB);
        for fn = remainingA'
            B(1).(fn{1}) = []; end
        B = orderfields(B,A);
    end
    % that being established, convert them to cell arrays, merge their cell arrays, and convert back to structs
    merged = merge(struct2cell(A),struct2cell(B));
    A = cell2struct(merged,fieldnames(A),1);
elseif isstruct(A) && ~isstruct(B)
    if ~isempty(B)
        error('One of the sub-items is a struct, and the other one is of a non-struct type.');
    else
        % we retain A
    end
elseif isstruct(B) && ~isstruct(A)
    if ~isempty(A)
        error('One of the sub-items is a struct, and the other one is of a non-struct type.');
    else
        % we retain B
        A = B;
    end
elseif iscell(A) && ~iscell(B)
    if ~isempty(B)
        error('One of the sub-items is a cell array, and the other one is of a non-cell type.');
    else
        % we retain A
    end
elseif iscell(B) && ~iscell(A)
    if ~isempty(A)
        error('One of the sub-items is a cell array, and the other one is of a non-cell type.');
    else
        % we retain B
        A = B;
    end
elseif isempty(A) && ~isempty(B)
    % we retain B
    A = B;
elseif isempty(B) && ~isempty(A)
    % we retain A
elseif ~isequal(A,B)
    % we retain A and warn about dropping B
    disp('Two non-empty (and non-identical) sub-elements occupied the same index; one was dropped. This warning will only be displayed once.');
end
end

function C = grow_cell(C,idx)
% grow a cell array to accomodate a particular index
% (assuming that this index is not contained in the cell array yet)
tmp = sprintf('%i,',idx);
eval(['C{' tmp(1:end-1) '} = [];']);
end

repo: REST


text:
%RPY2R Roll-pitch-yaw angles to rotation matrix
%
% R = RPY2R(ROLL, PITCH, YAW, OPTIONS) is an SO(3) orthonornal rotation
% matrix (3x3) equivalent to the specified roll, pitch, yaw angles angles.
% These correspond to rotations about the X, Y, Z axes respectively. If
% ROLL, PITCH, YAW are column vectors (Nx1) then they are assumed to
% represent a trajectory and R is a three-dimensional matrix (3x3xN), where
% the last index corresponds to rows of ROLL, PITCH, YAW.
%
% R = RPY2R(RPY, OPTIONS) as above but the roll, pitch, yaw angles angles
% angles are taken from consecutive columns of the passed matrix RPY =
% [ROLL, PITCH, YAW].  If RPY is a matrix (Nx3) then they are assumed to
% represent a trajectory and R is a three-dimensional matrix (3x3xN), where
% the last index corresponds to rows of RPY which are assumed to be [ROLL,
% PITCH, YAW].
%
% Options::
%  'deg'   Compute angles in degrees (radians default)
%  'zyx'   Return solution for sequential rotations about Z, Y, X axes (Paul book)
%
% Note::
% - In previous releases (<8) the angles corresponded to rotations about ZYX. Many
%   texts (Paul, Spong) use the rotation order ZYX. This old behaviour can be enabled
%   by passing the option 'zyx'
%
% See also TR2RPY, EUL2TR.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function R = rpy2r(roll, varargin)
    opt.zyx = false;
    opt.deg = false;
    [opt,args] = tb_optparse(opt, varargin);

    % unpack the arguments
    if numcols(roll) == 3
		pitch = roll(:,2);
		yaw = roll(:,3);
		roll = roll(:,1);
	elseif nargin >= 3
        pitch = args{1};
        yaw = args{2};
    else
        error('RTB:rpy2r:badarg', 'bad arguments')
    end

    % optionally convert from degrees
    if opt.deg
        d2r = pi/180.0;
        roll = roll * d2r;
        pitch = pitch * d2r;
        yaw = yaw * d2r;
    end

    if ~opt.zyx
        % XYZ order
        if numrows(roll) == 1
            R = rotx(roll) * roty(pitch) * rotz(yaw);
        else
            R = zeros(3,3,numrows(roll));
            for i=1:numrows(roll)
                R(:,:,i) = rotx(roll(i)) * roty(pitch(i)) * rotz(yaw(i));
            end
        end
    else
        % old ZYX order (as per Paul book)
        if numrows(roll) == 1
            R = rotz(roll) * roty(pitch) * rotx(yaw);
        else
            R = zeros(3,3,numrows(roll));
            for i=1:numrows(roll)
                R(:,:,i) = rotz(roll(i)) * roty(pitch(i)) * rotx(yaw(i));
            end
        end
    end

author: Allopart

repo: rbpf-gmapping


text:
function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)
% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
% Pearl p183 eq 4.52

switch msg_type
 case 'd',
  error('gaussian_CPD can''t create discrete msgs')
 case 'g',
  cps = ps(CPD.cps);
  cpsizes = CPD.sizes(CPD.cps);
  self_size = CPD.sizes(end);
  i = find_equiv_posns(p, cps); % p is n's i'th cts parent
  psz = cpsizes(i);
  if all(msg{n}.lambda.precision == 0) % no info to send on
    lam_msg.precision = zeros(psz, psz);
    lam_msg.info_state = zeros(psz, 1);
    return;
  end
  [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence);
  Bmu = m;
  BSigma = Q;
  for k=1:length(cps) % only get pi msgs from cts parents
    pk = cps(k);
    if pk ~= p
      %bk = block(k, cpsizes);
      bk = CPD.cps_block_ndx{k};
      Bk = W(:, bk);
      m = msg{n}.pi_from_parent{k};
      BSigma = BSigma + Bk * m.Sigma * Bk';
      Bmu = Bmu + Bk * m.mu;
    end
  end
  % BSigma = Q + sum_{k \neq i} B_k Sigma_k B_k'
  %bi = block(i, cpsizes);
  bi = CPD.cps_block_ndx{i};
  Bi = W(:,bi);
  P = msg{n}.lambda.precision;
  if (rcond(P) > 1e-3) | isinf(P)
    if isinf(P) % Y is observed
      Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance
      mu_lambda = msg{n}.lambda.mu; % observed_value;
    else
      Sigma_lambda = inv(P);
      mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;
    end
    C = inv(Sigma_lambda + BSigma);
    lam_msg.precision = Bi' * C * Bi;
    lam_msg.info_state = Bi' * C * (mu_lambda - Bmu);
  else
    % method that uses matrix inversion lemma to avoid inverting P
    A = inv(P + inv(BSigma));
    C = P - P*A*P;
    lam_msg.precision = Bi' * C * Bi;
    D = eye(self_size) - P*A;
    z = msg{n}.lambda.info_state;
    lam_msg.info_state = Bi' * (D*z - D*P*Bmu);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function check = triangular_check ( a, b )

%*****************************************************************************80
%
%% TRIANGULAR_CHECK checks the parameters of the Triangular CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    A < B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b <= a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULAR_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= A.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function t = isascii(c)
%ISASCII True for decimal digits.
%
%   For a string C, ISASCII(C) is 1 for any ASCII character code between 0
%   and 127, inclusive, and 0 otherwise.
%
%   See also: ISALNUM, ISALPHA, ISDIGIT, ISLOWER, ISPRTCHR, ISUPPER,
%   ISXDIGIT.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:51:07 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 1, nargin));

   t = ischar(c) & ( 0 <= c ) & ( c <= 127 );

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

%update preview in panels
function MU_update_preview(axes_handle,TMatrix,V)

    axes(axes_handle);
    cla(axes_handle);
    imagesc(TMatrix(:, :, V.Slice),[V.C_lower V.C_upper]);
    colormap(V.Color_map);
    if V.Color_bar==1
        colorbar;
    end
    axis image;
    axis off;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% read_coef.m
% function for reading *.coef file and saving the coefficients to *_des.mat
% file

% INPUT: filename, which is the name of the input coef file.
% If not specifying a filename, a dialogue window will be activated for choosing a
% gipl file.
%
% Wan Jing
% 11/06/2008 - create
%============================================

function fvec = read_coef(filename)

fid = fopen(filename,'r');
if(fid<0)
    fprintf('could not open file %s\n',filename);
    return
end

Count = fscanf(fid, '{ %d,', 1);
if(isempty(Count))
    fprintf('The type of this file is not coef. \n');
    return;
end

% read the coefficients to C
C = fscanf(fid, '{%lf, %lf, %lf},\n', [3 inf]);
C = C';
fprintf('The number of coefficients are: %d .\n ',length(C));

if (Count ~= length(C))
    fprintf('The file has error about the number of coefficients');
end

fvec = [];
i = 1;

md = sqrt(length(C))-1;
for L = 0:md
    for M=-L:L
%         disp(sprintf('i=%d, L=%d, M=%d',i,L,M));
        switch sign(M)
            case -1
                fvec(end+1,:) = [0 0 0];
            case 0
                fvec(end+1,:) = C(i,:); i=i+1;
            case 1
                fvec(end+1,:) = complex(C(i,:),C(i+1,:)); i=i+2;
        end
    end
end

fclose('all');

return

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function pass = test_times( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e2*pref.techPrefs.chebfuneps;

f = ballfun(@(x,y,z)x.*y);
Z = ballfun(@(x,y,z)0);

% Example 1:
F = ballfunv(f,Z,Z);
G = ballfunv(Z,f,Z);
H = times(F,G);
Hexact = ballfunv(Z,Z,Z);
pass(1) = norm(H - Hexact) < tol;

% Example 2:
F = ballfunv(f,Z,Z);
G = ballfunv(-f,Z,Z);
H = times(F,G);
Hexact = ballfunv(-power(f,2),Z,Z);
pass(2) = norm(H - Hexact) < tol;

% Example 3:
F = ballfunv(f,2*f,3*f);
G = ballfunv(-2*f,f,-f);
H = times(F,G);
Hexact = ballfunv(-2*power(f,2),2*power(f,2),-3*power(f,2));
pass(3) = norm(H - Hexact) < tol;

% Example 4:
% Multiply ballfunv by ballfun
V = ballfunv(@(x,y,z)x,@(x,y,z)y,@(x,y,z)z);
f = ballfun(@(x,y,z)cos(y));
exact = ballfunv(@(x,y,z)x.*cos(y),@(x,y,z)y.*cos(y),@(x,y,z)z.*cos(y));
pass(3) = norm(V*f-exact) < tol;
pass(4) = norm(f*V-exact) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef VehicleState < Simulink.IntEnumType
    % MATLAB enumeration class definition generated from template
    %   to track the active leaf state of FMS/FMS State Machine/Vehicle.

    enumeration
        None(0),
		Disarm(1),
		Standby(2),
		Offboard(3),
		Mission(4),
		InvalidAutoMode(5),
		Hold(6),
		Acro(7),
		Stabilize(8),
		Altitude(9),
		Position(10),
		InvalidAssistMode(11),
		Manual(12),
		InValidManualMode(13),
		InvalidArmMode(14),
		Land(15),
		Return(16),
		Takeoff(17)
    end

    methods (Static)

        function defaultValue = getDefaultValue()
            % GETDEFAULTVALUE  Returns the default enumerated value.
            %   If this method is not defined, the first enumeration is used.
            defaultValue = VehicleState.None;
        end

        function dScope = getDataScope()
            % GETDATASCOPE  Specifies whether the data type definition should be imported from,
            %               or exported to, a header file during code generation.
            dScope = 'Auto';
        end

        function desc = getDescription()
            % GETDESCRIPTION  Returns a description of the enumeration.
            desc = 'enumeration to track active leaf state of FMS/FMS State Machine/Vehicle';
        end

        function fileName = getHeaderFile()
            % GETHEADERFILE  Returns path to header file if non-empty.
            fileName = '';
        end

        function flag = addClassNameToEnumNames()
            % ADDCLASSNAMETOENUMNAMES  Indicate whether code generator applies the class name as a prefix
            %                          to the enumeration.
            flag = true;
        end

    end

end

author: Firmament-Autopilot

repo: FMT-Model


text:
function M = mrvMinmax(Y)
%
% M = mrvMinmax(Y);
%
% Return a 2-vector with the minimum and maximum values in a matrix.
%
% Example: Return the vector [1 100]:
%    M = mrvMinmax(magic(10));
%
% ras, 06/16/2008.
% JW: 11/30/2016 Renamed from minmax to mrvMinmax to avoid conflict with
%                the Matlab function minmax in the nnet toolbox
%

if isequal(class(Y), 'single')
	Y = double(Y);
elseif isequal(class(Y), 'cell')
    Y = double( [Y{:}] );
end

Y = Y(~isinf(Y) & ~isnan(Y));

M = [min(Y(:)) max(Y(:))];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function CreateImageSeq(inData, filePath, format)
[height, width, layer] = size(inData{1});
for n = 1:layer
    Img(:,:,1) = inData{1}(:,:,n);
    Img(:,:,2) = inData{2}(:,:,n);
    Img(:,:,3) = inData{3}(:,:,n);
    fileName = [filePath, int2str(n), format];
    imwrite(Img, fileName);
end

author: andrewssobral


text:
function v = norm(F)
%NORM   Frobenius norm of a SPHEREFUNV.
%   V = sqrt(norm(F1).^2 + norm(F2).^2 + norm(F3).^2).

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty(F) )
    v = [];
    return
end

v = 0;
for jj = 1:3
    v = v + sum(norm(F.components{jj}, 2).^2);
end
v = sqrt(v);

end


text:
% Test file for trigtech/min.m

function pass = test_min(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

%%
% Spot-check the extrema for a few functions.

pass(1) = test_spotcheck_min(testclass, @(x) -exp(-cos(2*pi*x)), -exp(1), pref);
pass(2) = test_spotcheck_min(testclass, @(x) -sin(10*pi*x), -1, pref);

pass(3) = test_spotcheck_min(testclass, @(x) -exp(sin(pi*x).^100), -exp(1), pref);
pass(4) = test_spotcheck_min(testclass, @(x) -exp(-sin(pi*x).^100), -1, pref);

% Approx to sign function
pass(5) = test_spotcheck_min(testclass, @(x) -4/pi*(sin(pi*x) + ...
    1/3*sin(3*pi*x) + 1/5*sin(5*pi*x) + 1/7*sin(7*pi*x) + 1/9*sin(9*pi*x)), -1.182328208857607, pref);

%%
% Check operation for array-valued inputs.
fun_op = @(x) -[exp(-cos(2*pi*x)) sin(10*pi*x) exp(-sin(pi*(x-0.32)).^100)];
f = testclass.make(fun_op, [], pref);
[y, x] = min(f);
exact_min = -[exp(1) 1 1];
fx = -[exp(-cos(2*pi*x(1))) sin(10*pi*x(2)) exp(-sin(pi*(x(3)-0.32)).^100)];
pass(6) = (all(abs(y - exact_min) < 100*eps) && ...
           all(abs(fx - exact_min) < 10*eps));

%%
% Test for complex-valued trigtech objects.
f = testclass.make(@(x) cos(pi*x) + exp(1i*pi*x), [], pref);
[y, x] = min(f);
exact_min = 1i; % Could be +/- 1i depending on machine.
pass(7) = ( (abs(y - exact_min) < 1e2*vscale(f)*eps) || ...
            (abs(y + exact_min) < 1e2*vscale(f)*eps) );

end

% Spot-check the results for a given function.
function result = test_spotcheck_min(testclass, fun_op, exact_min, pref)

f = testclass.make(fun_op, [], pref);
[y, x] = min(f);
fx = fun_op(x);
result = ((abs(y - exact_min) < 100*vscale(f)*eps) && ...
          (abs(fx - exact_min) < 10*vscale(f)*eps));

end


text:
function bezier_surface_rectangle_write ( rectangle_file_name, ...
  rectangle_num, rectangle_node )

%*****************************************************************************80
%
%% BEZIER_SURFACE_RECTANGLE_WRITE writes rectangle to a rectangle file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string RECTANGLE_FILE_NAME, the name of the
%    rectangle file.
%
%    Input, integer RECTANGLE_NUM, the number of rectangle.
%
%    Input, integer RECTANGLE_NODE(16,RECTANGLE_NUM), the nodes that
%    make up each rectangle.
%
  dim_num = 16;
  header = 1;

  i4mat_write ( rectangle_file_name, dim_num, rectangle_num, ...
    rectangle_node, header );

  return
end


text:
function [f_output f_votes]= eval_Stochastic_Bosque(Data,Random_Forest,varargin)

%Returns the output of the ensemble (f_output) as well
%as a [num_treesXnum_samples] matrix (f_votes) containing
%the outputs of the individual trees.
%
%The 'oobe' flag allows the out-of-bag error to be used to
%weight the final response (only for classification).

okargs =   {'oobe'};
defaults = {0};
[eid,emsg,oobe_flag] = internal.stats.getargs(okargs,defaults,varargin{:});

f_votes = zeros(numel(Random_Forest),size(Data,1));
oobe = zeros(numel(Random_Forest),1);

for i = 1 : numel(Random_Forest)
    f_votes(i,:) = eval_cartree(Data,Random_Forest(i));
    oobe(i) = Random_Forest(i).oobe;
end

switch lower(Random_Forest(1).method)
    case {'c','g'}
        [unique_labels,~,f_votes]= unique(f_votes);
        f_votes = reshape(f_votes,numel(Random_Forest),size(Data,1));
        [~, f_output] = max(weighted_hist(f_votes,~oobe_flag+oobe_flag*oobe,numel(unique_labels)),[],1); %#ok<ASGLU>
        f_output = unique_labels(f_output);
    case 'r'
        f_output = mean(f_votes,1);
    otherwise
        error('No idea how to evaluate this method');
end

author: Sable

repo: mcbench-benchmarks


text:
function A2 = project_to_simplex(A)
%PROJECT_TO_SIMPLEX Summary of this function goes here
%   Detailed explanation goes here
[N,M] = size(A);
[A1,IX] = sort(A',1,'descend');
J = ones(N,1)*(1:M); J = J';
A1cum = cumsum(A1,1)-1*ones(M,N);
A1minus = (A1 - (A1cum ./ J)) > 0;

rho = max(double(A1minus).*J,[],1);
theta = A1cum(sub2ind([M N],rho,(1:N)))./rho;

A1 = A1 - ones(M,1)*theta;
A1(A1<0) = 0;

% reorder
A2 = zeros(N,M);
for i = 1:M
    A2(:,i) = A1(IX==i);
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
classdef PTKCoordinateSystem
    % PTKCoordinateSystem. Legacy support class for backwards compatibility. Replaced by MimCoordinateSystem
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Constant)
        PTK = MimCoordinateSystem.PTK;                 % Coordinates in mm relatve to the top-left corner of the image volume
        Dicom = MimCoordinateSystem.Dicom              % Coordinates in mm relative to the scanner origin
        DicomUntranslated = MimCoordinateSystem.DicomUntranslated  % Origin in the centre of the first voxel of the image
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef nnsigmoid < nntest
  methods (Test)
   function basic(test)
      x = test.randn(5,5,1,1)/test.range ;
      y = vl_nnsigmoid(x) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nnsigmoid(x,dzdy) ;
      test.der(@(x) vl_nnsigmoid(x), x, dzdy, dzdx, 1e-3) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function plotStream(sx,sy,samplesNum,startFrom)
% PLOTSTREAM  Plots a pair of streams (x,y) for a certain number of samples
if nargin<4 startFrom=1; end
plot(cell2mat(streamRange(sx,startFrom,samplesNum)), ...
     cell2mat(streamRange(sy,startFrom,samplesNum)));


text:
% Copyright 2013 Oliver Johnson, Srikanth Patala
%
% This file is part of MisorientationMaps.
%
%     MisorientationMaps 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.
%
%     MisorientationMaps is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with MisorientationMaps.  If not, see <http://www.gnu.org/licenses/>.

function sphcoord = polarcoord(v)

%%%%% Cartesian to spherical coordinates
%%%%% v = n x 3 array of pts
%%%%% output is n x 3 array of spherical coordinates
%%%%% [r theta phi]

r = sqrt(v(:,1).^2 + v(:,2).^2 + v(:,3).^2);
ind1 = find(r == 0);
ind2 = find(r ~= 0);

if(size(ind1,2) > 0)
theta(ind1)=0; phi(ind1)=0;
end

if(size(ind2,2) > 0)
    theta(ind2) = acos(v(ind2,3)./r(ind2));
    phi(ind2) = atan2(v(ind2,2),v(ind2,1));
end

sphcoord = [r theta' phi'];
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [t,r] = db_index(D, cl, C, p, q)

% DB_INDEX Davies-Bouldin clustering evaluation index.
%
% [t,r] = db_index(D, cl, C, p, q)
%
%  Input and output arguments ([]'s are optional):
%    D     (matrix) data (n x dim)
%          (struct) map or data struct
%    cl    (vector) cluster numbers corresponding to data samples (n x 1)
%    [C]   (matrix) prototype vectors (c x dim) (default = cluster means)
%    [p]   (scalar) norm used in the computation (default == 2)
%    [q]   (scalar) moment used to calculate cluster dispersions (default = 2)
%
%    t     (scalar) Davies-Bouldin index for the clustering (=mean(r))
%    r     (vector) maximum DB index for each cluster (size c x 1)
%
% See also  KMEANS, KMEANS_CLUSTERS, SOM_GAPINDEX.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% input arguments

if isstruct(D),
    switch D.type,
    case 'som_map', D = D.codebook;
    case 'som_data', D = D.data;
    end
end

% cluster centroids
[~, dim] = size(D);
u = unique(cl);
c = length(u);
if nargin <3,
  C = zeros(c,dim);
  for i=1:c,
      me = nanstats(D(cl==u(i),:));
      C(i,:) = me';
  end
end

u2i = zeros(max(u),1); u2i(u) = 1:c;
D = som_fillnans(D,C,u2i(cl)); % replace NaN's with cluster centroid values

if nargin <4, p = 2; end % euclidian distance between cluster centers
if nargin <5, q = 2; end % dispersion = standard deviation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% action

% dispersion in each cluster
S = zeros(1, c);
for i = 1:c
  ind = find(cl==u(i)); % points in this cluster
  l   = length(ind);
  if l > 0
    S(i) = (mean(sqrt(sum((D(ind,:) - ones(l,1) * C(i,:)).^2,2)).^q))^(1/q);
  else
    S(i) = NaN;
  end
end

% distances between clusters
%for i = 1:c
%  for j = i+1:c
%    M(i,j) = sum(abs(C(i,:) - C(j,:)).^p)^(1/p);
%  end
%end
M = som_mdist(C,p);

% Davies-Bouldin index
R = NaN * zeros(c);
r = NaN * zeros(c,1);
for i = 1:c
  for j = i+1:c
    R(i,j) = (S(i) + S(j))/M(i,j);
  end
  r(i) = max(R(i,:));
end

t = mean(r(isfinite(r)));

return;

author: ilarinieminen

repo: SOM-Toolbox


text:
 function mass = xray_read_atomic_mass(mtype, varargin)
%function mass = xray_read_atomic_mass(mtype)
%|
%| Read atomic mass for a given material type.
%|
%| in
%|	mtype			'aluminum', 'copper', 2, '2', '02-helium', ...
%|				See xray_material_file_name.m
%|		(Optionally mtype can be a cell array of several materials.
%|		If so, kev_in is mandatory and the output will be an array.)
%|
%| out
%|	mass	[M 1]	atomic mass
%|
%| Copyright 2008-04-27, Jeff Fessler, The University of Michigan

% default is to show example
if nargin < 1, ir_usage, end
if nargin == 1 && streq(mtype, 'test'), xray_read_atomic_mass_test, return, end

if isnumeric(mtype) && length(mtype) > 1
	mtype = num2cell(mtype);
end

if iscell(mtype)
	mass = zeros(size(mtype)); % array
	for ll=1:length(mtype)
		mass(ll) = xray_read_atomic_mass(mtype{ll}, varargin{:});
	end
return
end

% defaults
arg.units = '';
arg = vararg_pair(arg, varargin);

% trick: mass file near is near material file, and name parsing
if isnumeric(mtype)
	file = xray_material_file_name(1); % trick to allow z > 92
	zz = mtype; % atomic number
else
	file = xray_material_file_name(mtype);

	% atomic number:
	zz = regexprep(file, '/.*\/', '');
	zz = regexprep(zz, '-.*', '');
	zz = str2num(zz);
end

% atomic mass file
file = regexprep(file, 'element.*', 'atomic-mass/atomic-mass-list');

tmp = load_ascii_skip_header(file); % read uncommented lines
if zz < 1 || zz > length(tmp)
	fail('bad element %d', zz)
end
tmp = tmp(zz,:);
if tmp(1) ~= zz, fail('bug %d', zz), end

mass = tmp(2);

% xray_read_atomic_mass_test()
% example usage
function xray_read_atomic_mass_test
mtypes = {'lead', 2, '02-helium'};
mass = xray_read_atomic_mass(mtypes);

zz = 1:118;
mass = xray_read_atomic_mass(zz);
if im
	plot(zz, zz ./ mass, '-o')
	xlabel 'Z'
	ylabel 'Z/A'
%	ir_savefig fig-z-a
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
clear
    % Generarea unei celule goale de dimensiunea 2x2
Q=cell(2,2)
    % Incarcarea celulei
Q={ones(3), [1 3 6 7 6]; ['TN ';'SL ';'DMI';'NDD'], logical(eye(3))};
    % Afisarea informatiilor aferente celulei
whos Q
    % Afisarea continutului inregei celule i a unui element a acesteia
disp('Afisarea intregei celule:')
disp(Q)
disp('Afisarea unui elelement din celula:')
disp(Q{2,1})


text:
clc
clear all
addpath(genpath('./'));

tol = 5e-6; % optimality tolerance for stopping_type 1
DB=20; % SNR of the video
frame_array=[35,100,125]; % frames that will be shown at the end

seed = 602;
randn('state',seed); rand('state',seed);
if exist('Hall_airport_1000_1496_497_144_176_gray.mat')==0
    error('Video file cannot be found. Please download it from: http://www2.ie.psu.edu/aybat/Hall_airport_1000_1496_497_144_176_gray.mat')
else
    load Hall_airport_1000_1496_497_144_176_gray.mat;
end
D = images(:,1:201);
n1=144*176; n2=201;
stdev = norm(D,'fro')/(sqrt(144*176*201)*10^(DB/20));
D = D+stdev*randn(144*176,201);
[X,S]=nsa(D,stdev,tol,1);
figure
plot_data(frame_array,D,X,S,144,176)

repo: lrslibrary


text:
 function y = mtimes_block(ob, x, istart, nblock)
%function y = mtimes_block(ob, x, istart, nblock)
% y = G(i'th block) * x	or y = G'(i'th block) * x
% in either case the projection data will be "small"
% istart is 1,...,nblock

% support 'exists' option for seeing if this routine is available
if nargin == 2 && ischar(x) && streq(x, 'exists')
	y = ~isempty(ob.handle_mtimes_block);
	return
end

if nargin ~= 4, ir_usage, end

if isempty(ob.handle_mtimes_block)
	error 'bug: no mtimes_block() method for this object'
end

if ob.is_transpose
	y = Fatrix_mtimes_block_back(ob, x, istart, nblock);
else
	y = Fatrix_mtimes_block_forw(ob, x, istart, nblock);
end

% caution: cascade_* may not work except for scalars here

function y = Fatrix_mtimes_block_forw(ob, x, istart, nblock);

x = do_cascade(ob.cascade_before, x, false, istart, nblock, true);
y = ob.handle_mtimes_block(ob.arg, ob.is_transpose, x, istart, nblock);
y = do_cascade(ob.cascade_after, y, false, istart, nblock, false);

function x = Fatrix_mtimes_block_back(ob, y, istart, nblock);

y = do_cascade(ob.cascade_after, y, true, istart, nblock, false);
x = ob.handle_mtimes_block(ob.arg, ob.is_transpose, y, istart, nblock);
x = do_cascade(ob.cascade_before, x, true, istart, nblock, true);


text:
classdef CEC2010_F16 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{16};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 10;
            obj.upper    = zeros(1,obj.D) + 10;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum((Z.^2/4000),2) - prod((cos(Z./repmat(sqrt(1:size(Z,2)),size(Z,1),1))),2) + 1;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = sum((Z.^2-100*cos(pi*Z)+10),2);
            PopCon(:,2) = prod(Z,2);
            PopCon(:,3) = abs(sum((Z.*sin(sqrt(abs(Z)))),2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ s, varargout ] = size( x, dim )

%   Disciplined convex/geometric programming information for SIZE:
%       SIZE imposes no convexity restrictions on its first argument.
%       The second argument, if supplied, must be a positive integer.

s = x.size_;
if nargin > 1,
    if nargout > 1,
        cvx_throw( 'Too many output arguments.' );
    elseif ~isnumeric( dim ) || length( dim ) ~= 1 || dim <= 0 || dim ~= floor( dim ),
        cvx_throw( 'Dimension argument must be a positive integer scalar.' );
    elseif dim > length( s ),
        s = 1;
    else
        s = s(dim);
    end
elseif nargout > 1,
    no = nargout;
    s( end + 1 : no ) = 1;
    s( no ) = prod( s( no : end ) );
    for k = 2 : no,
        varargout{k-1} = s(k); %#ok
    end
    s = s(1);
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% pop_runscript() - Run Matlab script
%
% Usage: >> pop_runscript;
%        >> pop_runscript( filename );
%
% Input:
%   filename - [string] name of the file.
%
% Author: Arnaud Delorme, SCCN / INC / UCSD, August 2009

% Copyright (C) Arnaud Delorme, August 2009
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function com = pop_runscript(filename);

com = [];
if nargin <1
    [filename filepath] = uigetfile('*.*', 'Please select input script -- pop_runscript()');

    if filename(1) == 0, return; end;

    filename = fullfile(filepath, filename);
end;

str = readtxtfile(filename);
try
    evalin('base', str);
catch
    lasterr
end;
com = sprintf('pop_runscript(''%s'');', filename);

author: PatternRecognition

repo: OpenBMI


text:
function [columnNames, data] = tableRead(fileName, separator)

% TABLEREAD Read in data which has column titles in the first line and separated values in each other line.
% FORMAT
% DESC reads in data from a file that has column titles in the first line and
% separated values in every other line.
% ARG fileName : file name in which the data is stored.
% ARG separator : separator between the columns (default ',').
% RETURN columnNames : the names of the columns taken from the
% first line.
% RETURN data : the data, taken from the remaining lines.
%
% SEEALSO : fopen, fgetl
%
% COPYRIGHT : Neil D. Lawrence, 2004

% NDLUTIL

if nargin < 2
  separator = ',';
end

fid = fopen(fileName);
i = 0;
while 1
  i = i + 1;
  lin=fgetl(fid);
  if ~ischar(lin), break, end
end
numLines = i;

fid = fopen(fileName);
lin = fgetl(fid);
columnNames = stringSplit(lin, separator);
numCol = length(columnNames);
i = 0;
data = zeros(numLines - 1, numCol);
while 1
  i = i+1;
  lin=fgetl(fid);
  if ~ischar(lin), break, end
  split = stringSplit(lin, separator);
  if length(split) ~= numCol
    error(['Error at line ' num2str(i) ' of file ' fileName ': wrong ' ...
                        'number of columns'])
  end
  for j = 1:length(split);
    data(i, j) = num2str(split{j});
  end
end
fclose(fid);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% GENERATE_RANDOM_OPTIONS Create random options according to some
% delimiters
%
% Usage
%	options = GENERATE_RANDOM_OPTIONS(white_list,type,values)
%
% Input
%    white_list (cell): the white-list corresponding to the name of the
%    field of the options
%    type (cell of string): the type of the variable that one should
%    expect. 'c' continuous, 'i' integer, 'd' discrete
%    values (cell of cell): the corresponding values to a given type.
%
% Output
%    filt_for_disp (numerical): displayed (real) image.
%
% Description
%    For continous and integer, one excepts some interval given by
%    {{s},{e}} where s is the first extremity of the interval and e the
%    end.
%    When the value is discrete, a value of the corresponding cell of
%    values is chosen.
%
% Example: If white_list = {'filter_type', 'precision', 'Q', 'J', 'L',
%   'sigma_phi','sigma_psi','xi_psi', 'slant_psi'};
%   type = {'d', 'i', 'i', 'i', 'i', 'c', 'c', 'c', 'c'}
%   values = { {'morlet','gabor'}, {0,1}, {1,10}, {1,10},
%   {1,10}, {0,50}, {0,30}, {0,30}, {0,10} };
%
%   then GENERATE_RANDOM_OPTIONS(white_list,type,values) is for instance
%
%     filter_type: 'morlet'
%       precision: 0
%               Q: 3
%               J: 8
%               L: 9
%       sigma_phi: 25.7629
%       sigma_psi: 26.1973
%          xi_psi: 1.8308
%       slant_psi: 8.2580
%
% See also
%   CHECK_OPTIONS_WHITE_LIST

function options=generate_random_options(white_list,type,values)

assert(length(white_list)==length(type));
assert(length(white_list)==length(values));

options=struct;

for i=1:length(white_list)
   switch type{i}
       case 'c'  % continuous
           r=values{i}{1}+rand*(values{i}{2}-values{i}{1});
           options=setfield(options,white_list{i},r);
       case 'i'  % integer
           r=floor(values{i}{1}+rand*(1+values{i}{2}-values{i}{1}-1));
           options=setfield(options,white_list{i},r);
       case 'd'  % discrete
           s=ceil(rand*(length(values{i})));
           options=setfield(options,white_list{i},values{i}{s});
   end
end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function out = spm_run_dcm_bms_vis(job)
% Review BMS results
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% CC Chen and Maria Joao Rosa
% $Id: spm_run_dcm_bms_vis.m 6005 2014-05-21 16:46:26Z guillaume $

if ~nargin || isempty(job.bmsmat{1})
    [bmsmat, sts] = spm_select(1,'^BMS\.mat$','Select BMS.mat');
    if ~sts, out = []; return; end
else
    bmsmat = job.bmsmat{1};
end

try
    load(bmsmat);
catch
    error('Cannot load file: %s', bmsmat);
end

lab = {};
val = {};
if isfield(BMS.DCM,'ffx')
    lab = [lab 'Fixed'];
    val = [val 'ffx'];
end
if isfield(BMS.DCM,'rfx')
    lab = [lab 'Random'];
    val = [val 'rfx'];
end
if numel(val) == 2
    lab = [lab 'Bye'];
    val = [val {''}];
end

method = '?';

while ~isempty(method)

    method = spm_input('Inference method','0','b',lab,char(val),numel(val));

    switch method

        case 'ffx'

            N   = size(BMS.DCM.ffx.F,2);
            N   = 1:N;
            out = spm_api_bmc(BMS.DCM.ffx.SF,N);

            if ~isfield(BMS.DCM,'rfx'), method = ''; end

        case 'rfx'

            N   = size(BMS.DCM.rfx.F,2);
            N   = 1:N;
            if isfield(BMS.DCM.rfx,'model')
                out = spm_api_bmc(BMS.DCM.rfx.SF,...
                    N,BMS.DCM.rfx.model.exp_r,BMS.DCM.rfx.model.xp);
            else % Older version (prior to family level)
                out = spm_api_bmc(BMS.DCM.rfx.SF,...
                    N,BMS.DCM.rfx.exp_r,BMS.DCM.rfx.xp);
            end

            if ~isfield(BMS.DCM,'ffx'), method = ''; end

    end

    spm_input('Thank you',1,'d');

end

author: spm


text:
classdef DataGen
  %DATAGEN Various artificial data generation methods.

  methods (Static=true)

    function [input output] = gaussian_single(means, stds, label)
      input = normrnd(means, stds);
      output = label;
    end

    function [inputs outputs] = gaussians(varargin)
      %GAUSSIANS Creates any number of multi-dimensional gaussian clusters,
      %and labels each cluster from 1 to K.
      %
      % GAUSSIANS(examples, mean1, std1, mean2, std2,...)
      %   Examples is the number of examples to generate (the number of
      %   classes must divide this exactly). The number of classes is
      %   implied by the number of following arguments. Each mean argument
      %   is a 1xM array of means, and each std argument is a 1xM array of
      %   standard deviations. The data is sampled using normrnd, and
      %   returned as an NxM matrix of inputs and an Nx1 vector of labels
      %   (labels are processed by DATAGEN/TIDYCLASSES)
      %
      % See also DATAGEN

      examples = varargin{1};
      classes = (nargin - 1) / 2;
      if (mod(examples, classes) ~= 0)
        disp('Error, examples must be a multiple of classes');
      end
      features = numel(varargin{2});
      examples_per_class = floor(examples / classes);

      inputs = reshape(permute(reshape(cell2mat(arrayfun(@(x)normrnd(repmat(varargin{x * 2}, examples_per_class, 1), ...
        repmat(varargin{(x * 2) + 1}, examples_per_class, 1), ...
        examples_per_class, features), 1:classes, 'UniformOutput', false)), ...
        examples_per_class, features, classes), [3 1 2]), examples, features);
      outputs = repmat((1:classes)', examples_per_class, 1);

      outputs = DataGen.tidy_classes(outputs);

    end

    function [inputs outputs] = checkerboard(count, box_width, tile_width, rotation)
      %A rotated checkerboard dataset

      inputs = reshape(arrayfun(@(x)(rand * box_width) - (box_width/2), 1:(count * 2)), count, 2);
      rot_mat = [cos(rotation) -sin(rotation); sin(rotation) cos(rotation)];
      outputs = ssign(double(feval(@(x)x(:, 1) == x(:, 2), mod(floor((rot_mat * inputs')' ./ tile_width), 2))));

    end

    function [inputs outputs] = circle(count, radius, class_prop)
      %A circle dataset

      rots = rand(count, 1) * 2 * pi;
      inputs = repmat(sqrt(rand(count, 1)) .* radius, 1, 2) .* [cos(rots) sin(rots)];
      outputs = ssign((inputs(:, 1).^2 + inputs(:, 2).^2) < class_prop*radius.^2);

    end

    function outputs = tidy_classes(outputs)
      %Class labels should be sequential natural numbers, or (for 2-class)
      %1 and -1.

      if (numel(unique(outputs)) == 2)
        outputs(outputs == min(outputs)) = -1;
        outputs(outputs ~= -1) = 1;
      end

    end

    function inputs = normalize(inputs)
      %Make the input data be 0-mean with a standard deviation of 1.

      inputs = inputs - repmat(mean(inputs), size(inputs, 1), 1);
      inputs = inputs ./ repmat(std(inputs), size(inputs, 1), 1);

    end

  end

end

author: epfl-lasa


text:
classdef MQTTAgent < RealTimeAgents.RealTimeAgentBase

end

author: CelsoReyes


text:
function x = clearconicinfo(x)
x.conicinfo = [0 0];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% The COBRAToolbox: testdynamicRFBA.m
%
% Purpose:
%     - testrFBA tests the dynamicRFBA function and its different outputs
%
% Author:
%     - Marouen BEN GUEBILA - 31/01/2017

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testdynamicRFBA'));
cd(fileDir);

%load model and test data
modelReg = getDistributedModel('modelReg.mat');
load('refData_dynamicRFBA.mat');

%Solver packages
solverPkgs = {'tomlab_cplex','ibm_cplex'};
%QP solvers
QPsolverPkgs = {'tomlab_cplex','ibm_cplex'};

for k =1:length(solverPkgs)

    solverLPOK = changeCobraSolver(solverPkgs{k},'LP', 0);

    for j=1:length(QPsolverPkgs)%QP solvers

        solverQPOK = changeCobraSolver(QPsolverPkgs{j}, 'QP', 0);

        if solverLPOK && solverQPOK
            substrateRxns ={'EX_glc(e)' 'EX_ac(e)'};
            initConcentrations = [10.8; 1];
            initBiomass = 0.001;
            timeStep = 1/100;
            nSteps = 200;

             %Assert
            [concentrationMatrixtest,excRxnNamestest,timeVectest,biomassVectest,drGenestest,constrainedRxnstest,statestest] = ...
                dynamicRFBA(modelReg,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps);

            tol = eps(0.5);%set tolerance
            %assertions
            assert(any(any(abs(concentrationMatrixtest-concentrationMatrix) < tol)))
            assert(isequal(excRxnNamestest,excRxnNames))
            assert(isequal(timeVectest,timeVec))
            assert(any(abs(biomassVectest-biomassVec) < tol))
            assert(isequal(drGenestest,drGenes))
            assert(isequal(constrainedRxnstest,constrainedRxns))
            assert(isequal(statestest,states))
        end
    end

    fprintf('Done.\n');
end

 % change the directory
 cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function out = domainCheck(f, g)
%DOMAINCHECK   True if the domains of two CHEBFUN3T objects are the same.
%   DOMAINCHECK(F, G) returns TRUE if the domains of the two CHEBFUN3T
%   objects F and G coincide up to a tolerance depending on their
%   horizontal scales or if both F and G are empty CHEBFUN objects.
%
% See also CHEBFUN/DOMAINCHECK.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The CHEBFUN3T class uses this function internally to compare the domains of
% CHEBFUN3T objects before attempting to perform operations on multiple
% CHEBFUN3T objects that require the CHEBFUN3T objects to reside on the same interval.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Empty check:
if ( isempty(f) && isempty(g) )
    out = true;
    return
elseif ( xor(isempty(f), isempty(g) ) )
    out = false;
    return
end

% Compute INF norm of the domains:
hscaleF = norm(f.domain, inf);
hscaleG = norm(g.domain, inf);
hscale = max(hscaleF, hscaleG);

% Compare the domains:
err = f.domain - g.domain;

% Should be less than tolerance.
out = all(abs(err) < 1e-15*hscale);

end

author: chebfun

repo: chebfun


text:
function op = linop_vec( sz )
%LINOP_VEC Matrix to vector reshape operator
%OP = LINOP_VEC( SZ )
%    Constructs a TFOCS-compatible linear operator that reduces a matrix
%    variable to a vector version using column-major order.
%    This is equivalent to X(:)
%    The transpose operator will reshape a vector into a matrix.
%
%    The input SZ should of the form [M,N] where [M,N] describe the
%    size of the matrix variable.  The ouput vector will be of
%    length M*N.  If SZ is a single entry, then M = N is assumed.
%    For advanced usage with multidimensional arrays,
%    use the more general linop_reshape function instead.
%
%   To do the reverse operation (from vector to matrix),
%   use this function together with linop_adjoint.
%
%   See also linop_reshape

if numel(sz) > 2, error('must supply a 2-entry vector'); end
if numel(sz) == 1, sz = [sz(1),sz(1)]; end

% Switch conventions of the size variable:
sz = { [sz(1),sz(2)], [sz(1)*sz(2),1] };

op = @(x,mode)linop_handles_vec( sz, x, mode );

function y = linop_handles_vec(sz, x, mode )
switch mode,
    case 0, y = sz;
    case 1, y = x(:);
    case 2,
        MN = sz{1}; M = MN(1); N = MN(2);
        y = reshape( x, M, N );

end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

% % SynPerData
% hdcTrajFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\SynPerData\06_hdc_trajectory.txt';
%
%
% % SynPanData
% hdcTrajFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\SynPanData\06_hdc_trajectory.txt';

% QUTCarparkData
hdcTrajFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\QUTCarparData\06_hdc_trajectory.txt';

plot_hdc_traj_main(hdcTrajFile);

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function r8lib_test119 ( )

%*****************************************************************************80
%
%% R8LIB_TEST119 tests R8VEC_EVEN3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  nold = 4;
  nval = 12;

  xold = [ 0.0, 5.1, 7.0, 10.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8LIB_TEST119\n' );
  fprintf ( 1, '  R8VEC_EVEN3 tries to evenly interpolate new data\n' );
  fprintf ( 1, '  between old values.\n' );

  r8vec_print ( nold, xold, '  Original vector:' );

  xval = r8vec_even3 ( nold, nval, xold );

  r8vec_print ( nval, xval, '  New vector:' );

  return
end

author: johannesgerer

repo: jburkardt-m