text:
function read(obj,filename,varargin)
%READ Load sequence from file.
%   READ(seqObj, filename, ...) Read the given filename and load sequence
%   data into sequence object.
%
%   optional parwameter 'detectRFuse' can be given to let the function
%   infer the currently missing flags concerning the intended use of the RF
%   pulses (excitation, refocusing, etc). These are important for the
%   k-space trajectory calculation
%
%   Examples:
%   Load the sequence defined in gre.seq in my_sequences directory
%
%       read(seqObj,'my_sequences/gre.seq')
%
% See also  write

detectRFuse=false;
if ~isempty(varargin) && ~isempty(strfind(varargin{:},'detectRFuse'))
    detectRFuse=true;
end

fid = fopen(filename);

if fid<0
    error('filed to open file ''%s''', filename);
end

% Clear sequence data
%obj.blockEvents = [];
obj.blockEvents = {};
obj.definitions = containers.Map();
obj.gradLibrary = mr.EventLibrary();
obj.shapeLibrary = mr.EventLibrary();
obj.rfLibrary = mr.EventLibrary();
obj.adcLibrary = mr.EventLibrary();
%obj.delayLibrary = mr.EventLibrary();
obj.trigLibrary = mr.EventLibrary();
obj.labelsetLibrary = mr.EventLibrary();
obj.labelincLibrary = mr.EventLibrary();
obj.extensionStringIDs={};
obj.extensionNumericIDs=[];

version_combined=0;

% Load data from file
while true
    section = skipComments(fid);
    if section == -1
        break
    end

    switch section
        case '[DEFINITIONS]'
            obj.definitions = readDefinitions(fid);
            v=obj.getDefinition('GradientRasterTime');
            if ~isempty(v)
                obj.gradRasterTime=v;
            end
            v=obj.getDefinition('RadiofrequencyRasterTime');
            if ~isempty(v)
                obj.rfRasterTime=v;
            end
            v=obj.getDefinition('AdcRasterTime');
            if ~isempty(v)
                obj.adcRasterTime=v;
            end
            v=obj.getDefinition('BlockDurationRaster');
            if ~isempty(v)
                obj.blockDurationRaster=v;
            end
        case '[SIGNATURE]'
            tmpSignDefs = readDefinitions(fid);
            if isKey(tmpSignDefs,'Type')
                obj.signatureType=tmpSignDefs('Type');
            end
            if isKey(tmpSignDefs,'Hash')
                obj.signatureValue=tmpSignDefs('Hash');
                obj.signatureFile='Text'; % we are reading a text file, so much is known for sure
            end
        case '[VERSION]'
            [version_major, ...
             version_minor, ...
             version_revision] = readVersion(fid);
            assert(version_major == obj.version_major, ...
                    'Unsupported version_major %d', version_major)
            %
            version_combined=1000000*version_major+1000*version_minor+version_revision;
            %
            if version_combined < 1002000
                error('Unsupported version %d.%d.%d, only file format revision 1.2.0 and above are supported', version_major, version_minor, version_revision);
            end
            if version_combined < 1003001
                warning('Loading older Pulseq format file (version %d.%d.%d) some code may function not as expected', version_major, version_minor, version_revision);
            end
        case '[BLOCKS]'
            if ~exist('version_major')
                error('Pulseq file MUST include [VERSION] section prior to [BLOCKS] section');
            end
            [obj.blockEvents,obj.blockDurations,delayInd_tmp] = readBlocks(fid, obj.blockDurationRaster, version_combined);
        case '[RF]'
            if version_combined >= 1004000
                obj.rfLibrary = readEvents(fid, [1 1 1 1 1e-6 1 1]); % this is 1.4.x format
            else
                obj.rfLibrary = readEvents(fid, [1 1 1 1e-6 1 1]); % this is 1.3.x and below
                % we will have to scan through the library later after all the shapes have been loaded
            end
        case '[GRADIENTS]'
            if version_combined >= 1004000
                obj.gradLibrary = readEvents(fid, [1 1 1 1e-6], 'g' ,obj.gradLibrary); % this is 1.4.x format
            else
                obj.gradLibrary = readEvents(fid, [1 1 1e-6], 'g' ,obj.gradLibrary); % this is 1.3.x and below
            end
        case '[TRAP]'
            obj.gradLibrary = readEvents(fid, [1 1e-6 1e-6 1e-6 1e-6], 't', obj.gradLibrary);
        case '[ADC]'
            obj.adcLibrary = readEvents(fid, [1 1e-9 1e-6 1 1]);
        case '[DELAYS]'
            if version_combined >= 1004000
                error('Pulseq file revision 1.4.0 and above MUST NOT contain [DELAYS] section');
            end
            tmp_delayLibrary = readEvents(fid, 1e-6);
        case '[SHAPES]'
            obj.shapeLibrary = readShapes(fid, (version_major==1 && version_minor<4));
        case '[EXTENSIONS]'
            obj.extensionLibrary = readEvents(fid);
        otherwise
            if     strncmp('extension TRIGGERS', section, 18)
                id=str2num(section(19:end));
                obj.setExtensionStringAndID('TRIGGERS',id);
                obj.trigLibrary = readEvents(fid, [1 1 1e-6 1e-6]);
            elseif strncmp('extension LABELSET', section, 18)
                id=str2num(section(19:end));
                obj.setExtensionStringAndID('LABELSET',id);
                obj.labelsetLibrary = readAndParseEvents(fid,@str2num,@(s)find(strcmp(mr.getSupportedLabels,s)));
            elseif strncmp('extension LABELINC', section, 18)
                id=str2num(section(19:end));
                obj.setExtensionStringAndID('LABELINC',id);
                obj.labelincLibrary = readAndParseEvents(fid,@str2num,@(s)find(strcmp(mr.getSupportedLabels,s)));
            else
                error('Unknown section code: %s', section);
            end
    end
end
fclose(fid);

%
if version_combined < 1002000
    error('Unsupported version %07d, only file format revision 1.2.0 (1002000) and above are supported', version_combined);
end
% fix blocks, gradients and RF objects imported from older versions
if version_combined < 1004000  % MZ: FIXME: the code below does not update key-to-id mapping, so the libraries are not fully functional... these are only partially updated by the shaped gradients patching below (first/last detection)
    % scan through the RF objects
    for i=1:length(obj.rfLibrary.data)
        % % need to (partially) decode the magnitude shape to find out the pulse duration
        %magSamples = obj.shapeLibrary.data(obj.rfLibrary.data(i).array(2)).array(1);
        % % create time shape
        %timeShape = mr.compressShape((1:magSamples)-0.5); % time shape is stored in units of RF raster
        %data = [timeShape.num_samples timeShape.data];
        %timeID = obj.shapeLibrary.find_or_insert(data);
        obj.rfLibrary.data(i).array = [obj.rfLibrary.data(i).array(1:3), 0, obj.rfLibrary.data(i).array(4:end)];
        obj.rfLibrary.lengths(i) = obj.rfLibrary.lengths(i) + 1;
    end

    % scan through the gradient objects and update 't'-s (trapezoids) und 'g'-s (free-shape gradients)
    for i=1:length(obj.gradLibrary.data)
        if obj.gradLibrary.type(i)=='t'
            if obj.gradLibrary.data(i).array(2)==0
                if abs(obj.gradLibrary.data(i).array(1))==0 && obj.gradLibrary.data(i).array(3) > 0
                    obj.gradLibrary.data(i).array(3)=obj.gradLibrary.data(i).array(3)-obj.gradRasterTime;
                    obj.gradLibrary.data(i).array(2)=obj.gradRasterTime;
                end
            end
            if obj.gradLibrary.data(i).array(4)==0
                if abs(obj.gradLibrary.data(i).array(1))==0 && obj.gradLibrary.data(i).array(3) > 0
                    obj.gradLibrary.data(i).array(3)=obj.gradLibrary.data(i).array(3)-obj.gradRasterTime;
                    obj.gradLibrary.data(i).array(4)=obj.gradRasterTime;
                end
            end
        end
        if obj.gradLibrary.type(i)=='g'
            % % need to (partially) decode the shape to find out the duration
            %nSamples = obj.shapeLibrary.data(obj.gradLibrary.data(i).array(2)).array(1);
            % % create time shape
            %timeShape = mr.compressShape((1:nSamples)-0.5); % time shape is stored in units of grad raster
            %data = [timeShape.num_samples timeShape.data];
            %timeID = obj.shapeLibrary.find_or_insert(data);
            obj.gradLibrary.data(i).array = [obj.gradLibrary.data(i).array(1:2), 0, obj.gradLibrary.data(i).array(3:end)];
            obj.gradLibrary.lengths(i) = obj.gradLibrary.lengths(i) + 1;
        end
    end

    % for versions prior to 1.4.0 blockDurations have not been initialized
    obj.blockDurations=zeros(1,length(obj.blockEvents));
    % scan trhough blocks and calculate durations
    for iB = 1:length(obj.blockEvents)
        b=obj.getBlock(iB);
        if delayInd_tmp(iB) > 0
            b.delay.type = 'delay';
            b.delay.delay = tmp_delayLibrary.data(delayInd_tmp(iB)).array;
        end
        obj.blockDurations(iB)=mr.calcDuration(b);
    end
end

gradChannels={'gx','gy','gz'};
gradPrevLast=zeros(1,length(gradChannels));
for iB = 1:length(obj.blockEvents)
    b=obj.getBlock(iB);
    block_duration=obj.blockDurations(iB);
    %obj.blockDurations(iB)=block_duration;
    % we also need to keep track of the event IDs because some Pulseq files written by external software may contain repeated entries so searching by content will fail
    eventIDs=obj.blockEvents{iB};
    % update the objects by filling in the fields not contained in the
    % pulseq file
    for j=1:length(gradChannels)
        grad=b.(gradChannels{j});
        if isempty(grad)
            gradPrevLast(j)=0;
            continue;
        end
        if strcmp(grad.type,'grad')
            if grad.delay>0
                gradPrevLast(j)=0;
            end
            if isfield(grad,'first')
                continue;
            end
            grad.first = gradPrevLast(j);
            % is this an extended trapezoid?
            if grad.time_id~=0
                grad.last=grad.waveform(end);
                grad_duration=grad.delay+grad.tt(end);
            else
                % restore samples on the edges of the gradient raster intervals
                % for that we need the first sample
                odd_step1=[grad.first 2*grad.waveform'];
                odd_step2=odd_step1.*(mod(1:length(odd_step1),2)*2-1);
                waveform_odd_rest=(cumsum(odd_step2).*(mod(1:length(odd_step2),2)*2-1))';
                grad.last = waveform_odd_rest(end);
                grad_duration=grad.delay+length(grad.waveform)*obj.gradRasterTime;
            end
            % bookkeeping
            gradPrevLast(j) = grad.last;
            if grad_duration+eps<block_duration
                gradPrevLast(j)=0;
            end
            %b.(gradChannels{j})=grad;
            % update library object
% this does not work s we don't know how the amplitude was defined
%             amplitude = max(abs(grad.waveform));
%             if amplitude>0
%                 [~,~,fnz]=find(grad.waveform,1); % find the first non-zero value and make it positive
%                 amplitude=amplitude*sign(fnz);
%             end
            % need to recover the amplidute from the library data directly...
            id=eventIDs(j+2);
            amplitude=obj.gradLibrary.data(id).array(1);
            %
            if version_combined>=1004000
                old_data = [amplitude grad.shape_id grad.time_id grad.delay];
            else
                old_data = [amplitude grad.shape_id grad.delay];
            end
            new_data = [amplitude grad.shape_id grad.time_id grad.delay grad.first grad.last];
            update_data(obj.gradLibrary, id, old_data, new_data,'g');
        else
            gradPrevLast(j)=0;
        end
    end
    %% copy updated objects back into the event library
    %obj.setBlock(iB,b);

%for iB=1:size(obj.blockEvents,1)
%     % update the objects by filling in the fields not contained in the
%     % pulseq file
%     for j=1:length(gradChannels)
%         grad=b.(gradChannels{j});
%         if isempty(grad)
%             continue;
%         end
%         if strcmp(grad.type,'grad')
%             grad.first = grad.waveform(1); % MZ: eventually we should use extrapolation by 1/2 gradient rasters here
%             grad.last = grad.waveform(end);
%             b.(gradChannels{j})=grad;
%         end;
%     end
%     % copy updated objects back into the event library
%     obj.setBlock(iB,b);
end

if detectRFuse
    % find the RF pulses, list flip angles
    % and work around the current (rev 1.2.0) Pulseq file format limitation
    % that the RF pulse use is not stored in the file
    for k=obj.rfLibrary.keys
        libData=obj.rfLibrary.data(k).array;
        rf=obj.rfFromLibData(libData);
        %flipAngleDeg=abs(sum(rf.signal))*rf.t(1)*360; %we use rfex.t(1) in place of opt.system.rfRasterTime
        flipAngleDeg=abs(sum(rf.signal(1:end-1).*(rf.t(2:end)-rf.t(1:end-1))))*360;
        offresonance_ppm=1e6*rf.freqOffset/obj.sys.B0/obj.sys.gamma;
        % fix library %%%% if length(obj.rfLibrary.type)>=eventInd(2)
        if flipAngleDeg < 90.01 % we add 0.01 degree to account for rounding errors which we've experienced for very short RF pulses
            obj.rfLibrary.type(k) = 'e';
        else
            if rf.shape_dur > 6e-3 && offresonance_ppm >= -3.5 && offresonance_ppm <= -3.4 % approx -3.45 ppm
                obj.rfLibrary.type(k) = 's'; % saturation (fat-sat)
            else
                obj.rfLibrary.type(k) = 'r';
            end
        end
%         % fix libData
%         if length(libData) < 9
%             if flipAngleDeg < 90.01 % we add 0.01 degree to account for rounding errors which we've experienced for very short RF pulses
%                 libData(9) = 0; % or 1 ?
%             else
%                 libData(9) = 2; % or 1 ?
%             end
%             obj.rfLibrary.data(k).array=libData;
%         end
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Helper functions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function def = readDefinitions(fid)
        %readDefinitions Read the [DEFINITIONS] section of a sequence file.
        %   defs=readDefinitions(fid) Read user definitions from file
        %   identifier of an open MR sequence file and return a map of
        %   key/value entries.

        def = containers.Map;
        %line = strip(fgetl(fid));
        line = strip(skipComments(fid));
        while ischar(line) && ~(isempty(line) || line(1) == '#')
            tok = textscan(line, '%s');
            def(tok{1}{1}) = str2double(tok{1}(2:end));
            if ~all(isfinite(def(tok{1}{1})))
                def(tok{1}{1}) = strtrim(line((length(tok{1}{1})+2):end));
            end
            line = fgetl(fid);
        end
    end

    function [major, minor, revision] = readVersion(fid)
        %readVersion Read the [VERSION] section of a sequence file.
        %   defs=readVersion(fid) Read Pulseq version from file
        %   identifier of an open MR sequence file and return it

        major = [];
        minor = [];
        revision = [];
        line = fgetl(fid);
        while ischar(line) && ~(isempty(line) || line(1)=='#')
            tok = textscan(line,'%s');
            switch tok{1}{1}
                case 'major'
                    major = str2double(tok{1}(2:end));
                case 'minor'
                    minor = str2double(tok{1}(2:end));
                case 'revision'
                    revision = str2double(tok{1}(2:end));
            end
            line = fgetl(fid);
        end
    end

    function [eventTable,blockDurations,delayIDs_tmp] = readBlocks(fid, blockDurationRaster, version_combined)
        %readBlocks Read the [BLOCKS] section of a sequence file.
        %   library=readBlocks(fid) Read blocks from file identifier of an
        %   open MR sequence file and return the event table.

        eventTable = {};
        blockDurations = [];
        delayIDs_tmp = [];
        line = fgetl(fid);
        while ischar(line) && ~(isempty(line) || line(1) == '#')
            blockEvents = sscanf(line, '%f')';
            %eventTable = [eventTable; blockEvents(2:end)];
            if version_combined<=1002001
                eventTable{blockEvents(1)} = [0 blockEvents(3:end) 0];
            else
                eventTable{blockEvents(1)} = [0 blockEvents(3:end)];
            end
            if version_combined>=1004000
                blockDurations(blockEvents(1)) = blockEvents(2)*blockDurationRaster;
            else
                delayIDs_tmp(blockEvents(1)) = blockEvents(2);
            end
            line = fgetl(fid);
        end
    end

    function eventLibrary = readEvents(fid, scale, type, eventLibrary)
        %readEvents Read an event section of a sequence file.
        %   library=readEvents(fid) Read event data from file identifier of
        %   an open MR sequence file and return a library of events.
        %
        %   library=readEvents(fid,scale) Read event data and scale
        %   elements according to column vector scale.
        %
        %   library=readEvents(fid,scale,type) Attach the type string to
        %   elements of the library.
        %
        %   library=readEvents(...,library) Append new events to the given
        %   library.
        if nargin < 2
            scale = 1;
        end
        if nargin < 4
            eventLibrary = mr.EventLibrary();
        end
        line = fgetl(fid);
        while ischar(line) && ~(isempty(line) || line(1) == '#')
            data = sscanf(line,'%f')';
            id = data(1);
            data = scale.*data(2:end);
            if nargin < 3
                eventLibrary.insert(id, data);
            else
                eventLibrary.insert(id, data, type);
            end

            line=fgetl(fid);
        end
    end

    function eventLibrary = readAndParseEvents(fid, varargin)
        %readAndParseEvents Read an event section of a sequence file.
        %   library=readAndParseEvents(fid) Read event data from file
        %   identifier of an open MR sequence file and return a library of
        %   events.
        %
        %   library=readAndParseEvents(fid,parser1,parser2,...) Read event
        %   data and convert the elements using to the provided parser.
        %   Default parser is str2num()
        %
        eventLibrary = mr.EventLibrary();
        line = fgetl(fid);
        while ischar(line) && ~(isempty(line) || line(1) == '#')
            datas=regexp(line, '(\s+)','split');
            data=zeros(1,length(datas)-1);
            id = str2num(datas{1});
            for i=2:length(datas)
                if i>nargin
                    data(i-1) = str2num(datas{i});
                else
                    data(i-1) = varargin{i-1}(datas{i});
                end
            end
            eventLibrary.insert(id, data);

            line=fgetl(fid);
        end
    end

    function shapeLibrary = readShapes(fid, forceConvertUncompressed)
        %readShapes Read the [SHAPES] section of a sequence file.
        %   library=readShapes(fid) Read shapes from file identifier of an
        %   open MR sequence file and return a library of shapes.

        shapeLibrary = mr.EventLibrary();
        line = skipComments(fid);
        while ~(~ischar(line) || isempty(line) || ~strcmp(line(1:8), 'shape_id'))
            tok = textscan(line, '%s');
            id = str2double(tok{1}(2));
            line = skipComments(fid);
            tok = textscan(line, '%s');
            num_samples = str2double(tok{1}(2));
            data = [];
            line = skipComments(fid);   % first sample
            while ischar(line) && ~(isempty(line) || line(1) == '#')
                data = [data sscanf(line, '%f')];
                %data = [data single(sscanf(line, '%f'))]; % C-code uses single precision and we had problems already due to the rounding during reading in of the shapes...
                line = fgetl(fid);
            end

            line = skipComments(fid, true); % MZ: second parameter to prevent readShapes from reading into the next section (long-standing bug)

            % check if conversion is needed: in v1.4.x we use length(data)==num_samples
            % as a marker for the uncompressed (stored) data. In older versions this condition could occur by chance
            if forceConvertUncompressed && length(data)==num_samples
                shape.data=data;
                shape.num_samples=num_samples;
                shape = mr.compressShape(mr.decompressShape(shape,true));
                data = [shape.num_samples shape.data];
            else
                data = [num_samples data];
            end
            shapeLibrary.insert(id, data);
        end
    end

    function nextLine = skipComments(fid, stopBeforeSection)
        %skipComments Read lines of skipping blank lines and comments.
        %   line=skipComments(fid) Read lines from valid file identifer and
        %   return the next non-comment line.

        if (nargin<2)
            stopBeforeSection=false;
        end

        tmpPos=ftell(fid);
        line = fgetl(fid);
        while ischar(line) && (isempty(line) || line(1) == '#')
            tmpPos=ftell(fid);
            line = fgetl(fid);
        end
        if ischar(line)
            if stopBeforeSection && line(1)=='['
                fseek(fid,tmpPos,-1); % restore the file position
                nextLine = ''; % feasible (non-error) dummy return
            else
                nextLine = line;
            end
        else
            nextLine = -1;
        end
    end
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function varargout = acos(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('decreasing','callback');
        operator.derivative = @(x)real(-((1 - x.^2).^-0.5));
        operator.inverse = @(x)cos(x);
        operator.inflection = [-inf 1 0 -1];
        operator.range = [0 pi];
        operator.domain = [-1 1];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function cnn_imagenet(varargin)
%CNN_IMAGENET   Demonstrates training a CNN on ImageNet
%  This demo demonstrates training the AlexNet, VGG-F, VGG-S, VGG-M,
%  VGG-VD-16, and VGG-VD-19 architectures on ImageNet data.

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.dataDir = fullfile('data','ILSVRC2012') ;
opts.modelType = 'alexnet' ;
opts.networkType = 'simplenn' ;
opts.batchNormalization = true ;
opts.weightInitMethod = 'gaussian' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.modelType ;
if opts.batchNormalization, sfx = [sfx '-bnorm'] ; end
sfx = [sfx '-' opts.networkType] ;
opts.expDir = fullfile('data', ['imagenet12-' sfx]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.numFetchThreads = 12 ;
opts.lite = false ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.train = struct([]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                             Prepare model
% -------------------------------------------------------------------------

net = cnn_imagenet_init('model', opts.modelType, ...
                        'batchNormalization', opts.batchNormalization, ...
                        'weightInitMethod', opts.weightInitMethod, ...
                        'networkType', opts.networkType) ;

% -------------------------------------------------------------------------
%                                                              Prepare data
% -------------------------------------------------------------------------

if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
else
  imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% Set the class names in the network
net.meta.classes.name = imdb.classes.name ;
net.meta.classes.description = imdb.classes.description ;

% Compute image statistics (mean, RGB covariances, etc.)
imageStatsPath = fullfile(opts.expDir, 'imageStats.mat') ;
if exist(imageStatsPath)
  load(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
else
  [averageImage, rgbMean, rgbCovariance] = getImageStats(opts, net.meta, imdb) ;
  save(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
end

% Set the image average (use either an image or a color)
%net.meta.normalization.averageImage = averageImage ;
net.meta.normalization.averageImage = rgbMean ;

% Set data augmentation statistics
[v,d] = eig(rgbCovariance) ;
net.meta.augmentation.rgbVariance = 0.1*sqrt(d)*v' ;
clear v d ;

% -------------------------------------------------------------------------
%                                                                     Learn
% -------------------------------------------------------------------------

switch opts.networkType
  case 'simplenn', trainFn = @cnn_train ;
  case 'dagnn', trainFn = @cnn_train_dag ;
end

[net, info] = trainFn(net, imdb, getBatchFn(opts, net.meta), ...
                      'expDir', opts.expDir, ...
                      net.meta.trainOpts, ...
                      opts.train) ;

% -------------------------------------------------------------------------
%                                                                    Deploy
% -------------------------------------------------------------------------

net = cnn_imagenet_deploy(net) ;
modelPath = fullfile(opts.expDir, 'net-deployed.mat')

switch opts.networkType
  case 'simplenn'
    save(modelPath, '-struct', 'net') ;
  case 'dagnn'
    net_ = net.saveobj() ;
    save(modelPath, '-struct', 'net_') ;
    clear net_ ;
end

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
useGpu = numel(opts.train.gpus) > 0 ;

bopts.numThreads = opts.numFetchThreads ;
bopts.imageSize = meta.normalization.imageSize ;
bopts.border = meta.normalization.border ;
bopts.averageImage = meta.normalization.averageImage ;
bopts.rgbVariance = meta.augmentation.rgbVariance ;
bopts.transformation = meta.augmentation.transformation ;

switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(bopts,x,y) ;
  case 'dagnn'
    fn = @(x,y) getDagNNBatch(bopts,useGpu,x,y) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getSimpleNNBatch(opts, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

if ~isVal
  % training
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0) ;
else
  % validation: disable data augmentation
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0, ...
                              'transformation', 'none') ;
end

if nargout > 0
  labels = imdb.images.label(batch) ;
end

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, useGpu, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

if ~isVal
  % training
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0) ;
else
  % validation: disable data augmentation
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0, ...
                              'transformation', 'none') ;
end

if nargout > 0
  if useGpu
    im = gpuArray(im) ;
  end
  labels = imdb.images.label(batch) ;
  inputs = {'input', im, 'label', imdb.images.label(batch)} ;
end

% -------------------------------------------------------------------------
function [averageImage, rgbMean, rgbCovariance] = getImageStats(opts, meta, imdb)
% -------------------------------------------------------------------------
train = find(imdb.images.set == 1) ;
train = train(1: 101: end);
bs = 256 ;
opts.networkType = 'simplenn' ;
fn = getBatchFn(opts, meta) ;

for t=1:bs:numel(train)
  batch_time = tic ;
  batch = train(t:min(t+bs-1, numel(train))) ;
  fprintf('collecting image stats: batch starting with image %d ...', batch(1)) ;
  temp = fn(imdb, batch) ;
  z = reshape(permute(temp,[3 1 2 4]),3,[]) ;
  n = size(z,2) ;
  avg{t} = mean(temp, 4) ;
  rgbm1{t} = sum(z,2)/n ;
  rgbm2{t} = z*z'/n ;
  batch_time = toc(batch_time) ;
  fprintf(' %.2f s (%.1f images/s)\n', batch_time, numel(batch)/ batch_time) ;
end
averageImage = mean(cat(4,avg{:}),4) ;
rgbm1 = mean(cat(2,rgbm1{:}),2) ;
rgbm2 = mean(cat(3,rgbm2{:}),3) ;
rgbMean = rgbm1 ;
rgbCovariance = rgbm2 - rgbm1*rgbm1' ;

author: willard-yuan

repo: cnn-for-image-retrieval


function Svals = Simulate_RegimeSwitching_Diffusion_Unbiased( N_sim, T, S_0, drift_vec, sigma_vec, Q, initial_state)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Simulates Terminal Value (S_T) of Regime Switching Diffusion Models
%        This scheme is unbiased, but only generates the terminal S values
% Returns: vector of size N_sim
% Author: Justin Lars Kirkby
%
% -----------------
% Params
% -----------------
% N_sim = # samples
% T = time to maturity, ie path is on [0,T]
% S_0 = initial underlying value (e.g. S_0=100)
% drift_vec = vector of drift coefficient by regime state, e.g. r_i - q_i, where r is interest rate in state i, and q is div yield
% sigma_vec = diffusion coefficients in each state

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

Svals = zeros(N_sim,1);
lambdas = -diag(Q);  % Transition rates
P = (Q + diag(lambdas))./lambdas;  % Remove diagonal, no self-transition
cdfs = cumsum(P, 2);

for n = 1:N_sim
   S = S_0;
   t_last = 0;
   state = initial_state;
   while t_last < T
       tau = exprnd(lambdas(state));
       t = min(T, t_last + tau);
       tau = t - t_last;

       Sigsqdt = sigma_vec(state)*sqrt(tau);
       drift = (drift_vec(state) - 0.5*sigma_vec(state)^2)*tau;

       S = S.*exp(drift + Sigsqdt*randn());  %log scheme
       state = next_state(cdfs(state, :));
       t_last = t;
   end
   Svals(n) = S;
end

end

function j = next_state(cdf)
u = rand();
j = find(u <= cdf, 1);
end


text:
function [xi,w]=transformSimplexTriPts(xi,w,v)
%%TRANSFORMSIMPLEXTRIPTSTOTRI Given cubature points and weights for
%   integration over a the standard simplex triangle with vertices (1,0),
%   (0,1), and (0,0), transform the points and weights for integration over
%   the triangle with vertices given in v.
%
%INPUTS: xi The 2XnumCubPoints set of cubature points for the standard 2D
%           simplex triangle.
%         w The numCubPointsX1 associated cubature weights for the standard
%           2D simplex triangle.
%         v The 2X3 set of vertices of the new triangle.
%
%OUTPUTS: xi The 2XnumCubPoints set of cubature points for the specified
%            triangle.
%          w The 2XnumCubPoints set of transformed weights.
%
%With cubature formulae as in [1], the weights must sum to the area of the
%triangle. Thus, the first step here is to scale the weights by 2*the
%triangle area (the 2 is because the area of the basic simplex triangle is
%1/6). An affine transformation, as described in Chapter 5 of [2], is then
%performed to move the cubature weights into the new triangle.
%
%EXAMPLE:
%In this example, we find the integral of x^3*y^2 using fifth-order simplex
%points. The exact solution is 128/5 and one can see that the relative
%error between the solution found and the exact solution is within what one
%would expect due to finite precision limits.
% vertices=[0,2,4;
%           0,2,0];
% [xi,w]=fifthOrderSimplexCubPoints(2);
% [xi,w]=transformSimplexTriPts(xi,w,vertices);
% alpha=[3;2];
% cubSol=findMomentFromSamp(alpha,xi,w)
% exactSol=128/5
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%[2] K. Anjyo and H. Ochiai, Mathematical Basics of Motion and Deformation
%    in Computer Graphics, 2nd ed. Morgan and Claypool Publishers, 2017.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Adjust w for the triangle area.
T=triangleArea(v(:,1),v(:,2),v(:,3),true);
w=w*(2*T);

%Transforms the points into the new triangle.
A=[v(:,1)-v(:,3),v(:,2)-v(:,3)];
xi=A*xi+v(:,3);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function F = asin(F, varargin)
%ASIN   Inverse sine of a CHEBFUN.
%   ASIN(F) computes the inverse sine of the CHEBFUN F.
%
%   ASIN(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also SIN, ASIND.

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

% Call the compose method:
F = compose(F, @asin, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function update_background_parallel(obj, use_parallel)
%% update the background related variables in CNMF framework
% input:
%   use_parallel: boolean, do initialization in patch mode or not.
%       default(true); we recommend you to set it false only when you want to debug the code.

%% Author: Pengcheng Zhou, Columbia University, 2017
%% email: zhoupc1988@gmail.com

%% process parameters

try
    % map data
    mat_data = obj.P.mat_data;

    % folders and files for saving the results
    log_file =  obj.P.log_file;
    flog = fopen(log_file, 'a');
    log_data = matfile(obj.P.log_data, 'Writable', true); %#ok<NASGU>

    % dimension of data
    dims = mat_data.dims;
    d1 = dims(1);
    d2 = dims(2);
    T = dims(3);
    obj.options.d1 = d1;
    obj.options.d2 = d2;

    % parameters for patching information
    patch_pos = mat_data.patch_pos;
    block_pos = mat_data.block_pos;

    % number of patches
    [nr_patch, nc_patch] = size(patch_pos);
catch
    error('No data file selected');
end
fprintf('\n-----------------UPDATE BACKGROUND---------------------------\n');

% frames to be loaded for initialization
frame_range = obj.frame_range;
T = diff(frame_range) + 1;

% threshold for detecting large residuals
thresh_outlier = obj.options.thresh_outlier;

% use parallel or not
if ~exist('use_parallel', 'var')||isempty(use_parallel)
    use_parallel = true; %don't save initialization procedure
end

% options
options = obj.options;
nb = options.nb;
bg_ssub = options.bg_ssub;
bg_model = options.background_model;
with_projection = options.bg_acceleration;

% previous estimation
A = cell(nr_patch, nc_patch);
C = cell(nr_patch, nc_patch);
sn = cell(nr_patch, nc_patch);
W = obj.W;
b0 = obj.b0;
b = obj.b;
f = obj.f;
RSS = cell(nr_patch, nc_patch);

%% check whether the bg_ssub was changed
if strcmpi(bg_model, 'ring')
    tmp_block = block_pos{1};    % block position
    nr_block = diff(tmp_block(1:2))+1;
    nc_block = diff(tmp_block(3:4))+1;
    [~, temp] = size(W{1});
    [d1s, d2s] = size(imresize(zeros(nr_block, nc_block), 1/bg_ssub));
    if temp~=d1s*d2s
        rr = ceil(obj.options.ring_radius/bg_ssub);    % radius of the ring
        [r_shift, c_shift] = get_nhood(rr, obj.options.num_neighbors);    % shifts used for acquiring the neighboring pixels on the ring
        parfor mpatch=1:(nr_patch*nc_patch)
            tmp_patch = patch_pos{mpatch};    % patch position
            tmp_block = block_pos{mpatch};    % block position
            nr = diff(tmp_patch(1:2)) + 1;
            nc = diff(tmp_patch(3:4)) + 1;
            nr_block = diff(tmp_block(1:2))+1;
            nc_block = diff(tmp_block(3:4))+1;
            b0{mpatch} = zeros(nr*nc, 1);

            if bg_ssub==1
                [csub, rsub] = meshgrid(tmp_patch(3):tmp_patch(4), tmp_patch(1):tmp_patch(2));
                csub = reshape(csub, [], 1);
                rsub = reshape(rsub, [], 1);
                ii = repmat((1:numel(csub))', [1, length(r_shift)]);
                csub = bsxfun(@plus, csub, c_shift);
                rsub = bsxfun(@plus, rsub, r_shift);
                ind = and(and(csub>=1, csub<=d2), and(rsub>=1, rsub<=d1));
                jj = (csub-tmp_block(3)) * (diff(tmp_block(1:2))+1) + (rsub-tmp_block(1)+1);

                temp = sparse(ii(ind), jj(ind), 1, nr*nc, nr_block*nc_block);
                W{mpatch} = bsxfun(@times, temp, 1./sum(temp, 2));
            else
                d1s = ceil(nr_block/bg_ssub);
                d2s = ceil(nc_block/bg_ssub);

                [csub, rsub] = meshgrid(1:d2s, 1:d1s);
                csub = reshape(csub, [], 1);
                rsub = reshape(rsub, [], 1);
                ii = repmat((1:numel(csub))', [1, length(r_shift)]);
                csub = bsxfun(@plus, csub, c_shift);
                rsub = bsxfun(@plus, rsub, r_shift);
                jj = (csub-1) * d1s + rsub;
                % remove neighbors that are out of boundary
                ind = and(and(csub>=1, csub<=d2s), and(rsub>=1, rsub<=d1s));
                temp = sparse(ii(ind), jj(ind), 1, d1s*d2s, d1s*d2s);
                W{mpatch} = bsxfun(@times, temp, 1./sum(temp, 2));
            end
        end
    end
end

%% start updating the background
for mpatch=1:(nr_patch*nc_patch)
    tmp_block = block_pos{mpatch};

    % find the neurons that are within the block
    mask = zeros(d1, d2);
    mask(tmp_block(1):tmp_block(2), tmp_block(3):tmp_block(4)) = 1;

    ind = (reshape(mask(:), 1, [])* obj.A>0);
    A{mpatch}= obj.A(logical(mask), ind);
    C{mpatch} = obj.C(ind, :);
    temp = obj.P.sn(logical(mask));
    if bg_ssub==1
        sn{mpatch} = temp;
    else
        nr_block = diff(tmp_block(1:2))+1;
        nc_block = diff(tmp_block(3:4))+1;
        sn{mpatch} = imresize(reshape(temp, nr_block, nc_block), 1/bg_ssub, 'nearest')*bg_ssub;
    end
end

% check whether this is the first run of updating background components
if strcmpi(bg_model, 'ring')
    flag_first = (length(unique(W{1}(1, :)))==2);
else
    flag_first = (mean2(b{1})==0);
end

if use_parallel
    % load data before running parfor
    %     data_patch = cell(nr_patch, nc_patch);
    %     flag_ignore = cell(nr_patch, nc_patch);
    %     for mpatch=1:(nr_patch*nc_patch)
    %         tmp_patch = patch_pos{mpatch};
    %         A_block = A{mpatch};
    %
    %         % stop updating B because A&C doesn't change in this area
    %         if isempty(A_block) && (~flag_first)
    %             flag_ignore{mpatch} = true;
    %             [r, c] = ind2sub([nr_patch, nc_patch], mpatch);
    %
    %             % keep the current results. this step looks rediculous, but it
    %             % is needed for some computer/matlab. very weird.
    %             W{mpatch} = W{mpatch};
    %             b0{mpatch} = b0{mpatch};
    %             b{mpatch} = b{mpatch};
    %             f{mpatch} = f{mpatch};
    %             fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
    %             continue;
    %         else
    %             flag_ignore{mpatch} = false;
    %             % pull data
    %             data_patch{mpatch} = get_patch_data(mat_data, tmp_patch, frame_range, true);
    %         end
    %     end
    % do the actual computation
    parfor mpatch=1:(nr_patch*nc_patch)
        %         if flag_ignore{mpatch}
        %             continue;
        %         end
        tmp_patch = patch_pos{mpatch};
        tmp_block = block_pos{mpatch};

        A_block = A{mpatch};
        sn_block = sn{mpatch};
        C_block = C{mpatch};

        % stop updating B because A&C doesn't change in this area
        if isempty(A_block) && (~flag_first)
            [r, c] = ind2sub([nr_patch, nc_patch], mpatch);

            % keep the current results. this step looks rediculous, but it
            % is needed for some computer/matlab. very weird.
            W{mpatch} = W{mpatch};
            b0{mpatch} = b0{mpatch};
            b{mpatch} = b{mpatch};
            f{mpatch} = f{mpatch};
            fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
            continue;
        end

        % use ind_patch to indicate pixels within the patch and only
        % update (W, b0) corresponding to these pixels
        ind_patch = false(diff(tmp_block(1:2))+1, diff(tmp_block(3:4))+1);
        ind_patch((tmp_patch(1):tmp_patch(2))-tmp_block(1)+1, (tmp_patch(3):tmp_patch(4))-tmp_block(3)+1) = true;

        % pull data
        %         Ypatch = data_patch{mpatch};
        Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
        [nr_block, nc_block, T_block] = size(Ypatch);
        if strcmpi(bg_model, 'ring')
            % get the previous estimation
            W_old = W{mpatch};
            Ypatch = reshape(Ypatch, [], T);

            % run regression to get A, C, and W, b0
            if bg_ssub==1
                sn_patch = sn_block(ind_patch);
                [W{mpatch}, b0{mpatch}] = fit_ring_model(Ypatch, A_block, C_block, W_old, thresh_outlier, sn_patch, ind_patch, with_projection);
            else
                % downsapmle data first
                temp = reshape(double(Ypatch)-A_block*C_block, nr_block, nc_block, T_block);
                tmp_b0 = mean(temp, 3);
                b0{mpatch} = tmp_b0(ind_patch);
                Ypatch = imresize(temp, 1./bg_ssub, 'nearest');
                Ypatch = reshape(Ypatch, [], T_block);

                [W{mpatch}, ~] = fit_ring_model(Ypatch, [], [], W_old, thresh_outlier, sn_block(:), [],  with_projection);
                %                 tmp_b0 = imresize(reshape(tmp_b0, size(sn_block)), [nr_block, nc_block]);
                %                 b0{mpatch} = tmp_b0(ind_patch(:));
            end
        elseif strcmpi(bg_model, 'nmf')
            b_old = b{mpatch};
            f_old = f{mpatch};
            Ypatch = reshape(Ypatch, [], T);
            sn_patch = sn_block(ind_patch);
            [b{mpatch}, f{mpatch}] = fit_nmf_model(Ypatch, nb, A_block, C_block, b_old, f_old, thresh_outlier, sn_patch, ind_patch);
        else
            b_old = b{mpatch};
            f_old = f{mpatch};
            Ypatch = reshape(Ypatch, [], T);
            sn_patch = sn_block(ind_patch);
            [b{mpatch}, f{mpatch}, b0{mpatch}] = fit_svd_model(Ypatch, nb, A_block, C_block, b_old, f_old, thresh_outlier, sn_patch, ind_patch);
        end
        [r, c] = ind2sub([nr_patch, nc_patch], mpatch);
        fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);

    end
else
    for mpatch=1:(nr_patch*nc_patch)
        tmp_patch = patch_pos{mpatch};
        tmp_block = block_pos{mpatch};

        A_block = A{mpatch};
        sn_block = sn{mpatch};
        C_block = C{mpatch};

        % stop the updating B because A&C doesn't change in this area
        if isempty(A_block) && (~flag_first)
            [r, c] = ind2sub([nr_patch, nc_patch], mpatch);
            fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
            continue;
        end

        % use ind_patch to indicate pixels within the patch and only
        % update (W, b0) corresponding to these pixels
        ind_patch = false(diff(tmp_block(1:2))+1, diff(tmp_block(3:4))+1);
        ind_patch((tmp_patch(1):tmp_patch(2))-tmp_block(1)+1, (tmp_patch(3):tmp_patch(4))-tmp_block(3)+1) = true;

        % pull data
        Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
        [nr_block, nc_block, T_block] = size(Ypatch);
        if strcmpi(bg_model, 'ring')
            % get the previous estimation
            W_old = W{mpatch};
            Ypatch = reshape(Ypatch, [], T_block);

            % run regression to get A, C, and W, b0
            if bg_ssub==1
                sn_patch = sn_block(ind_patch);
                [W{mpatch}, b0{mpatch}] = fit_ring_model(Ypatch, A_block, C_block, W_old, thresh_outlier, sn_patch, ind_patch, with_projection);
            else
                % downsapmle data first
                temp = reshape(double(Ypatch)-A_block*C_block, nr_block, nc_block, T_block);
                tmp_b0 = mean(temp, 3);
                b0{mpatch} = tmp_b0(ind_patch);
                Ypatch = imresize(temp, 1./bg_ssub, 'nearest');
                Ypatch = reshape(Ypatch, [], T_block);

                [W{mpatch}, ~] = fit_ring_model(Ypatch, [], [], W_old, thresh_outlier, sn_block(:), [],  with_projection);
                %                 tmp_b0 = imresize(reshape(tmp_b0, size(sn_block)), [nr_block, nc_block]);
                %                 b0{mpatch} = tmp_b0(ind_patch(:));
            end
        elseif strcmpi(bg_model, 'nmf')
            b_old = b{mpatch};
            f_old = f{mpatch};
            Ypatch = reshape(Ypatch, [], T);
            sn_patch = sn_block(ind_patch);
            [b{mpatch}, f{mpatch}] = fit_nmf_model(Ypatch, nb, A_block, C_block, b_old, f_old, thresh_outlier, sn_patch, ind_patch);
        else
            b_old = b{mpatch};
            f_old = f{mpatch};
            Ypatch = reshape(Ypatch, [], T);
            sn_patch = sn_block(ind_patch);
            [b{mpatch}, f{mpatch}, b0{mpatch}] = fit_svd_model(Ypatch, nb, A_block, C_block, b_old, f_old, thresh_outlier, sn_patch, ind_patch);
        end
        [r, c] = ind2sub([nr_patch, nc_patch], mpatch);
        fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);

    end
end
obj.b = b;
obj.f = f;
obj.b0 = b0;
obj.W = W;
obj.b0_new = obj.reconstruct_b0();
obj.A_prev = obj.A;
obj.C_prev = obj.C;

%% save the results to log
fprintf('Finished updating background using %s model.\n', bg_model);

fprintf(flog, '[%s]\b', get_minute());
fprintf(flog, 'Finished updating background using %s model.\n', bg_model);
if obj.options.save_intermediate
    bg.b = obj.b;
    bg.f = obj.f;
    bg.b0 = obj.b0;
    bg.W = obj.W;
    tmp_str = get_date();
    tmp_str=strrep(tmp_str, '-', '_');
    eval(sprintf('log_data.bg_%s = bg;', tmp_str));
    fprintf(flog, '\tThe results were saved as intermediate_results.bg_%s\n\n', tmp_str);
end
fclose(flog);

author: zhoupc


text:
function niederreiter2_test ( )

%*****************************************************************************80
%
%% NIEDERREITER2_TEST tests the NIEDERREITER2 library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NIEDERREITER2_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the NIEDERREITER2 library.\n' );

  niederreiter2_test01 ( );
  niederreiter2_test02 ( );
  niederreiter2_test03 ( );
  niederreiter2_test04 ( );
  niederreiter2_test05 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NIEDERREITER2_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function res = spm_eeg_regressors_movement_ctf(S)
% Generate movement regressors for CTF MEG data
% S                     - input structure
% fields of S:
%    S.D                - M/EEG object
%
%    Additional parameters can be defined specific for each plugin
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided the plugin returns
%______________________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak

SVNrev = '$Rev: 6186 $';

if nargin == 0

    Dmov        = cfg_files;
    Dmov.tag    = 'Dmov';
    Dmov.name   = 'Movement dataset name';
    Dmov.filter = 'mat';
    Dmov.num    = [1 1];
    Dmov.help   = {'Select the M/EEG mat file containing continuous head localisation data.',...
        'This might or might not be the same as the dataset being analysed.'};

    movement_ctf = cfg_branch;
    movement_ctf.tag = 'movement_ctf';
    movement_ctf.name = 'CTF head movements';
    movement_ctf.val = {Dmov};

    res = movement_ctf;

    return
end

%read HLC-channels
%HLC0011 HLC0012 HLC0013 x, y, z coordinates of nasion-coil in m.
%HLC0021 HLC0022 HLC0023 x, y, z coordinates of lpa-coil in m.
%HLC0031 HLC0032 HLC0033 x, y, z coordinates of rpa-coil in m.
hlc_chan_label = {'HLC0011' 'HLC0012' 'HLC0013'...
    'HLC0021' 'HLC0022' 'HLC0023'...
    'HLC0031' 'HLC0032' 'HLC0033'};

%-Startup
%--------------------------------------------------------------------------
spm('sFnBanner', mfilename, SVNrev);
spm('FigName','CTF movement regressors');

if iscell(S.Dmov)
    S.Dmov = char(S.Dmov);
end

Dmov = spm_eeg_load(S.Dmov);
D    = spm_eeg_load(S.D);

hlc_chan_ind = Dmov.indchannel(hlc_chan_label);

if length(hlc_chan_ind) ~= 9
    error('Nine CTF HLC channels should be present to define movement regressors.');
end

if ~isequal(Dmov.transformtype, 'time')
    error('The movement dataset should be in the time domain.');
end

if isequal(D.type, 'continuous')

    if ~isequal(Dmov.type, 'continuous') || (D.time(1) < Dmov.time(1)) || (D.time(end)>Dmov.time(end))
        error('All times of the input dataset should be within the movement dataset.');
    end

    data = Dmov(hlc_chan_ind, :);

    if D.fsample ~= Dmov.fsample
        [data, alpha] = spm_timeseries_resample(data, D.fsample/Dmov.fsample);
    else
        alpha = 1;
    end

    start = round(alpha*Dmov.indsample(D.time(1)));

    data = data(:, start:(start+D.nsamples-1));

else
    if D.ntrials ~= Dmov.ntrials
        error('Trial numbers should be equal between input and movement dataset.');
    end

    data = Dmov(hlc_chan_ind, :, :);

    if S.summarise
        data = spm_squeeze(mean(data, 2), 2);
    else
        data = reshape(data, size(data, 1), []);
    end

end

res.R     = hpi2mov(data);

res.names = {'x', 'y', 'z', 'pitch', 'roll', 'yaw'};

spm('FigName','CTF movement regressors: done');

function P = hpi2mov(data)

% John Ashburner

ref   = reshape(data(:,1),3,3)*1000;
muRef = mean(ref,2);
ref0  = ref - repmat(muRef, 1, 3);
P     = zeros(size(data,2),6);

ns   = size(data, 2);

spm_progress_bar('Init', ns, 'Computing movement parameters'); drawnow;
if ns > 100, Ibar = floor(linspace(1, ns,100));
else Ibar = 1:ns; end

for i=1:ns
    cur     = reshape(data(:,i),3,3)*1000;
    muCur   = mean(cur,2);
    cur0    = cur - repmat(muCur, 1, 3);
    [U,SS,V] = svd(cur0*ref0');
    Ri      = U*V';
    if det(Ri)<0, V(:,end)=-V(:,end); Ri=U*V'; end
    p       = spm_imatrix([Ri muCur-Ri*muRef; 0 0 0 1]);
    P(i,:)  = p(1:6);

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

P = detrend(P, 'constant');
P = P./repmat(std(P), size(P, 1), 1);

author: spm


text:
function EEG = hlp_blockAverageConnectivity(EEG,numblocks)
    % INPUTS:
    %
    %   EEG         - EEG structure containing EEG.CAT.Conn
    %
    % Optional:
    %
    %   numblocks   - number of blocks to average over
    %
    % OUTPUTS:
    %
    %   EEG     - contained modifed Conn structure

    % average connectivity across trials or blocks
    % trials/blocks are assumed to be concatenated columnwise in last
    % dimension of EEG.CAT.Conn.<connmethod>
    % e.g. EEG.CAT.Conn.<connmethod> is assumed to be [nchs x nchs x nfreqs x numblocks*blocklength]
    % numblocks must exactly devide the number of time points in Conn matrices

    % References:
    %
    % [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
    %     Theoretical Handbook and User Manual. Chapter 4
    %     Available at: http://www.sccn.ucsd.edu/wiki/Sift
    %
    % Author: Tim Mullen 2011, SCCN/INC, UCSD.
    % Email:  tim@sccn.ucsd.edu

    % This function is part of the Source Information Flow Toolbox (SIFT)
    %
    % 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, write to the Free Software
    % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    connmethod = hlp_getConnMethodNames(EEG.CAT.Conn);

    if nargin < 2
        % default number of blocks is number of trials
        numblocks = size(EEG.CAT.srcdata,3);
    end

    [nchs nchs nfreq npoints] = size(EEG.CAT.Conn.(connmethod{1}));

    for m =1:length(connmethod)
        C = EEG.CAT.Conn.(connmethod{m});
        C = reshape(C,[nchs nchs nfreq npoints/numblocks numblocks]);
        EEG.CAT.Conn.(connmethod{m}) = mean(C,5);

    end

    EEG.CAT.Conn.erWinCenterTimes = EEG.CAT.Conn.erWinCenterTimes(1:npoints/numblocks);

repo: REST


text:
function full_forward()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.activations{curr_layer_idx} = bsxfun(@plus, config.weights{curr_layer_idx} * mem.layer_inputs{curr_layer_idx}, config.weights{config.layer_num+curr_layer_idx});
    config.misc.current_layer = curr_layer_idx + 1;
end

author: jimmy-ren


text:
function bvr_sendcommand(fcn, varargin)
%BVR_SENDCOMMAND - Control the BrainVision Recorder software
%
%Synopsis:
% bvr_sendcommand(FCN, <ARGS>)
%
%Arguements:
% FCN: Name of the function to be executed. See in the bvr_*.vbs files in
%     acquistion_tools folder for a list of options. Here are some:
%     'loadworkspace' - load BV workspace into the recorder; ARG: name of
%        the workspace (extension '.rwsp' may be left out)
%     'startrecording' - Start EEG recording; ARG: name of the file with
%        full path, without extension.
%     'startimprecording' - Make Impedance measurement first and start
%        recording afterward (impedance values are saved into the EEG
%        file); ARG as above.
%     'stoprecording' - Stops the recording.
%     'viewsignals' - Switch to monitoring state
%     'viewsignalsandwait' - Switch to monitoring mode and wait (unless monitoring
%        mode is already active). Example: bvr_sendcommand('viewsignalsandwait','3000');

mypath= fileparts(which(mfilename));
vbs_function= fullfile(mypath, ['bvr_' fcn '.vbs']);

system(['"' vbs_function '"' sprintf(' %s',varargin{:})]);

author: bbci


text:
function im = bbox_draw(im,boxes,c,t)

% copied from Ross Girshick
% Fast R-CNN
% Copyright (c) 2015 Microsoft
% Licensed under The MIT License [see LICENSE for details]
% Written by Ross Girshick
% --------------------------------------------------------
% source: https://github.com/rbgirshick/fast-rcnn/blob/master/matlab/showboxes.m
%
%
% Fast R-CNN
%
% Copyright (c) Microsoft Corporation
%
% All rights reserved.
%
% MIT License
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED *AS IS*, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
% OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
% ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
% OTHER DEALINGS IN THE SOFTWARE.

image(im);
axis image;
axis off;
set(gcf, 'Color', 'white');

if nargin<3
  c = 'r';
  t = 2;
end

s = '-';
if ~isempty(boxes)
    x1 = boxes(:, 1);
    y1 = boxes(:, 2);
    x2 = boxes(:, 3);
    y2 = boxes(:, 4);
    line([x1 x1 x2 x2 x1]', [y1 y2 y2 y1 y1]', ...
        'color', c, 'linewidth', t, 'linestyle', s);
    for i = 1:size(boxes, 1)
        text(double(x1(i)), double(y1(i)) - 2, ...
            sprintf('%.4f', boxes(i, end)), ...
            'backgroundcolor', 'b', 'color', 'w', 'FontSize', 10);
    end
end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function spm_MDP_da(MDP)
% Simulated histograms of dopamine firing
% FORMAT spm_MDP_da(MDP)
%
% See also: spm_MDP_game, which generalises this scheme and replaces prior
% beliefs about KL control with minimisation of expected free energy.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_da.m 6451 2015-05-26 09:26:03Z karl $

% deconvolve to simulate dopamine responses
%--------------------------------------------------------------------------
wn   = MDP.d;
da   = pinv( tril(toeplitz(exp(-((1:length(wn)) - 1)'/8))) )*wn;
da   = da(4:end);

% peristimulus time histogram
%--------------------------------------------------------------------------
pst  = (1:length(da))*256;
r    = 256;
psth = r*da + randn(size(da)).*sqrt(r*da);
bar(pst,psth,1)
title('Simulated dopamine responses','FontSize',16)
xlabel('Peristimulus time (ms)','FontSize',12)
ylabel('Spikes per bin','FontSize',12)
axis([pst(1) pst(end) 0 max(psth)*(1 + 1/4)])
axis square


text:
function dlap_rhs_write ( iunit, n, rhs )

%*****************************************************************************80
%
%% DLAP_RHS_WRITE writes a right hand side vector to a DLAP file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Mark Seager,
%    A SLAP for the Masses,
%    Lawrence Livermore National Laboratory,
%    Technical Report UCRL-100267, December 1988.
%
%  Parameters:
%
%    Input, integer IUNIT, the device to which
%    the matrix information is to be written.
%
%    Input, integer N, the order of the matrix.
%
%    Input, real RHS(N), the right hand side vector to be written.
%
  fprintf ( iunit, ' %16.7e\n', rhs(1:n) );

  return
end


text:
%Program for Fusing 2 images

%Author : Athi Narayanan S
%M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://sites.google.com/site/athisnarayanan/
%s_athi1983@yahoo.co.in

%Program Description
%This program is the main entry of the application.
%This program fuses/combines 2 images
%It supports both Gray & Color Images
%Alpha Factor can be varied to vary the proportion of mixing of each image.
%With Alpha Factor = 0.5, the two images mixed equally.
%With Alpha Facotr < 0.5, the contribution of background image will be more.
%With Alpha Facotr > 0.5, the contribution of foreground image will be more.

function fusedImg = FuseImages(bgImg, fgImg, alphaFactor)

bgImg = double(bgImg);
fgImg = double(fgImg);

fgImgAlpha = alphaFactor .* fgImg;
bgImgAlpha = (1 - alphaFactor) .* bgImg;

fusedImg = fgImgAlpha + bgImgAlpha;

author: Sable

repo: mcbench-benchmarks


text:
function Grasping_and_twist

%% Grasping and twisting Function with EEG recording
clear; close; clc;
%% Initializing EEG recording device
bbci_acquire_bv('close');
startup_bbci;

%% RoboticArm communication channel opened
jc = JacoComm;
connect(jc);
calibrateFingers(jc);

%% Query individual object properties
jc.JointPos
%%
jc.JointVel
%%
jc.JointTorque
%%
jc.JointTemp
%%
jc.FingerPos
%%
jc.FingerVel
%%
jc.FingerTorque
%%
jc.FingerTemp
%%
jc.EndEffectorPose
%%
jc.EndEffectorWrench
%%
jc.ProtectionZone
%%
jc.EndEffectorOffset
%%
jc.DOF
%%
jc.TrajectoryInfo

%% Methods to query joint and finger values all at once
pos = getJointAndFingerPos(jc);
vel = getJointAndFingerVel(jc);

torque = getJointAndFingerTorque(jc);
temp = getJointAndFingerTemp(jc);

setPositionControlMode(jc);
goToHomePosition(jc);

current_pos=jc.EndEffectorPose;
home_pos=jc.EndEffectorPose;
previous_pos=current_pos;

setPositionControlMode(jc);
fCmd = 0*ones(3,1);
sendFingerPositionCommand(jc,fCmd);

%% Offline model training
global EEG_MAT_DIR
EEG_MAT_DIR = '';

dd = 'MotorImagery Converted Data\';
filelist= {'20191113_demo_bhkwon'};

%% Preprocessing Offline Data
Bandpass_Filter = [4 40];

% Offline Classifier Training
[cnt, mrk, mnt] = eegfile_loadMatlab([dd filelist{1}]);

ival = [0 3000]; % 3s interval

cnt_filt = proc_filtButter(cnt, 5, Bandpass_Filter);
epo = cntToEpo(cnt_filt,mrk,ival);

[fv, Out.csp_w]=  proc_multicsp(epo, 3);
fv = proc_variance(fv); fv= proc_logarithm(fv);

fv.classifier_param = {'RLDAshrink','prior', nan, 'store_prior', 1, 'store_means', 1, ...
    'store_cov', 1, 'store_invcov', 1, 'scaling', 1};

proc = {'wr_multiClass','policy','one-vs-all','coding','hamming'}; % one-vs-all all-pairs

Out.C = trainClassifier(fv, proc);
Out.out_eeg = applyClassifier(fv, 'wr_multiClass', Out.C);

%% Onling Initialization
params = struct;
state = bbci_acquire_bv('init', params);
EEG_data = [];
mnt = getElectrodePositions(state.clab);

epo.clab = state.clab;
epo.fs = state.fs;
epo.title = filelist{1};
Ans = zeros(1,3);
%% Get EEG data
py.SharedDemo_6_Vision_LBH.getPerspectiveMine(); % vision model upload

coord=py.SharedDemo_6_Vision_LBH.mainRunning(1);
coord=double(py.array.array('d',py.numpy.nditer(coord)));
coord=uint8(coord/5);
pause(1);

home_pos=jc.EndEffectorPose;
current_pos=jc.EndEffectorPose;
prev_pos=current_pos;

targetX = coord(1);
targetY = coord(2);
targetX = double(targetX);
targetY = double(targetY);
x_origin = (102 - targetX)/100
y_origin = (5+targetY)/100

desired_pos=[x_origin; -y_origin; 0.1; home_pos(4); home_pos(5); home_pos(6)];
moveToCP(jc,desired_pos);
% YOLO to RoboticArm coordinate axis
while 1
    %% Sound cue
    [A,AFs] = audioread('grasp.mp3');
    sound(A,AFs);
    pause(2);

    eog_ch = [1 ,31];
    eog_th = 18;
    time_window = 20;

    eog_test(eog_ch, eog_th, time_window);
    disp('Receiving brain signal')

    pause(2);

    %% Initializing EEG recording
    % Initialization communication channel
    bbci_acquire_bv('close');
    EEG_MAT_DIR = '';
    params=struct;
    state=bbci_acquire_bv('init',params);
    EEGData=[];
    mnt=getElectrodePositions(state.clab);
    epo.clab = state.clab;
    epo.fs = state.fs;
    epo.title = filelist{1};
    pause('on')

    [B,BFs] = audioread('censor-beep-4.wav');
    sound(B,BFs);
    pause(4);

    data = bbci_acquire_bv(state);
    EEG_data = [EEG_data; data];

    % Sampling Rate: 250Hz, get data every 3 sec
    if size(EEG_data,1) >= 750
        epo.x = EEG_data;

        % EEG filtering
        Wps= [42 49]/epo.fs*2;
        [n, Ws]= cheb2ord(Wps(1), Wps(2),3, 40);
        [filt.b, filt.a]= cheby2(n, 50, Ws);
        epo = proc_filt(epo, filt.b, filt.a);

        % Feature Extraction and Classification
        Classification_Result = MotorImagery_Online_Fn(epo, Bandpass_Filter, Out);

        disp('Signal processing');
        pause(2);
        disp('Decoding signal');
        pause(2);

        disp('Grasp');

        disp('robotic arm activation');
        pause(2)

        setPositionControlMode(jc);
        fCmd = 4000*ones(3,1);
        sendFingerPositionCommand(jc,fCmd);

        pause(1);

        desired_pos=[0.6; -0.2; 0.15; home_pos(4); home_pos(5); home_pos(6)];
        moveToCP(jc,desired_pos);

        pause(1);

    end

    [B,BFs] = audioread('drinkwater.mp3');
    sound(B,BFs);
    pause(2);

    eog_ch = [1 ,31];
    eog_th = 18;
    time_window = 20;

    eog_test(eog_ch, eog_th, time_window);
    disp('Receiving brain signal')

    bbci_acquire_bv('close');
    EEG_MAT_DIR = '';
    params=struct;
    state=bbci_acquire_bv('init',params);
    EEGData=[];
    data=[];
    mnt=getElectrodePositions(state.clab);
    epo.clab = state.clab;
    epo.fs = state.fs;
    epo.title = filelist{1};

    [B,BFs] = audioread('censor-beep-4.wav');
    sound(B,BFs);
    pause(4);

    data = bbci_acquire_bv(state);
    EEG_data = [EEG_data; data];

    if size(EEG_data,1) >= 750
        epo.x = EEG_data;

        % EEG filtering
        Wps= [42 49]/epo.fs*2;
        [n, Ws]= cheb2ord(Wps(1), Wps(2),3, 40);
        [filt.b, filt.a]= cheby2(n, 50, Ws);
        epo = proc_filt(epo, filt.b, filt.a);

        % Feature Extraction and Classification
        Classification_Result = MotorImagery_Online_Fn(epo, Bandpass_Filter, Out);

        disp('Signal processing');
        pause(2);
        disp('Decoding signal');
        pause(2);

        desired_pos=[0.6; 0.15; 0.15; home_pos(4); home_pos(5); home_pos(6)];
        moveToCP(jc,desired_pos);

        jntVelCmd = [0;0;0;0;0;0;0.2]; %7DOF
        for i=1:300
            sendJointVelocityCommand(jc,jntVelCmd);
        end

        pause(5);

        jntVelCmd = [0;0;0;0;0;0;-0.2]; %7DOF
        for i=1:260
            sendJointVelocityCommand(jc,jntVelCmd);
        end

        pause(1);

        desired_pos=[x_origin; -y_origin; 0.08; home_pos(4); home_pos(5); home_pos(6)];
        moveToCP(jc,desired_pos);

        pause(1);

        setPositionControlMode(jc);
        fCmd = 0*ones(3,1);
        sendFingerPositionCommand(jc,fCmd);

        desired_pos=[x_origin; -y_origin; 0.2; home_pos(4); home_pos(5); home_pos(6)];
        moveToCP(jc,desired_pos);

        pause(1);

        goToHomePosition(jc);

        break;
    end

    disp('Decoding Error!');
    disp('Receiving brain signal again....');

    pause(2);

    desired_pos=[x_origin; -y_origin; 0.07; home_pos(4); home_pos(5); home_pos(6)];
    moveToCP(jc,desired_pos);

    pause(1);

    setPositionControlMode(jc);
    fCmd = 0*ones(3,1);
    sendFingerPositionCommand(jc,fCmd);

    pause(1);

    desired_pos=[x_origin; -y_origin; 0.2; home_pos(4); home_pos(5); home_pos(6)];
    moveToCP(jc,desired_pos);

    goToHomePosition(jc);

    pause(1);

    break

    bbci_acquire_bv('close');
end
end

author: PatternRecognition

repo: OpenBMI


text:
%% QUAD_POOL uses the MATLABPOOL command to run the QUAD code.
%
%  Discussion:
%
%    Output printed by the function appears directly on the screen.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
  matlabpool open local 4

  n = 100000;
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_POOL\n' );
  fprintf ( 1, '  Estimate integral using %d points\n', n );

  value = quad_fun ( n );

  matlabpool close

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integral estimate is %f\n', value );
  fprintf ( 1, '  Exact integral    is %f\n', pi );
  fprintf ( 1, '  Error is             %e\n', abs ( value - pi ) );


text:
function test_suite = test_regression_sparse2

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_REGRESSION_SPARSE2
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

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

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName(), { ...
    {'Eft_full', @(x)x(1:100)}, {'Varft_full', @(x)x(1:100)}, ...
    {'Eft_var', @(x)x(1:100)}, {'Varft_var', @(x)x(1:100)} })
end

function testPredictionsFull(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_full', 'Varft_full'}, @mean, ...
    'RelTolElement', 0.1)
end

function testPredictionsVar(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_var', 'Varft_var'}, ...
    'RelTolElement', 0.1, 'RelTolRange', 0.02)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function weight = subset_weight ( n, t )

%*****************************************************************************80
%
%% SUBSET_WEIGHT computes the Hamming weight of a set.
%
%  Discussion:
%
%    The Hamming weight is simply the number of elements in the set.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the order of the master set, of which T
%    is a subset.  N must be positive.
%
%    Input, integer T(N), defines the subset T.
%    T(I) is 1 if I is an element of T, and 0 otherwise.
%
%    Output, integer WEIGHT, the Hamming weight of the subset T.
%

%
%  Check.
%
  subset_check ( n, t );

  weight = sum ( t(1:n) );

  return
end


text:
function fh = isolate_axes(ah, vis)
%ISOLATE_AXES Isolate the specified axes in a figure on their own
%
% Examples:
%   fh = isolate_axes(ah)
%   fh = isolate_axes(ah, vis)
%
% This function will create a new figure containing the axes/uipanels
% specified, and also their associated legends and colorbars. The objects
% specified must all be in the same figure, but they will generally only be
% a subset of the objects in the figure.
%
% IN:
%    ah - An array of axes and uipanel handles, which must come from the
%         same figure.
%    vis - A boolean indicating whether the new figure should be visible.
%          Default: false.
%
% OUT:
%    fh - The handle of the created figure.

% Copyright (C) Oliver Woodford 2011-2014, Yair Altman 2015-
%{
% Thank you to Rosella Blatt for reporting a bug to do with axes in GUIs
% 16/03/12: Moved copyfig to its own function. Thanks to Bob Fratantonio
%           for pointing out that the function is also used in export_fig.m
% 12/12/12: Add support for isolating uipanels. Thanks to michael for suggesting it
% 08/10/13: Bug fix to allchildren suggested by Will Grant (many thanks!)
% 05/12/13: Bug fix to axes having different units. Thanks to Remington Reid for reporting
% 21/04/15: Bug fix for exporting uipanels with legend/colorbar on HG1 (reported by Alvaro
%           on FEX page as a comment on 24-Apr-2014); standardized indentation & help section
% 22/04/15: Bug fix: legends and colorbars were not exported when exporting axes handle in HG2
% 02/02/21: Fix axes, figure size to preserve input axes image resolution (thanks @Optecks)
% 25/10/21: Bug fix: subplots were not isolated properly leading to print error (issue #347)
% 24/03/23: Remove any non-legendable objects from the legend (workaround for copyobj bug)
%}

    % Make sure we have an array of handles
    if ~all(ishandle(ah))
        error('ah must be an array of handles');
    end

    % Check that the handles are all for axes or uipanels, and are all in the same figure
    fh = ancestor(ah(1), 'figure');
    nAx = numel(ah);
    for a = 1:nAx
        if ~ismember(get(ah(a), 'Type'), {'axes', 'uipanel'})
            error('All handles must be axes or uipanel handles.');
        end
        if ~isequal(ancestor(ah(a), 'figure'), fh)
            error('Axes must all come from the same figure.');
        end
    end

    % Tag the objects so we can find them in the copy
    old_tag = get(ah, 'Tag');
    if nAx == 1
        old_tag = {old_tag};
    end
    set(ah, 'Tag', 'ObjectToCopy');

    % Mark all non-legendable objects in the original figure (if any legend shown)
    % a workaround to the copyobj bug of not copying the ApplicationData property
    lh = findall(fh, 'Tag','legend');
    if ~isempty(lh)
        ax = arrayfun(@(h) get(h,'Axes'), lh);
        c = findall(ax);
        hb = ~arrayfun(@(h)hasbehavior(h,'legend'), c);
        c = c(hb);
        oldTags = get(c,'Tag');
        set(c,'Tag','non-legendable!')
    end

    % Create a new figure exactly the same as the old one
    fh = copyfig(fh); %copyobj(fh, 0);

    % Fix Axes & Figure size for image catpuring to have almost exact resolution
    % of the Input Axes (thanks @Optecks)
    allaxes = findall(fh, 'type', 'axes');
    if ~isempty(ah)
        sz = get(ah(1), 'OuterPosition');
        un = get(ah(1), 'Units');
        set(allaxes(1), 'Units',un, 'OuterPosition', [0 0 sz(3) sz(4)]);
        set(allaxes(1), 'Units','pixels');
        sz = get(allaxes(1), 'OuterPosition');
        set(fh, 'Units','pixels', 'Position',[0 0 sz(3) sz(4)]+1);
    end

    if nargin < 2 || ~vis
        set(fh, 'Visible', 'off');
    end

    % Restore the axes tags in the original figure back to their original values
    for a = 1:nAx
        set(ah(a), 'Tag', old_tag{a});
    end

    % Find the objects to save
    ah = findall(fh, 'Tag', 'ObjectToCopy');
    if numel(ah) ~= nAx
        close(fh);
        error('Incorrect number of objects found.');
    end

    % Set the axes tags in the new figure to what they should be
    for a = 1:nAx
        set(ah(a), 'Tag', old_tag{a});
    end

    % Keep any legends and colorbars which overlap the subplots
    % Note: in HG1 these are axes objects; in HG2 they are separate objects, therefore we
    %       don't test for the type, only the tag (hopefully nobody but Matlab uses them!)
    lh = findall(fh, 'Tag','legend', '-or', 'Tag','Colorbar');
    nLeg = numel(lh);
    if nLeg > 0
        % Remove any non-legendable objects from the legend (if previously set)
        c2 = findall(ah,'Tag','non-legendable!');
        if ~isempty(c2)
            arrayfun(@(h)hasbehavior(h,'legend',false), c2);
            arrayfun(@(h,t)set(h,'Tag',t{1}),c,oldTags);  %restore object tags to orig values
        end

        set([ah(:); lh(:)], 'Units', 'normalized');
        try
            ax_pos = get(ah, 'OuterPosition'); % axes and figures have the OuterPosition property
        catch
            ax_pos = get(ah, 'Position'); % uipanels only have Position, not OuterPosition
        end
        if nAx > 1
            ax_pos = cell2mat(ax_pos(:));
        end
        ax_pos(:,3:4) = ax_pos(:,3:4) + ax_pos(:,1:2);
        try
            leg_pos = get(lh, 'OuterPosition');
        catch
            leg_pos = get(lh, 'Position');  % No OuterPosition in HG2, only in HG1
        end
        if nLeg > 1
            leg_pos = cell2mat(leg_pos);
        end
        leg_pos(:,3:4) = leg_pos(:,3:4) + leg_pos(:,1:2);
        ax_pos = shiftdim(ax_pos, -1);
        % Overlap test
        M = bsxfun(@lt, leg_pos(:,1), ax_pos(:,:,3)) & ...
            bsxfun(@lt, leg_pos(:,2), ax_pos(:,:,4)) & ...
            bsxfun(@gt, leg_pos(:,3), ax_pos(:,:,1)) & ...
            bsxfun(@gt, leg_pos(:,4), ax_pos(:,:,2));
        ah = [ah; lh(any(M, 2))];
    end

    % Get all the objects in the figure
    axs = findall(fh);

    % Delete everything except for the input objects and associated items
    delete(axs(~ismember(axs, [ah; allchildren(ah); allancestors(ah)])));
end

function ah = allchildren(ah)
    ah = findall(ah);
    if iscell(ah)
        ah = cell2mat(ah);
    end
    ah = ah(:);
end

function ph = allancestors(ah)
    ph = [];
    for a = 1:numel(ah)
        h = get(ah(a), 'parent');
        while h ~= 0
            ph = [ph; h]; %#ok<AGROW>
            h = get(h, 'parent');
        end
    end
end

author: altmany


text:
function R = compute_residuals(Y,A,b,C,f)

% compute residual traces for each component

% INPUTS
% Y     matrix of data in 2D format of pointer to file
% A     set of spatial components
% b     set of spatial background components
% C     set of temporal components
% f     set of temporal background component

AA = A'*A;
AY = mm_fun(A,Y);
nA2 = sum(A.^2,1);
R = bsxfun(@times, AY - AA*C - (A'*b)*f,1./nA2(:));

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [U,T] = autoGen_cartPoleEnergy(x,q,dx,dq,m1,m2,g,l)
%AUTOGEN_CARTPOLEENERGY
%    [U,T] = AUTOGEN_CARTPOLEENERGY(X,Q,DX,DQ,M1,M2,G,L)

%    This function was generated by the Symbolic Math Toolbox version 6.3.
%    30-Nov-2015 20:13:07

t2 = cos(q);
U = -g.*l.*m2.*t2;
if nargout > 1
    t3 = dx+dq.*l.*t2;
    t4 = sin(q);
    T = dx.^2.*m1.*(1.0./2.0)+m2.*(t3.^2+dq.^2.*l.^2.*t4.^2).*(1.0./2.0);
end

author: MatthewPeterKelly

repo: dscTutorials


function [price, opt, L] = PROJ_GMDB_DCA_Fast(proj_params, S_0, gmdb_params, r, q, modelInput)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for DCA-Style Garuanteed Minimum Withdraw Benefit (GMWB) using PROJ method
%        This version is based on a dollar cost average style investment account (see reference below)
%
% Terminal Payoff:  Payoff(tau) = L*exp(g*tau) + (Gam(tau) - L*exp(g*tau))^+
%                      Gam(tau) = S_M * sum_{m=0}^M(alpha*gamma / S_m)
%                          tau  = time of death (discrete periods)
%                            M  = number of periods until time of death (each period length dt)
%
% Models Supported: Levy Processes, including jump diffusions and Black-Scholes model
% Returns: price of contract
%
% Author: Justin Lars Kirkby
% References: 1) Equity-Linked  Guaranteed Minimum Death Benefits with Dollar Cost Averaging, J.L.Kirkby & D.Nguyen, 2021
%
% ----------------------
% Contract/Model Params
% ----------------------
% S_0 = initial stock price (e.g. 100)
% r   = interest rate (e.g. 0.05)
% q   = dividend yield (e.g. 0.05)
% M   = number of subintervals of [0,T] (total of M+1 monitoring points in time grid, including S_0)
% gmdb_params = container of GMDB contract params, see below
% modelInput =  model inputs, see below
%
% ----------------------
% Numerical (PROJ) Params
% ----------------------
% proj_params = numerical params
%   proj_params.N = number of basis elements, e.g. N = 2^10
%   proj_params.L1 = gridwidth param, e.g. L1 = 8
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ------------------
% GMDB Contract Params
% ------------------
L = gmdb_params.L;   % Guarantee Level: Set L = -1 to use ATMF value for L
alpha = gmdb_params.alpha;  % Period premium payment, paid every dt time units
gamma = gmdb_params.gamma;  % Proportion of investment retained by policyholder (fee is 1-gamma)
contract_type = gmdb_params.contract_type;  % Contract type: 1 = GMDB, 2 = GMDB-RS (Ratchet strike)
p = gmdb_params.death_prob;  % death probability distribution, must be consistent with dt
g = gmdb_params.g;

% ------------------
% Model Inputs
% ------------------
dt = modelInput.dt;  % Time increment, premiums paid / underlying S is monitored every dt
phiR = modelInput.rnCHF;  % Risk neutral CHF for time period dt

call = 1;
ER = 0;

Z = gen_func(-r, dt, p);
if g == 0
    Zrg = Z;
else
    Zrg = gen_func(-(r-g), dt, p);
end

if L == -1
    MF = gen_func(r - q - g, dt, p);
    Zg = gen_func(-g, dt, p);
    L = alpha * gamma * (exp((r-q)*dt)*MF - Zg) / (exp((r-q)*dt) - 1);
end

Mmax = length(p); Tmax = Mmax*dt;
pr_alpha = getTruncationAlpha(Tmax, proj_params.L1, modelInput, modelInput.model);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = proj_params.N;
dx = 2*pr_alpha/(N-1); a = 1/dx;
A    = 32*a^4;
C_aN = A/N;

%%% SHIFTS
x1    = zeros(1,Mmax);

if contract_type == 1
    strikes = S_0*L*exp(g*dt*(1:Mmax)) ./ (alpha*gamma*(2:Mmax + 1));
else
    strikes = S_0*ones(1,Mmax);
end

for m=1:Mmax
    if m == 1
        x1(m) = ER;
    else
        x1(m) = ER + log(1+exp(x1(m-1)));  %%BENHAMOU SHIFT
    end
    %x1(m) = log(m) + .5*(m+1)*ER;    %% LOWER BOUND SHIFT derived in APROJ paper
end

Nm   = floor(a*(x1-ER));
x1   = ER + (1-N/2)*dx + Nm*dx;
NNM  = N + Nm(Mmax-1);   %Number of columns of PSI

%%% Now check that we wont fall off the grid later
for m=1:Mmax
    ystar = log((m+1)*strikes(m)/S_0 -1);
    nbar  = floor((ystar-x1(m))*a+1);
    % fprintf('%.0f\n', nbar);
    if nbar + 1 > N
        proj_params.L1 = proj_params.L1*1.25;
        [price, opt, L] = PROJ_GMDB_DCA_Fast(proj_params, S_0, gmdb_params, r, q, modelInput);
        return;
    end
end

dxi   = 2*pi*a/N;
xi    = dxi*(1:(N-1))';
PhiR = [1; phiR(xi)];

b0    = 1208/2520; b1 = 1191/2520; b2 = 120/2520; b3 = 1/2520;
zeta  = (sin(xi/(2*a))./xi).^4./(b0 + b1*cos(xi/a) +b2*cos(2*xi/a) +b3*cos(3*xi/a));
AA = 1/A;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PSI Matrix: 5-Point GAUSSIAN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI = make_PSI(N,NNM,x1(1),dx,dxi);

%%%%%%%%%%%%%%
% STEP 1) Value the European!!!
%%%%%%%%%%%%%%
beta  = [AA; zeta.*PhiR(2:N).*exp(-1i*x1(1)*xi)];   %grand(end)=.5*grand(end);
beta  = real(fft(beta));

s = european_price(zeta,PhiR,xi,dx, r, q, dt, strikes(1), S_0, ER, N, a, call);
s = p(1) * 2 * s;

%%%%%%%%%%%%%%
% STEP 2) Value the rest
%%%%%%%%%%%%%%

PhiR  = C_aN*PhiR;
beta  = PSI(:,1:N)*beta.*PhiR;  %Nm(1)=0

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

%%%%% Loop to find PSI_M
for n=2:Mmax
    opt_v = intermediate_asian_price(N, dx, dt, xi, zeta, beta, x1(n), n, r, q, strikes(n), S_0);
    % fprintf('%.12f \n',opt_v);

    beta(2:N) = zeta.*beta(2:N).*exp(-1i*x1(n)*xi); beta(1) = AA;
    beta      = real(fft(beta));
    if n < Mmax
        beta      = PSI(:,Nm(n)+1:Nm(n)+N)*beta.*PhiR;
    end

    s = s + p(n) * (n + 1) * opt_v;
end

opt = s * alpha * gamma / S_0;
price = L*Zrg - alpha * (exp(r*dt) - Z) / (exp(r*dt) - 1) + opt;

end

function PSI = make_PSI(N,NNM,x_1,dx,dxi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PSI Matrix: 5-Point GAUSSIAN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI     = zeros(N,NNM);    %The first row will remain ones
PSI(1,:) = ones(1,NNM);

%%%% Sample
Neta  = 5*(NNM) + 15;   %sample size
Neta5 = (NNM) + 3;
g2    = sqrt(5-2*sqrt(10/7))/6;
g3    = sqrt(5+2*sqrt(10/7))/6;
v1    = .5*128/225;
v2    = .5*(322+13*sqrt(70))/900;
v3    = .5*(322 - 13*sqrt(70))/900;

thet                 = zeros(1,Neta);   %sample initialized
thet(5*(1:Neta5)-2)  = x_1 -1.5*dx + dx*(0:Neta5-1);
thet(5*(1:Neta5)-4)  = x_1 -1.5*dx + dx*(0:Neta5-1) - dx*g3;
thet(5*(1:Neta5)-3)  = x_1 -1.5*dx + dx*(0:Neta5-1) - dx*g2;
thet(5*(1:Neta5)-1)  = x_1 -1.5*dx + dx*(0:Neta5-1) + dx*g2;
thet(5*(1:Neta5))    = x_1 -1.5*dx + dx*(0:Neta5-1) + dx*g3;

%%%% Weights
sig      = [-1.5-g3, -1.5-g2, -1.5, -1.5+g2, -1.5+g3, -.5-g3, -.5-g2, -.5, -.5+g2, -.5+g3,];
sig(1:5) = (sig(1:5) + 2).^3/6;
sig(6:10) = 2/3 - .5*(sig(6:10)).^3 - (sig(6:10)).^2;

sig([1 5 6 10]) = v3*sig([1 5 6 10]);
sig([2 4 7 9]) = v2*sig([2 4 7 9]);
sig([3 8]) = v1*sig([3 8]);

%%%% Fill Matrix
zz  = exp(1i*dxi*log(1+exp(thet)));
thet   = zz;

for j=2:N-1
    PSI(j,:) =  sig(1)*(thet(1:5:Neta-19) + thet(20:5:Neta)) ...
              + sig(2)*(thet(2:5:Neta-18) + thet(19:5:Neta-1)) ...
              + sig(3)*(thet(3:5:Neta-17)  + thet(18:5:Neta-2)) ...
              + sig(4)*(thet(4:5:Neta-16)  + thet( 17:5:Neta-3)) ...
              + sig(5)*(thet(5:5:Neta-15)  + thet( 16:5:Neta-4)) ...
              + sig(6)*(thet(6:5:Neta-14)  + thet( 15:5:Neta-5)) ...
              + sig(7)*(thet(7:5:Neta-13)  + thet( 14:5:Neta-6)) ...
              + sig(8)*(thet(8:5:Neta-12)  + thet( 13:5:Neta-7)) ...
              + sig(9)*(thet(9:5:Neta-11)  + thet( 12:5:Neta-8)) ...
              + sig(10)*(thet(10:5:Neta-10)  + thet( 11:5:Neta-9));

    thet = thet.*zz;
end
end

function Val = intermediate_asian_price(N, dx, dt, xi, zeta, chf, x_1, M, r, q, W, S_0)
call = 1;
a = 1/dx;
A    = 32*a^4;
AA = 1/A;
C_aN = A/N;
T = M*dt;

%%%%% FINAL VALUE
ystar = log((M+1)*W/S_0 -1);
nbar  = floor((ystar-x_1)*a+1);
C     = S_0/(M+1);
D     = W - C;
x_1   = ystar- (nbar-1)*dx;

beta(2:N) = zeta.*chf(2:N).*exp(-1i*x_1*xi); beta(1)=AA;
beta      = real(fft(beta)).';

Cc1 = C*( exp(-7/4*dx)/54 + exp(-1.5*dx)/18 + exp(-1.25*dx)/2 + 7*exp(-dx)/27 )/20;

Cc2 = C*.05*(28/27 + exp(-7/4*dx)/54 + exp(-1.5*dx)/18 + exp(-1.25*dx)/2 + 14*exp(-dx)/27 ...
              + 121/54*exp(-.75*dx) + 23/18*exp(-.5*dx) + 235/54*exp(-.25*dx));

Cc3 = C*( (28 + 7*exp(-dx))/3 ...
              + ( 14*exp(dx) + exp(-7/4*dx) + 242*cosh(.75*dx) + 470*cosh(.25*dx))/12 ...
              +.25*(exp(-1.5*dx) + 9*exp(-1.25*dx) + 46*cosh(.5*dx)))/90;

Cc4 = C*( 14/3*(2+cosh(dx)) ...
              + .5*(cosh(1.5*dx) + 9*cosh(1.25*dx) +23*cosh(.5*dx))...
              +  1/6*(cosh(7/4*dx) + 121*cosh(.75*dx) +235*cosh(.25*dx)))/90;

G           = zeros(nbar+1,1);
E           = exp(ystar-(nbar-1)*dx+dx*(0:nbar));

G(nbar+1)   = D/24    - Cc1*E(nbar+1);
G(nbar)     = .5*D    - Cc2*E(nbar);
G(nbar-1)   = 23*D/24 - Cc3*E(nbar-1);
G(1:nbar-2) = D       - Cc4*E(1:nbar-2);

Val = C_aN*exp(-r*T)*sum(beta(1:nbar+1).*G);
if call==1  %Call Option, use Put-Call-Parity
    if r - q == 0
        mult = M + 1;
    else
        mult = (exp((r-q)*T*(1+1/M))-1)/(exp((r-q)*dt)-1);
    end
    Val = Val + C*exp(-r*T)*mult - W*exp(-r*T);
end
Val = max(0, Val);

end

function price = european_price(zeta,PhiR,xi,dx, r, q, T, W, S_0, c1, N, a, call)

W = 2*W - S_0;  % we shift the strike then divide by 2 at the end, this is asian with M=2

lws = log(W/S_0);
lam = c1 -(N/2 -1)*dx;
nbar = floor(a*(lws-lam)+1);
if nbar>=N
    nbar = N-1;
end
xmin = lws - (nbar-1)*dx;

Cons = 32*a^4;

beta  = [1/Cons; zeta.*PhiR(2:N).*exp(-1i*xmin*xi)];   %grand(end)=.5*grand(end);
beta  = real(fft(beta)).';

G = zeros(1,nbar +1);
G(nbar +1) = W*(1/24 - 1/20*exp(dx)*(exp(-7/4*dx)/54 + exp(-1.5*dx)/18 + exp(-1.25*dx)/2 + 7*exp(-dx)/27));

G(nbar )   =  W*(.5 -.05*(28/27 + exp(-7/4*dx)/54 + exp(-1.5*dx)/18 + exp(-1.25*dx)/2 + 14*exp(-dx)/27 ...
            + 121/54*exp(-.75*dx) + 23/18*exp(-.5*dx) + 235/54*exp(-.25*dx)));

G(nbar -1) = W*( 23/24 - exp(-dx)/90*( (28 + 7*exp(-dx))/3 ...
            + ( 14*exp(dx) + exp(-7/4*dx) + 242*cosh(.75*dx) + 470*cosh(.25*dx))/12 ...
            +.25*(exp(-1.5*dx) + 9*exp(-1.25*dx) + 46*cosh(.5*dx))) );

G(1: nbar -2) = W - S_0*exp(xmin +dx*(0:nbar-3))/90*( 14/3*(2+cosh(dx)) ...
                + .5*(cosh(1.5*dx) + 9*cosh(1.25*dx) +23*cosh(.5*dx))...
                +  1/6*(cosh(7/4*dx) + 121*cosh(.75*dx) +235*cosh(.25*dx)));

if call == 1  % Use put-call parity
    price = Cons*exp(-r*T)/N*G*(beta(1:length(G))') + S_0*exp(-q*T) - W*exp(-r*T);
else
    price = Cons*exp(-r*T)/N*G*(beta(1:length(G))');
end

price = 0.5*max(price, 0);  % Protect against deep out of money case

end


text:
function  varargout = block_interface(varargin)
%BLOCK_INTERFACE Common block processing backend
%
%  Object-like interface for sharing data between block handling
%  functions.

if nargin<1
     error('%s: Too few input parameters.',upper(mfilename));
end

% Persistent data
persistent pos;
persistent datapos; % Second data pointer
persistent pageNo;
persistent sourceName;
persistent maxBufCount;
persistent pageList;
persistent playChanList;
persistent recChanList;
persistent skipCounter;

persistent classid;
persistent toPlayBlock;
persistent anaOverlap;
persistent synOverlap;
persistent dispLoad;

persistent Ls;
persistent Fs;

persistent loop;

% DEFAULTS
persistent bufLen;

persistent outFile;

persistent offline;

command = varargin{1};

switch command
   case {'reset','clearAll'}
      pos = 0;
      datapos = 0;
      sourceName = [];
      maxBufCount = 3;
      pageList = [];
      playChanList = [];
      recChanList = [];
      pageNo = 0;
      skipCounter = 0;
      classid = 'double';
      anaOverlap = [];
      synOverlap = [];
      toPlayBlock = [];
      dispLoad = 1;
      loop = 0;
      Ls = -1;
      bufLen = 1024;
      outFile = [];
      offline = 0;
      Fs = 0;
%% SETTERS %%%
   case 'setLs'
      Ls = varargin{2};
   case 'setFs'
      Fs = varargin{2};
   case 'setOutFile'
      outFile = varargin{2};
   case 'setOffline'
      offline = varargin{2};
   case 'setPos'
      pos = varargin{2};
   case 'setDatapos'
      datapos = varargin{2};
   case 'setBufCount'
      maxBufCount = varargin{2};
   case 'setPlayChanList'
      playChanList = varargin{2};
   case 'setRecChanList'
      recChanList = varargin{2};
   case 'setPageNo'
      pos = varargin{2};
   case 'setSkipped'
      skipCounter = varargin{2};
   case 'setBufLen'
      bufLen = varargin{2};
   case 'setClassId'
      classid = varargin{2};
   case 'setAnaOverlap'
      anaOverlap = varargin{2};
   case 'setSynOverlap'
      synOverlap = varargin{2};
   case 'setDispLoad'
      dispLoad = varargin{2};
   case 'setToPlay'
      toPlayBlock=varargin{2};
   case 'setIsLoop'
      loop = varargin{2};
   case 'setSource'
      sourceName = varargin{2};
%% GETTERS %%%
   case 'getLs'
      varargout{1}=Ls;
   case 'getFs'
      varargout{1}=Fs;
   case 'getOutFile'
      varargout{1}=outFile;
   case 'getOffline'
      varargout{1}=offline;
   case 'getPos'
      varargout{1}=pos;
   case 'getDatapos'
      varargout{1}=datapos;
   case 'getBufCount'
      varargout{1}= maxBufCount;
   case 'getPlayChanList'
      varargout{1}=playChanList;
   case 'getRecChanList'
      varargout{1}=recChanList;
   case 'getPageList'
      varargout{1}=pageList;
   case 'getPageNo'
      varargout{1}=pageNo;
   case 'getSkipped'
      varargout{1}=skipCounter;
   case 'getBufLen'
      varargout{1}=bufLen;
   case 'getClassId'
      varargout{1}=classid;
   case 'getAnaOverlap'
      varargout{1}=anaOverlap;
   case 'getSynOverlap'
      varargout{1}=synOverlap;
   case 'getDispLoad'
      varargout{1}=dispLoad;
   case 'getToPlay'
      varargout{1}=toPlayBlock;
      toPlayBlock = [];
   case 'getIsLoop'
      varargout{1}=loop;
   case 'getSource'
      if isnumeric(sourceName)
         varargout{1}='numeric';
      else
         varargout{1}=sourceName;
      end
   case 'getEnqBufCount'
      varargout{1}= numel(pageList);

%% OTHER %%%
   case 'incPageNo'
      pageNo = pageNo +1;
   case 'flushBuffers'
      anaOverlap = [];
      synOverlap = [];
   case 'popPage'
      varargout{1}=pageList(1);
      pageList = pageList(2:end);
   case 'pushPage'
      pageList = [pageList, varargin{2}];

%   case 'incSkipped'
%      skipCounter = skipCounter + varargin{2};
%   case 'readNumericBlock'
%      L = varargin{2};
%      varargout{1}=sourceName(pos+1:pos+1+L,:);
   otherwise
      error('%s: Unrecognized command.',upper(mfilename));
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function plotFeatureImportanceRF(pathRF)

startpath = pwd;

cd(pathRF), load('testingVariableImportance')
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
% fSetNames = fieldnames(training.textures.(nameOutcomes{1})); nFset = numel(fSetNames);
fSetNames = {'PET','CT','PETCT'}; nFset = numel(fSetNames);

for o = 1:nOutcomes
    for f = 1:nFset
        percentAUCdecrease = variableImportance.(nameOutcomes{o}).(fSetNames{f}).percentAUCdecrease;
        varNames = variableImportance.(nameOutcomes{o}).(fSetNames{f}).varNames; nVar = numel(varNames);
        figure
        barh(1:nVar,percentAUCdecrease*100)
        for i = 1:nVar
            ind = strfind(varNames{i},'_'); nInd = numel(ind);
            if ~isempty(ind)
                for n = 1:nInd
                    varNames{i} = [varNames{i}(1:(ind(n)-1+(n-1))),'\',varNames{i}((ind(n)+(n-1)):end)];
                end
            end
        end
        set(gca,'yticklabel',varNames)
        title(['RANDOM PERMUTATIONS (RF):',nameOutcomes{o},' -- ',fSetNames{f}])
        xlabel('Percent AUC decrease')
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function[]=makefigs_wavetrans
%MAKEFIGS_WAVETRANS  Makes a sample figure for WAVETRANS.

load bravo94
%not using 'use' as it's not compatible in scripts with some systems
num=bravo94.rcm.num;
cv=bravo94.rcm.cv;

vindex(cv,2,2);

gamma=3;beta=2;
fs=morsespace(gamma,beta,{0.05,pi},pi/1000,4);
[wp,wn]=wavetrans(cv,conj(cv),{gamma,beta,fs,'bandpass'});
h=wavespecplot(yearfrac(num),vfilt(cv,24),1./fs,sqrt(squared(wp)+squared(wn)));
colormap lansey
axes(h),ylim([-35 55]),hlines(0,'k:'),caxis([0.25 21])

%To print
if 0
    currentdir=pwd;
    cd([whichdir('jlab_license') '/figures'])
    print -dpng wavetrans
    crop wavetrans.png
    cd(currentdir)
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function x=demean(x,dim)

% DEMEAN(X)
% Removes the Average or mean value.
%
% DEMEAN(X,DIM)
% Removes the mean along the dimension DIM of X.

if(nargin==1),
   dim = 1;
   if(size(x,1) > 1)
      dim = 1;
   elseif(size(x,2) > 1)
      dim = 2;
   end;
end;

dims = size(x);
dimsize = size(x,dim);
dimrep = ones(1,length(dims));
dimrep(dim) = dimsize;

x = x - repmat(mean(x,dim),dimrep);

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
% given an ND tensor of size DxN1xN2xN3.., compute the outer product of the
% first dimension independently to return DxDxN1xN2xN3...
%
function output = outProdND(data)

[D,N1,N2,N3] = size(data);
A = reshape(data,D, N1*N2*N3);

% B = permute(bsxfun(@times, A, conj(permute(A,[3 2 1]))), [1 3 2]); % slower
B = bsxfun(@times, permute(A, [1 3 2]), permute(conj(A), [3 1 2]));

output = reshape(B, D,D,N1,N2,N3);

end

author: snsun


text:
function cvx_clearspath

%CVX_CLEARSPATH   Clears the cvx solver path.
%   CVX_CLEARSPATH removes the internal solver directories to the Matlab
%   path. CVX automatically clears the solver path after completion of a
%   model computation, so calling this function should not be necessary.
%   Nevertheless, we provide it for completeness and debugging.

global cvx___
osolv = cvx___.solvers.active;
if osolv,
    tstr = cvx___.solvers.list(osolv).path;
    if ~isempty( tstr ),
        path(strrep(path,tstr,''));
    end
    cvx___.solvers.active = 0;
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:
%  Insert text from a file
%
%  ::
%
%     document.include('option_name',option_value);
%
%  Args:
%
%     varargin: arguments need to come in pairs
%
%        - **filename** : name of the file to be included
%
%

author: jmaih