text:
% POP_SPECPARAMS - Set plotting and statistics parameters for computing
%                    STUDY component spectra.
% Usage:
%       >> STUDY = pop_specparams(STUDY, 'key', 'val');
%
% Inputs:
%   STUDY        - EEGLAB STUDY set
%
% Plot options:
%   'topofreq'   - [real] Plot Spectrum scalp maps at one specific freq. (Hz).
%                  A frequency range [min max] may also be defined (the
%                  spectrum is then averaged over the interval) {default: []}
%   'freqrange'  - [min max] spectral frequency range (in Hz) to plot.
%                  {default: whole frequency range} .
%   'ylim'       - [mindB maxdB] spectral plotting limits in dB
%                  {default: from data}
%   'plotgroups' - ['together'|'apart'] 'together' -> plot subject groups
%                  on the same figure in different colors, else ('apart') on
%                  different figures {default: 'apart'}
%   'plotconditions' - ['together'|'apart'] 'together' -> plot conditions
%                  on the same figure in different colors, else ('apart')
%                  on different figures. Note: keywords 'plotgroups' and
%                  'plotconditions' cannot both be set to 'together'.
%                  {default: 'apart'}
%   'subtractsubjectmean' - ['on'|'off'] subtract individual subject mean
%                  from each spectrum before plotting and computing
%                  statistics. Default is 'off'.
%   'averagechan' - ['rms'|'on'|'off'] average data channels when several are
%                  selected ('on') or compute root mean square ('rms').
%
% See also: STD_SPECPLOT
%
% Authors: Arnaud Delorme, CERCO, CNRS, 2006-

% Copyright (C) Arnaud Delorme, CERCO
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [ STUDY, com ] = pop_specparams(STUDY, varargin);

STUDY = default_params(STUDY);
TMPSTUDY = STUDY;
com = '';
if isempty(varargin)

    enablecond  = 'off';
    enablegroup = 'off';
    if length(STUDY.design(STUDY.currentdesign).variable) > 0 && length(STUDY.design(STUDY.currentdesign).variable(1).value)>1, enablecond  = 'on'; end
    if length(STUDY.design(STUDY.currentdesign).variable) > 1 && length(STUDY.design(STUDY.currentdesign).variable(2).value)>1, enablegroup = 'on'; end;
    plotconditions     = fastif(strcmpi(STUDY.etc.specparams.plotconditions, 'together'), 1, 0);
    plotgroups         = fastif(strcmpi(STUDY.etc.specparams.plotgroups,'together'), 1, 0);
    submean            = fastif(strcmpi(STUDY.etc.specparams.subtractsubjectmean,'on'), 1, 0);
    tmpFreqRange = STUDY.etc.specparams.freqrange;
    if strcmpi(STUDY.etc.specparams.averagechan,'off')
        if isempty(STUDY.etc.specparams.topofreq) || any(isnan(STUDY.etc.specparams.topofreq))
            multipleChansVal   = 1; % scalp array
        else
            multipleChansVal   = 2; % scalp topo
            tmpFreqRange = STUDY.etc.specparams.topofreq;
        end
    else
        if strcmpi(STUDY.etc.specparams.averagechan,'on')
            multipleChansVal   = 3; % average channels
        else
            multipleChansVal   = 4; % root mean square
        end
    end

    cb_radio = [ 'set(findobj(gcbf, ''userdata'', ''radio''), ''value'', 0);' ...
                 'set(gcbo, ''value'', 1);' ...
                 'set(findobj(gcbf, ''tag'', ''topofreq''), ''string'', '''');' ];
    cb_edit  = [ 'set(findobj(gcbf, ''userdata'', ''radio''), ''value'', 0);' ...
                 'set(findobj(gcbf, ''tag'', ''scalptopotext''), ''value'', 1);' ];
    cb_multiplechan    = [ 'if get(gcbo, ''value'') == 2 ' ...
                         '    if isempty(get(findobj(gcbf, ''tag'', ''freqrange''), ''string'')),' ...
                         '       set(gcbo, ''value'', 1);' ...
                         '       warndlg2([ ''Select frequency range first for plotting topography,'' 10 ''then select that setting again.'' ]);' ...
                         '    end;' ...
                         'else,' ...
                         '    if ~isempty(get(findobj(gcbf, ''tag'', ''freqrange''), ''string'')) && length(unique(str2num(get(findobj(gcbf, ''tag'', ''freqrange''), ''string'')))) ==1,' ...
                         '       set(findobj(gcbf, ''tag'', ''freqrange''), ''string'', '''');' ...
                         '    end;' ...
                         'end;' ];

    uilist = { ...
           {'style' 'text'       'string' 'Spectrum plotting options' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'text'       'string' 'Frequency [low_Hz high_Hz]' } ...
           {'style' 'edit'       'string' num2str(tmpFreqRange) 'tag' 'freqrange' } ...
        {} {'style' 'text'       'string' 'Plot limits [low high]'} ...
           {'style' 'edit'       'string' num2str(STUDY.etc.specparams.ylim) 'tag' 'ylim' } ...
        {} {'style' 'checkbox'   'string' 'Subtract individual subject mean spectrum' 'value' submean 'tag' 'submean' } ...
        {} ...
           {'style' 'text'       'string' 'Spectrum plotting format' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'checkbox'   'string' 'Plot first variable on the same panel' 'value' plotconditions 'enable' enablecond  'tag' 'plotconditions' } ...
        {} {'style' 'checkbox'   'string' 'Plot second variable on the same panel' 'value' plotgroups 'enable' enablegroup 'tag' 'plotgroups' } ...
        {} ...
           {'style' 'text'       'string' 'Multiple channels selection' 'fontweight' 'bold' 'tag', 'spec' 'fontsize', 12} ...
        {} {'style' 'popupmenu'  'string' { 'Plot channels individually' 'Plot averaged topography over frequency range' 'Average power of selected channels' 'Compute RMS power of selected channels' } 'value' multipleChansVal 'tag' 'multiplechan' 'callback' cb_multiplechan } };
    cbline = [0.07 1.1];
    otherline = [ 0.07 0.6 .3];
    chanline  = [ 0.07 0.8];
    geometry = { 1 otherline otherline cbline 1 1 cbline cbline 1 1 chanline };
    geomvert = [1.2 1 1 1 0.5 1.2 1 1 0.5 1.2 1 ];

    % component plotting
    % ------------------
    if isnan(STUDY.etc.specparams.topofreq)
        geometry(end-2:end) = [];
        geomvert(end-2:end) = [];
        uilist(end-3:end) = [];
    end

    [out_param userdat tmp res] = inputgui( 'geometry' , geometry, 'uilist', uilist, 'geomvert', geomvert, ...
                                            'title', 'Spectrum plotting options -- pop_specparams()');
    if isempty(res), return; end
    if ~isfield(res, 'multiplechan'), res.multiplechan = 0; end

    % decode inputs
    % -------------
    %if res.plotgroups && res.plotconditions, warndlg2('Both conditions and group cannot be plotted on the same panel'); return; end
    if res.submean   , res.submean    = 'on'; else res.submean    = 'off'; end
    if res.plotgroups, res.plotgroups = 'together'; else res.plotgroups = 'apart'; end
    if res.plotconditions , res.plotconditions  = 'together'; else res.plotconditions  = 'apart'; end
    res.freqrange = str2num( res.freqrange );
    res.ylim      = str2num( res.ylim );

    % build command call
    % ------------------
    options = {};
    if ~strcmpi( res.plotgroups, STUDY.etc.specparams.plotgroups), options = { options{:} 'plotgroups' res.plotgroups }; end
    if ~strcmpi( res.plotconditions , STUDY.etc.specparams.plotconditions ), options = { options{:} 'plotconditions'  res.plotconditions  }; end
    if ~strcmpi( res.submean   , STUDY.etc.specparams.subtractsubjectmean ), options = { options{:} 'subtractsubjectmean'  res.submean  }; end
    if ~isequal(res.ylim, STUDY.etc.specparams.ylim),               options = { options{:} 'ylim' res.ylim      }; end
    if ~isequal(res.freqrange, STUDY.etc.specparams.freqrange) && res.multiplechan ~= 2,     options = { options{:} 'freqrange' res.freqrange }; end

    % multiple channel option
    % -----------------------
    if res.multiplechan == 1
        if ~isequal('off', STUDY.etc.specparams.averagechan), options = { options{:} 'averagechan' 'off' }; end
        if ~isempty(       STUDY.etc.specparams.topofreq),    options = { options{:} 'topofreq' [] }; end
    elseif res.multiplechan == 2
        if ~isequal('off', STUDY.etc.specparams.averagechan), options = { options{:} 'averagechan' 'off' }; end
        if ~isequal(res.freqrange, STUDY.etc.specparams.topofreq) options = { options{:} 'topofreq' res.freqrange }; end
        if ~isequal([], STUDY.etc.specparams.freqrange) options = { options{:} 'freqrange' [] }; end
        if isempty(res.freqrange)
            disp('Warning: you must select a frequency range to plot scalp topographies, plotting individual channels instead');
        end
    elseif res.multiplechan > 2
        if ~isempty(       STUDY.etc.specparams.topofreq),    options = { options{:} 'topofreq' [] }; end
        if res.multiplechan == 3
            if ~isequal('on', STUDY.etc.specparams.averagechan), options = { options{:} 'averagechan' 'on' }; end
        else
            if ~isequal('rms', STUDY.etc.specparams.averagechan), options = { options{:} 'averagechan' 'rms' }; end
        end
    end

    % execute option
    % --------------
    if ~isempty(options)
        STUDY = pop_specparams(STUDY, options{:});
        com = sprintf('STUDY = pop_specparams(STUDY, %s);', vararg2str( options ));
    end
else

    if strcmpi(varargin{1}, 'default')
        STUDY = default_params(STUDY);
    else
        for index = 1:2:length(varargin)
            if ~isempty(strmatch(varargin{index}, fieldnames(STUDY.etc.specparams), 'exact'))
                STUDY.etc.specparams = setfield(STUDY.etc.specparams, varargin{index}, varargin{index+1});
            end
        end
    end
end

% scan clusters and channels to remove specdata info if freqrange has changed
% ----------------------------------------------------------
if ~isequal(STUDY.etc.specparams.freqrange, TMPSTUDY.etc.specparams.freqrange) || ...
    ~isequal(STUDY.etc.specparams.subtractsubjectmean, TMPSTUDY.etc.specparams.subtractsubjectmean)
    rmfields = { 'specdata' 'specfreqs' };
    for iField = 1:length(rmfields)
        if isfield(STUDY.cluster, rmfields{iField})
            STUDY.cluster = rmfield(STUDY.cluster, rmfields{iField});
        end
        if isfield(STUDY.changrp, rmfields{iField})
            STUDY.changrp = rmfield(STUDY.changrp, rmfields{iField});
        end
    end
end

function STUDY = default_params(STUDY)
    if ~isfield(STUDY.etc, 'specparams'), STUDY.etc.specparams = []; end
    if ~isfield(STUDY.etc.specparams, 'topofreq'),             STUDY.etc.specparams.topofreq = []; end
    if ~isfield(STUDY.etc.specparams, 'freqrange'),            STUDY.etc.specparams.freqrange = []; end
    if ~isfield(STUDY.etc.specparams, 'ylim'     ),            STUDY.etc.specparams.ylim      = []; end
    if ~isfield(STUDY.etc.specparams, 'subtractsubjectmean' ), STUDY.etc.specparams.subtractsubjectmean  = 'off'; end
    if ~isfield(STUDY.etc.specparams, 'plotgroups'),           STUDY.etc.specparams.plotgroups = 'apart'; end
    if ~isfield(STUDY.etc.specparams, 'plotconditions'),       STUDY.etc.specparams.plotconditions  = 'apart'; end
    if ~isfield(STUDY.etc.specparams, 'averagechan') ,         STUDY.etc.specparams.averagechan  = 'off'; end
    if ~isfield(STUDY.etc.specparams, 'detachplots') ,         STUDY.etc.specparams.detachplots  = 'on'; end % deprecated

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% example_radar.m
% ====================================================>
% The following seem to be descent parameters:
%   - ProbOfDetection = 0.8
%   - ProbOfDeath = 0.2
%   - VelocityErrVariance = 2
%   - ObsErrVariance = 50
%   - ProbOfConfirm = 0.9
%   - PHD.BirthScheme = {'Expansion', 5000}
%

%% Radar specific settings
RadarName = 'staddon';
switch(RadarName)
    case('staddon')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    -212.289393440959;... % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
%         RadarCoords.lat = 50.33933333;
%         RadarCoords.lon = -4.12527778;
        RadarCoords.lat = 50.346069;
        RadarCoords.lon = -4.113670;
    case('longroom')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    3000;...              % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
        RadarCoords.lat = 50.36286670714617;
        RadarCoords.lon = -4.156833300366998;
end

% Load dataset
load(strcat(RadarName,'.mat'));
N = size(DataList,2); % Simulation length

%% Focus on region of high clutter
% lon_lim = [-4.195 -4.11];
% lat_lim = [50.31 50.375];
% lon_lim = [-4.15 -4.13];
% lat_lim = [50.335 50.35];
% lon_lim = [-4.1538   -4.1152];
% lat_lim = [50.3097   50.35];

lon_lim = [-4.1538   -4.1152];
lat_lim = [50.3097   50.3373];

% lon_lim = [-4.173972 -4.145608];
% lat_lim = [50.351393 50.362318];

% Convert LLA limits to NED
[y_lim, x_lim] = geodetic2ned(lat_lim,...
                              lon_lim, ...
                              0,...
                              RadarCoords.lat,...
                              RadarCoords.lon,...
                              0,...
                              referenceEllipsoid('wgs84'));
V_bounds = [x_lim(1);... % X-min |
            x_lim(2);... % X-max | Surveillance region bounding
            y_lim(1);... % Y-min | box coordinates (m)
            y_lim(2)]';  % Y-max |
for i = 1:N
    j = 1;
    while j <= size(DataList{i},2)
        if(DataList{i}(3,j)>lat_lim(2) || DataList{i}(3,j)<lat_lim(1) || DataList{i}(4,j)<lon_lim(1) || DataList{i}(4,j)>lon_lim(2))
            DataList{i}(:,j) = [];
            j = j -1;
        end
        j = j + 1;
    end
end

radar_meas_convert;
%% Plot & Recording settings
% Plot settings
ShowPlots = 1;              % Set to 0 to prevent showing any plots
ShowUpdate = 1;             % Set to 0 to skip showing update plots
ShowTrackInfo = 0;
NumPersistFrames = 10000;
PlotLatLon = 1;

% Recording settings
Record = 1;                 % Set to (0|1) to turn video recording (off|on)
FrameRate = 0.5;            % Number of frames per second
VideoQuality = 100;         % Set to desired quality percentage
VideoPathName = strcat(RadarName,'_heavy_clutter_staddon.avi'); % Set to the desired path and name of produced recording

% Model parameter shortcuts
lambdaV = 2; % Expected number of clutter measurements over entire surveillance region
V = (abs(V_bounds(2)-V_bounds(1))*abs(V_bounds(4)-V_bounds(3))); % Total area of surveillance region
P_D = 0.5;    % Probability of detection
timestep_duration = duration(0,0,2);

%% Instantiation of necessary components

% Instantiate a Dynamic model (CV model with q = 1 m/s^2)
% transition_model = ConstantVelocityX('VelocityErrVariance', 0.1,...
%                                      'NumDims', 2,...
%                                      'TimestepDuration', timestep_duration);
transition_model = ConstantHeadingX('VelocityErrVariance',1,...
                      'HeadingErrVariance',0.03,...
                      'TimestepDuration',timestep_duration);
% Instantiate an Observation model (Variance of 50m^2 on each coordinate)
measurement_model = LinearGaussianX('NumMeasDims', 2,...
                                    'NumStateDims', 4,...
                                    'MeasurementErrVariance', 10^2,...
                                    'Mapping', [1 2]);
clutter_model = PoissonRateUniformPosition2X('ClutterRate',lambdaV,...
                                            'Limits',[V_bounds(1:2);...
                                                      V_bounds(3:4)]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                V_bounds(3:4);...
                                                                                [-7 7 ];...
                                                                                [0 2*pi ]]),...
                                           'BirthIntensity', 1);
% Instantiate birth model
numBirthComponents = 9;
%numBirthComponents = round(sqrt(numBirthComponents))^2; % Adjust to the nearest square-rootable number

x_birth = linspace(V_bounds(1), V_bounds(2), sqrt(numBirthComponents));
y_birth = linspace(V_bounds(4), V_bounds(3), sqrt(numBirthComponents));
[X,Y] = meshgrid(x_birth, y_birth);
X = reshape(X, [1,numel(X)]);
Y = reshape(Y, [1,numel(Y)]);
BirthComponents.Means = [X;
                         mvnrnd(zeros(numel(X),1), 0.5)';
                         Y;
                         mvnrnd(zeros(numel(Y),1), 0.5)'];
obs_covar= measurement_model.covar();
BirthComponents.Covars = repmat(50*transition_model.covar() + 10*blkdiag(obs_covar(1,1), obs_covar(2,2), 0, 0),1,1,numBirthComponents);
BirthComponents.Weights = 0.01*ones(1,numBirthComponents);
birth_distribution = GaussianMixtureX(BirthComponents.Means,BirthComponents.Covars, BirthComponents.Weights);
birth_model2 = DistributionBasedBirthModelX('Distribution', birth_distribution,...
                                           'BirthIntensity', 1);
% Compile the State-Space model
model = StateSpaceModelX(transition_model, measurement_model,...
                        'Clutter',clutter_model,...
                        'Detection', detection_model,...
                        'Birth', birth_model);
model2 = StateSpaceModelX(transition_model, measurement_model,...
                        'Clutter',clutter_model,...
                        'Detection', detection_model,...
                        'Birth', birth_model2);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), obs_covar(2,2), 0, 0));
base_filter = ExtendedKalmanFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionModel = detection_model;
jpdaf = JointIntegratedProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_phd.Model = model;
[priorParticles, priorWeights] = model.Birth.random(50000);
config_phd.StatePrior = ParticleStateX(priorParticles,10*priorWeights);
config_phd.BirthScheme = {'Expansion', 5000};
config_phd.SurvivalProbability = 0.995;

config_phd2.Model = model2;
[priorParticles, priorWeights] = model2.Birth.random(50000);
config_phd2.Filter = KalmanFilterX('Model',model2);
config_phd2.StatePrior = GaussianMixtureStateX(birth_distribution);
config_phd2.StatePrior.Weights = config_phd2.StatePrior.Weights*100;
config_phd2.SurvivalProbability = 0.99;
config_phd2.DetectionProbability = P_D;

% Instantiate PHD filter
myphd = SMC_PHDFilterX(config_phd);
%myphd = GM_PHDFilterX(config_phd2);
% Initiate Tag Generator
tag_gen = UuidTagGeneratorX();

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.PhdFilter = myphd;
config_ti.ConfirmThreshold = 0.8;

% Create the track initiator
myti = PhdTrackInitiatorX(config_ti);

%% Track Deleter
config_td = struct('Fieldname', 'ExistenceProbability',...
                   'ReferenceValue', 0.1,...
                   'ReferenceOperand', 'lt');
mytd = FieldBasedDeleterX(config_td);

TrackList = [];

%% Create plot windows
if(ShowPlots)

    % Map plot
    figure('units','normalized','outerposition',[0 0 .5 0.9])
    ax(1) = gca;
    if PlotLatLon
        plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
        axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
    end
%     axis(ax(1),[-4.195 -4.11 50.31 50.375])

%     % PHD Intensity plot
%     figure('units','normalized','outerposition',[.5 0 .5 1])
%     ax(2) = gca;

    plots = [];
end

%% START OF SIMULATION
% ===================>

for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = MeasurementScans(k);
    timestamp_km1 = MeasurementScans(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate(TrackList, MeasurementList);
    jpdaf.updateTracks();

    % Delete all plots (other than map)
    for i = 1:numel(plots)
        delete(plots(i))
    end
    plots = [];
    hold on;

%     if(MeasurementList.NumMeasurements>0)
%         % Convert measurements to LLA and plot them
%         [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
%                                    MeasurementList.Vectors(1,:),...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%         plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
%
%         for i=1:numel(lat)
%             plots(end+1) = text(ax(1), lon(:,i)+0.0001,lat(:,i)+0.00027,num2str(i),'FontSize',8,'Color','w');
%         end
% %             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
% %                                '-s','MarkerSize',20,...
% %                                'MarkerEdgeColor','red',...
% %                                'MarkerFaceColor',[1 .6 .6]);
%     end
%     pause(0.01);
     %% Perform Track initiation
    [TrackList] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix);

    %% Perform Track deletion
    TrackList = mytd.deleteTracks(TrackList);

    % Plot update step results
    if(ShowPlots && ShowUpdate)

        % Delete all plots (other than map)
        for i = 1:numel(plots)
            delete(plots(i))
        end
        plots = [];
        hold on;

        if(MeasurementList.NumMeasurements>0)
            if PlotLatLon
                % Convert measurements to LLA and plot them
                [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
                                           MeasurementList.Vectors(1,:),...
                                           0,...
                                           RadarCoords.lat,...
                                           RadarCoords.lon,...
                                           0,...
                                           referenceEllipsoid('wgs84'));
                plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
            else
                plots(end+1) = plot(ax(1), MeasurementList.Vectors(1,:),MeasurementList.Vectors(2,:),'b*','MarkerSize', 10);
            end
%             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
%                                '-s','MarkerSize',20,...
%                                'MarkerEdgeColor','red',...
%                                'MarkerFaceColor',[1 .6 .6]);

        end
        if PlotLatLon
%             x = [-3560.53119790131, -3560.53119790131,  -3080.76025084906,  -3080.76025084906,  -3367.20269177091,  -3367.20269177091,  -3560.53119790131];
%             y = [1349.54778996881,  1590.19881386064,   1590.19881386064,   1108.03444256592,   1108.03444256592,   1349.54778996881,   1349.54778996881];
%             [y,x,~] = ned2geodetic(y,...
%                                    x,...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%             plots(end+1) = plot(ax(1), x, y, 'm-');
%             x = [-3367.20269177091,-3367.00101398088,-3083.23234936543,-3083.41702988195,-3367.20269177091];
%             y = [1108.20727096473,1425.35687149648,1425.18405344934,1108.03444256592,1108.20727096473];
%             [y,x,~] = ned2geodetic(y,...
%                                    x,...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%             plots(end+1) = plot(ax(1), x, y, 'm-');
            x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
            y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
            [y,x,~] = ned2geodetic(y,...
                                   x,...
                                   0,...
                                   RadarCoords.lat,...
                                   RadarCoords.lon,...
                                   0,...
                                   referenceEllipsoid('wgs84'));
            plots(end+1) = plot(ax(1), x, y, 'r-');
        else
%             x = [-3560.53119790131, -3560.53119790131,  -3080.76025084906,  -3080.76025084906,  -3367.20269177091,  -3367.20269177091,  -3560.53119790131];
%             y = [1349.54778996881,  1590.19881386064,   1590.19881386064,   1108.03444256592,   1108.03444256592,   1349.54778996881,   1349.54778996881];
%             plots(end+1) = plot(ax(1), x, y, 'm-');
%             x = [-3367.20269177091,-3367.00101398088,-3083.23234936543,-3083.41702988195,-3367.20269177091];
%             y = [1108.20727096473,1425.35687149648,1425.18405344934,1108.03444256592,1108.20727096473];
%             plots(end+1) = plot(ax(1), x, y, 'm-');
            x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
            y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
            plots(end+1) = plot(ax(1), x, y, 'r-');
        end
        % Plot all existing tracks
        for j=1:numel(TrackList)
            track = TrackList{j};
            % Convert track trajectory to LLA and plot it
            means = [track.Trajectory.Mean];
            if PlotLatLon
                [lat,lon,~] = ned2geodetic(means(2,:),...
                                           means(1,:),...
                                           0,...
                                           RadarCoords.lat,...
                                           RadarCoords.lon,...
                                           0,...
                                           referenceEllipsoid('wgs84'));
                traj_length = size(lon,2);
                if(traj_length>NumPersistFrames)
                    start = traj_length-NumPersistFrames;
                else
                    start = 1;
                end
                plots(end+1) = plot(ax(1), lon(:,start:end),lat(:,start:end),'-.w','LineWidth',2);
                plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'ws','MarkerSize',15);

                % Convert track velocity to LLA and plot it
                [lat_vel,lon_vel,~] = ned2geodetic(track.State.Mean(3,end)*sin(track.State.Mean(4,end)),...
                                                   track.State.Mean(3,end)*cos(track.State.Mean(4,end)),...
                                                   0,...
                                                   lat(:,end),...
                                                   lon(:,end),...
                                                   0,...
                                                   referenceEllipsoid('wgs84'));
                lat_vel = lat_vel-lat(:,end);
                lon_vel = lon_vel-lon(:,end);
                plots(end+1) = quiver(ax(1), lon(:,end),lat(:,end),20*lon_vel,20*lat_vel,'r','LineWidth',1.5);
                if ShowTrackInfo
                    speed_kmph = sqrt(track.State.Mean(4,end)^2+track.State.Mean(2,end)^2)*3.6;
                    speed_knot = speed_kmph/1.852;
                    plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)+0.00027,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','w');
                    plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)-0.00027,strcat("PoE:",num2str(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','w');
                end
            else
                traj_length = size(lon,2);
                if(traj_length>NumPersistFrames)
                    start = traj_length-NumPersistFrames;
                else
                    start = 1;
                end
                plots(end+1) = plot(ax(1), means(1,start:end),means(2,start:end),'-.k','LineWidth',2);
                plots(end+1) = plot(ax(1), means(1,end),means(2,end),'ks','MarkerSize',15);
%                 lat_vel = lat_vel-lat(:,end);
%                 lon_vel = lon_vel-lon(:,end);
%
                x_vel = track.State.Mean(3,end)*cos(track.State.Mean(4,end));
                y_vel = track.State.Mean(3,end)*sin(track.State.Mean(4,end));
                plots(end+1) = quiver(ax(1), means(1,end),means(2,end),20*x_vel,20*y_vel,'r','LineWidth',1.5);

                if ShowTrackInfo
                    speed_kmph = track.State.Mean(3,end)*3.6;
                    speed_knot = speed_kmph/1.852;
                    plots(end+1) = text(ax(1), means(1,end)+60,means(2,end)+50,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','k');
                    plots(end+1) = text(ax(1), means(1,end)+60,means(2,end)-50,strcat("PoE:",num2str(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','k');
                end
            end

            % TODO: Convert heading and state covariance to LLA and plot them
            % [lat,lon,h] = ned2geodetic(North,East,0,50.346069,-4.113670,0,referenceEllipsoid('wgs84'));
            % h2 = plot_gaussian_ellipsoid(TrackList{j}.Filter.StateMean([1 3]), TrackList{j}.Filter.StateCovar([1 3],[1 3]),1,20,ax(1));
            % plots(end+1) = h2;
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-5,int2str(TrackList{j}.TrackID));
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-130,num2str(TrackList{j}.ProbOfExist,2));
        end

        % Add axis labels
%         xlabel('Longitude')
%         ylabel('Latitude')
set(gca,'visible','off')
        pause(0.0001)
        % Store video frame
        if(Record)
            F(k) = getframe(ax(1));
        end
    end

end

% Create video file and write to it
if(Record)
    F = F(2:end);
    vidObj = VideoWriter(char(VideoPathName));
    vidObj.Quality = VideoQuality;
    vidObj.FrameRate = FrameRate;
    open(vidObj);
    writeVideo(vidObj, F);
    close(vidObj);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function SRhdplot(HD,VARopt)
% =======================================================================
% Plot the HD shocks computed with SR (sign restriction procedure)
% =======================================================================
% SRhdplot(HD,VARopt)
% -----------------------------------------------------------------------
% INPUT
%   - HD: structure from SR
%   - VARopt: options of the VAR (from VARmodel and SR)
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

%% Check inputs
%================================================
if ~exist('VARopt','var')
    error('You need to provide VAR options (VARopt from VARmodel)');
end
% If there is VARopt check that vnames and snames are not empty
vnames = VARopt.vnames;
snames = VARopt.snames;
if isempty(vnames)
    error('You need to add label for endogenous variables in VARopt');
end
if isempty(snames)
    error('You need to add label for shocks in VARopt');
end

%% Check inputs Define some parameters
%===============================================
filename = [VARopt.figname 'HD_SR_'];
quality = VARopt.quality;
suptitle = VARopt.suptitle;
pick = VARopt.pick;

% Initialize HD matrix
nshocks = length(snames); [nsteps, nvars, ~] = size(HD);

% If one shock is chosen, set the right value for nshocks
if pick<0 || pick>nvars
    error('The selected shock is non valid')
else
    if pick==0
        pick=1;
    else
        nshocks = pick;
    end
end

%% Plot
%================================================
FigSize
for ii=pick:nvars
    colormap(parula)
    BarPlot(HD(:,1:nshocks,ii));
    xlim([1 nsteps]);
    title([vnames{ii}], 'FontWeight','bold','FontSize',10);
    % Save
    FigName = [filename num2str(ii)];
    if quality
        if suptitle==1
            Alphabet = char('a'+(1:nvars)-1);
            SupTitle([Alphabet(jj) ') HD of '  vnames{ii}])
        end
        opt = LegOption; LegSubplot(vnames,opt);
        set(gcf, 'Color', 'w');
        export_fig(FigName,'-pdf','-png','-painters')
    else
        legend(snames)
        print('-dpng','-r100',FigName);
        print('-dpdf','-r100',FigName);
    end
    clf('reset');
end

close all

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
classdef PTKSegmentalBronchi < PTKPlugin
    % PTKSegmentalBronchi. Plugin to find the segmental bronchi in each lobe
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Segmental<br>bronchi'
        ToolTip = 'Finds the segmental bronchi in each lobe'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'DoNothing'
        HidePluginInDisplay = true
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = false
        Visibility = 'Developer'
end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            lobes = dataset.GetResult('PTKLobes');
            centreline_tree = dataset.GetResult('PTKAirwayCentreline');
            results_image = dataset.GetTemplateImage(PTKContext.LungROI);
            centreline_tree.AirwayCentrelineTree.GenerateBranchParameters;
            new_centreline_tree = PTKGetSegmentalBronchiCentrelinesForEachLobe(centreline_tree.AirwayCentrelineTree, lobes, results_image, reporting);

            results = [];
            results.StartBranches = new_centreline_tree;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [A, b, x, ProbInfo] = PRblurgauss(varargin)
% PRblurgauss Image deblurring problem with a Gaussian point spread function
%
% [A, b, x, ProbInfo] = PRblurgauss
% [A, b, x, ProbInfo] = PRblurgauss(n)
% [A, b, x, ProbInfo] = PRblurgauss(n, options)
% [A, b, x, ProbInfo] = PRblurgauss(options)
%
% Generates data for use in image deblurring problems with a Gaussian
% point spread function.
%
% Input:
%  n      -  size of the image. Can be either a scalar n (in this case
%            the size is n x n) or a vector [nrow, ncol] (in this case the
%            size is (nrow x ncol).
%            Default: n = 256.
%  options - Structure containing the following optional fields:
%    trueImage  : test image of size n, of type numeric, 2-D only,
%                 or character string indicating
%                 'pattern1'  : geometrical image
%                 'pattern2'  : geometrical image
%                 'sppattern' : sparse (edges of ) a geometrical image
%                 'ppower'    : random image with patterns of nonzero pixels
%                 'smooth'    : very smooth image
%                 'dot2'      : two small Gaussian shaped dots, e.g., a
%                               binary star
%                 'dotk'      : n/2 small Gaussian shaped dots, e.g., stars
%                               (placement is random, reset using rng(0))
%                 'satellite' : satellite test image
%                 'hst'       : image of the Hubble space telescope
%                 Default: 'hst'.
%                 This image is then stored in the output vector x.
%    BlurLevel  : If choosing one of the built-in PSFs, this sets the
%                 severity of the blur to one of the following:
%                 'mild'
%                 'medium'
%                 'severe'
%                 Default is 'medium'
%    BC         : Specify boundary condition:
%                 'zero'
%                 'periodic'
%                 'reflective' (or 'neumann' or 'reflexive')
%                 Default: 'reflective'
%                 Note that in this case an extended (or padded) test image
%                 is blurred using 'zero' boundary conditions, and then the
%                 central subimage of size n is extracted from the exact and
%                 the blurred image. No inverse crime is committed,
%                 i.e., A*x ~= b.
%    CommitCrime: To get an exact system Ax = b (i.e., commit the inverse
%                 crime), set this to:
%                 'on'
%                 Default is 'off' (do not commit the inverse crime).
%
% Output:
%  A        - blurring matrix, which is a either a psfMatrix object in
%             the case of spatically invariant blur, or a sparse matrix
%             in the case of spatially variant blur
%  b        - blurred vector (i.e., blurred image with stacked columns)
%  x        - image vector, i.e., exact (unknown) image with stacked columns
%  ProbInfo - structure whose fields contain information about problem:
%               problemType : kind of test problem generated
%                             (in this case: 'deblurring')
%               xType       : solution type (in this case 'image2D')
%               bType       : data type (in this case 'image2D')
%               xSize       : size of image x
%               bSize       : size of image b
%               psf         : point spread function
%
% See also: PRblur, PRblurdefocus, PRblurmotion, PRblurrotation,
% PRblurshake, PRblurspeckle, PRdiffusion, PRinvinterp2, PRnmr,
% PRseismic, PRspherical, PRtomo, PRnoise, PRshowb, PRshowx, fspecial

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('trueImage', 'hst', 'BlurLevel', 'medium', ...
    'BC', 'reflective', 'CommitCrime', 'off');

% If input is 'defaults,' return the default options in X
if nargin == 1 && nargout <= 1 && strcmp(varargin,'defaults')
    A = defaultopt;
    return;
end

% Check for acceptable number of optional input arguments
switch length(varargin)
    case 0
        n = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            n = varargin{1}; options = [];
        else
            n = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            n = varargin{1}; options = varargin{2};
        else
            n = varargin{2}; options = varargin{1};
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = PRset(defaultopt, options);
options = PRset(options, 'PSF', 'gauss');
[A, b, x, ProbInfo] = PRblur(n, options);

repo: IRtools


text:
function [T1 Mls]=cp_Tparam(M, R)
[u,s,v]=svd(M');
T1=v(:,size(M,2)+1:size(M,1))';

%ls matrix
Mls=lsmat_v000(M,R);

repo: ekfmonoslam


text:
%% CALIBRATION - multiClass grasping actions (realMove data only, EMG pattern classification approach)
clear all; close all; clc;
tic

%% Load Data
dd = 'C:\Users\Doyeunlee\Desktop\Analysis\01_Raw data\';
filelist = 'session1_sub1_multigrasp_MI_EMG'; % this code is only for realMove data

% for image dataset building
% order_count = 0; % from 0 ~

%% Parameter setting for EEG
filt_band_eeg = [0.1 40]; % [0.3 3]Hz, [0.3 35]Hz -> reference: Decoding natural reach-and-grasp actions from human EEG
% ival_rest = [-7000 -3000];
ival = [0 4000]; % for MOTOR IMAGERY % 250Hz sampling? Should besame for EEG and EMG
% subChannel_eeg = [11 12 13 15 16 17 18 20 21 22 23 43 44 45 47 48 49 51 52 53]; %20 out of 60 channels
subChannel_eeg = [6 7 8 9 11 12 13 15 16 17 18 20 21 22 23 25 26 27 28 39 40 10 ...
    43 44 45 47 48 49 51 52 53 55 56 57]; %34 out of 60 channels
main_class = {'Cylindrical','Spherical','Lumbrical','Rest'};

% Parameter setting for EMG
filt_band_emg = [10 225]; % EMG [10 225]Hz
subChannel_emg = [65:70]; % 4 EMG corresponding to hand and wrist movements only (grasp actions)
% Total 71ch(60 EEG, 4 EOG, 7 EMG) <- 6 EMG channels [65:70] = 6 EMG

%% load thres.mat; % threshold data of EEG and EMG signals, respectively.
ival_thres = [-6000 -3000];
bline_list = [0.75 1.0 1.25 1.5];

for bline_count = 1:size(bline_list, 2)
    bline_val = bline_list(bline_count);
    [thres_EEG, thres_EMG] = thresCalib(dd, filelist, ival_thres, filt_band_eeg, filt_band_emg, ...
        subChannel_eeg, subChannel_emg, main_class);

    %% converted data file to cnt - EMG, EEG all same
    [cnt, mrk, mnt]=eegfile_loadMatlab([dd filelist], 'fs', 500); % Load cnt, mrk, mnt variables to Matlab (grasp data)
    % separate cnt into EEG's and EMG's
    %for EEG
    cnt_eeg = proc_commonAverageReference(cnt, subChannel_eeg);
    cnt_eeg = proc_selectChannels(cnt_eeg, subChannel_eeg); % Channel Selection
    % for EMG
%     cnt_emg = proc_commonAverageReference(cnt, subChannel_emg);
    cnt_emg = proc_selectChannels(cnt_emg, subChannel_emg);

    %% IIR Filter (4th butterworth filtter -> reference: Decoding natural reach-and-grasp actions)
    %for EEG
    cnt_eeg = proc_filtButter(cnt_eeg, 3, filt_band_eeg);  % Band-pass filtering

    % FIR FIlter - Option 1
%     [y1, D1] = bandpass_computeGPU(cnt.x, filt_band, cnt.fs, 'ImpulseResponse', 'fir');
%     cnt.x = gather(y1);

    % FIR Filter - Option 2
%     Fs = cnt_eeg.fs;  % Sampling Frequency
%     N    = 16;     % Order
%     Fc1  = filt_band_eeg(1);    % First Cutoff Frequency
%     Fc2  = filt_band_eeg(2);    % Second C utoff Frequency
%     flag = 'scale';  % Sampling Flag
%     % Create the window vector for the design algorithm.
%     win = blackman(N+1);
%     % Calculate the coefficients using the FIR1 function.
%     b  = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);
%     Hd = dfilt.dffir(b);
%     cnt_eeg.x = filter(Hd, cnt_eeg.x);

    %for EMG
%     cnt_emg = proc_filtButter(cnt_emg, 5, filt_band_emg);  % Band-pass filtering

    %% MNT_ADAPTMONTAGE - Adapts an electrode montage to another electrode set
    mnt = mnt_adaptMontage(mnt, cnt); %MNT_ADAPTMONTAGE - Adapts an electrode montage to another electrode set

    %% Noise elimination
    % cnt_eeg.x = medfilt1(cnt_eeg.x); % median values, to remove the spike noise on data
    % cnt_emg.x = medfilt1(cnt_emg.x);

    %% Normalization, rescale, smoothing filter
%     cnt_eeg.x = normalize(cnt_eeg.x);
%     cnt_emg.x = normalize(cnt_emg.x);
%     cnt_eeg.x = hampel(cnt_eeg.x, 2000);
%     cnt_emg.x = hampel(cnt_emg.x, 2000);

    %% cnt to epoch
%     w_size = 1010; %(ms), window size 400, step size 150, overlap size = 250(ms) -> 25 segments
%     step_size = 130;
%     overlap_size = w_size - step_size; % 250ms
%     org_ival = ival; % [0 4000]= 4000ms
%     seg_length = abs(org_ival(2) - org_ival(1));
%     n_segment = ((seg_length - w_size)/(w_size - overlap_size)) + 1; % 25 segments
%
%     h = 0;
%     for h = 1:ch

    k = 0;
    for k = 1:size(main_class, 2) % repeat 3 classes
            i = 0;
            for i = 1:n_segment % repeat 25 segments
                s_point = ((i-1)*step_size) + org_ival(1);
                end_point = s_point + w_size;
                sub_ival = [s_point end_point];

                % make epoch for EMG
                epo_emg = cntToEpo(cnt_emg, mrk, sub_ival);
                epo_emg = proc_selectClasses(epo_emg, main_class{k});
                % make epoch for EEG
                epo_eeg = cntToEpo(cnt_eeg, mrk, sub_ival);
                epo_eeg = proc_selectClasses(epo_eeg, main_class{k});

                for j = 1:50 % number of trials
                    % create muscle patterns
                    temp.muscle_fv = epo_emg.x(:, :, j);
                    temp.ch1 = temp.muscle_fv(:, 1); temp.ch2 = temp.muscle_fv(:, 2);
                    temp.ch3 = temp.muscle_fv(:, 3); temp.ch4 = temp.muscle_fv(:, 4);
                    temp.ch5 = temp.muscle_fv(:, 5); temp.ch6 = temp.muscle_fv(:, 6);

                    temp.rms_fv(1) = rms(temp.ch1); temp.rms_fv(2) = rms(temp.ch2);
                    temp.rms_fv(3) = rms(temp.ch3); temp.rms_fv(4) = rms(temp.ch4);
                    temp.rms_fv(5) = rms(temp.ch5); temp.rms_fv(6) = rms(temp.ch6);

                    % build EEG epoches per each class
                    max_size = size(epo_eeg.x(:, :, j), 1);
                    temp.x{j}(1:max_size, :, i) = epo_eeg.x(:, :, j);

                    ch_num = 0;
                    for ch_num = 1:6 % number of channels
                        if temp.rms_fv(ch_num) > bline_val * thres_EMG{ch_num} % ???% of rest state rms(EMG value)
                            temp.emg_ptrn(i, ch_num, j) = 1;
                        else
                            temp.emg_ptrn(i, ch_num, j) = 0;

                        end
                    end
                end % for j
            end % for i
        emg_ptrn{k} = temp.emg_ptrn;
        temp.x_all{k} = temp.x;

    end % for k
end

X = abs(temp.muscle_fv);
X_1 = movmean(X, 30);
plot(X_1);
hold on;
xlim([0 4000]);

%     k = 0;
%     for k = 1:3
%         i = 0; temp.epo_x{k} = temp.x_all{k}{1};
%         temp.y{k} = emg_ptrn{k}(:, :, 1)'; % create label
%         for i = 2:50
%             % epoched EEG data
%             temp.epo_x{k} = cat(3, temp.epo_x{k}, temp.x_all{k}{i});
%             % labeling from EMG pattern 0, 1
%             temp.y{k} = cat(2, temp.y{k}, emg_ptrn{k}(:, :, i)');
%
%         end
%
%     end
%
%     epo_all.x = cat(3, temp.epo_x{1}, temp.epo_x{2}, temp.epo_x{3});
%     epo_all.y = cat(2, temp.y{1}, temp.y{2}, temp.y{3});
%
%     epo_all.className = {'CH1', 'CH2', 'CH3', 'CH4', 'CH5', 'CH6'};
%     epo_all.clab = epo_eeg.clab; epo_all.fs = epo_eeg.fs; epo_all.title = epo_eeg.title; epo_all.file = epo_eeg.file;
%
%     % build epo
%     epo.x = epo_all.x;
%     epo.clab = epo_all.clab; epo.fs = epo_all.fs; epo.title = epo_all.title; epo.file = epo_all.file;
%
%     %% repeat cross-validation for six times (classification results by channel)
%     for n = 1:6
%         %% create epo by each class
%         % One class vs. others (binary approach 6 times)
%         epo.y = epo_all.y(n, :); epo.className = {'1', '0'}; % one of six channels
%
%         temp.epo_y = zeros(2, size(epo.y, 2));
%         temp.epo_y(1, :) = epo.y;
%          for i = 1:size(epo.y, 2)
%             if temp.epo_y(1, i) == 1
%                 temp.epo_y(2, i) = 0;
%             else
%                 temp.epo_y(2, i) = 1;
%
%             end
%
%          end
%         epo.y = temp.epo_y;
%
%         %% eliminating NaN values
%         % epo.x = standardizeMissing(epo.x, -99);
%         % epo.x = fillmissing(epo.x, 'previous');
%
%         %% CSP feature extraction
%         [csp_fv, csp_w{bline_count}{n}] = proc_csp3(epo); % apply CSP for binary classification
% %         [csp_fv, csp_w{bline_count}{n}] = proc_csp_regularised(epo, 17, 0.7); % regulized CSP
%         % % [csp_fv, csp_w] = proc_cspscp(epo_all, 2, 1); %CSP slow cortical potential variations
% %         [csp_fv, csp_w{bline_count}{n}] = proc_csssp(epo, 17); % Common Sparse Spectrum Spatial Pattern
% %         [csp_fv, csp_w] = proc_cspp_auto(epo); %auto csp patches, only for binary-class
%         csp_fv = proc_variance(csp_fv);
%         csp_fv = proc_logarithm(csp_fv);
%
%         %% all feature vectors
%     %     all_csp_fv{n} = csp_fv; all_csp_w{n} = csp_w;
%
%         %% train classifier ver1
%         % N = length(epo.y);
%         % idxTr = 1:N/2;
%         % C = trainClassifier(csp_fv, 'LDA', idxTr);
%
%         %% train classifier ver2 (Recommended)
%         csp_fv.classifier_param = {'LDA', 'prior', nan, 'store_prior', 1, 'store_means', 1, ...
%             'store_cov', 1, 'store_invcov', 1, 'scaling', 1};
%
%         proc = {'wr_multiClass','policy','one-vs-all','coding','hamming'}; % one-vs-all all-pairs
%         c_out{bline_count}{n}.C = trainClassifier(csp_fv, proc);
%         c_out{bline_count}{n}.out_eeg = applyClassifier(csp_fv, 'wr_multiClass', c_out{bline_count}{n}.C);
%
%         %% apply classifier
%         [c_result{bline_count}{n}] = tnsre_applyClassifer(epo, c_out{bline_count}{n}, csp_w{bline_count}{n});
%
%         %% check for classification accuracy each channel
%         temp.true_label = epo.y(1, :); temp.test_label = c_result{bline_count}{n};
%         temp.score = abs(temp.true_label - temp.test_label);
%         c_score{bline_count}{n} = 1 - (sum(temp.score)/size(epo.y, 2));
%
%         %% Cross-validation for performance evaluation (Including CSP feature extraction process)
%     %     proc = struct('memo', 'csp_w');
%     %     proc.train= ['[fv, csp_w, csp_eig]= proc_csp3(fv, 3); ' ...
%     %         'fv= proc_variance(fv); ' ...
%     %         'fv= proc_logarithm(fv);'];
%     %     proc.apply= ['fv= proc_linearDerivation(fv, csp_w); ' ...
%     %         'fv= proc_variance(fv); ' ...
%     %         'fv= proc_logarithm(fv);'];
%     %     [C_eeg, loss_eeg_std, out_eeg.out, memo] = xvalidation(epo, 'LDA', 'proc', proc, 'kfold', 5);
%     %
%     %     % Result after cross validation = 1-error rate
%     %     Result{n} = 1 - C_eeg;
%     %     Result_Std{n} = loss_eeg_std;
%
%     end
%     % % classification result from cross-validation test for each channel
%     % Result
%     % Result_Std
%
%     % estimted EMG pattern from EEG data decoding
%     esti_ptrn = vertcat(c_result{bline_count}{1}, c_result{bline_count}{2}, c_result{bline_count}{3}, c_result{bline_count}{4}, ...
%         c_result{bline_count}{5}, c_result{bline_count}{6});
%
%     %% True label (true_ptrn) and test label (test_ptrn)
%     i = 0; s1 = 1; s2 = n_segment*50*1+1; s3 = n_segment*50*2+1;
%     true_l{bline_count} = epo_all.y; test_l{bline_count} = esti_ptrn;
%
% end
% %% save the trained classifier
% % % save C and csp_w (for version 1)
% % save('classifierVer1_out.mat', 'C', 'csp_w', 'filt_band_eeg', 'ival', 'subChannel_eeg', 'epo');
% % disp('Saving classifier is done!');
%
% % save C, out_eeg (for version 2)
% save classifier6emg_out.mat c_out csp_w;
% % disp('Calibration work is done!');
%
% %%
% i = 0; true_label = true_l{1}; test_label = test_l{1};
% for i = 2:bline_count
%     true_label = vertcat(true_label, true_l{i});
%     test_label = vertcat(test_label, test_l{i});
%
% end
%
% for i = 1:50
%     nd1 = s1+(n_segment-1); nd2 = s2+(n_segment-1); nd3 = s3+(n_segment-1);
%     % _ptrn{1} is for class 1 ... to 3
%     true_ptrn{1}(:, :, i) = true_label(:, s1:nd1);
%     true_ptrn{2}(:, :, i) = true_label(:, s2:nd2);
%     true_ptrn{3}(:, :, i) = true_label(:, s3:nd3);
%
%     test_ptrn{1}(:, :, i) = test_label(:, s1:nd1);
%     test_ptrn{2}(:, :, i) = test_label(:, s2:nd2);
%     test_ptrn{3}(:, :, i) = test_label(:, s3:nd3);
%
%     s1 = nd1+1; s2 = nd2+1; s3 = nd3+1;
%
% end
%
% %% save the dataset
% save truelabel6emg.mat true_label;
%
% %% test_ptrn to epo_ptrn dataset
% i = 0; temp.epo_ptrn = cat(3, test_ptrn{1}, test_ptrn{2}, test_ptrn{3});
% for i = 1:150
%     epo_ptrn.x(:, :, i) = temp.epo_ptrn(:, :, i)';
%
% end
% epo_ptrn.y = zeros(3, 150);
% epo_ptrn.y(1, 1:50) = 1; epo_ptrn.y(2, 51:100) = 1; epo_ptrn.y(3, 101:150) = 1;
% epo_ptrn.className = main_class;
%
% %% ptrn matching test
% class_num = 0; %1: cylindrical, 2: spherical, 3: lumbrical
% for class_num = 1:3
%     j = 0; n = 0; i = 0;
%     match_result = zeros(50, 3); match_val = zeros(1, 50);
%     for j = 1:50 % test_data trials
%         for n = 1:3
%             for i = 1:50 % trials
%                 test_data = test_ptrn{class_num}(:, :, j);
%                 ref_data = true_ptrn{n}(:, :, i);
%                 % matching method (ssim, immse, psnr, ...)
%                 error_val = immse(test_data, ref_data);
% %                 match_val(i) = ssim(test_data, ref_data);
%                 match_val(i) = 1 - error_val;
%
%             end
%             match_result(j, n) = max(match_val);
%
%         end
%     end
%
%     [M, I] = max(match_result');
%     R_percent(class_num) = sum(I==class_num)/50; R(class_num) = mode(I);
%
% end
%
% R
% R_percent
% fianl_accuracy = sum(R_percent)/3
%
%
% %% save img_data (50 images for each class)
% % create 36 by 36 size images
% % n = 0; i = 0; s2 = '.jpg';
% % dir1 = '/Users/JEONG/Documents/MATLAB/robotArm_code/multiGraspDecoding/Online/temp_imgData/';
% % for n = 1:3
% %     dir2 = num2str(n-1); dir2 = strcat(dir2, '/'); % label folder starts from 0
% %     for i = 1:50
% %         s1 = num2str(i+(50*order_count));
% %         filename = strcat(s1, s2);
% %         dir = strcat(dir1, dir2);
% %         data = test_ptrn{n}(:, :, i);
% %         imwrite(data, [dir filename]);
% %
% %     end
% % end
% %
% % disp('Image saving is completed!');
%
% toc

author: PatternRecognition

repo: OpenBMI


text:
function run_test_single_model
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% Given a model, evaluate the performance.
    baseDir = '../../../';
    addpath([baseDir, filesep, 'codes']);
    addpath([baseDir, filesep, 'codes', filesep, 'denoising']);

    addpath([baseDir, filesep, 'tools', filesep,'bss_eval']);
    addpath([baseDir, filesep, 'tools', filesep,'bss_eval_3']);
    addpath([baseDir, filesep, 'tools', filesep,'labrosa']);

    ModelPath=[baseDir, filesep, 'codes',filesep,'denoising', filesep, 'demo'];

    global SDR;
    global SDR_bss3;

    SDR.deviter=0;   SDR.devmax=0;   SDR.testmax=0;
    SDR.devsar=0; SDR.devsir=0; SDR.testsar=0; SDR.testsir=0;
    SDR_bss3.deviter=0;   SDR_bss3.devmax=0;   SDR_bss3.testmax=0;
    SDR_bss3.devsar=0; SDR_bss3.devsir=0; SDR_bss3.testsar=0; SDR_bss3.testsir=0;

    j=870;

    % Load model
    load([ModelPath, filesep, 'denoising_model_', num2str(j),'.mat']);
    eI.saveDir = [baseDir, filesep, 'codes', filesep, 'denoising', ...
        filesep, 'demo', filesep, 'results', filesep];
    %%
    index = 2;
    [speech, fs] = audioread(['wav', filesep, 'original_speech', num2str(index), '.wav']);
    [noise, fs] = audioread(['wav', filesep, 'original_noise',num2str(index),'.wav']);

    x = speech + noise;
    eI.fs = fs;
    %%
    output = test_denoising_general_kl_bss3(x', theta, eI, 'testall', 0);
    %%
  	sz = 1024.*[1 1/4];
    wn = sqrt( hann( sz(1), 'periodic')); % hann window
    wav_singal = stft2( output.source_signal, sz(1), sz(2), 0, wn);
    wav_noise = stft2( output.source_noise, sz(1), sz(2), 0, wn);
    wav_singal = wav_singal./max(abs(wav_singal));
    wav_noise = wav_noise./max(abs(wav_noise));

    audiowrite([eI.saveDir, filesep,'separated_speech',num2str(index),'.wav'], wav_singal, fs);
    audiowrite([eI.saveDir, filesep,'separated_noise',num2str(index),'.wav'], wav_noise, fs);

    % Get separation stats
    [sdr,sir,sar,stoi] = sep_perf(wav_singal, [speech'; noise'], fs);
end

author: posenhuang

repo: deeplearningsourceseparation


text: %TIMES Mapping overload

author: marianux

repo: ecg-kit

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


text:
function showrateh3(h1,err1,k1,opt1,str1,h2,err2,k2,opt2,str2,h3,err3,k3,opt3,str3)
%% SHOWRATEH3 rate of two error sequences
%
% showrate3(N1,err1,k1,opt1,str1,N2,err2,k2,opt2,str2,N3,err3,k3,opt3,str3)
% plots the err1 vs N1, err2 vs N2, and err3 vs N3 in the loglog scale.
% Additional input
%
%   - k1, k2, k3: specify the starting indices; see showrate
%   - opt1, opt2, opt3: the line color and style
%   - str1, str2, str3: strings used in legend
%
% Example
%
% showrate2(N,energyErr,1,'r-+','||u-u_h||_A',...
%           N,L2Err,1,'b-+','||u-u_h||');
%
% See also showrate, showresult, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N1 = 1./h1; N2 = 1./h2; N3 =1./h3;
if (nargin<=2)
    k1 = 1; opt1 = '-*';
end
r1 = showrate(N1,err1,k1,opt1);
hold on
r2 = showrate(N2,err2,k2,opt2);
r3 = showrate(N3,err3,k3,opt3);
h_legend = legend(str1,['C_1h^{' num2str(-r1) '}'],...
                  str2,['C_2h^{' num2str(-r2) '}'],...
                  str3,['C_3h^{' num2str(-r3) '}'],'LOCATION','Best');
set(h_legend,'FontSize',12);
xlabel('log(1/h)');

repo: Aero-matlab


text:
function writePajekNet(model)
% Builds a metabolite centric directed graph from a COBRA model
% and outputs a graph in a .net format ready
% to use for most graph analysis software e.g. Pajek, it does one fba to
% set the link width equal to reaction fluxes.
%
% USAGE:
%
%    writePajekNet(model)
%
% INPUT:
%    model:    a COBRA structured model
%
% OUTPUT:
%    .net:     file containing the graph
%
% Ex: `A + B -> C` (hypergraph) with `v = 0` => no output (empty line)
%
% if `v > 0` then it becomes `A -> C`; `B -> C` (graph),
%
% if `v < 0` then the order is reversed
%
% .. Author: - Marouen BEN GUEBILA 20/01/2016

m = length(model.mets); %reaction and metabolite number
n = length(model.rxns);

%performs one FBA
FBA = solveCobraLP(model);

%creats a .net file
fileID = fopen('COBRAmodel.net', 'w');
l = 0;  % initialize edge number

% write node id
fprintf(fileID, '*Vertices %d\n', m);

for i = 1:m
    fprintf(fileID, [num2str(i) ' "' model.mets{i} '"\n']);
end

for i = 1:n
    % cleans graph from biomass
    biomassRxn = strfind(model.rxns{i}, 'biomass');

    % cleans graph from objective functions
    objRxn = strfind(model.rxns{i}, 'objective');
    if isequal(biomassRxn, []) && isequal(objRxn, [])
        biomassRxn = 0;
    else
        biomassRxn = 1;
    end

    if (FBA.full(i) == 0 || biomassRxn)
        continue %controls for active reactions
    else
        metPos = find(model.S(:, i) > 0);
        metNeg = find(model.S(:, i) < 0);

        %cleans graph from demand and sink reactions
        if isequal(metPos, zeros(0, 1)) || isequal(metNeg, zeros(0, 1))
            continue
        end
        l = l + (length(metNeg) * length(metPos));
    end
end

%write edge number
fprintf(fileID, '*Edges %d\n', l);

for i = 1:n
    biomassRxn = strfind(model.rxns{i}, 'biomass');
    objRxn = strfind(model.rxns{i}, 'objective');

    if isequal(biomassRxn, []) && isequal(objRxn, [])
        biomassRxn = 0;
    else
        biomassRxn = 1;
    end

    if (FBA.full(i) == 0 || biomassRxn)
        continue
    else
        metPos = find(model.S(:, i) > 0);
        metNeg = find(model.S(:, i) < 0);

        if isequal(metPos, zeros(0, 1)) || isequal(metNeg, zeros(0, 1))
            continue
        end

        for j = 1:length(metNeg)
            for k = 1:length(metPos)
                %take into account flux value and directionality
                if FBA.full(i) > 0
                    fprintf(fileID, [num2str(metNeg(j)) ' ' num2str(metPos(k)) ' ' num2str(abs(FBA.full(i))) '\n']);  % write edges
                else
                    fprintf(fileID, [num2str(metPos(k)) ' ' num2str(metNeg(j)) ' ' num2str(abs(FBA.full(i))) '\n']);
                end

                l = l + 1;
            end
        end
    end
end

%close file
fclose(fileID);

end

author: opencobra

repo: cobratoolbox


text:
% VL_IMINTEGRAL  Compute integral image
%   J = VL_IMINTEGRAL(I) calculates the integral image J of the image
%   I.  I must a matrix with DOUBLE, SINGLE, UINT32, or INT32 storage
%   class. J is given by
%
%    J(i,j) = sum(I(1:i,1:j)).
%
%   J has the same size as I and the same storage class.
%
%   Example::
%     The following identity holds:
%       VL_IMINTEGRAL(ONES(3)) = [ 1 2 3 ;
%                                  2 4 6 ;
%                                  3 6 9 ]
%
%   See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [qd] = stateToQd(x)
%Converts qd struct used in hardware to x vector used in simulation
% x is 1 x 13 vector of state variables [pos vel quat omega]
% qd is a struct including the fields pos, vel, euler, and omega

%current state
qd.pos = x(1:3);
qd.vel = x(4:6);

Rot = QuatToRot(x(7:10)');
[phi, theta, yaw] = RotToRPY_ZXY(Rot);

qd.euler = [phi; theta; yaw];
qd.omega = x(11:13);

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function transformationMatrix = create_transformation_matrix3d(varargin)
% CREATE_TRANSFORMATION_MATRIX3D Create a 3D transformation matrix
%
% [transformationMatrix rotationMatrix shearMatrix scaleMatrix] = ...
%   create_transformation_matrix3d()
%
% INPUT ARGUMENTS
% N/A
%
% Optional input arguments
% 'phiX'            - Rotation angle around x-axis
% 'phiY'            - Rotation angle around y-axis
% 'phiZ'            - Rotation angle around z-axis
% 'shearXY'         - Shear XY
% 'shearXZ'         - Shear XZ
% 'shearYX'         - Shear YX
% 'shearYZ'         - Shear YZ
% 'shearZX'         - Shear ZX
% 'shearZY'         - Shear ZY
% 'scaleX'          - Scale along x-axis
% 'scaleY'          - Scale along y-axis
% 'scaleZ'          - Scale along z-axis
% 'translationX'    - Translation along x-axis
% 'translationY'    - Translation along y-axis
% 'translationZ'    - Translation along z-axis
%
% OUTPUT ARGUMENTS
% transformationMatrix  - Total transformation matrix
% rotationMatrix        - Rotation matrix
% shearMatrix           - Shear matrix
% scaleMatrix           - Scale matrix

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

phiX = 0;
phiY = 0;
phiZ = 0;
shearXY = 0;
shearXZ = 0;
shearYX = 0;
shearYZ = 0;
shearZX = 0;
shearZY = 0;
scaleX = 1;
scaleY = 1;
scaleZ = 1;
translationX = 0;
translationY = 0;
translationZ = 0;

for k=1:2:length(varargin),         % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

rotationX = [1 0         0;
             0 cos(phiX) -sin(phiX);
             0 sin(phiX) cos(phiX)];
rotationY = [cos(phiY)  0 sin(phiY);
             0          1 0;
             -sin(phiY) 0 cos(phiY)];
rotationZ = [cos(phiZ) -sin(phiZ) 0;
             sin(phiZ) cos(phiZ)  0;
             0         0          1];

shearMatrix = [1       shearXY shearXZ;
               shearYX 1       shearYZ;
               shearZX shearZY 1];

scaleMatrix = [scaleX 0      0;
               0      scaleY 0;
               0      0      scaleZ];

transformationMatrix = zeros(3,4);
transformationMatrix(1:3,1:3) = rotationX*rotationY*rotationZ*shearMatrix*scaleMatrix;
transformationMatrix(:,4) = [translationX translationY translationZ]';

author: fordanic

repo: image-registration


text:
function [F,S] = bodyForceSSD(A,B,BGrad,U,X,NumRows,NumCols,NumPages,VoxSize,RegDim,HX,HY,HZ);
% bodyForceSSD: compute SSD body force
%
% arguments:
%   A - reference image
%   B - floating image
%   BGrad - gradient of floating image
%   U - displacement field
%   X - field of grid positions
%   NumRows - number of rows in image
%   NumCols - number of columns in image
%   NumPages (optional) - number of pages in (3D) images
%   F - body force field
%   S - value of similarity measure
%
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3 licence.
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

if RegDim == 2 % 2-D

    % evaluate floating image and gradient at deformation grid
    % x and y index are flipped!
    XNew = X-U;

    BWarp = interp2(X(:,:,2),X(:,:,1),B,XNew(:,:,2),XNew(:,:,1),'*linear',0);
    BWarpGrad = cat(3,imfilter(BWarp,HX,'replicate','same'),...
        imfilter(BWarp,HY,'replicate','same'));

    % compute difference image
    DiffImage = A - BWarp;

    % compute SSD value
    S = sum(DiffImage(:).^2);

    % finally, compute body force
    F = repmat(DiffImage,[1 1 2]).*BWarpGrad;

else % 3-D

    % evaluate floating image and gradient at deformation grid
    % x and y index are flipped!
    XNew = X-U;

    BWarp = interp3(X(:,:,:,2),X(:,:,:,1),X(:,:,:,3),B,XNew(:,:,:,2),XNew(:,:,:,1),XNew(:,:,:,3),'*linear',0);
    BWarpGrad = cat(4,imfilter(BWarp,HX,'replicate','same'),...
        imfilter(BWarp,HY,'replicate','same'),...
        imfilter(BWarp,HZ,'replicate','same'));

    % compute difference image
    DiffImage = A - BWarp;

    % compute SSD value
    S = sum(DiffImage(:).^2);

    % finally, compute body force
    F = repmat(DiffImage,[1 1 1 3]).*BWarpGrad;

end


text:
function tr = traces(a)

%function tr = traces(a)
%
% Provides a vector of traces of the vector LTI a
%

%S. de Waele, March 2003.

a_len = kingsize(a,3);
tr = zeros(a_len,1);

for t = 1:a_len,
   tr(t) = trace(a(:,:,t));
end

author: Sable

repo: mcbench-benchmarks


text:
function t_opf_dc_mips(quiet)
%T_OPF_DC_MIPS  Tests for DC optimal power flow using MIPS solver.

%   MATPOWER
%   Copyright (c) 2004-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

num_tests = 43;

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
    sing_matrix_warn_id = 'Octave:singular-matrix';
else
    sing_matrix_warn_id = 'MATLAB:singularMatrix';
    near_sing_matrix_warn_id = 'MATLAB:nearlySingularMatrix';
    s3 = warning('query', near_sing_matrix_warn_id);
    warning('off', near_sing_matrix_warn_id);
end
s2 = warning('query', sing_matrix_warn_id);

t0 = 'DC OPF (MIPS): ';
mpopt = mpoption('out.all', 0, 'verbose', verbose);
mpopt = mpoption(mpopt, 'opf.dc.solver', 'MIPS', 'mips.comptol', 1e-9);

%% reveals a bug introduced with summer 2017 OPF refactorization
if have_feature('pdipmopf')
    mpopt = mpoption(mpopt, 'opf.ac.solver', 'PDIPM');
end

    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved DC power flow case from MAT-file
    load soln9_dcopf;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    t = [t0 'w/angle diff lims : '];
    mpc = loadcase(casefile);
    mpc.branch(4, ANGMAX) = 3;
    mpc.branch(7, ANGMIN) = -4.5;
    r = rundcopf(mpc, mpopt);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(   f, 6456.7213, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,PG        ),    [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    e = zeros(size(branch, 1), 1);
    e(4) = 297.83776;
    e(7) = -26.94788;
    t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']);

    t = [t0 'w/ignored angle diff lims : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    r = rundcopf(mpc, mpopt1);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% two new z variables
    %%      0 <= z1, P3 - P1 <= z1
    %%      0 <= z2, P3 - P2 <= z2
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 4, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 4, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 3, [t 'user costs']);

    %% with A and N sized for AC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 2 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 4, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 4, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 3, [t 'user costs']);

    t = [t0 'infeasible : '];
    warning('off', sing_matrix_warn_id);
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);   %% Pg1 + Pg2
    mpc.u = Inf;
    mpc.l = 600;
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(~success, [t 'no success']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc = loadcase(casefile);
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = rundcopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
else
    warning(s3.state, near_sing_matrix_warn_id);
end
warning(s2.state, sing_matrix_warn_id);

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [v,t]=histndim(x,b,mode)
%HISTNDIM - generates and/or plots an n-dimensional histogram
%
%  Inputs:  X(m,d)   is the input data: each row is one d-dimensiona data point
%           B(3,d)   specifies the histogram bins.
%                         B(1,:) gives the number of bins in each dmension [default 10]
%                         B(2,:) gives the minimum of the first bin in each dimension [default min(X)]
%                         B(3,:) gives the maximum of the last bin in each dimension [default max(X)]
%                    If B has only one column, the same values are use for al dimensions
%                    If B(1,i)=0 then that dimension will be ignored (and excluded from V)
%           MODE     is a character string containing a combination of the following:
%                        'z' for zero base in the 2D plot [default base = min(V)]
%                        'p' to scale V as probabilities [default actual counts]
%                        'h' to plot a histogam even if output arguments are present
%
% Outputs:  V        d-dimensional array containing the histogram counts
%           T        d-element cell array. d{i} contains the bin boundary values for
%                    the i'th dimension. The length of d{i} is one more than the number of bins
%                    in that dimension.
%
%                    Note that if any of B(1,:) are zero then the number of dimensions in V and elements
%                    of T will be correspondingly reduced.
%
% Example: histndim(randn(100000,2),[20 -3 3]','pz');

%	   Copyright (C) Mike Brookes 2004
%      Version: $Id: histndim.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,d]=size(x);
if nargin<3
    mode=' ';
    if(nargin<2)
        b=repmat(10,1,d);
    end
end

if size(b,2)==1
    b=repmat(b,1,d);
end
if size(b,1)<3
    mi=min(x,[],1);
    ma=max(x,[],1);
    w=(ma-mi)./(b(1,:)-0.001);  % nudge slightly to make sure al points included
    b(3,:)=ma+0.0005*w;
    b(2,:)=mi-0.0005*w;
end

acd=find(b(1,:)>0);
sv=b(1,acd);
nbt=prod(sv);
t=cell(length(acd),1);

% loop through each dimension
k=1;        % indexing factor
ok=repmat(1>0,n,1);
ix=repmat(nbt-sum(cumprod(sv)),n,1);
for i=1:length(acd)
    j=acd(i);
    bw=b(1,j)/(b(3,j)-b(2,j));
    bi=ceil((x(:,j)-b(2,j))*bw);
    ok=ok & (bi>0) & (bi<=b(1,j));
    ix(ok)=ix(ok)+k*bi(ok);
    k=k*b(1,j);
    t{i}=b(2,j)+(0:b(1,j))/bw;
end
v=full(sparse(ix(ok),1,1,nbt,1));
if length(sv)>1
    v=reshape(v,sv);
end
if any(mode=='p')
    v=v/n;
end

if ~nargout | any(mode=='h')
    svg=find(sv>1);
    if length(svg)==1
        j=acd(svg);
        bar(b(2,j)+(0.5:sv(svg)-0.5)*(b(3,j)-b(2,j))/b(1,j),v(:));
    elseif length(svg)==2
        j=acd(svg(1));
        k=acd(svg(2));
        bj=b(1,j);
        bk=b(1,k);
        %     imagesc(b(2:3,k),b(2:3,j),reshape(v,b(1,j),b(1,k)));
        vda=kron(reshape(v,bj,bk),[1 1;1 1]);
        if any(mode=='z')
            ba=0;
        else
            ba=min(vda(:));
        end
        vda=[repmat(ba,1,2*bk+2);repmat(ba,2*bj,1) vda repmat(ba,2*bj,1);repmat(ba,1,2*bk+2)];
        jda=kron(t{svg(1)},[1 1]);
        jda=jda-(jda(3)-jda(2))*0.01*[-0.5 (-1).^(1:2*bj) 0.5]; % nudge slightly to avoid MATLAB plotting bug
        kda=kron(t{svg(2)},[1 1]);
        kda=kda-(kda(3)-kda(2))*0.01*[-0.5 (-1).^(1:2*bk) 0.5];
        surf(jda,kda,vda');
        ylabel(sprintf('Axis %d',k));
        xlabel(sprintf('Axis %d',j));
        colorbar;
    else
        fprintf(2,'Error in %s: Cannot plot 3-D histogram\n',mfilename);
    end
end

author: decouples


text:
function pass = test_max3(pref)
% Test @chebfun3/max3 command.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 100*pref.cheb3Prefs.chebfun3eps;

% Check the MAX3 function for a cosine function
f = chebfun3(@(x,y,z) cos(x.*y.*z));
g = chebfun(@(x,y,z) 1 + 0*x);
h1 = max3(f);
pass(1) = norm(h1-g) < tol;

% Check the MAX function for a sine function
f = chebfun3(@(x,y,z) sin(x+y+z));
h2 = max3(f);
pass(2) = norm(h2-g) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function z = g_three ( x, y )

%*****************************************************************************80
%
%% G_THREE evaluates a function which has 6 minima, 2 global, in [-3,3]x[-2,2].
%
%  Discussion:
%
%    Actually, F_THREE evaluates that function; this function returns the
%    log of that function plus 2, which is useful when trying to draw contours.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the arguments.
%
%    Output, real Z, the function value.
%
  z = ( 4.0 - 2.1 * x.^2 + x.^4 / 3.0 ) .* x.^2 ...
    + x .* y + ( - 4.0 + 4.0 * y.^2 ) .* y.^2;

  z = log ( 2 + z );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function nc2010_plot_compexperiment(flatW, datatype, convergence)

n = 200;
m = 50;
d = 10;

model = 'vbpca';
if flatW
  stringW = 'flatW';
else
  stringW = 'hierW';
end

% SELECT THE SEED
seed = 4132
datasets = 10
% Initialise variables
norottime = zeros(m,datasets);
norotloglike = zeros(m,datasets);
rottime = zeros(m,datasets);
rotloglike = zeros(m,datasets);
for ncomps = 1:m
  for ind = 1:datasets
    rotstring = 'rotate=0';
    norotfilename = sprintf(['/home/jluttine/matlab/neurocomputing2010/' ...
                        'nc2010_%s_compexperiment_%s_subspace=%s_n=%d_m=' ...
                        '%d_d=%d_ncomps=%d_%s_seed=%d_datasetindex=%d.mat'], ...
                            model,stringW,datatype,n,m,d, ncomps, rotstring, ...
                            seed, ind);

    rotstring = 'rotate=1';
    rotfilename = sprintf(['/home/jluttine/matlab/neurocomputing2010/' ...
                        'nc2010_%s_compexperiment_%s_subspace=%s_n=%d_m=' ...
                        '%d_d=%d_ncomps=%d_%s_seed=%d_datasetindex=%d.mat'], ...
                            model,stringW,datatype,n,m,d, ncomps, rotstring, ...
                            seed, ind);
    norot = load(norotfilename, 'duration', 'cost', 'rmse', 'rmse_test');
    rot = load(rotfilename, 'duration', 'cost', 'rmse', 'rmse_test');
    %plot_results(norot, rot,ncomps);

    % Find the point of convergence
    switch 2
     case 1 % by relative change in loglikelihood
      len = length(norot.cost);
      dcost = diff(norot.cost);
      dstep = abs(dcost ./ norot.cost(1:(len-1)));
      norotind = find(dstep < convergence, 1, 'first');
      len = length(rot.cost);
      dcost = diff(rot.cost);
      dstep = abs(dcost ./ rot.cost(1:(len-1)));
      rotind = find(dstep < convergence, 1, 'first');
     case 2 % being inside a threshold compared to converged loglikelihood
      convind = find(~isnan(norot.cost), 1, 'last');
      dif = abs((norot.cost - norot.cost(convind))/norot.cost(convind));
      norotind = find(dif < convergence, 1, 'first');
      convind = find(~isnan(rot.cost), 1, 'last');
      dif = abs((rot.cost - rot.cost(convind))/rot.cost(convind));
      rotind = find(dif < convergence, 1, 'first');
    end

    %  norotind = find(~isnan(norot.cost), 1, 'last');
    %  rotind = find(~isnan(rot.cost), 1, 'last');
    norottime(ncomps,ind) = norot.duration(norotind);
    norotloglike(ncomps,ind) = norot.cost(norotind);
    rottime(ncomps,ind) = rot.duration(rotind);
    rotloglike(ncomps,ind) = rot.cost(rotind);
  end
end

figure
norotmedian = prctile(norottime,50,2);
norot90 = prctile(norottime,100,2);
norot10 = prctile(norottime,0,2);
%norotmean = mean(norottime,2);
%rotmean = mean(rottime,2);
rotmedian = prctile(rottime,50,2);
rot90 = prctile(rottime,100,2);
rot10 = prctile(rottime,0,2);
semilogy(1:ncomps, norotmedian(:), 'r-');
hold on
x = [1:ncomps, ncomps:-1:1]';
y = [rot10(:); rot90(end:-1:1)];
c = [0.75 0.75 0.75];
fill(x,y,c,'EdgeColor',c);
%patch(x,y, [0.6 0.6 0.6]);
% $$$ semilogy(1:ncomps, rot90(:), 'k:', 'linewidth', 1);
% $$$ semilogy(1:ncomps, rot10(:), 'k:', 'linewidth', 1);
semilogy(1:ncomps, norotmedian(:), 'r-');
semilogy(1:ncomps, norot90(:), 'r-');
semilogy(1:ncomps, norot10(:), 'r-');
semilogy(1:ncomps, rotmedian(:), 'k:', 'linewidth', 2);
%semilogy(1:ncomps, rot90(:), 'k:', 'linewidth', 1);
%semilogy(1:ncomps, rot10(:), 'k:', 'linewidth', 1);
%[norotloglike(:) - rotloglike(:)];
xlabel('number of components, D')
ylabel('time (seconds)')
xlim([1 m]);

return

% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function plot_results(norot, rot,ncomps)
% $$$
% $$$ % $$$ norot.cost(norot.cost==-inf) = -1e10;
% $$$ % $$$ rot.cost(rot.cost==-inf) = -1e10;
% $$$ %suffix = sprintf(' (n=%d m=%d d=%d ncomps=%d)',n,m,d,ncomps);
% $$$
% $$$ figure
% $$$ subplot(1,3,1)
% $$$ semilogx([norot.duration(:), rot.duration(:)], -[norot.cost(:), rot.cost(:)])
% $$$ %title(['Negative loglikelihood', suffix]);
% $$$ %filename = sprintf('fig_cost_n=%d_m=%d_d=%d_ncomps=%d', n,m,d,ncomps);
% $$$ % $$$ set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
% $$$ % $$$ pos = get(gcf, 'position');
% $$$ % $$$ set(gcf, 'position', [pos(1:2), 8,8]);
% $$$ % $$$ pos = get(gcf, 'paperposition');
% $$$ % $$$ set(gcf, 'paperposition', [pos(1:2),8,8])
% $$$
% $$$ %figure
% $$$ subplot(1,3,2)
% $$$ semilogx([norot.duration(:), rot.duration(:)], [norot.rmse(:), rot.rmse(:)])
% $$$ % $$$ %title(['Training set RMSE', suffix]);
% $$$ % $$$ %filename = sprintf('fig_rmse_n=%d_m=%d_d=%d_ncomps=%d', n,m,d,ncomps);
% $$$ % $$$ set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
% $$$ % $$$ pos = get(gcf, 'position');
% $$$ % $$$ set(gcf, 'position', [pos(1:2), 8,8]);
% $$$ % $$$ pos = get(gcf, 'paperposition');
% $$$ % $$$ set(gcf, 'paperposition', [pos(1:2),8,8])
% $$$
% $$$
% $$$ %figure
% $$$ subplot(1,3,3)
% $$$ semilogx([norot.duration(:), rot.duration(:)], [norot.rmse_test(:), rot.rmse_test(:)])
% $$$ % $$$ %title(['Test set RMSE', suffix]);
% $$$ % $$$ %filename = sprintf('fig_rmsetest_n=%d_m=%d_d=%d_ncomps=%d', n,m,d,ncomps);
% $$$ % $$$ set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
% $$$ % $$$ pos = get(gcf, 'position');
% $$$ % $$$ set(gcf, 'position', [pos(1:2), 8,8]);
% $$$ % $$$ pos = get(gcf, 'paperposition');
% $$$ % $$$ set(gcf, 'paperposition', [pos(1:2),8,8])
% $$$
% $$$

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%% SVD_filter
% Below is a demonstration of the features of the |SVD_filter| function

%%
clear; close all; clc;

%% Syntax
% |[Zm]=SVD_filter(Z,P,T);|

%% Description
% This funciton uses singular value decomposition to smoothen 2D matrix
% data

%% Examples
%

%%
% Plot settings
font_size=20;
cmap=resampleColormap([1 0 0; 1 1 1; 0 0 1],250);

%%
% Create example data

%Create clean data
n=35;
s=5;
[X,Y]=ndgrid(linspace(-4*s,4*s,n));
Z=n*exp( -0.5.*((X./s).^2+(Y./s).^2));
Z(X<0)=-Z(X<0); %Add sharp feature

%Create noise eroded data
Zn=Z+n/30*randn(size(Z));

%% Using |SVD_filter|
P=[1-1e-4 1e-4];
T=1;
[Zm]=SVD_filter(Zn,P,T);

[F,V]=surf2patch(Z);
[~,Vn]=surf2patch(Zn);
[~,Vm]=surf2patch(Zm);
Cn=Vn(:,3)-V(:,3);
Cm=Vm(:,3)-V(:,3);

c=max(abs(Cn(:)));

%%
% Visualize

cFigure;
subplot(1,3,1); hold on;
title('Clean','FontSize',font_size);
hp=gpatch(F,V,V(:,3),'k');
hp.FaceColor='interp';
axisGeom(gca,font_size);
colormap(gca,gjet(250));
camlight headlight;

subplot(1,3,2); hold on;
title('Raw','FontSize',font_size);
hp=gpatch(F,Vn,Cn,'k');
hp.FaceColor='interp';
axisGeom(gca,font_size);
colormap(gca,cmap); colorbar;
caxis([-c,c]);
camlight headlight;

subplot(1,3,3); hold on
title('SVD filtered','FontSize',font_size);
% gpatch(F,V,'g','none',0.5);
hp=gpatch(F,Vm,Cm,'k');
hp.FaceColor='interp';
axisGeom(gca,font_size);
colormap(gca,cmap); colorbar;
camlight headlight;
caxis([-c,c]);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [Y,W,SetupStruc] = Process_AuxIVA(s,Transfer,SetupStruc)
K = SetupStruc.AuxIVA.K;
hop = SetupStruc.AuxIVA.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.AuxIVA.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(Num,frame_N,K);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
X_sp = zeros(Num,frame_N,K_m);
W_IVA = zeros(Num,Num,K_m);
V_sp = zeros(Num,N,K_m);
theta = 10^-6;
for i = 1:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Initialize W by PCA
    [E,D] = PCA(X_f,1,Num);
    V = sqrt(D)\E';
    V_sp(:,:,i) = V;
    X_sp(:,:,i) = V*X_f;
    %%%%%%%%%%%% Adjust amplitude of 'w'
    W_o = eye(Num);
    y_f = W_o*V*X_f;
%     norm = max(abs(y_f),[],2);
%     if(norm>10)
%         norm = repmat(norm,1,Num);
%         W_o = W_o./norm;
%         y_f = W_o*V*X_f;
%     end
    W_IVA(:,:,i) = W_o;
    Y_f(:,:,i) = y_f;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% IVA iterations
max_iteration = 200;
Y_k = zeros(Num,frame_N);
epsi = 1e-6;
pObj = inf;
A = zeros(1001,2)-1; %%%% Show the decrease of non-linear correlation, IVA max iterations 1000
for iteration = 1:max_iteration
    for i = 1:Num
        y_temp = permute(Y_f(i,:,:),[3 2 1]);
        Y_k(i,:) = sqrt(sum(abs(y_temp(1:K_m,:)).^2))+epsi;
    end
    dlw = 0;
    for i = 1:K_m
        W = W_IVA(:,:,i);
        X_f = X_sp(:,:,i);
        dlw = dlw +log(abs(det(W))+epsi);
        for i_n = 1:Num
            G_ = Y_k(i_n,:).^-1;
            G_ = repmat(G_,Num,1);
            Vk = (G_.*X_f)*X_f'/frame_N;
            if rcond(Vk)<theta
                Vk = Vk+eye(Num)*max(eig(Vk))*theta;
            end
            wk = inv(W*Vk);
            wk = wk(:,i_n);
            wk = wk/(sqrt(wk'*Vk*wk)+epsi);
            W(i_n,:) = wk';
        end
        W_IVA(:,:,i) = W;
        Y_f(:,:,i) = W*X_f;
    end
    Obj = (sum(sum(Y_k))/frame_N-2*dlw)/(Num*K_m);
    dObj = pObj-Obj;
    pObj = Obj;
    A(iteration,:) = [Obj,abs(dObj)/abs(Obj)];
    if(abs(dObj)/abs(Obj)<theta)
        break;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Post processing
W = zeros(Num,N,K_m);
Y_f(:,:,1) = zeros(Num,frame_N);
for i = 2:K_m
    W_inv = pinv(W_IVA(:,:,i)*V_sp(:,:,i));
    for ii = 1:Num
        Y_f(ii,:,i) = Y_f(ii,:,i)*W_inv(1,ii);
        W_IVA(ii,:,i) = W_IVA(ii,:,i)*W_inv(1,ii);
    end
    W(:,:,i) = W_IVA(:,:,i)*V_sp(:,:,i);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if(i~=K_m)
        Y_f(:,:,K+2-i) = conj(Y_f(:,:,i));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    y_temp = permute(Y_f(i,:,:),[3 2 1]);
    Y(:,i) = overlapadd(real(ifft(y_temp))',win,hop);
end
return;

repo: Voice-Separation-and-Enhancement


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

figure; % open a new figure window
plot(x, y, 'rx','MarkerSize', 10);
ylabel('Profit in $10,000s');
xlabel('Population of City in 10,000s');

% ============================================================

end

author: anirudhjayaraman

repo: Machine-Learning


text:
function add_to_buffer
global po pso pco
global Fs s
global ao
global pla pau
global hdls
global vol
global hd1 hd2
if pla
    i1=1+pso*(pco-1);
    i2=pso*pco;
    L=size(s,1);
    if (i1>L)||(i2>L)
        pla=false;
        stop(ao);
        delete(ao);
    else
        if pau
            d1=0*s(i1:i2,:);

            putdata(ao,d1);
        else
            d1=vol*s(i1:i2,:);
            d1=[filter(hd1,d1(:,1))  filter(hd2,d1(:,2)) ];

            putdata(ao,d1);
            display_time_script;
            pco=pco+1;
        end
    end
else
    stop(ao);
    delete(ao);
end


% read bboxes from txt file given by the MOTChallege
function bboxes = read_bboxes(file_name)

fid = fopen(file_name, 'r');
% <frame>, <id>, <bb_left>, <bb_top>, <bb_width>, <bb_height>, <conf>, <x>, <y>, <z>
line = textscan(fid, '%d %d %f %f %f %f %f %f %f %f', 'Delimiter', ',');
fclose(fid);

% build the bboxes structure for detections
bboxes.fr = line{1};
bboxes.id = line{2};
bboxes.x = line{3};
bboxes.y = line{4};
bboxes.w = line{5};
bboxes.h = line{6};
bboxes.r = line{7};


text:
%% Wigner-D functions
%
%% Theorie
%
% The Wigner-D functions are special functions on the rotation group
% $SO(3)$.
%
% In terms of Matthies (ZYZ-convention) Euler angles ${\bf R} = {\bf
% R}(\alpha,\beta,\gamma)$ the $L_2$-normalized Wigner-D function of degree
% $n$ and orders $k,l \in \{-n,\dots,n\}$ is defined by
%
% $$ D_n^{k,l}({\bf R}) = \sqrt{2n+1} \, \mathrm e^{-\mathrm i k\gamma}
% \mathrm d_n^{k,l}(\cos\beta) \,e^{-\mathrm i l\alpha} $$
%
% where $d_n^{k,l}$, denote the real valued Wigner-d functions, which are
% defined in terms of Jacobi polynomial $P_s^{a,b}$  by
%
% $$ d_n^{k,l}(x) = (-1)^{\nu} \binom{2n-s}{s+a}^{\frac12}
% \binom{s+b}{b}^{-\frac12} \left(\frac{1-x}{2}\right)^{\frac{a}{2}}
% \left(\frac{1+x}{2}\right)^{\frac{b}2} P_s^{a,b}(x)$$
%
% using the constants $a =|k-l|$, $b =|k+l|$, $s = n - \max\{|k|,|l|\}$ and
% $\nu = \min\{0,k\}+\min\{0,l\}$ if $l \geq k$; $\nu =
% \min\{0,k\}+\min\{0,l\} + k+l$ otherwise.
%
% This definition is slightly different to other well known definitions,
% because the Wigner-D functions are defined compatible to the
% <SphericalHarmonics.html spherical harmonics> which form an orthonormal
% basis on the 2-sphere.
%
%%
% In MTEX the Wigner-D and Wigner-d functions are available through the
% command <Wigner_D.html |Wigner_D|>

% the Wigner-d function of degree 1
beta = 0.5;
d = Wigner_D(1,beta)

% the Wigner-D function of degree 1
R = rotation.rand;
D = sqrt(3) * Wigner_D(1,R)

%%
% Here the orders $k$, $l$ work as row and column indices.
%
%% Series Expansion
%
% The Wigner-D functions form an orthonormal basis in $L_2(SO(3))$. Hence,
% we can describe functions on the rotation group $SO(3)$ by there harmonic
% representation using the class <SO3FunHarmonicRepresentation.html |SO3FunHarmonic|>.
%
% Hence we define the Wigner-D function $D_1^{1,-1}$ by

D = SO3FunHarmonic([0;0;0;1])
D.eval(R)

%%
% Various normalizations for the Wigner-D functions are common in the
% literature.
%
% Here we define the $L_2$-norm by
%
% $$ \| f \|_2 = \left(\frac1{8\pi^2}\,\int_{SO(3)} \lvert f(\bf R) \rvert^2 \,\mathrm d \bf R \right)^{1/2} $$
%
% such that the norm of the constant function $f=1$ is $1$. Take a look on the section
% <SO3FunOperations.html#6 Integration of SO3Fun's>.
%
% Using that definition the Wigner-D functions in MTEX are normalized, i.e. $\|
% D_n^{k,l} \|_2 = 1$ for all $n,k,l$.

norm(D)

%% Some important formulas for Wigner-D functions
%
% The Wigner-D functions are the matrix elements of the representations
% $D_n \colon SO(3) \to \mathbb C^{(2n+1)\times(2n+1)}$ on $SO(3)$.
% Since representations are group homomorphisms, we have
% $D_n(\bf{R} \, \bf{Q}) = \frac1{\sqrt{2n+1}} \, D_n(\bf{Q}) \, D_n(\bf{R}).$
% Hence we get
%
% $$ D_n^{k,l}(\bf{R}\,\bf{Q}) = \frac1{2n+1} \sum_{j=-n}^n D_n^{k,j}(\bf{Q})\,D_n^{j,l}(\bf{R}). $$
%
%%
% Some symmetry properties of Wigner-D functions yields
%
% $$ D_n^{k,l}(\bf{R}) = \overline{D_n^{l,k}(\bf{R}^{-1})}. $$
%

%% Symmetry properties of Wigner-d functions
%
% The Wigner-d functions by construction fulfill a lot of symmetry
% properties. Some importants are
%
% $$ d_n^{k,l}(x) = d_n^{-k,-l}(x) = (-1)^{k+l}\, d_n^{l,k}(x) = (-1)^{k+l}\, d_n^{-l,-k}(x)$$
%
% $$ d_n^{k,l}(x) = (-1)^{n+k+l}\,d_n^{-k,l}(-x) = (-1)^{n+k+l}\,d_n^{k,-l}(-x) $$
%
% $$d_n^{k,l}(\cos\beta) = (-1)^{k+l}\,d_n^{k,l}(\cos(-\beta))$$
%

author: mtex-toolbox

repo: mtex

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


text:
function grad = f_grad(XX,XY,w)
%f_grad
%
grad = XX*w-XY;
end

author: xiubooth