text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REAL-TIME BCILAB/SIFT INTERFACE DEMO SCRIPT
%
% This script demonstrates a beta interface between BCILAB
% (Christian Kothe: www.sccn.ucsd.edu/wiki/BCILAB) and SIFT (Tim Mullen:
% www.sccn.ucsd.edu/wiki/SIFT) to perform real-time data extraction,
% cleaning, connectivity modeling and visualization.
%
% This script was prepared for the following whitepaper/demonstration:
% [1] Tim Mullen, Christian Kothe, Yu Chi, Tzyy-Ping Jung,
%     and Scott Makeig. "Real-Time Estimation and 3D Visualization of Sparse
%     Multivariate Effective Connectivity Using Wearable EEG." July 15th,
%     IEEE EMB/CAS/SMC Workshop on Brain-Machine-Body Interfaces. 34th
%     Annual International IEEE EMBS Conference of the IEEE Engineering in
%     Medicine and Biology Society (IEEE EMBS).
%
% Please refer to [1] when using any functionality from this demo.
%
% Copyright Tim Mullen, Christian Kothe, July 2012, SCCN/INC/UCSD
%
% Author: Tim Mullen,       July 2012, SCCN/INC/UCSD
%         Christian Kothe,  July 2012, SCCN/INC/UCSD
%
% 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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% SET UP CONFIGURATION OPTIONS
% NOTE: all paths are relative to 'datapath' which is a
% platform-independent path which itself can be relative to bcilab root
% folder (i.e. data:/ is the userdata folder in the bcilab root dir)

RUN_LSL          = false;   % if 'online' streaming, RUN_LSL = true, otherwise playback TestingDataFile (below)
COLOR_SOURCE_ROI = false;   % ONLY if using source reconstruction, this will color sources according to ROI
HEAD_MODEL_NAME  = '';      % path to head model object for source reconstruction (relative to 'datapath'). Leave empty if you aren't doing source reconstruction
GUI_CONFIG_NAME  = 'DARPA_DEMO_METACP_CFG_FEWCHANS.mat';    % relative path to a default pipeline configuration
GUI_BRAINMOVIE_CONFIG_NAME = 'DARPA_DEMO_BM_CFG.mat';       % relative path to BrainMovie configuration

% establish where we will read/write prefs, etc from
datapath         = 'data:/';  % this is relative to the BCILAB root dir
TrainingDataFile = 'calibration_mindo.xdf';     % this is the relative path to the
TestingDataFile  = 'calibration_mindo.xdf';

% Set up the name of the stream we will write to in the workspace
streamname              = 'EEGDATA';
LSL_SelectionProperty   = 'type';  % can also be 'name'
LSL_SelectionValue      = 'EEG';

%% load head model object
if ~isempty(HEAD_MODEL_NAME)
    hmObj = headModel.loadFromFile(env_translatepath([datapath HEAD_MODEL_NAME]));
    % load meshes for visualization
    surfData    = load(hmObj.surfaces);
    surfData    = surfData.surfData;
    scalpMesh   = surfData(1);
    csfMesh     = surfData(2);
    cortexMesh  = surfData(3);
end

%% load calibration recording
calibData = exp_eval(io_loadset([datapath TrainingDataFile],'markerchannel',{'remove_eventchns',false}));
flt_pipeline('update');

%% start LSL streaming
if RUN_LSL == true
    run_readlsl('MatlabStream',streamname,'SelectionProperty',LSL_SelectionProperty,'SelectionValue',LSL_SelectionValue);
end
%% ... OR read from datafile
if RUN_LSL == false
    run_readdataset('MatlabStream',streamname,'Dataset',io_loadset([datapath TestingDataFile],'markerchannel',{'remove_eventchns',false}));
end

%% initialize some variables
[benchmarking.preproc     ...
 benchmarking.modeling    ...
 benchmarking.brainmovie] = deal(NaN);

newPipeline     = true;
newBrainmovie   = true;
specOpts        = [];

% create a minimal moving average pipeline...
% noflats = exp_eval_optimized(flt_movavg(flt_selchans(calibData,{'FC3', 'FC4', 'FC2', 'FC1', 'C3', 'C4', 'C2', 'C1', 'CP3', 'CP4', 'CP2', 'CP1', 'P3', 'P4', 'P2', 'P1'})));
noflats = flt_seltypes(calibData,'EEG');
%noflats = exp_eval_optimized(flt_movavg(flt_seltypes(calibData,'EEG')));
highpassed   = exp_eval_optimized(flt_fir(noflats,[0.5 1],'highpass'));
goodchannels = exp_eval_optimized(flt_clean_channels('signal',highpassed,'min_corr',0.35,'ignored_quantile',0.3));

%% initialize options
if ~isempty(GUI_CONFIG_NAME) && exist(env_translatepath([datapath GUI_CONFIG_NAME]),'file')
    % try to load a configuration file
    io_load([datapath GUI_CONFIG_NAME]);
else
    % manually set some defaults
    opts.miscOptCfg.dispBenchmark  = true;
    opts.miscOptCfg.doSIFT         = [];
    opts.holdPipeline              = false;
    opts.exitPipeline              = false;

    opts.siftPipCfg.preproc.verb = 0;

    % set some defaults
    opts.fltPipCfg = ...
        struct('DataCleaning', ...
                struct('RetainPhases',true, ...
                       'HaveBursts', true, ...
                       'HaveChannelDropouts', false, ...
                       'DataSetting', 'drycap'), ...
                'Rereferencing', {{}});

    opts.siftPipCfg = ...
        struct('preproc', ...
                    struct('normalize', ...
                        struct('verb', 0, ...
                               'method', {{'time'}})));

end

% additional options regarding brainmovie node/edge adaptation
opts.adaptation.adaptationHL    = 10;  % HL of exp. win. MA in frames
opts.adaptation.updateInterval  = 1;
opts.adaptation.bufferTime      = 1;

%% initialize the Meta Control Panel
if ~isempty(HEAD_MODEL_NAME)
    calibData.srcpot = 1; % allow sources
end
figHandles.MetaControlPanel = gui_metaControlPanel(calibData,opts);
waitfor(figHandles.MetaControlPanel,'UserData','initialized');

%% Initialize figures
figHandles.specDisplay      = [];
figHandles.BMDisplay        = [];
figHandles.BMControlPanel   = [];

%% try to load a BrainMovie configuration
try
    io_load([datapath GUI_BRAINMOVIE_CONFIG_NAME]);
    BMCFG.BMopts.bmopts_suppl = {'title',{'Multivariate ','Granger Causality  '}};
    BMCFG.BMopts.caption = true;
catch
    BMCFG = struct([]);
end

BMRenderMode = 'init_and_render';

%% Main loop
% -------------------------------------------------------------------------
while ~opts.exitPipeline

    try

        if newPipeline
            % create a new pipeline on training data
            fprintf('Pipeline changed\n');
            cleaned_data = exp_eval(flt_pipeline('signal',calibData,opts.fltPipCfg));
            pipeline     = onl_newpipeline(cleaned_data,streamname);
            newPipeline  = false;
        end

        if opts.holdPipeline
            % pause the pipeline
            pause(1);
            continue;
        end

        % grab a chunk of data from the stream and preprocess it
        % ---------------------------------------------------------------------
        prepbench = tic;
        [eeg_chunk,pipeline] = onl_filtered(pipeline, ...
                                            round(opts.miscOptCfg.winLenSec*cleaned_data.srate), ...
                                            opts.miscOptCfg.suppress_console_output, ...
                                            fastif(opts.fltPipCfg.psourceLocalize.arg_selection,{'srcpot'},{}));
        benchmarking.preproc = toc(prepbench);

        % model the chunk via sift
        % ---------------------------------------------------------------------
        if opts.miscOptCfg.doSIFT.arg_selection
            try
                % force the modeling window length to agree with chunk length
                opts.siftPipCfg.modeling.winlen = opts.miscOptCfg.winLenSec;

                % import ROI Names
                if strcmpi(opts.siftPipCfg.preproc.sigtype.arg_selection,'sources')
                    opts.siftPipCfg.preproc.varnames = eeg_chunk.roiLabels;
                end

                % select a subset of channels
                if ~isempty(opts.miscOptCfg.doSIFT.channelSubset)
                    [dummy eeg_chunk] = evalc('hlp_scope({''disable_expressions'',1},@flt_selchans,''signal'',eeg_chunk,''channels'',opts.miscOptCfg.doSIFT.channelSubset,''arg_direct'',true)');
                end

                % run the SIFT modeling pipeline
                siftbench = tic;
                [eeg_chunk cfg] = hlp_scope({'is_online',1},@onl_siftpipeline,'EEG',eeg_chunk,opts.siftPipCfg);
                benchmarking.modeling = toc(siftbench);

                if eeg_chunk.pnts == 0
                    fprintf('No data!\n');
                    pause(0.01);
                    continue;
                end

            catch err
                hlp_handleerror(err);
                pause(0.1);
                continue;
            end

            % visualize the brainmovie
            % -------------------------------------------------------------
            if opts.miscOptCfg.doSIFT.dispBrainMovie.arg_selection

                if newBrainmovie ...
                    | isempty(ishandle(figHandles.BMDisplay)) ...
                    | ~ishandle(figHandles.BMDisplay)

                    % initialize brainmovie
                    if ~isempty(BMCFG)
                        BMCFG.showNodeLabels.nodelabels = eeg_chunk.CAT.curComponentNames;
                        figHandles.BMControlPanel = gui_causalBrainMovie3D_online(eeg_chunk,eeg_chunk.CAT.Conn,struct('arg_direct',0),BMCFG,'timeRange',[]);
                    else
                        figHandles.BMControlPanel = gui_causalBrainMovie3D_online(eeg_chunk,eeg_chunk.CAT.Conn,struct('arg_direct',0),'showNodeLabels',{'nodelabels',eeg_chunk.CAT.curComponentNames});
                    end
                    waitfor(figHandles.BMControlPanel,'UserData','newBrainmovie');

                    % set the closing behavior
                    set(figHandles.BMControlPanel,'CloseRequestFcn','evalin(''base'',''figHandles.BMControlPanel = [];''); delete(gcbf);');

                    BM_CFG_CHANGED       = true;
                    newBrainmovie        = false;
                    figHandles.BMDisplay = [];
                    BMStateVars          = [];
                end

                if BM_CFG_CHANGED
                    BMRenderMode    = 'init_and_render';
                    BMStateVars     = [];

                    % determine whether to adapt limits
                    adapt_nodeSizeDataRange  = isempty(BMCFG.BMopts.graphColorAndScaling.nodeSizeDataRange)  && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;
                    adapt_edgeSizeDataRange  = isempty(BMCFG.BMopts.graphColorAndScaling.edgeSizeDataRange)  && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;
                    adapt_nodeColorDataRange = isempty(BMCFG.BMopts.graphColorAndScaling.nodeColorDataRange) && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;
                    adapt_edgeColorDataRange = isempty(BMCFG.BMopts.graphColorAndScaling.edgeColorDataRange) && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;

                    % clear state
                    nodeSizeState   = [];
                    nodeColorState  = [];
                    edgeSizeState   = [];
                    edgeColorState  = [];

                    % reset flag
                    BM_CFG_CHANGED = false;
                end

                % Handle special rendering of source meshes and colors
                % ---------------------------------------------------------
                if COLOR_SOURCE_ROI
                    BG_COLOR = [0.1 0.1 0.1];

                    if BMCFG.BMopts.Layers.custom.arg_selection
                        BMCFG.BMopts.Layers.custom.volumefile = eeg_chunk.dipfit.surfmesh;
                        BMCFG.BMopts.Layers.custom.meshcolor  = hlp_getROIVertexColorTable( ...
                                size(eeg_chunk.dipfit.surfmesh.vertices,1), ...
                                eeg_chunk.roiVertices,BG_COLOR,@(x)distinguishable_colors(x,BG_COLOR));
                    end

                    % if the cortex mesh is a custom mesh with 'constant'
                    % coloring...
                    if strcmpi(BMCFG.BMopts.Layers.cortex.cortexres,'custommesh') ...
                       && strcmp(BMCFG.BMopts.Layers.cortex.cortexcolor.arg_selection,'Constant')
                            % ... replace the colors with ROI coloring
                            BMCFG.BMopts.Layers.cortex.cortexcolor.colormapping = ...
                                hlp_getROIVertexColorTable( ...
                                    size(eeg_chunk.dipfit.surfmesh.vertices,1), ...
                                    eeg_chunk.roiVertices,BG_COLOR,@(x)distinguishable_colors(x,BG_COLOR));
%                         {hmObj.atlas.label, hlp_getROIColorTable(hmObj.atlas.label,eeg_chunk.roiLabels,[0.5 0.5 0.5],[1 0 0])};
                    end
                end

                % Render the Brain Movie
                % ---------------------------------------------------------
                bmbench = tic;
                [BMcfg_tmp BMhandles BMout] = ...
                    vis_causalBrainMovie3D(eeg_chunk,eeg_chunk.CAT.Conn,BMCFG, ...
                        'timeRange',[], ...
                        'BMopts',hlp_mergeVarargin(BMCFG.BMopts, ...
                                    'figurehandle',figHandles.BMDisplay, ...
                                    'RenderMode',BMRenderMode, ...
                                    'speedy',true, ...
                                    'InternalStateVariables',BMStateVars));
                benchmarking.brainmovie = toc(bmbench);

                if isempty(BMStateVars)
                    set(figHandles.BMDisplay, 'MenuBar','figure', ...
                                'CloseRequestFcn','evalin(''base'',''figHandles.BMDisplay = [];''); delete(gcbf);', ...
                                'ToolBar','none', 'Name','BrainMovie3D');
                end

                % save figure handle and state
                figHandles.BMDisplay = BMhandles.figurehandle;
                BMStateVars          = BMcfg_tmp.BMopts.vars;
                BMRenderMode         = 'render';

                % adapt the node and edge color/size limits
                % ---------------------------------------------------------
                if opts.miscOptCfg.doSIFT.dispBrainMovie.arg_selection ...
                   && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits

                    % adapt BrainMovie node/edge limits
                    if adapt_nodeSizeDataRange
                        [BMCFG.BMopts.graphColorAndScaling.nodeSizeDataRange, nodeSizeState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.NodeSize,     ...
                                                    'state',    nodeSizeState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end
                    if adapt_nodeColorDataRange
                        [BMCFG.BMopts.graphColorAndScaling.nodeColorDataRange, nodeColorState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.NodeColor,     ...
                                                    'state',    nodeColorState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end
                    if adapt_edgeSizeDataRange
                        [BMCFG.BMopts.graphColorAndScaling.edgeSizeDataRange, edgeSizeState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.EdgeSize,     ...
                                                    'state',    edgeSizeState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end
                    if adapt_edgeColorDataRange
                        [BMCFG.BMopts.graphColorAndScaling.edgeColorDataRange, edgeColorState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.EdgeColor,     ...
                                                    'state',    edgeColorState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end

                end % adaptBrainMovieLimits

            else % dispBrainMovie == false
                % close the brainmovie figure if control panel is closed
                if isempty(figHandles.BMControlPanel) ...
                        || ~ishandle(figHandles.BMControlPanel) ...
                        && ishandle(figHandles.BMDisplay)
                    close(figHandles.BMDisplay);
                end
                benchmarking.brainmovie = NaN;
            end

            % display the VAR spectrum
            % -------------------------------------------------------------
            if opts.miscOptCfg.doSIFT.dispSpectrum
                try
                    if ishandle(figHandles.specDisplay)
                        % update the SpecViewer
                        vis_autospectrum(specOpts,'FigureHandle',figHandles.specDisplay,'stream',eeg_chunk);
                    else
                        % generate a GUI and initialize SpecViewer
                        [figHandles.specDisplay specOpts] = arg_guidialog(@vis_autospectrum,'Parameters',{'stream',eeg_chunk});

                        set(figHandles.specDisplay,'CloseRequestFcn','evalin(''base'',''opts.miscOptCfg.dispSpectrum=false; figHandles.specDisplay = [];''); delete(gcbf);');
                        specOpts.arg_direct = 0;
                    end
                catch e
                    disp(e.message);
                    opts.miscOptCfg.doSIFT.dispSpectrum = false;
                end
            else
                if ishandle(figHandles.specDisplay)
                    close(figHandles.specDisplay);
                end
            end

        else
            benchmarking.modeling = NaN;
        end

        if opts.miscOptCfg.dispBenchmark
            % TODO: implement this as a figure
            vis_benchmark(benchmarking);
        end

    catch err
         hlp_handleerror(err);
    end

    pause(0.005);

end

%% all done, close all open figures
fnames = fieldnames(figHandles);
for f=1:length(fnames)
    if ishandle(figHandles.(fnames{f}))
        close(figHandles.(fnames{f}));
    end
end

repo: REST


text:
function tapas_linear_validate_data(y, u)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

if ~iscell(y)
    error('tapas:linear:validate_data', 'y should be a cell array');
end

if ~iscell(u)
    error('tapas:linear:validate_data', 'y should be a cell array');
end

if ~ (size(y, 2) == 1)
     error('tapas:linear:validate_data', 'y should be Nx1');
end

if ~ (size(u, 2) == 1)
     error('tapas:linear:validate_data', 'u should be Nx1');
end

if ~ (size(y, 1) == size(u, 1))
    error('tapas:linear:validate_data', 'y and u should have the same size');
end

np = size(y, 1);

[nt, d2] = size(y{1});

if d2 ~= 1
    error('tapas:linear:validate_data', 'y{1} should be Nx1');
end

[nu, nb] = size(u{1});

if nu ~= nt
    error('tapas:linear:validate_data', ...
        'u{1} first dimension should be %d', nt);
end

for i = 1:np
    if ~ (size(y{i}, 2) == 1)
        error('tapas:linear:validate_data', ...
            'y{%d} first dimension should be 1', i);
    end

    if ~ (size(u{i}, 1) == size(y{i}, 1))
         error('tapas:linear:validate_data', ...
            'y{%d} and u{%d} first dimension should be equal', i, i);
    end

    if ~(size(u{i}, 2) == nb)
        error('tapas:linear:validate_data',  ...
            'u{%d} second dimension is inconsistent', i);
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
function islands=bwislands(img)
%
% islands=bwislands(img)
%
% return the indices of non-zero elements in a 2D or 3D image
% grouped by connected regions in a cell array
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%	 img: a 2D or 3D array
% output:
%	 islands: a cell array, each cell records the indices
%		  of the non-zero elements in img for a connected
%		  region (or an island)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

img=logical(1-img);
idx=find(1-img(:));
islands={};

count=1;
out=cell(1,ndims(img));
while(~isempty(idx))
    [out{:}]=ind2sub(size(img),idx(1));
	imgnew=imfill(img,cell2mat(out));
	islands{count}=find(imgnew~=img);
	count=count+1;
	img=imgnew;
	idx=find(1-img(:));
end

author: fangq


text:
function nirfastmesh=readnirfast(filestub)
%
% nirfastmesh=readnirfast(v,f,filestub)
%
% load a group of mesh files saved in the NIRFAST format
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%      filestub: output file stub, output will include multiple files
%          filestub.node: node file
%          filestub.elem: element file to store the surface or tet mesh
%          filestub.param: parameter file
%          filestub.region: node label file
%          filestub.excoef: extinction coeff list
%
% output:
%      nirfastmesh.nodes: node list, 3 columns
%      nirfastmesh.elements: element list, 3 or 4 columns integers
%      nirfastmesh.bndvtx: boundary flag for each node, 1: on the boundary
%      nirfastmesh.region: node segmentation labels
%      nirfastmesh.dimension: dimension of the mesh
%      nirfastmesh.excoef: extinction coeff list
%      nirfastmesh.excoefheader: extinction coeff list field names
%      nirfastmesh.type: the header of the .param file
%      nirfastmesh.prop: optical property list (non-standard, need further processing)
%
%   format definition see http://www.dartmouth.edu/~nir/nirfast/tutorials/NIRFAST-Intro.pdf
%
% example:
%    [node,face,elem]=meshabox([0 0 0],[10 10 10],0.3,1);
%    savenirfast(node,elem,'test', [], ones(size(node)), 'user');
%    mymesh=readnirfast('test')
%    plotmesh([mymesh.nodes mymesh.bndvtx], mymesh.elements,'x>5')
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fname=[filestub,'.node'];
if(~exist(fullfile(pwd,fname),'file'))
    error([fname ' could not be found']);
end
nirfastmesh.nodes=load(fname);

nirfastmesh.bndvtx=nirfastmesh.nodes(:,1);
nirfastmesh.nodes(:,1)=[];

fname=[filestub,'.elem'];
if(~exist(fullfile(pwd,fname),'file'))
    error([fname ' could not be found']);
end
nirfastmesh.elements=load(fname);
nirfastmesh.dimension=size(nirfastmesh.elements,2)-1;

fname=[filestub,'.region'];
if(exist(fullfile(pwd,fname),'file'))
    nirfastmesh.region=load(fname);
end

fname=[filestub,'.excoef'];
fid=fopen(fname,'rt');
if(fid>=0)
    linenum=0;
    textheader={};
    while(~feof(fid))
        oneline=fgetl(fid);
        linenum=linenum+1;
        [data, count]=sscanf(oneline,'%f');
        if(count>1)
            params=fscanf(fid,repmat('%f ',1,count),inf);
            params=reshape(params,length(params)/count, count);
            params(2:end+1,:)=params;
            params(1,:)=data(:)';
            nirfastmesh.excoef=params;
            nirfastmesh.excoefheader=textheader;
            break;
        else
            textheader{end+1}=oneline;
        end
    end
    fclose(fid);
end

fname=[filestub,'.param'];
fid=fopen(fname,'rt');
if(fid>=0)
    linenum=0;
    params=[];
    while(~feof(fid))
        oneline=fgetl(fid);
        if(linenum==0)
            nirfastmesh.type=oneline;
        end
        linenum=linenum+1;
        [data, count]=sscanf(oneline,'%f');
        if(count>1)
            params=fscanf(fid,repmat('%f ',1,count),inf);
            params=reshape(params,length(params)/count, count);
            params(2:end+1,:)=params;
            params(1,:)=data(:)';
            nirfastmesh.prop=params;
            break;
        end
    end
    fclose(fid);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ChannelTag
   %ChannelTag Summary of this class goes here
   % ChannelTag Contains network, station, location and channel information
   %
   %
   % ChannelTag Properties:
   %   network  - network tag   (Generally 2 characters)
   %   station  - station tag   (Up to 5 alphanumeric characters)
   %   location - location tag  (2 characters)
   %   channel  - channel tag   (3 character channel tag)
   %
   % ChannelTag Methods:
   %   char - Retrieve properties as 'network.station.location.channel'
   %   string - Retrieve properties as 'network.station.location.channel'
   %   scn - Retrieve properties as 'station_channel_network'
   %   fixedlengthstrings -  return strings of a fixed length
   %   eq - == for channeltags, with simple '*' wildcard support
   %   ne - ~= for channeltag (no wilcard support)
   %
   %   sort - Sort in assending order by net.cha.loc.sta
   %   matching - Field-by-field comparison
   %
   % Examples of ChannelTag construction:
   %  chaTag = ChannelTag('NW.STA.LO.CHA'); % create from text
   %  chaTag = ChannelTag(network, station, location, channel);
   %  chaTag = ChannelTag(); %default, blank nscl
   %
   %  Example using cells of stationID strings
   %    x = ChannelTag.array({'NW.STA.LO.CHA','NW.STA.LO.CHA'[,...]});
   %
   %  Example using individual parts
   %    cha_tags = ChannelTag.array({'IU','ANMO','00',{'BHZ','BH1','BH2'});
   %
   %  note, an array of cha_tags may be created by including cell arrays of
   %  choices. However, all cells that contain multiple strings must be the
   %  same size and shape.
   %
   %  OK: x = ChannelTag.array('IU',{'ANMO','ANTO},'00',{'BHE','BH1'}
   %        results in IU.ANMO.00.BHE & IU.ANTO.00.BH1
   %
   %  NOTOK: x = ChannelTag.array('IU',{'ANMO','ANTO},'00',{'BHE';'BH1'}
   %        error because stations (1x2) not same size as channels (2x1)
   %

   %{
   %  ----------- WILDCARDS ---------------------
   %  ChannelTag is merely a storage unit for net-sta-chan-loc information. It
   %  is blind to wildcards. However, when used as an argument in WAVEFORM,
   %  wildcards in the ChannelTag take on meaning that depends somewhat on the
   %  DATASOURCE type. In most cases, * wildcards are understood without
   %  issue. Note that '*' differs from ''. The latter excludes this term from
   %  the search altogether. Station and channel cannot be excluded from the
   %  search.
   %
   %  ChannelTag('R*','BHZ','XE','') match all stations begining with R
   %  ChannelTag('MCK','*','XE','')  match all channels for MCK
   %  ChannelTag('R*','*','XE','')   all R* stations. All channels
   %
   %  The uses above have been tested against Antelope data sources. Note that
   %  in Antelope it would be more common to use the literal wildcards 'R.*'
   %  or '.*'. When waveform interprets ChannelTags, it considers these use
   %  the same.
   %
   % EXAMPLE:
   %   s={'ANMO','ANTO'} , n='IU' , L='00' , c={'BHZ,BHE'}
   %   chaTag=ChannelTag(s,c,n,L); %  returns a 1x2 ChannelTag:
   %    chaTag(1) contains ANMO, IU, 00, BHZ
   %    chaTag(2) contains ANTO, IU, 00, BHE
   %}

   properties
      network  = '';
      station  = '';
      location = '';
      channel  = '';
   end

   methods
      function obj = ChannelTag(varargin)
         %ChannelTag    Construct a ChannelTag from a string or 4 strings.
         %  chaTag = ChannelTag('IU', 'ANMO', '00', 'BHZ'); % net, sta, loc, cha
         %  chaTag = ChannelTag('IU.ANMO.00.BHZ');
         %
         %  any of these may be empty ([] or '').
         %  to create multiple ChannelTags at once, use ChannelTag.array()
         %
         %  See also: array
         % TODO: Update the details. ChannelTag.array is not necessary anymore

         switch nargin
            case 4
               % inputsThatAreCells = cellfun(@iscell,{N,S,L,C});
               % if any(cellfun(@iscell,{N,S,L,C}))
               if any(cellfun(@iscell,varargin))
                  error('ChannelTag:InvalidConversion',...
                     ['expected strings but received cells.\n', ...
                     'To create multiple ChannelTags, use ChannelTag.array()']);
               end
               varargin = strtrim(varargin);
               [obj.network, obj.station, obj.location, obj.channel] = deal(varargin{:});
            case 1
               inObj = varargin{1};
               switch class(inObj)
                  case 'ChannelTag'
                     obj = inObj;
                  case {'char','struct'}
                     if size(inObj,1) == 1
                     [obj.network, obj.station, obj.location, obj.channel] = ...
                        ChannelTag.parse(inObj);
                     else
                        obj = ChannelTag.array(varargin{1});
                        %error('ChannelTag:InvalidConversion',...
                        %   'To create multiple ChannelTags, use ChannelTag.array()');
                     end

                  case 'cell'
                     if numel(inObj) == 1
                        [obj.network, obj.station, obj.location, obj.channel] = ...
                           ChannelTag.parse(inObj{1});
                     else
                        obj = ChannelTag.array(varargin{1});
                        %error('ChannelTag:InvalidConversion',...
                        %   'To create multiple ChannelTags, use ChannelTag.array()');
                     end
                  otherwise
                     error ('ChannelTag:InvalidConversion','Invalid number of arguments');
               end

            case 0
               return
            otherwise
               error ('ChannelTag:InvalidConversion','Invalid number of arguments');
         end
      end

      function [IDX, objs] = matching(objs, N, S, L, C)
         %matching   Field-by-field comparison
         %
         % idx = cha_tags.matching(Net, Sta, Loc, Cha);
         % idx = cha_tags.matching('net.sta.loc.cha');
         % returns an index vector the same size as cha_tags
         % [idx, matched] = cha_tags.matching(...)
         %   additionally returns the actual matching cha_tags
         %
         % so, cha_tags(matches) will return the actual values
         % use [] or '' as wildcard
         %  either:
         %    cha_tags.matching('IU','ANMO',[],'BHZ')   or
         %    cha_tags.matching('IU.ANMO..BHZ')
         % will match BOTH
         %   'IU.ANMO.00.BHZ' and 'IU.ANMO.01.BHZ'
         if nargin == 2
            if (ischar(N) && any(N == '.')) || isa(N,'ChannelTag')
               % query used single string : 'net.sta.loc.cha'
               [N, S, L, C] = ChannelTag.parse(N);
            else
               error('CHANNELTAG:matching:UnknownMatchType','don''t know how to match it');
            end
         end
         IDX = true(size(objs));
         if ~isempty(N) && all(N ~= '*')
            IDX = strcmp({objs.network}, N);
         end
         if exist('S', 'var') && ~isempty(S) && all(S ~= '*')
            IDX = IDX & strcmp({objs.station}, S);
         end
         if exist('L', 'var')  && ~isempty(L) && all(L ~= '*')
            IDX = IDX & strcmp({objs.location}, L);
         end
         if exist('C', 'var') && ~isempty(C) && all(C ~= '*')
            IDX = IDX & strcmp({objs.channel}, C);
         end
         if nargout == 2
            objs = objs(IDX);
         end
      end

      function [Y, I] = sort(cha_tags)
         %sort   Sort channels assending name order
         %   sorted = sort(channeltags) will sort channels in assending
         %   order using the precidence: Network, Station, Location,
         %   Channel. The result will be an array of sorted ChannelTags.
         %
         %   [sorted, I] = sort(chantags) will additionally return the
         %   index, so that chantags(I) = sorted
         %
         %   sort is designed to handle SEED-like names, so it will pad each
         %   value to 5 places before comparing.
         %      Example:
         %        c = ChannelTag({'IU.ANMO.00.BHZ','AV.OKTU..EHZ'})
         %        sortedC = c.sort()  % or sort(c)
         %
         %  The above example would compare:
         %    'IU    .ANMO  .00    .BHZ   ' vs
         %    'AV    .OKTU  .      .EHZ   '
         %
         [~,I] = sort(cha_tags.fixedlengthstrings(5,5,5,5));
         Y = cha_tags(I);
      end
      result = eq(A, B)

      function result = ne(A, B)
         %~=   for ChannelTag.
         %   wilcards aren't considered
         result = ~strcmp(A.string,B.string);
      end

      function val = get(obj, prop)
         %GET for the chaTag object
         %  result = get(ChannelTag, property), where PROPERTY is one of the
         %  following:
         %    STATION, CHANNEL, LOCATION, NETWORK, chaTag_STRING
         %
         %  *Deprecated*
         %
         % If the results of a single chaTag are requested, then a string is returned.
         % Otherwise, a cell of values will be returned.

         %warning('Old Usage. instead of "get", ,use the field name directly.  eg. x = ch.station');

         prop = lower(prop);

         switch prop
            case{'station','channel','network','location'}
               val = {obj.(prop)};
            case {'nscl_string'}
               %stuff = chaTag.string;
               warning('CHANNELTAG:get:obsoleteUsage',...
                  ['get(''nsclstring'') obsolete. Use chaTag.string().\n'...
                  '  notice it will return N.S.L.C instead of N_S_C_L']);
               val=strcat({obj.network},'_',{obj.station},'_',{obj.channel},'_',{obj.location});
            otherwise
               error('CHANNELTAG:get:UnrecognizedProperty',...
                  'Unrecognized property name : %s',  upper(prop));
         end

         %if a single chaTag, then return the string representation, else return a
         %cell of strings.
         if numel(val) == 1
            val = val{1};
         else
            val = reshape(val,size(obj));
         end
      end%get

      function obj = set(obj, varargin)
         %SET  Set properties for ChannelTag
         %    s = Set(s,prop_name, val, ...)
         %    Valid property names:
         %    STATION, LOCATION, NETWORK, CHANNEL
         %
         %    * Deprecated*
         %    See also: get.channel, get.location, get.network, get.station

         %warning('Old usage: using set for a ChannelTag. use field notation instead. eg. ch.station = X');
         idx = [1 2];
         while numel(varargin) > 1 % for each property
            [prop, val] = deal(varargin{idx});
            varargin(idx) = []; % remove the prop and val
            switch upper(prop)
               case {'STATION','LOCATION','NETWORK','CHANNEL'}
                  if iscell(val) && numel(val) > 1
                     val = val{1};
                     warning('ChannelTag:set:TooManyValues',...
                        'Too many property values, only the first will be used');
                  end
                  [obj.(lower(prop))] = deal(val);
               otherwise
                  error('ChannelTag:set:UnrecognizedProperty',...
                     'Unrecognized property name : %s',  upper(prop));
            end; %switch
         end; %each property
         if ~isempty(varargin)
            error('ChannelTag:set:PropertyValueMismatch',...
               'Each property should have a value. [%s] has no matching value', varargin{1});
         end
      end

      function s = fixedlengthstrings(objs, netLen, staLen, locLen, chaLen)
         %fixedlengthstrings   return strings of a fixed length
         %    s = fixedlengthstrings(objs, netLen, staLen, locLen, chaLen)

         strformat = sprintf('%%-%ds.%%-%ds.%%-%ds.%%-%ds',netLen,staLen,locLen,chaLen);
         for n = numel(objs) : -1 : 1
            chaTag = objs(n);
            s(n) = {sprintf(strformat,chaTag.network,chaTag.station,chaTag.location,chaTag.channel)};
         end
      end

      function disp(obj)
         %disp   display the ChannelTag details
         if numel(obj) == 1
            disp('ChannelTag with network.station.location.channel:')
            fprintf('   network: ''%s''\n   station: ''%s''\n  location: ''%s''\n   channel: ''%s''\n',...
               obj.network,obj.station,obj.location, obj.channel);
         else
            disp('ChannelTag array (network.station.location.channel):')
            for n=1:numel(obj)
               if strcmp(obj(n).string,'...')
                  disp('  ... (empty ChannelTag)');
               else
                  disp(['  ', obj(n).string]);
               end
            end
         end
      end

      function c = char(obj)
         %char   convert a ChannelTag to its string representation
         c = obj.string([],'nocell');
      end

      function s = string(obj, delim, option)
         %string returns string representation of the channelTag(s)
         % s = chaTag.string()  will return the string representation
         %      (1xn char) in the format NET.STA.LOC.CHA
         %
         % s = chaTag.string(DELIM) will use DELIM to separate
         %     fields.
         %     ex.  chaTag.string('-'); will return NET-STA-LOC-CHA
         %
         % If chaTag is an array, then results are returned as a
         % a cell of strings of same shape as chaTag will be returned.
         %
         % s = chaTag.string(DELIM,'nocell') to overrides functionality
         % to return a padded NxM char array
         % if DELIM is empty, then '.' will be used.
         if ~exist('delim','var') || isempty(delim)
            delim = '.';
         end
         if numel(obj) == 1
            s = getDelimitedString(obj, delim);
         else
            if exist('option','var') && strcmpi(option,'nocell')
               s = '';
               for n=1 : numel(obj)
                  tmp = getDelimitedString(obj(n), delim);
                  s(n,1:numel(tmp)) = tmp;
               end
            else
               s = cell(size(obj));
               for n=1 : numel(obj)
                  s(n) = {getDelimitedString(obj(n), delim)};
               end
            end
         end

          function s = getDelimitedString(obj, delim)
            %getDelimitedString   return a string with delimited fields
            %   getDelimitedString(chtag, delim) returns the fields in
            %   network,station, location, channel order, with the
            %   specified delimiter in between.
            %
            %   example:
            %     >> ch = ChannelTag('NT.STA.00.BHZ');
            %     >> s = getDelimitedString(ch, ':');
            %     s =
            %        'NT:STA:00:BHZ'
            %
            s = [obj.network, delim, obj.station, delim,...
                     obj.location, delim, obj.channel];
          end
      end

      function s = scn(obj, delim, option)
         %scn returns scn representation of the channelTag(s)
         % s = chaTag.scn()  will return the string representation
         %      (1xn char) in the format STA_CHA_NET
         %
         % s = chaTag.string(DELIM) will use DELIM to separate
         %     fields.
         %     ex.  chaTag.string('-'); will return NET-STA-LOC-CHA
         %
         % If chaTag is an array, then results are returned as a
         % a cell of strings of same shape as chaTag will be returned.
         %
         % s = chaTag.string(DELIM,'nocell') to overrides functionality
         % to return a padded NxM char array
         % if DELIM is empty, then '.' will be used.
         if ~exist('delim','var') || isempty(delim)
            delim = '_';
         end
         if numel(obj) == 1
            s = [obj.station, delim, obj.channel, delim, obj.network];
         else
            if exist('option','var') && strcmpi(option,'nocell')
               s = '';
               for n=1 : numel(obj)
                  tmp = [obj.station, delim, obj.channel, delim, obj.network];
                  s(n,1:numel(tmp)) = tmp;
               end
            else
               s = cell(size(obj));
               for n=1 : numel(obj)
                  s(n) = {[obj.station, delim, obj.channel, delim, obj.network]};
               end
            end
         end

      end

   end%methods

   methods(Static)
      function  [N, S, L, C] = parse(val)
         %parse   parse a period-delimeted string into NET,STA,LOC,CHA
         %     [N, S, L, C] = ChannelTag.parse('net.sta.loc.cha')

         if isstruct(val) || isa(val,'ChannelTag')
            N = val.network;
            S = val.station;
            L = val.location;
            C = val.channel;
         elseif ischar(val)
            % test validity, then split.  in R2015a+ strsplit() could do
            % this. But, this is being designed on R2012a.
            parts =  strsplit(val, '.');
            if numel(parts) ~= 4
               error('ChannelTag:parse:InvalidFieldCount',...
                  'Expected ''A.B.C.D'' (4 fields), but received %d.', numel(parts));
            end
            if any(val == ' '),
               parts = strtrim(parts);
            end
            [N, S, L, C] = deal(parts{:});
         else
            error('ChannelTag:parse:UnknownClass',...
               'unknown how to parse: %s', class(val));
         end

         function  C = strsplit(s, delim)
            %strsplit   return a cell created from a string.
            % C = strsplit(str, delim) returns cell created from string
            %  ex. C = strsplit('A.B..D','.') -> {'A', 'B', '', 'D'}
            % can (probably) be commented out or removed in r2015a+
            splitpts = find(s == delim);
            starts = [1, splitpts+1];
            ends = [splitpts-1, numel(s)];
            C = cell(1,numel(starts));
            for n = 1:numel(starts)
               C(n) = {s(starts(n):ends(n))};
            end
         end
      end

      function cha_tags = array(varargin)
         % chaTag.array creates multiple ChannelTags from the input.
         % cha_tags = chaTag.array({'NET.STA.LOC.CHA','NET2.STA2..CH2'[,...]})
         % cha_tags = chaTag.array(Net/list, Sta/list, Loc/list, Cha/list)
         switch nargin
            case 0
               % return an empty nscltag array
               cha_tags = ChannelTag('...');
               cha_tags(:)=[];
            case 1
               % parse a cell array of stations. OR parse an NxM char, where
               % each row is NET.STA.LOC.CHA
               if iscell(varargin{1})
                  for n = numel(varargin{1}) : -1 : 1
                     cha_tags(n) = ChannelTag(varargin{1}{n});
                  end
                  cha_tags = reshape(cha_tags,size(varargin{1}));
               elseif ischar(varargin{1})
                  char_array = varargin{1};
                  nRows = size(char_array,1);
                  for n = nRows : -1 : 1
                     cha_tags(1,n) = ChannelTag(strtrim(char_array(n,:)));
                  end
               else
                  error('ChannelTag:array:UnknownInput','expected cell of strings or a char array');
               end
            case 4
               % expect 1xN char OR arbitrary-sized CELL of 1xN char
               % if multiple cells, then they should all be same size
               inputsThatAreCells = cellfun(@iscell,varargin);
               if any(inputsThatAreCells)
                  %make a cell array of N.S.L.C
                  A = strcat(varargin{1},'.',varargin{2},'.',varargin{3},'.',varargin{4});
               end
               for n = numel(A): -1 : 1
                  cha_tags(n) = ChannelTag(A{n});
               end
         end
      end

   end %static methods
end %classdef

author: geoscience-community-codes

repo: GISMO


text:
function updateLayout(mtexFig)
% resize figure and reorder subfigs

if isempty(mtexFig.children), return;end

% store old units and perform all calculations in pixel
old_units = get(mtexFig.parent,'Units');
set(mtexFig.parent,'Units','pixels');

figSize = get(mtexFig.parent,'Position');
figSize = figSize(3:4) - sum(reshape(mtexFig.figTightInset,2,2),2).';

% compute layout
[mtexFig.ncols,mtexFig.nrows] = calcPartition(mtexFig,figSize);
[mtexFig.axisWidth,mtexFig.axisHeight] = calcAxesSize(mtexFig,figSize);

% align axes according to layout
for i = 1:length(mtexFig.children)

  % compute position in raster
  [col,row] = ind2sub([mtexFig.ncols mtexFig.nrows],i);

  aw = mtexFig.axisWidth + mtexFig.innerPlotSpacing + sum(mtexFig.tightInset([1,3]));
  ah = mtexFig.axisHeight + mtexFig.innerPlotSpacing + sum(mtexFig.tightInset([2,4]));

  axisPos = [mtexFig.figTightInset(1:2),0,0] + ...
  [1 + (col-1)*aw + mtexFig.tightInset(1),...
    1 + figSize(2) - row * ah ...
    + mtexFig.innerPlotSpacing + mtexFig.tightInset(2),...
    mtexFig.axisWidth,mtexFig.axisHeight];
  axisPos(axisPos<0)=0;
  set(mtexFig.children(i),'Units','Pixel','Position',axisPos);

  % position the colorbars
  if numel(mtexFig.cBarAxis) == numel(mtexFig.children)

    resizeColorBar(mtexFig.cBarAxis(i))

  end
end

if length(mtexFig.cBarAxis)==1 && i>1
  pos = get(mtexFig.cBarAxis,'position');

  if pos(4)>pos(3) %Vertical bar

    pos(4) = mtexFig.nrows * mtexFig.axisHeight + ...
      (mtexFig.nrows-1) * (mtexFig.innerPlotSpacing + sum(mtexFig.tightInset([2,4])));
    pos(2) = axisPos(2)+1;
    pos(1) = mtexFig.ncols*(mtexFig.axisWidth + mtexFig.innerPlotSpacing + ...
      sum(mtexFig.tightInset([1,3]))) + mtexFig.outerPlotSpacing;

  else  %Horizontal bar

    pos(3)=mtexFig.ncols*(mtexFig.axisWidth) + ...
      (mtexFig.ncols-1) * (mtexFig.innerPlotSpacing + sum(mtexFig.tightInset([1,3]))); %c_bar width
    pos(2) = mtexFig.outerPlotSpacing + 2*pos(4);%axisPos(2) - 2*pos(4);
    pos(1) = mtexFig.outerPlotSpacing + mtexFig.tightInset(1); %c_bar left

  end
  set(mtexFig.cBarAxis,'position',pos);
end

% revert figure units
set(mtexFig.parent,'Units',old_units);

  function resizeColorBar(cBar)

    pos = get(cBar,'position');
    if pos(4) > pos(3) % vertical
      set(cBar,'position',...
        [axisPos(1)+mtexFig.axisWidth+10,...
        axisPos(2)+1,...
        pos(3),mtexFig.axisHeight-1]);
    else % horizonal
      set(cBar,'position',...
        [axisPos(1),...
        axisPos(2)-mtexFig.tightInset(2) + mtexFig.cBarShift,...
        mtexFig.axisWidth-1,pos(4)]);
    end
  end

function testit

close all
mtexFig = mtexFigure;
mtexFig.gca
rectangle('position',[0,0,1,1])
axis equal  tight
title('asdsa')
xlabel('asd')
mtexFig.nextAxis;
rectangle('position',[0,0,1,1])
axis equal tight
xlabel('asd')

title('asdasd2')
axis(mtexFig.children(1),'off')
axis(mtexFig.children(2),'off')

mtexFig.drawNow

end

end

author: mtex-toolbox

repo: mtex

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


text:
function [output_canvas] = image_blending_average(warped_img1,warped_img2)
    w1 = imfill(im2bw(uint8(warped_img1), 0),'holes');
    w2 = imfill(im2bw(uint8(warped_img2), 0),'holes');

    w1 = mat2gray(w1);
    w2 = mat2gray(w2);

    warped_img1 = double(warped_img1);
    warped_img2 = double(warped_img2);
    output_canvas(:,:,1) = ((warped_img1(:,:,1).*w1)+(warped_img2(:,:,1).*w2))./(w1+w2);
    output_canvas(:,:,2) = ((warped_img1(:,:,2).*w1)+(warped_img2(:,:,2).*w2))./(w1+w2);
    output_canvas(:,:,3) = ((warped_img1(:,:,3).*w1)+(warped_img2(:,:,3).*w2))./(w1+w2);
    output_canvas = uint8(output_canvas);
end

author: YaqiLYU

repo: AANAP

save_path: github-repos/MATLAB/YaqiLYU-AANAP


text:
function varargout = weights_int(mode, n, t0, Hdyn, Zdyn, Tdyn, Rdyn, Qdyn, Hmat, Zmat, Tmat, Rmat, Qmat, P, tol)
% mode:
%   0 - all primitive outputs.
%   1 - calculate Kalman filter weights.
%   2 - calculate Kalman smoother weights.

%% Initialization %%
m       = size(P, 1);
D       = (P == Inf);
init    = any(any(D)); % use exact initialization if init = true.
if init, d = n+1; P(D) = 0; P_inf = double(D);
else d = 0; end
stationary  = ~Hdyn && ~Zdyn && ~Tdyn && ~Rdyn && ~Qdyn;
converged   = false;
RQdyn       = Rdyn || Qdyn;
if ~Hdyn, H = Hmat; end
if ~Zdyn, Z = Zmat; end
if ~Tdyn, T = Tmat; end
if ~Rdyn, R = Rmat; end
if ~Qdyn, Q = Qmat; end
if ~RQdyn, RQRt = R*Q*R'; end

%% Preallocate Output Results %%
switch mode
    case 0 % all primitive outputs
        Result_P    = cell(1, n);
        Result_P{1} = P;
        Result_invF = cell(1, n);
        Result_K    = cell(1, n);
        Result_L    = cell(1, n);
        Result_N    = cell(1, n);
        Result_W    = cell(1, n);
    case 1 % calculate Kalman filter weights
        Result_K        = cell(1, t0-1);
        Result_L        = cell(1, t0-1);
        Result_omega    = cell(1, t0-1);
    case 2 % calculate Kalman smoother weights
        Result_invF         = cell(1, n);
        Result_K            = cell(1, n);
        Result_L            = cell(1, n);
        Result_LLasc        = cell(1, n);
        if t0 > 1, Result_LLasc{t0-1} = eye(m);
        else Result_P = P; end % t0 == 1
        Result_omega        = cell(1, t0-1);
        Result_omegaalpha   = cell(1, n);
end
Fns     = true(1, n); % Is F_inf nonsingular for each iteration

%% Kalman filter loop %%
for t = 1 : n
    if Hdyn, H = Hmat{t}; end
    if Zdyn, Z = Zmat{t}; end
    if Tdyn, T = Tmat{t}; end
    if Rdyn, R = Rmat{t}; end
    if Qdyn, Q = Qmat{t}; end
    if RQdyn, RQRt = R*Q*R'; end
    if ~converged
        if init
            %% Exact initial Kalman filter %%
            M       = P*Z';
            M_inf   = P_inf*Z';
            A_inf   = T*P_inf;
            if abs(M_inf) < tol % F_inf is zero
                Fns(t)  = false;
                invF    = inv(Z*M + H);
                K       = T*M*invF;
                L       = T - K*Z;
                P       = T*P*L' + RQRt;
                P_inf   = A_inf*T';
            else % F_inf is assumed to be nonsingular
                invF    = inv(Z*M_inf); % This is actually invF1
                F2      = -invF*(Z*M + H)*invF;
                K       = T*M_inf*invF;
                K1      = T*(M*invF + M_inf*F2);
                L       = T - K*Z;
                L1      = -K1*Z;
                P       = A_inf*L1' + T*P*L' + RQRt;
                P_inf   = A_inf*L';
            end
            if abs(P_inf) < tol, d=t; init=false; end
        else
            %% Normal Kalman filter %%
            M       = P*Z';
            invF    = inv(Z*M + H);
            K       = T*M*invF;
            L       = T - K*Z;
            prevP   = P;
            P       = T*P*L' + RQRt;
            if stationary, if abs(P-prevP) < tol, converged = true; end, end
        end
    end
    %% Store results for this time point %%
    switch mode
        case 0 % all primitive outputs
            Result_P{t+1}   = P;
            Result_invF{t}  = invF;
            Result_K{t}     = K;
            Result_L{t}     = L;
        case 1 % calculate Kalman filter weights
            Result_K{t}     = K;
            Result_L{t}     = L;
            if t >= t0 - 1, break; end
        case 2 % calculate Kalman smoother weights
            if t == t0 - 1, Result_P    = P; end
            if t == t0, Result_LLasc{t} = L'; end
            Result_invF{t}  = invF;
            Result_K{t}     = K;
            Result_L{t}     = L;
            if t > t0
                Result_LLasc{t} = Result_LLasc{t-1}*L';
            end
    end
end

%% Backwards recursion %%
switch mode
    case 0 % all primitive outputs
        N       = zeros(m);
        for t = n : -1 : 1
            if Hdyn, H = Hmat{t}; end
            if Zdyn, Z = Zmat{t}; end
            invF        = Result_invF{t};
            K           = Result_K{t};
            L           = Result_L{t};
            Result_W{t} = H*(invF*Z - K'*N*L);
            Result_N{t} = N;
            if t <= d && Fns(t), N = L'*N*L; else N = Z'*invF*Z + L'*N*L; end
        end
    case 1 % calculate Kalman filter weights
        LLdsc   = eye(m);
        for t = t0-1 : -1 : 1
            Result_omega{t} = LLdsc*Result_K{t};
            LLdsc           = LLdsc*Result_L{t};
        end
    case 2 % calculate Kalman smoother weights
        N       = zeros(m);
        LLdsc   = eye(m);
        for t = n : -1 : 1
            if t >= t0
                if t <= d && Fns(t)
                    K       = Result_K{t};
                    L       = Result_L{t};
                    W2      = -K'*N*L;
                    N       = L'*N*L;
                else
                    if Zdyn, Z = Zmat{t}; end
                    invF    = Result_invF{t};
                    K       = Result_K{t};
                    L       = Result_L{t};
                    W2      = invF*Z - K'*N*L;
                    N       = Z'*invF*Z + L'*N*L;
                end
                if t == t0
                    Result_omegaalpha{t} = Result_P*W2';
                    IPN = eye(m) - Result_P*N;
                else Result_omegaalpha{t} = Result_P*Result_LLasc{t-1}*W2'; end
            else
                Result_omega{t}         = LLdsc*Result_K{t};
                LLdsc                   = LLdsc*Result_L{t};
                Result_omegaalpha{t}    = IPN*Result_omega{t};
            end
        end
end

%% Output Results %%
switch mode
    case 0 % all primitive outputs
        varargout = {Result_P Result_invF Result_K Result_L Result_W Result_N};
    case 1 % calculate Kalman filter weights
        varargout = {Result_omega};
    case 2 % calculate Kalman smoother weights
        varargout = {Result_omega Result_omegaalpha};
end


text:
% Functions specific for simulating scenarios.
%
% 3D objects (with edges information)
%   obj2pnts           - Get points matrix from 3D object.
%   obj2segs           - Get segments matrix from 3D object.
%
% 3D patches (objects without edges information)
%   camGraphics        - Create a camera graphics 3D object.
%   camSimpleGraphics  - Create a simple camera graphics 3D object.
%   swing              - Create a graphics structure with a wing shape.
%   thickVehicle       - Create a vehicle graphics 3D object.
%   wing               - Create a wing graphics 3D object.
%
% 3D points (3-by-N matrices)
%   obj2pnts           - Get points matrix from 3D object.
%   cloister           - Generates features in a 2D cloister shape.
%   thickCloister      - Generates features in a 3D cloister shape.
%
% 3D segments (6-by-N matrices)
%   obj2segs           - Get segments matrix from 3D object.
%   makeSegment        - Make a segment out of two endpoints
%   xSegment           - Horizontal segment in the X direction
%   ySegment           - Horizontal segment in the Y direction
%   zSegment           - Vertical segment
%   XYRectangle        - Rectangle in the XY plane.
%   XZRectangle        - Rectangle in the XZ plane.
%   YZRectangle        - Rectangle in the YZ plane.
%   house              - House made out of segments.
%
% Work with pre-processed data
%   readControlSignal  - Read robot-s control signal from file.
%   readProcessedImg   - Read processed image data from file.
%   writeControlSignal - Write robot's control signal to file.
%   writeProcessedImg  - Write processed image data to file.

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function s = weekday_to_name_islamic ( w )

%*****************************************************************************80
%
%% WEEKDAY_TO_NAME_ISLAMIC returns the name of an Islamic weekday.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer W, the weekday index.
%
%    Output, string S, the weekday name.
%
  weekday_name = [ ...
    'Yom ilHadd   '; ...
    'Yom litneen  '; ...
    'Yom ittalat  '; ...
    'Yom larba    '; ...
    'Yom ilkhamiis'; ...
    'Yom ilguma   '; ...
    'Yom issabt   '];
%
%  Make a local copy of the weekday number.
%
  w2 = w;
%
%  Check the weekday number.
%
  weekday_check_common ( w2 );
%
%  Return the weekday name.
%
  s = weekday_name(w2,:);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function R = dark_channel_prior_dehaze(I)
%DARK_CHANNEL_PRIOR_DEHAZE  Dehaze input image following the pipeline proposed
%in Single Image Haze Removal Using Dark Channel Prior.
%
%   INPUTS:
%
%   -|I|: H-by-W-by-|image_channels| hazy image.
%
%   OUTPUTS:
%
%   -|R|: output dehazed image representing the estimate for the true radiance
%   of the scene, with same size as |I|.

% Add required paths.
current_script_full_name = mfilename('fullpath');
current_script_directory = fileparts(current_script_full_name);
addpath(fullfile(current_script_directory, '..', '..', 'utilities'));
addpath_relative_to_caller(current_script_full_name,...
    fullfile('..', '..', 'Fog_simulation'));
addpath_relative_to_caller(current_script_full_name,...
    fullfile('..', '..', 'Dehazing'));

% Set parameters.
neighborhood_size_dark_channel = 15;
t_thresh = 0.1;
window_size_guided_filter = 41;
epsilon = 1e-3;

% Dehaze using Dark Channel Prior, following the paper of He, Sun and Tang (CVPR
% 2009). The refinement of the raw transmission map is performed according to
% Guided Image Filtering by He, Sun and Tang (ECCV 2010), in order to gain
% speed.
[I_dark, I_eroded] = get_dark_channel(I, neighborhood_size_dark_channel);
L = estimate_atmospheric_light_dcp(I_dark, I);
t_initial = transmission_initial(I_eroded, L);
t = transmission_guided_filtering(t_initial, I,...
    window_size_guided_filter, epsilon);
t = clip_to_unit_range(t);
R = inverse_haze_linear(I, t, L, t_thresh);
R = clip_to_unit_range(R);

end

author: sakaridis


text:
function SO3F = approximation(nodes, y, varargin)
% computes a least square problem to get an approximation
%
% Syntax
%   SO3F = SO3FunHarmonic.approximation(SO3Grid, f)
%   SO3F = SO3FunHarmonic.approximation(SO3Grid, f,'constantWeights')
%   SO3F = SO3FunHarmonic.approximation(SO3Grid, f, 'bandwidth', bandwidth, 'tol', TOL, 'maxit', MAXIT, 'weights', W)
%
% Input
%  SO3Grid - rotational grid
%  f       - function values on the grid (maybe multidimensional)
%
% Options
%  constantWeights  - uses constant normalized weights (for example if the nodes are constructed by equispacedSO3Grid)
%  bandwidth        - maximum degree of the Wigner-D functions used to approximate the function
%  tol              - tolerance for lsqr
%  maxit            - maximum number of iterations for lsqr
%  weights          - weight w_n for the nodes (default: Voronoi weights)
%

if isa(nodes,'orientation')
  SRight = nodes.CS; SLeft = nodes.SS;
  if nodes.antipodal
    nodes.antipodal = 0;
    varargin{end+1} = 'antipodal';
  end
else
  [SRight,SLeft] = extractSym(varargin);
  nodes = orientation(nodes,SRight,SLeft);
end

% make points unique
s = size(y);
y = reshape(y, length(nodes), []);
[nodes,~,ind] = unique(nodes(:),'tolerance',0.02);

% take the mean over duplicated nodes
for k = 1:size(y,2)
  yy(:,k) = accumarray(ind,y(:,k),[],@nanmean); %#ok<AGROW>
end

y = reshape(yy, [length(nodes) s(2:end)]);

tol = get_option(varargin, 'tol', 1e-6);
maxit = get_option(varargin, 'maxit', 50);

% TODO: What bandwidth?
bw = get_option(varargin, 'bandwidth', min(dim2deg(length(nodes)*2),getMTEXpref('maxSO3Bandwidth')));

W = get_option(varargin, 'weights');
if check_option(varargin,'constantWeights')
  W = 1/length(nodes);
elseif isempty(W)
  % TODO: calcVoronoiVolume is bad estimated
  [~,i,~] = unique(quaternion(nodes));
  nodes=nodes(i); y=y(i,:);
  W = calcVoronoiVolume(nodes);
else
  W = accumarray(ind,W);
end

W = sqrt(W(:));

b = W.*y;

% create plan
% SO3FunHarmonic([1;1]).eval(nodes,'createPlan','nfsoft')
% SO3FunHarmonic.quadrature(nodes,1,'createPlan','nfsoft','bandwidth',bw)

% least squares solution
for index = 1:size(y,2)
  [fhat(:, index),flag] = lsqr( ...
    @(x, transp_flag) afun(transp_flag, x, nodes, W,bw), b(:, index), tol, maxit);
end

% kill plan
% SO3FunHarmonic.quadrature(1,'killPlan','nfsoft')
% SO3FunHarmonic(1).eval(1,'killPlan','nfsoft')

SO3F = SO3FunHarmonic(fhat,SRight,SLeft);

end

function y = afun(transp_flag, x, nodes, W,bw)

if strcmp(transp_flag, 'transp')

  x = x.*W;
%   F = SO3FunHarmonic.quadrature(nodes,x,'keepPlan','nfsoft','bandwidth',bw);
  F = SO3FunHarmonic.quadrature(nodes,x,'bandwidth',bw);
  y = F.fhat;

elseif strcmp(transp_flag, 'notransp')

  F = SO3FunHarmonic(x,nodes.CS,nodes.SS);
  F.bandwidth = bw;
%   y = F.eval(nodes,'keepPlan','nfsoft');
  y = F.eval(nodes);
  y = y.*W;

end

end

% Possibly split the quadrature and eval functions. Therefore do the
% precomputations before the lsqr-Method.

% TODO: Try Cross-Vallidation, if there are to much points
% TODO: Possibly use the Round2ClenshawCurtis-function (see quadrature) or delete it


text:
function test_errors
%TEST_ERRORS tests error handling for the factorize object methods
%
% Example
%   test_errors
%
% See also test_all, factorize1, factorize.

% Copyright 2009, Timothy A. Davis, University of Florida

% the matrix A must be square for factorize1
fprintf ('\nTesting error-handling, error messages expected:\n') ;
ok = true ;
A = rand (3,2) ;
try
    F = factorize1 (A) ;                                               %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% the matrix A must be 2D for factorize
A = ones (2,2,2) ;
try
    F = factorize (A) ;                                                %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% the matrix A must be full-rank
A = ones (3) ;
try
    F = factorize (A) ;                                                %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

try
    F = factorize1 (A) ;                                               %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

A = ones (3,2) ;
try
    F = factorize (A) ;                                                %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

A = ones (2,3) ;
try
    F = factorize (A) ;                                                %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% cannot use cell indexing
A = rand (3,2) ;
F = factorize (A) ;
try
    C = F {1} ;                                                        %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% cannot use linear indexing of the inverse
S = inverse (A) ;
try
    C = S (1) ;                                                        %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% invalid indexing
try
    C = F (1,1).L ;                                                    %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% invalid indexing
try
    C = F.L (1,1).stuff ;                                              %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% non-existent field
try
    C = F.junk ;                                                       %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% can only update/downdate a dense Cholesky factorization
A = rand (2) ;
F = factorize (A) ;
w = rand (2,1) ;
try
    F = F + w ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end
try
    F = F - w ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

if (~ok)
    error ('error-handling failed') ;
end
fprintf ('\nAll error-handing tests passed\n') ;


text:
% .............................................................
    function [shape,naturalDerivatives]=shapeFunctionQ4(xi,eta)

    % shape function and derivatives for Q4 elements
    % shape : Shape functions
    % naturalDerivatives: derivatives w.r.t. xi and eta
    % xi, eta: natural coordinates (-1 ... +1)

    shape=1/4*[ (1-xi)*(1-eta);
                (1+xi)*(1-eta);
                (1+xi)*(1+eta);
                (1-xi)*(1+eta)];
    naturalDerivatives=...
           1/4*[-(1-eta), -(1-xi);
                   1-eta, -(1+xi);
		           1+eta,    1+xi;
                -(1+eta),    1-xi];

    end % end function shapeFunctionQ4

repo: Aero-matlab


text:
function k = biasKernDiagCompute(kern, x)

% BIASKERNDIAGCOMPUTE Compute diagonal of BIAS kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the bias kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : biasKernParamInit, kernDiagCompute, kernCreate, biasKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

k = repmat(kern.variance, size(x, 1), 1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function d = distancePointLine3d(point, line)
%DISTANCEPOINTLINE3D Euclidean distance between 3D point and line.
%
%   D = distancePointLine3d(POINT, LINE);
%   Returns the distance between point POINT and the line LINE, given as:
%   POINT : [x0 y0 z0]
%   LINE  : [x0 y0 z0 dx dy dz]
%   D     : (positive) scalar
%
%   See also:
%   lines3d, isPointOnLine3d, distancePointEdge3d, projPointOnLine3d,
%
%
%   References
%   http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 23/05/2005.
%

%   HISTORY
%   15/01/2007 unify size of input data
%   31/01/2007 typo in data formatting, and replace norm by vecnorm3d
%   12/12/2010 changed to bsxfun implementation - Sven Holcombe

% cf. Mathworld (distance point line 3d)  for formula
d = bsxfun(@rdivide, vectorNorm3d( ...
        crossProduct3d(line(:,4:6), bsxfun(@minus, line(:,1:3), point)) ), ...
        vectorNorm3d(line(:,4:6)));

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [c1,c2] = match(c1,c2,toleranceInSeconds)

   %Retain matching elements of two correlation objects.
   %
   % [C1,C2] = MATCH(C1,C2) Returns those elements of c1 and c2 which have
   % matching trigger times. Typically this will be used to sync correlation
   % objects from two different stations covering the same time frames, or
   % sync'ing two channels from the same station.
   %
   % [C1,C2] = MATCH(C1,C2,TOLERANCE) allows traces to match when their
   % trigger times are within TOLERANCE seconds of each other. When this
   % parameter is not included, the default value is 1 second.
   %
   % CAVEATS: Note that if C2 has multiple traces which match a trigger time
   % in C1, there is no way of controling which single trace in C2 will be
   % considered a match. Similarly, if C1 contains multiple traces with the
   % same trigger time, a match for each one will be found in C2.

   % AUTHOR: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   narginchk(2,3);

   if ~exist('tolerance','var')
      toleranceInSeconds = 1;
   end

   assert(isnumeric(toleranceInSeconds), 'Tolerance must be a number in seconds')

   disp(['Tolerance for successful match is: ' num2str(toleranceInSeconds) ' seconds']);

   % FIND ELEMENTS THAT MATCH
   trig1 = c1.trig;
   trig2 = c2.trig;
   matchIndex = zeros(size(trig1));
   for n = 1:numel(trig1)
      [minVal,minIndex] = min(abs(trig2 - trig1(n))*86400);
      if minVal < toleranceInSeconds
         matchIndex(n) = minIndex;
      else
         matchIndex(n) = 0;
      end
   end

   % SUBSET FIRST CORRELATION OBJECT
   index = find(matchIndex);
   c1 = subset(c1,index);

   % SUBSET FIRST CORRELATION OBJECT
   index2 = matchIndex(index);
   c2 = subset(c2,index2);

end


text:
function objset=readAnnotation(Name_batch,theConf)
MaxObjNum=5000;

minArea=theConf.data.minArea;
objset(MaxObjNum).folder=[];
objset(MaxObjNum).filename=[];
objset(MaxObjNum).name=[];
objset(MaxObjNum).bndbox=[];
objset(MaxObjNum).ID=[];
load([theConf.data.imgdir,Name_batch,'_obj/img/data.mat'],'samples');
root=[theConf.data.imgdir,Name_batch,'_obj/img/img/'];
files=dir(root);
for i=1:length(files)-2
    folder=[Name_batch,'_obj/img/img/'];
    filename=sprintf('%05d.jpg',i);
    xmin=samples(i).obj.bndbox.xmin;
    ymin=samples(i).obj.bndbox.ymin;
    xmax=samples(i).obj.bndbox.xmax;
    ymax=samples(i).obj.bndbox.ymax;
    bndbox.xmin=int2str(xmin);
    bndbox.xmax=int2str(xmax);
    bndbox.ymin=int2str(ymin);
    bndbox.ymax=int2str(ymax);
    if(~IsAreaValid(bndbox,minArea))
        continue;
    end
    objset(i).folder=folder;
    objset(i).filename=filename;
    objset(i).name=filename;
    objset(i).bndbox=bndbox;
    objset(i).ID=i;
end
objset=objset(1:length(files)-2);
end

function pd=IsAreaValid(bndbox,minArea)
xmin=str2double(bndbox.xmin);
xmax=str2double(bndbox.xmax);
ymin=str2double(bndbox.ymin);
ymax=str2double(bndbox.ymax);
if((xmax-xmin+1)*(ymax-ymin+1)>=minArea)
    pd=true;
else
    pd=false;
end
end

repo: interpretableCNN


text:
path('toolbox',path);
%% read pcd and show it
clear;clc;close all;
filename = '../data/2D/n.mat';
load(filename);
P.pts = M.verts;
clear('M');
npts = size(P.pts, 1);
[P.bbox, P.diameter] = GS.compute_bbox(P.pts);

%% call farthest_sampling_by_sphere
P.sample_radius = P.diameter*0.02;
[P.spls,P.corresp] = farthest_sampling_by_sphere(P.pts, P.sample_radius);

%% call connect_by_inherit_neigh
k=5;
kdtree = kdtree_build(P.pts);
P.neigh = zeros(npts, k);
for i = 1:npts
%     P.neigh(i,:)  = flipud( kdtree_k_nearest_neighbors(kdtree, P.pts(i,:), k))';
    P.neigh(i,:) = kdtree_k_nearest_neighbors( kdtree, P.pts(i,:), k)';
end

P.spls_adj = connect_by_inherit_neigh(P.pts, P.spls, P.corresp, P.neigh);

figure; movegui('northeast');set(gcf,'color','white');hold on;
plot3( P.spls(:,1), P.spls(:,2), P.spls(:,3), '.r', 'markersize', 1);
axis off; axis equal;set(gcf,'Renderer','OpenGL');
GS.plot_connectivity(P.spls, P.spls_adj, 1)
view3d zoom;

kdtree_delete( kdtree );

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
close all;
clear all;
clc;
rng('default');

png_export = true;
pdf_export = false;

mf = spx.graphics.Figures();

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

mf.new_figure('Recovery probability with S for different MMV recovery algorithms');

load ('bin/success_with_s_comparison.mat');

hold all;
legends = cell(1, 4);

plot(Ss, success_with_s.ra_ormp, '-+');
legends{1} = 'RA-ORMP';
% plot(Ss, success_with_s.ra_omp, '-o');
% legends{2} = 'RA-OMP';
plot(Ss, success_with_s.somp, '-s');
legends{2} = 'SOMP';
plot(Ss, success_with_s.gomp_2_mmv, '-d');
legends{3} = 'GOMP-MMV (L=2)';
plot(Ss, success_with_s.gomp_4_mmv, '-x');
legends{4} = 'GOMP-MMV (L=4)';

grid on;
xlabel('Number of signals');
ylabel('Recovery Probability');
legend(legends, 'Location', 'southeast');
title('Comparison of recovery performance for MMV algorithms');

author: indigits

repo: sparse-plex


text:
function destPath = compute_visual_feats(destPath, varargin)

	opts.gpus = 2 ;
  opts.limit = inf ;
  opts.clobber = false ;
  opts.batchSize = 128 ;
  opts.numEmotions = 8 ;
  opts.targetDataset = 'rml' ;
  opts.modelName = 'senet50-ferplus' ;
  opts = vl_argparse(opts, varargin) ;
  opts.imdbPath = fullfile(vl_rootnn, ...
             sprintf('data/xEmo18/%s_imdb/imdb.mat', opts.targetDataset)) ;

  if exist(destPath, 'file') && ~opts.clobber
    fprintf('found features at %s... skipping\n', destPath) ;
    return
  end

  switch opts.targetDataset
    case 'afew'
      opts.dataDir = fullfile(vl_rootnn, 'data/datasets/emotiw2016') ;
      getImdb = @(x) getAfewImdb(x, 'dropTracksWithNoDets', 1) ;
    case 'afew-6' % extracted at a stride of 6
      opts.dataDir = fullfile(vl_rootnn, 'data/datasets/emotiw2016') ;
      getImdb = @(x) getAfewImdb(x, 'dropTracksWithNoDets', 1, ...
                                 'subsampleStride', 6) ;
    case 'enterface'
      opts.dataDir = fullfile(vl_rootnn, 'data/datasets/enterface') ;
      getImdb = @(x) getEnterfaceImdb(x) ;
    case 'rml'
      opts.dataDir = fullfile(vl_rootnn, 'data/datasets/rml') ;
      getImdb = @(x) getRmlImdb(x) ;
    otherwise, error('unknown dataset %s\n', opts.targetDataset) ;
  end

  % work relative to audio files
  if exist(opts.imdbPath, 'file')
    fprintf('loading imdb from memory...') ; tic ;
    imdb = load(opts.imdbPath) ;
    fprintf('done in %g s\n', toc) ;
  else
    fprintf('building imdb...') ; tic ;
    imdb = getImdb(opts) ;
    mkdir(fileparts(opts.imdbPath)) ;
    save(opts.imdbPath, '-struct', 'imdb') ; % avoid compression if poss
    fprintf('done in %g s\n', toc) ;
  end

  [dag, pretrained] = ferPlusZoo(opts.modelName) ;
  assert(pretrained, 'loaded incorrect model') ;

  % remove losses
  lossLayers = arrayfun(@(x) isa(x.block, 'dagnn.Loss'), dag.layers) ;
  removables = {dag.layers(lossLayers).name} ;
  for ii = 1:numel(removables)
    dag.removeLayer(removables{ii}) ;
  end
	dag.mode = 'test' ;
	if numel(opts.gpus), gpuDevice(opts.gpus) ; dag.move('gpu') ; end

	inVars = dag.getInputs() ;
	assert(numel(inVars) == 1, 'too many inputs') ;

  % flatten frames into list, and keep track of the mapping to cells
  frames = zs_flattenCell(imdb.tracks.paths) ;
  frameIdx = arrayfun(@(x) ...
                     {x*ones(1, numel(imdb.tracks.paths{x}))}, ...
                     1:numel(imdb.tracks.paths)) ;
  frameIdx = [frameIdx{:}] ;
  assert(numel(frames) == numel(frameIdx), 'incorrect idx pairing') ;

  % compute for first `limit` tracks
  firstId = frameIdx(1) ; % use first in partition as offset
  numKeep = sum(frameIdx <= firstId + opts.limit) ; % number of tracks
  numIms = min(numel(frames), numKeep) ;
  meta = dag.meta ;
	inVars = dag.getInputs() ;
	assert(numel(inVars) == 1, 'too many inputs') ;

  logits = zeros(numIms, opts.numEmotions, 'single') ;

	fprintf('processing images with %s\n', opts.modelName) ;
	for ii = 1:opts.batchSize:numIms
		tic ;
		batchStart = ii ;
		batchEnd = min(batchStart+opts.batchSize-1, numIms) ;
		batch = batchStart:batchEnd ;
    imPaths = fullfile(opts.dataDir, frames(batch)) ;
		data = getImageBatch(imPaths, dag, meta) ;
		dag.eval({inVars{1}, data}) ;
		out = gather(squeeze(dag.vars(end).value)) ; % risky use of end variable
		logits(batch, :) = out' ;
		rate = numel(batch) / toc ;
    etaStr = zs_eta(rate, ii, numIms) ;
		fprintf(['processed image %d/%d at (%.3f Hz) ' ...
             '(%.3f%% complete) (eta:%s)\n'], ...
             ii, numIms, rate, 100 * ii/numIms, etaStr) ;
	end

  % save the logits alongside the frames for easier processing during training
  keptIdx = frameIdx(1:numIms) ; msg = 'unexpected indices' ;
  assert(isequal(1:numel(unique(keptIdx)), unique(keptIdx)), msg) ;
  numKeptTracks = numel(unique(keptIdx)) ;

  faceLogits = cell(1, numKeptTracks) ;
  for ii = 1:numKeptTracks
    fprintf('splitting logits into cells (%d/%d)\n', ii, numKeptTracks) ;
    sel = frameIdx == ii ;
    faceLogits(ii) = {logits(sel,:)} ;
  end

	fprintf('saving features to %s...', destPath) ; tic ;
  imdb.faceLogits = faceLogits ;
	if ~exist(fileparts(destPath), 'dir')
		zs_mkdirRec(fileparts(destPath)) ;
	end
	save(destPath, '-struct', 'imdb') ;
	fprintf('done in %g (s)\n', toc) ;
	if numel(opts.gpus), dag.move('cpu') ; end
end

% --------------------------------------------------------
function data = getImageBatch(imagePaths, dag, meta)
% --------------------------------------------------------
  numThreads = 10 ;
  opts.prefetch = false ; % can optimise this
  imageSize = meta.normalization.imageSize(1:2) ;

  % Afew is already tightly cropped
  args{1} = {imagePaths, ...
             'NumThreads', numThreads, ...
             'Pack', ...
             'Interpolation', 'bilinear', ... % bilinear to reproduce trainig
             'Resize', imageSize, ...
             'CropLocation', 'center'} ; % centre crop for testing
  args = horzcat(args{:}) ;

  if opts.prefetch
    vl_imreadjpeg(args{:}, 'prefetch') ;
    data = [] ;
  else
    data = vl_imreadjpeg(args{:}) ;
  end

  % throw away color, since it was not used during training
	data = uint8(data{1}) ;
	data_ = permute(data, [1 2 4 3]) ; sz = size(data_) ;
	data_ = reshape(data_, sz(1), sz(2) * sz(3), 3) ;
	gray = rgb2gray(data_) ;
	gray = single(reshape(gray, sz(1), sz(2), 1, [])) ;
	data = normalizeFace(gray, dag) ;
	if strcmp(dag.device, 'gpu'), data = gpuArray(data) ; end
end

% --------------------------------------------------------------------
function face = normalizeFace(greyFace, dag)
% --------------------------------------------------------------------
  greyFace = imresize(greyFace, dag.meta.normalization.imageSize(1:2)) ;
  face = repmat(greyFace, [1 1 3]) ;

  % handle average image consisting of color channel averages
  avgIm = reshape(dag.meta.normalization.averageImage, 1, 1, 3) ;
  face = bsxfun(@minus, face, avgIm) ;
end

author: albanie

repo: mcnCrossModalEmotions


text:
function [cu, cv, cx, cz, footprint, feat] = poly2contacts(dt, u, v, imsize, v0, yc, f, minp, im)

if ~exist('minp', 'var')
    minp = 0.5;
end

imh = imsize(1);
imw = imsize(2);

if isempty(f)
    f = 1.38; %S*max(size(im)) / imh;
end

u = u(:)';
v = v(:)';
if exist('im', 'var')
    figure(1), hold off, imshow(im), hold on
    plot([u(:) ; u(1)]', [v(:) ; v(1)]', '-b');
end

[tmp, ind] = min(u);
u = [u(ind:end)  u(1:ind-1)];
v = [v(ind:end)  v(1:ind-1)];

u = (u - imw/2) ./ imh;
v = 1 - (v ./ imh);

if exist('im', 'var')
    figure(1), plot(u*imh+imw/2, (1-v)*imh, '-g')
end
cind = convhull(u*imh+imw/2, (1-v)*imh);
if exist('im', 'var')
    figure(1), plot(u(cind)*imh+imw/2, (1-v(cind))*imh, '-y')
end

[v1, ind1] = min(v); % lowest point on object in image

[footx, footz] = computeGroundPosition([min(u) max(u)], [v1 v1], v0, yc, f);
footz(2) = footz(1) + footx(2)-footx(1);

[cx, cz] = computeGroundPosition(u, v, v0, yc, f);

data.u = (u*imh+imw/2)';
data.v = ((1-v)*imh)';
data.x3d = cx;
data.z3d = cz;
data.foot = [footx footz];

feat = contactdata2features(data);

[b, nodes] = treeval(dt, feat);

b = b - 1;

ind = dt.classprob(nodes, 2)>minp;%(b==1);

cu = u(ind)*imh+imw/2;
cv = (1-v(ind))*imh;
cz = cz(ind);
cx = cx(ind);
footprint = data.foot;

if exist('im', 'var')
    figure(1), plot(cu, cv, '*r')
    figure(1), plot([1 1000], (1-v0)*imh*ones(1,2), '-b');
end
%keyboard
% get how many points are in the footprint
%figure(2), plot(cx, cz, '*'), axis equal

%disp(num2str([cx ; cz]))

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

function [x, z] = computeGroundPosition(u, v, v0, yc, f)

z = yc*f./max((v0-v), 0.001);
x = u.*z./f;

author: Cloud-CV

repo: object-proposals


text:
% CVX: Matrix structure definitions and utilities.
%    CVX provides a keyword-based method for definiting matrices
%    with one or more types of structure; e.g.
%        variable X(n,n) symmetric toeplitz tridiagonal;
%    CVX automatically computes an efficient basis for the requested
%    structure. The files in this directory implement those computations.
%
%    None of these files should be called directly---matrix structure is
%    selected in the VARIABLE declaration; see VARIABLE for more details.
%    Below are the keywords that are available, and the structures they
%    represent. Keywords can be freely combined (see the above example),
%    but of course some combinations are degenerate, yielding only the
%    all-zero matrix; e.g.,
%       variable X(n,n)
%
% Structures:
%   banded            - (U,L)-banded matrices.
%   complex           - Complex variables of all sizes.
%   diagonal          - Diagonal matrices.
%   hankel            - Hankel matrices.
%   hermitian         - Complex Hermitian matrices.
%   lower_bidiagonal  - Lower bidiagonal matrices.
%   lower_hessenberg  - Lower Hessenberg matrices.
%   lower_triangular  - Lower triangular matrices.
%   scaled_identity   - Scaled identity: t*eye(n).
%   skew_symmetric    - Skew-symmetric matrices.
%   sparse            - Matrices with a fixed sparsity pattern.
%   symmetric         - Symmetric matrices.
%   toeplitz          - Toeplitz matrices.
%   tridiagonal       - Tridiagional matrices.
%   upper_bidiagonal  - Upper bidiagonal matrices.
%   upper_hankel      - Upper Hankel matrices.
%   upper_hessenberg  - Upper Hessenberg matrices.
%   upper_triangular  - Upper triangular matrices.

% 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:
%  Figure 6.61      Feedback Control of Dynamic Systems, 6e
%                   Franklin, Powell, Emami
%

clear all
%close all
clf

num=9;
den=conv([1 0.5],[1 1]);
den=conv(den,[1 2]);
w=logspace(-1,1,500);
[mag,phas]=bode(num,den,w);
[OLGM,OLPM,OLWcg,OLWcp]=margin(mag,phas,w)

%Lead compensator, first guess
numl=[1 1];
denl=0.333*[1 3];
numc=conv(num,numl);
denc=conv(den,denl);
[magc,phasc]=bode(numc,denc,w);
[D1GM,D1PM,D1Wcg,D1Wcp]=margin(magc,phasc,w)
dencl=denc+[0 0 0 numc];
t=0:.1:20;
y=step(numc,dencl,t);
%Design Iteration, next guess
numl=[1 1.5];
denl=0.1*[1 15];
numcc=conv(num,numl);
dencc=conv(den,denl);
[magcc,phascc]=bode(numcc,dencc,w);
[D2GM,D2PM,D2Wcg,D2Wcp]=margin(magcc,phascc,w)
%subplot(2,1,1)
%loglog(w,mag,'-',w,magc,'--',w,magcc,'-.',w,ones(500,1),'-');
%grid;
%xlabel('w (rad/sec)');
%ylabel('Magnitude');
%title('Fig. 6.60 Bode Plot for lead-compensation design (a) magnitude');
%subplot(2,1,2)
%semilogx(w,phas,'-',w,phasc,'--',w,phascc,'-.',w,-180*ones(500,1));
%grid;
%xlabel('w (rad/sec)');
%ylabel('phase (deg)');
%title('Fig. 6.60 (b) phase');
dencl=dencc+[0 0 0 numcc];
t=0:.1:20;
yy=step(numcc,dencl,t);
plot(t,y,'--',t,yy,'-');
xlabel('Time (sec)');
ylabel('y');
title('Fig. 6.61  Step Response');
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
function g = sdrbfKernGradient(sdrbfKern, t1, t2, covGrad)

% SDRBFKERNGRADIENT Gradient of SDRBF kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the switching
% dynamical radial basis function kernel's parameters. As well as the
% kernel structure and the input positions, the user provides a matrix
% PARTIAL which gives the partial derivatives of the function with respect
% to the relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t : the input locations for which the gradients are being
% computed.
% ARG covGrad : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in t1.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t1 : the input locations associated with the rows of the
% kernel matrix.
% ARG t2 : the input locations associated with the columns of the
% kernel matrix.
% ARG covGrad : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as t1 and the same number of columns
% as t2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO sdrbfKernParamInit, kernGradient
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    covGrad = t2;
    t2 = t1;
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end

%Form the basic kernels
rbfKern = struct();

% Create structures that will make easy the computation of the kernels
% spVector = [sdrbfKern.switchingTimes t1(end)+0.1]; % For the last interval include the last point

spVector = [cumsum(sdrbfKern.switchingTimes) t1(end)+100];

dim1 = zeros(1, sdrbfKern.nIntervals);
dim2 = zeros(1, sdrbfKern.nIntervals);

for i =1:sdrbfKern.nIntervals
    for j=1:sdrbfKern.nlfPerInt
        rbfKern(i,j).variance = sdrbfKern.variance;
        rbfKern(i,j).inverseWidth = sdrbfKern.inverseWidth(j,i);
        rbfKern(i,j).isNormalised = sdrbfKern.isNormalised;
    end
    newt1 = t1(t1> spVector(i) & t1<spVector(i+1));
    newt2 = t2(t2> spVector(i) & t2<spVector(i+1));
    dim1(i) = length(newt1);
    dim2(i) = length(newt2);
end

gradIW = zeros(sdrbfKern.nlfPerInt, sdrbfKern.nIntervals);
gradSP = zeros(1, sdrbfKern.nIntervals);
%dimParam = sdrbfKern.nIntervals*ones(1, 2);
%gC = cell(1, sdrbfKern.nlfPerInt);
%indxIW = 1:dimParam(1);
%indxSP = dimParam(1)+1:sum(dimParam);

covGrad2 = covGrad;
for q=1:sdrbfKern.nlfPerInt
    startValOne = 1;
    endValOne   = 0;
    startValThree = 1;
    endValThree = 0;
    gradSPTemp = zeros(1, sdrbfKern.nIntervals);
    if iscell(covGrad2)
        covGrad = covGrad2{q};
    end
    for i=1:sdrbfKern.nIntervals
        endValOne = endValOne + dim1(i);
        endValThree = endValThree + dim2(i);
        gP = rbfKernGradient(rbfKern(i,q), t1(startValOne:endValOne) - spVector(i), ...
            t2(startValThree:endValThree) - spVector(i), ...
            covGrad(startValOne:endValOne, startValThree:endValThree));
        gradIW(q,i)= gP(1);
        gX = rbfKernGradX(rbfKern(i,q), t1(startValOne:endValOne) - spVector(i), ...
            t2(startValThree:endValThree) - spVector(i));
        gXR = reshape(gX, dim1(i), dim2(i))';
        gradSPTemp(i) = gradSPTemp(i) - 2*sum(sum(gXR.*covGrad(startValOne:endValOne, startValThree:endValThree)));
        startValThree = endValThree + 1;
        startValOne = endValOne + 1;
    end
    gradSP = gradSP + gradSPTemp;
    %gC{q}(indxIW) = gradIW(q,:);
    %gC{q}(indxSP) = gradSPTemp;
end

%g = [gradIW(:)' gradSP];

%%%%%%
tempGradSP = fliplr(gradSP);
tempGradSP = cumsum(tempGradSP);
g = [gradIW(:)' fliplr(tempGradSP)];
%%%%%%


text:
function [M,P] = ckf_predict(M,P,f,Q,f_param)
% CKF_PREDICT - Cubature Kalman filter prediction step
%
% Syntax:
%   [M,P] = CKF_PREDICT(M,P,[f,Q,f_param])
%
% In:
%   M - Nx1 mean state estimate of previous step
%   P - NxN state covariance of previous step
%   f - Dynamic model function as a matrix A defining
%       linear function f(x) = A*x, inline function,
%       function handle or name of function in
%       form f(x,param)                   (optional, default eye())
%   Q - Process noise of discrete model   (optional, default zero)
%   f_param - Parameters of f               (optional, default empty)
%
% Out:
%   M - Updated state mean
%   P - Updated state covariance
%
% Description:
%   Perform additive form spherical-radial cubature Kalman filter (CKF)
%   prediction step.
%
%   Function f(.) should be such that it can be given a
%   DxN matrix of N sigma Dx1 points and it returns
%   the corresponding predictions for each sigma
%   point.
%
% See also:
%   CKF_UPDATE, CRTS_SMOOTH, CKF_TRANSFORM, SPHERICALRADIAL
%
% References:
%   Arasaratnam and Haykin (2009). Cubature Kalman Filters.
%    IEEE Transactions on Automatic Control, vol. 54, no. 5, pp.1254-1269

% Copyright (c) 2010 Arno Solin
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
%%

  %
  % Check which arguments are there
  %
  if nargin < 2
    error('Too few arguments');
  end
  if nargin < 3
    f = [];
  end
  if nargin < 4
    Q = [];
  end

  %
  % Apply defaults
  %
  if isempty(f)
    f = eye(size(M,1));
  end
  if isempty(Q)
    Q = zeros(size(M,1));
  end

  %
  % Do transform and add process noise
  %
  if nargin < 5
    [M,P] = ckf_transform(M,P,f);
  else
    [M,P] = ckf_transform(M,P,f,f_param);
  end
  P = P + Q;

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
% magnitude of velocity of tail
function [data,units] = compute_velmag_tail(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  % location of tail
  tailx = trx(fly).x_mm + 2*cos(-trx(fly).theta).*trx(fly).a_mm;
  taily = trx(fly).y_mm + 2*sin(-trx(fly).theta).*trx(fly).a_mm;
  dx = diff(tailx);
  dy = diff(taily);

  if trx(fly).nframes < 2,
    data{i} = [];
  else
    % magnitude of velocity vector
    data{i} = sqrt(dx.^2 + dy.^2)./trx(fly).dt;
  end
end
units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
function coef = sgmga_vcn_coef_naive ( dim_num, level_weight, x_max, x, q_max )

%*****************************************************************************80
%
%% SGMGA_VCN_COEF_NAIVE returns the "next" constrained vector's coefficient.
%
%  Discussion:
%
%    This function uses a naive approach to the computation, resulting in
%    a set of 2^DIM_NUM tasks.  Hence it is not suitable for cases where
%    DIM_NUM is moderately large.  The function SGMGA_VCN_COEF carries out
%    a more complicated but more efficient algorithm for the same computation.
%
%    We are given a nonnegative vector X of dimension DIM_NUM which satisfies:
%
%      sum ( 1 <= I <= DIM_NUM ) LEVEL_WEIGHT(I) * X(I) <= Q_MAX
%
%    This routine computes the appropriate coefficient for X in the
%    anisotropic sparse grid scheme.
%
%    The coefficient is calculated as follows:
%
%      Let B be a binary vector of length N, and let ||B|| represent
%      the sum of the entries of B.
%
%      Coef = sum ( all B such that X+B satisfies constraints ) (-1)^||B||
%
%    Since X+0 satisfies the constraint, there is always at least one
%    summand.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    An Anisotropic Sparse Grid Stochastic Collocation Method for Partial
%    Differential Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2411-2442.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the vector.
%
%    Input, real LEVEL_WEIGHT(DIM_NUM), the anisotropic weights.
%
%    Input, integer X_MAX(DIM_NUM), the maximum
%    values allowed in each component.
%
%    Input, integer X(DIM_NUM), a point which satisifies the constraints.
%
%    Input, real Q_MAX, the upper limit on the sum.
%
%    Output, real COEF, the combinatorial coefficient.
%

%
%  Force column vectors.
%
  level_weight = level_weight ( : );
  x = x ( : );

  b = zeros ( dim_num, 1 );
  coef = 1.0;

  while ( 1 )
%
%  Generate the next binary perturbation.
%
    b = binary_vector_next ( dim_num, b );
    b_sum = sum ( b(1:dim_num) );
%
%  We're done if we've got back to 0.
%
    if ( b_sum == 0 )
      break
    end
%
%  Perturb the vector.
%
    x2(1:dim_num,1) = x(1:dim_num,1) + b(1:dim_num,1);
%
%  Does it satisfy the XMAX constraint?
%  (THIS CHECK IS SURPRISINGLY NECESSARY, PARTLY BECAUSE OF ZERO WEIGHT).
%
    if ( any ( x_max(1:dim_num) < x2(1:dim_num) ) )
      continue
    end
%
%  Does it satisfy the Q_MAX constraint?
%
    q = level_weight(1:dim_num)' * x2(1:dim_num);

    if ( q <= q_max )
      coef = coef + r8_mop ( b_sum );
    end

  end

  return
end


text:
function [FS,read] = spm_voice_FS(wfile)
% Sampling frequency and function handle for handling sound signals
% FORMAT [FS,read] = spm_voice_FS(wfile)
%
% wfile  - .wav file, audio object or (double) timeseries
%
% FS     - sampling frequency
% read   - function handle: Y = read(wfile);
%
%  This auxiliary routine finds the sampling frequency and returns a
%  function handle appropriate for the sound format in question.
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_FS.m 7750 2019-12-05 17:54:29Z spm $

% get timeseries from audio recorder(or from a file)
%--------------------------------------------------------------------------
global VOX

% default
%--------------------------------------------------------------------------
if ~nargin
    try
        FS = get(VOX.audio,'SampleRate');
    catch
        FS = 22050;
    end
    return
end

% get source (recorder) and FS
%--------------------------------------------------------------------------
if isa(wfile,'audiorecorder')

    FS     = get(wfile,'SampleRate');
    read   = @getaudiodata;

elseif isnumeric(wfile)

    % timeseries
    %----------------------------------------------------------------------
    try
        FS = get(VOX.audio,'SampleRate');
    catch
        try
            FS = VOX.FS;
        catch
            FS = 22050;
        end
    end
    read   = @(Y)Y;

else

    % sound file
    %----------------------------------------------------------------------
    try
        xI     = audioinfo(wfile);
        FS     = xI.SampleRate;
        read   = @audioread;
    catch
        [~,FS] = wavread(wfile,[1 1]);
        read   = @wavread;
    end

end

% place sampling frequency in global VOX structure
%----------------------------------------------------------------------
% VOX.FS = FS;

author: spm


text:
%%********************************************************************
%% This is an example for the interface
%% This example is to solve the following SDP
%%
%% min trace(X1)+trace(X2)+sum(X3)
%% s.t. -X1(1,2)+2X2(3,3)+2X3(2)=4
%%       2X1(2,3)+X2(4,2)-X3(4) = 3
%%       2<= -X1(1,2)-2X2(3,3)+2X3(2)<=7
%%       X1 \in S^{6}_{+}, X2 \in R^{5*5}, X3 \in R^{7}_{+}
%%       0<=X1<=10, 0<=X2<=8
%% SDPNAL+:
%% Copyright (c) 2017 by
%% Yancheng Yuan, Kim-Chuan Toh, Defeng Sun and Xinyuan Zhao
%%********************************************************************
clear all;
n1 = 6; n2 = 5; n3 = 7;
%% Initialize model
model = ccp_model('Example_simple');
    X1 = var_sdp(n1,n1);
    X2 = var_nn(n2, n2);
    X3 = var_nn(n3);
    model.add_variable(X1, X2, X3);
    model.minimize(trace(X1)+trace(X2)+sum(X3));
    options=3;
    if (options==1)
       model.add_affine_constraint(-X1(1,2)+2*X2(3,3)+2*X3(2)==4);
       model.add_affine_constraint(2*X1(2,3)+X2(4,2)-X3(4) == 3);
    else
       %%two other ways to add the above equality constraints
       A1 = {sparse(n1,n1), sparse(n1,n1)};
       A2 = {sparse(n2,n2), sparse(n2,n2)}; A3 = sparse(2,n3);
       A1{1}(1,2)= -1;  A1{2}(2,3)= 2;
       A2{1}(3,3)= 2;   A2{2}(4,2)= 1;
       A3(1,2) = 2;     A3(2,4) = -1;
       b = [4;3];
       if (options==2)
          model.add_affine_constraint(A1*X1+A2*X2+A3*X3 == b);
       elseif (options==3)
          blk{1,1} = 's'; blk{1,2} = n1;
          %A1vec = [svec(blk,A1{1}),svec(blk,A1{2})]';
          A1vec = [A1{1}(:), A1{2}(:)]';
          %model.add_affine_constraint(A1vec*map_svec(X1)+A2*X2+A3*X3 == b);
          model.add_affine_constraint(A1vec*map_vec(X1)+A2*X2+A3*X3 == b);
       end
    end
    model.add_affine_constraint(2<=-X1(1,2)-2*X2(3,3)+2*X3(2)<=7);
    model.add_affine_constraint(0<= X1 <=10);
    model.add_affine_constraint(X2<=8);
    model.add_affine_constraint(map_diag(X1)==ones(n1,1));
    model.setparameter('BBTsolve.method','direct');
model.solve;
%%*********************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC