text:
%IMM_FILTER  Interacting Multiple Model (IMM) Filter prediction and update steps
%
% Syntax:
%   [X_i,P_i,MU,X,P] = IMM_FILTER(X_ip,P_ip,MU_ip,p_ij,ind,dims,A,Q,Y,H,R)
%
% In:
%   X_ip  - Cell array containing N^j x 1 mean state estimate vector for
%           each model j after update step of previous time step
%   P_ip  - Cell array containing N^j x N^j state covariance matrix for
%           each model j after update step of previous time step
%   MU_ip - Vector containing the model probabilities at previous time step
%   p_ij  - Model transition matrix
%   ind   - Indices of state components for each model as a cell array
%   dims  - Total number of different state components in the combined system
%   A     - State transition matrices for each model as a cell array.
%   Q     - Process noise matrices for each model as a cell array.
%   Y    - Dx1 measurement vector.
%   H    - Measurement matrices for each model as a cell array.
%   R    - Measurement noise covariances for each model as a cell array.
%
%
% Out:
%   X_p  - Updated state mean for each model as a cell array
%   P_p  - Updated state covariance for each model as a cell array
%   MU   - Model probabilities as vector
%   X    - Combined state mean estimate
%   P    - Combined state covariance estimate
%
% Description:
%   IMM filter prediction and update steps. Use this instead
%   of separate prediction and update functions, if you don't need
%   the prediction estimates.
%
% See also:
%   IMM_UPDATE, IMM_SMOOTH, IMM_FILTER

% History:
%   01.11.2007 JH The first official version.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% $Id: imm_update.m 111 2007-11-01 12:09:23Z jmjharti $
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [X_i,P_i,MU,X,P] = imm_filter(X_ip,P_ip,MU_ip,p_ij,ind,dims,A,Q,Y,H,R)
    % Number of models
    m = length(X_ip);

    % Default values for state mean and covariance
    MM_def = zeros(dims,1);
    PP_def = diag(20*ones(dims,1));

    % Normalizing factors for mixing probabilities
    c_j = zeros(1,m);
    for j = 1:m
        for i = 1:m
            c_j(j) = c_j(j) + p_ij(i,j).*MU_ip(i);
        end
    end

    % Mixing probabilities
    MU_ij = zeros(m,m);
    for i = 1:m
        for j = 1:m
            MU_ij(i,j) = p_ij(i,j) * MU_ip(i) / c_j(j);
        end
    end

    % Calculate the mixed state mean for each filter
    X_0j = cell(1,m);
    for j = 1:m
        X_0j{j} = zeros(dims,1);
        for i = 1:m
            X_0j{j}(ind{i}) = X_0j{j}(ind{i}) + X_ip{i}*MU_ij(i,j);
        end
    end

    % Calculate the mixed state covariance for each filter
    P_0j = cell(1,m);
    for j = 1:m
        P_0j{j} = zeros(dims,dims);
        for i = 1:m
            P_0j{j}(ind{i},ind{i}) = P_0j{j}(ind{i},ind{i}) + MU_ij(i,j)*(P_ip{i} + (X_ip{i}-X_0j{j}(ind{i}))*(X_ip{i}-X_0j{j}(ind{i}))');
        end
    end

    % Space for estimates
    X_p = cell(1,m);
    P_p = cell(1,m);
    X_i = cell(1,m);
    P_i = cell(1,m);
    lambda = zeros(1,m);

    % Filter the estimates for each model
    for i = 1:m
        % Predict the estimates
        [X_p{i}, P_p{i}] = kf_predict(X_0j{i}(ind{i}),P_0j{i}(ind{i},ind{i}),A{i},Q{i});
        % Update the estimates
        [X_i{i}, P_i{i}, K, IM, IS] = kf_update(X_p{i},P_p{i},Y,H{i},R{i});

        % Calculate likelihoods
        lambda(i) = kf_lhood(X_p{i},P_p{i},Y,H{i},R{i});
    end

     % Calculate the model probabilities
    MU = zeros(1,m);
    c = sum(lambda.*c_j);
    MU = c_j.*lambda/c;

    % Output the combined updated state mean and covariance, if wanted.
    if nargout > 3
        % Space for estimates
        X = zeros(dims,1);
        P = zeros(dims,dims);
        % Updated state mean
        for i = 1:m
            X(ind{i}) = X(ind{i}) + MU(i)*X_i{i};
        end
        % Updated state covariance
        for i = 1:m
            P(ind{i},ind{i}) = P(ind{i},ind{i}) + MU(i)*(P_i{i} + (X_i{i}-X(ind{i}))*(X_i{i}-X(ind{i}))');
        end
    end

author: EEA-sensors

repo: ekfukf

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


text:
function bp=backprojection(c);
% function backp=backprojection(c);
% calculates backprojections from eigenvalues (e) and Group Space
%(c.GroupSpace) which are the outputs of cmdscale.

eigenvals=diag(c.eigenvalues(1:c.ndims),0);
A=c.GroupSpace(:,1:c.ndims)*sqrt(eigenvals);
for n=1:size(c.data,3);
    mdata=c.data(:,:,n);
    backproj(:,:,n)=A'*(mdata*inv(mdata'*mdata))';
end
bp=mean(backproj,3);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function spline_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests BEZ_VAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 2;
  a = 0.0E+00;
  b = 1.0E+00;
  nsample = 21;
  y = [ 1.0E+00, 0.0E+00, 1.0E+00 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  BEZ_VAL evaluates a Bezier function.\n' );
%
%  One point on the curve should be (0.75, 20/32).
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  T, X(T), Y(T)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : nsample

    x = ( ( nsample - i     ) * a   ...
        + (           i - 1 ) * b ) ...
        / ( nsample     - 1 );

    fprintf ( 1, '%8d  %14f  %14f\n', i, x, bez_val ( n, x, a, b, y ) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  When X = 0.75E+00\n' );
  fprintf ( 1, '  BEZ_VAL(X) should be 0.625E+00.\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


function [evects,evals] = pca(X)
% [evects,evals] = pca(X)
%
% finds principal components of
%
% input:
%  X  dxn matrix (each column is a dx1 input vector)
%
% output:
% evects  columns are principal components (leading from left->right)
% evals   corresponding eigenvalues
%
% See also applypca
%
% copyright by Kilian Q. Weinberger, 2006

[d,N]  = size(X);

mm = mean(X,2);
X = X - mm*ones(1,N); % remove mean from data

cc = cov(X',1); % compute covariance
[cvv,cdd] = eig(cc); % compute eignvectors
[zz,ii] = sort(diag(-cdd)); % sort according to eigenvalues
evects = cvv(:,ii); % pick leading eigenvectors
cdd = diag(cdd); % extract eigenvalues
evals = cdd(ii); % sort eigenvalues


text:
function varargout = kj_analyze(varargin)
% KJ_ANALYZE Application M-file for kj_analyze.fig
%    FIG = KJ_ANALYZE launch kj_analyze GUI.
%    KJ_ANALYZE('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.0 25-Oct-2001 09:38:02

if (nargin == 2) & ~ischar(varargin{1})  % LAUNCH GUI

  fig = openfig(mfilename,'reuse');

  % Use system color scheme for figure:
  set(fig,'Color',get(groot,'defaultUicontrolBackgroundColor'));

  % Generate a structure of handles to pass to callbacks, and store it.
  handles = guihandles(fig);
  % Store the data
  handles.params = varargin{1};
  handles.hParentFigure = varargin{2};

  % Set the used splitting time if available
  try
    fSplit = handles.params.fSplitTime;
  catch
    fSplit = ((max(handles.params.mCatalog.Date) - min(handles.params.mCatalog.Date)) / 2) + min(handles.params.mCatalog.Date);
  end
  set(handles.txtSplit, 'String', num2str(fSplit, 6));

  % Set the length of the learning period if available
  try
    set(handles.txtFirst, 'String', num2str(handles.params.fLearning));
  catch
    set(handles.txtFirst, 'String', '1');
  end

  % Set the length of the forecasting period if available
  try
    set(handles.txtSecond, 'String', num2str(handles.params.fForecast));
  catch
    set(handles.txtSecond, 'String', '1');
  end

  plotit(handles);

  guidata(fig, handles);

  if nargout > 0
    varargout{1} = fig;
  end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

  try
    [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
  catch ME
    disp(ME.message);
  end

end

% --------------------------------------------------------------------
function plotit(handles)

% Plot the frequency magnitude distribution
if get(handles.radEntireCatalog, 'Value') == 1
  mPlotCatalog = handles.params.mCatalog;
else
  mPlotCatalog = handles.mNodeCatalog;
end

axes(handles.axsFMD);
set(handles.axsFMD, 'NextPlot', 'replace');
set(handles.axsTimePlot, 'NextPlot', 'replace');
if get(handles.radEntireTime, 'Value') == 1
  % Plot the frequency magnitude distribution
  [fBValue, fAValue, fStdDev] = plot_FMD(mPlotCatalog,...
    'Axes',handles.axsFMD, 'CalcMc', handles.params.nCalculateMC, 'MarkerSize', 10);
  set(handles.lblBValue1, 'String', num2str(fBValue, 3));
  set(handles.lblAValue1, 'String', num2str(fAValue, 3));
  set(handles.lblStdValue1, 'String', num2str(fStdDev, 3));
  set(handles.lblBValue2, 'String', '');
  set(handles.lblAValue2, 'String', '');
  set(handles.lblStdValue2, 'String', '');
  % Plot the cumulative number
  [~, vIndices] = sort(mPlotCatalog(:,3));
  mSortedCatalog = mPlotCatalog(vIndices(:,1),:) ;
  axes(handles.axsTimePlot)
  plot(mSortedCatalog(:,3),(1:length(mSortedCatalog(:,3))),'k-');
else
  % Split the catalog
  fFirst = str2double(get(handles.txtFirst, 'String'));
  fSecond = str2double(get(handles.txtSecond, 'String'));
  fSplit = str2double(get(handles.txtSplit, 'String'));
  [mFirst, mSecond] = ex_SplitCatalog(mPlotCatalog, fSplit, 1, fFirst, 1, fSecond);
  % Plot the frequency magnitude distribution for the first period
  [fBValue, fAValue, fStdDev] = plot_FMD(mFirst, 'Axes', handles.axsFMD, ...
      'CalcMethod',handles.params.nCalculateMC, 'MarkerSize', 3);
  set(handles.lblBValue1, 'String', num2str(fBValue, 3));
  set(handles.lblAValue1, 'String', num2str(fAValue, 3));
  set(handles.lblStdValue1, 'String', num2str(fStdDev, 3));
  % Plot the frequency magnitude distribution for the second period
  set(handles.axsFMD, 'NextPlot', 'add');
  [fBValue, fAValue, fStdDev] = plot_FMD(mSecond, 'Axes', handles.axsFMD,...
      'Marker','^', 'Color','r', 'CalcMethod', handles.params.nCalculateMC,'MarkerSize',3);
  set(handles.lblBValue2, 'String', num2str(fBValue, 3));
  set(handles.lblAValue2, 'String', num2str(fAValue, 3));
  set(handles.lblStdValue2, 'String', num2str(fStdDev, 3));

  % Plot the Kagan & Jackson test
  plot_FMD(mSecond, 'ShowCumulative', false, 'Axes', handles.axsTest, 'Marker','o', 'Color', 'b', 'MarkerSize', 4);
  set(handles.axsTest, 'NextPlot', 'add');
  % Calculate Mc
  fMc = calc_Mc(mPlotCatalog, handles.params.nCalculateMC);
  if isnan(fMc)
    fMc = handles.params.fMcOverall;
  elseif isempty(fMc)
    fMc = handles.params.fMcOverall;
  end
  % Define magnitude range for testing
  if handles.params.bMinMagMc
    fMinMag = fMc;
  else
    fMinMag = handles.params.fMinMag;
  end
  % Perform the test
  [fDeltaProbability, fProbabilityK, fProbabilityO, fWeightK, fWeightO, fBValueO, vProbK, vProbO, mPredictionFMD] ...
    = kj_poissonian(mFirst, handles.params.fLearning, mSecond, handles.params.fForecast, handles.params.nTestMethod, ...
    handles.params.nMinimumNumber, fMc, handles.params.fBValueOverall, handles.params.fStdDevOverall, ...
    fMinMag, handles.params.fMaxMag, 1);
  set(handles.txtProbK, 'String', num2str(fProbabilityK));
  set(handles.txtProbO, 'String', num2str(fProbabilityO));
  set(handles.txtProbDiff, 'String', num2str(fDeltaProbability));
  % Plot the results
  vX = (fMinMag:0.1:handles.params.fMaxMag)';
  semilogy(vX, mPredictionFMD(:,1), 'r');
  semilogy(vX, mPredictionFMD(:,2), 'k');
  set(handles.axsTest, 'NextPlot', 'replace');
  axes(handles.axsLikelihood);
  pl1 = plot(vX,cumsum((vProbK)),'k','LineWidth',2.0);
  set(handles.axsLikelihood, 'NextPlot', 'add');
  pl1 = plot(vX,cumsum((vProbO)),'r','LineWidth',2.0);
  set(handles.axsLikelihood, 'NextPlot', 'replace');
  vXLim = get(handles.axsTest, 'XLim');
  set(handles.axsLikelihood, 'XLim', vXLim);

  % Normalize the two subcatalogs in time
  fMinTime = min(mFirst(:,3));
  mFirst(:,3) = mFirst(:,3) - fMinTime;
  fMinTime = min(mSecond(:,3));
  mSecond(:,3) = mSecond(:,3) - fMinTime;
  % Plot the cumulative number for the first period
  [~, vIndices] = sort(mFirst(:,3));
  mSortedCatalog = mFirst(vIndices(:,1),:) ;
  axes(handles.axsTimePlot)
  plot(mSortedCatalog(:,3),(1:length(mSortedCatalog(:,3))),'k-');
  % Plot the cumulative number for the second period
  [~, vIndices] = sort(mSecond(:,3));
  mSortedCatalog = mSecond(vIndices(:,1),:) ;
  axes(handles.axsTimePlot)
  set(handles.axsTimePlot, 'NextPlot', 'add');
  plot(mSortedCatalog(:,3),(1:length(mSortedCatalog(:,3))),'r:');
end

% --------------------------------------------------------------------
function varargout = radEntireCatalog_Callback(h, eventdata, handles, varargin)

if get(handles.radEntireCatalog, 'Value') == 1
  set(handles.radPickLocations, 'Value', 0);
  set(handles.btnPick, 'Enable', 'off');
else
  set(handles.radPickLocations, 'Value', 1);
  set(handles.btnPick, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = radPickLocations_Callback(h, eventdata, handles, varargin)

if get(handles.radPickLocations, 'Value') == 1
  set(handles.radEntireCatalog, 'Value', 0);
  set(handles.btnPick, 'Enable', 'on');
else
  set(handles.radEntireCatalog, 'Value', 1);
  set(handles.btnPick, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = radEntireTime_Callback(h, eventdata, handles, varargin)

if get(handles.radEntireTime, 'Value') == 1
  set(handles.radSplitPeriods, 'Value', 0);
  set(handles.txtFirst, 'Enable', 'off');
  set(handles.txtSecond, 'Enable', 'off');
  set(handles.txtSplit, 'Enable', 'off');
else
  set(handles.radSplitPeriods, 'Value', 1);
  set(handles.txtFirst, 'Enable', 'on');
  set(handles.txtSecond, 'Enable', 'on');
  set(handles.txtSplit, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = radSplitPeriods_Callback(h, eventdata, handles, varargin)

if get(handles.radSplitPeriods, 'Value') == 1
  set(handles.radEntireTime, 'Value', 0);
  set(handles.txtFirst, 'Enable', 'on');
  set(handles.txtSecond, 'Enable', 'on');
  set(handles.txtSplit, 'Enable', 'on');
else
  set(handles.radEntireTime, 'Value', 1);
  set(handles.txtFirst, 'Enable', 'off');
  set(handles.txtSecond, 'Enable', 'off');
  set(handles.txtSplit, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = btnRePlot_Callback(h, eventdata, handles, varargin)

plotit(handles);

% --------------------------------------------------------------------
function varargout = btnClose_Callback(h, eventdata, handles, varargin)

% Close the dialog
delete(handles.dlgAnalyze);

% --------------------------------------------------------------------
function varargout = btnPick_Callback(h, eventdata, handles, varargin)

% Get the axes handle of the plotwindow
axes(kj_result('GetAxesHandle', handles.hParentFigure, [], guidata(handles.hParentFigure)));

% Select a point in the plot window with the mouse
[fX, fY] = ginput(1);

% Prepare the data
% If cross-section calculate the lenght along cross-section
if ~handles.params.bMap
  [nRow, nColumn] = size(handles.params.mCatalog);
  xsecx2 = handles.params.mCatalog(:,nColumn);  % length along x-section
  xsecy2 = handles.params.mCatalog.Depth;        % depth of hypocenters
end

% Calculate distance from center point and sort with distance
if handles.params.bMap
  vDistances = sqrt(((handles.params.mCatalog.Longitude-fX)*cosd(fY)*111).^2 + ((handles.params.mCatalog.Latitude-fY)*111).^2);
else
  vDistances = sqrt(((xsecx2 - fX)).^2 + ((xsecy2 + fY)).^2);
end

% Create the subcatalog
if handles.params.bNumber
  % Use first nNumberEvents events
  [vTmp, vIndices] = sort(vDistances);
  mNodeCatalog = handles.params.mCatalog.subset(vIndices);
  mNodeCatalog = mNodeCatalog(1:handles.params.nNumberEvents);
else
  % Use all events within fRadius
  vSelection = (vDistances <= handles.params.fRadius);
  mNodeCatalog = handles.params.mCatalog.subset(vSelection);
end

% Store the subcatalog for later use
handles.mNodeCatalog = mNodeCatalog;
guidata(handles.dlgAnalyze, handles);

% Replot the frequency magnitude distribution and cumulative number
plotit(handles);

author: CelsoReyes


text:
function tagS = content_identification_macro_tags
%"content_identification_macro_tags"
%   Return the tags used to represent content as specified by C.10.9 in
%   PS3.3 of 2006 DICOM specification. Tags are returned in a struct array
%   with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%JRA 06/06/06
%YWU Modified 03/01/08
%
%Usage:
%   tagS = patient_module_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

tagC = {'00200013','00700080','00700081','00700084'};

tagDecC = {2097171, 7340160, 7340161, 7340164};

tagTypeC = {'1','1','2','2'};

childC = {[],[],[],[]};

tagS = struct('tag', tagC, 'tagdec', tagDecC, 'type', tagTypeC, 'children', childC);

%
% %Initialize the tagS structure.
% tagS = struct('tag', {}, 'type', {}, 'children', {});
%
% %Create an empty tagS template for sequence creation.
% template = tagS;
%
% %Add tags based on PS3.3 attribute lists.
%
% %Instance Number
% tagS(end+1) = struct('tag', ['00200013'], 'type', ['1'], 'children', []);
%
% %Content Label
% tagS(end+1) = struct('tag', ['00700080'], 'type', ['1'], 'children', []);
%
% %Content Description
% tagS(end+1) = struct('tag', ['00700081'], 'type', ['2'], 'children', []);
%
% %Content Creator's Name
% tagS(end+1) = struct('tag', ['00700084'], 'type', ['2'], 'children', []);
%

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function chrpak_test001 ( )

%*****************************************************************************80
%
%% TEST001 tests A_TO_I4 and I4_TO_A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 January 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST001\n' );
  fprintf ( 1, '  A_TO_I4: Alphabetic character => I\n' );
  fprintf ( 1, '  I4_TO_A: I => Alphabetic character\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  1:26 = A:Z\n' );
  fprintf ( 1, '  27:52 = a:z\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I  ==>  A  ==>  I\n' );
  fprintf ( 1, '\n' );

  for i4 = 0 : 3 : 55

    a = i4_to_a ( i4 );
    j4 = a_to_i4 ( a );

    fprintf ( 1, '  %8d  "%c"  %8d\n', i4, a, j4 );

  end

  return
end


text:
function xyl = Xl(s,row,b0,b1)

if nargin == 1
    row = 2. ;
    b0 = 1. ;
    b1 = 2. ;
end

x = 0 ;

y = b0+(b1-b0)*s ;

xyl = [x ; y] ;

author: Sable

repo: mcbench-benchmarks


text:
% Fig. 6.29   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=[1 1];
den=conv([1 0],[0.1 -1]);
rlocus(num,den);

axis([-8 12 -10 10]);
axis equal
grid;
title('Fig. 6.29 Root locus for Example 6.10');


text:
function spm_eeg_inv_checkforward(varargin)
% Check M/EEG forward model
% FORMAT spm_eeg_inv_checkforward(D, val, ind)
%__________________________________________________________________________
% Copyright (C) 2008-2018 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_inv_checkforward.m 7702 2019-11-22 11:32:26Z guillaume $

%-SPM data structure
%--------------------------------------------------------------------------
[D, val] = spm_eeg_inv_check(varargin{:});

forward = D.inv{val}.forward;

if nargin < 3
    str = sprintf('%s|', forward(:).modality);
    str = str(1:(end-1));

    ind = spm_input('What to display?','+1','b',str,1:numel(forward),1);
else
    ind = varargin{3};
end

try
    vol      = forward(ind).vol;
    modality = forward(ind).modality;
    if isfield(forward(ind), 'siunits') && forward(ind).siunits
        sens     = forward(ind).sensors;
    else %backward compatibility
        sens     = D.inv{val}.datareg(ind).sensors;
    end
    Mcortex  = forward(ind).mesh;
catch
    spm('alert!','Please coregister and create forward model',mfilename);
    return
end

%-Display
%--------------------------------------------------------------------------
Fgraph  = spm_figure('GetWin','Graphics');
spm_figure('Focus',Fgraph);
spm_figure('Clear',Fgraph);

spm('Pointer', 'Watch');

chanind = strmatch(modality, D.chantype);
chanind = setdiff(chanind, D.badchannels);
if isempty(chanind)
    error(['No good ' modality ' channels were found.']);
end

if ischar(vol)
    vol = ft_read_headmodel(vol);
end

face    = Mcortex.face;
vert    = Mcortex.vert;
h_ctx   = patch('vertices',vert,'faces',face,'EdgeColor','b','FaceColor','b');

hold on

[volp, sens] = ft_prepare_vol_sens(vol, sens, 'channel', D.chanlabels(chanind));

ft_plot_headmodel(vol, 'edgecolor', [0 0 0], 'facealpha', 0);

hold on

try
    ft_plot_sens(sens, 'style', '*', 'edgecolor', 'g', 'elecsize', 20, 'coil', ft_senstype(sens, 'eeg'));
catch
    ft_plot_sens(sens, 'edgecolor', 'g', 'coilshape', 'point', 'elecsize', 20, 'coil', true);
end

rotate3d on;

axis off
axis vis3d
axis equal

spm('Pointer', 'Arrow');

author: spm


text:
function [DEM] = spm_LAP(DEM)
% Laplacian model inversion (see also spm_LAPS)
% FORMAT DEM   = spm_LAP(DEM)
%
% DEM.M  - hierarchical model
% DEM.Y  - response variable, output or data
% DEM.U  - explanatory variables, inputs or prior expectation of causes
%__________________________________________________________________________
%
% generative model
%--------------------------------------------------------------------------
%   M(i).g  = v     =  g(x,v,P)   {inline function, string or m-file}
%   M(i).f  = dx/dt =  f(x,v,P)   {inline function, string or m-file}
%
%   M(i).ph = pi(v) = ph(x,v,h,M) {inline function, string or m-file}
%   M(i).pg = pi(x) = pg(x,v,g,M) {inline function, string or m-file}
%                                  (assumed to be linear in v and x)
%
%   pi(v,x) = vectors of log-precisions; (h,g) = precision parameters
%
%   M(i).pE = prior expectation of p model-parameters
%   M(i).pC = prior covariances of p model-parameters
%   M(i).hE = prior expectation of h log-precision (cause noise)
%   M(i).hC = prior covariances of h log-precision (cause noise)
%   M(i).gE = prior expectation of g log-precision (state noise)
%   M(i).gC = prior covariances of g log-precision (state noise)
%   M(i).xP = precision (states)
%   M(i).Q  = precision components (input noise)
%   M(i).R  = precision components (state noise)
%   M(i).V  = fixed precision (input noise)
%   M(i).W  = fixed precision (state noise)
%
%   M(i).P  = optional initial value for parameters (defaults to M(i).pE)
%
%   M(i).m  = number of inputs v(i + 1);
%   M(i).n  = number of states x(i);
%   M(i).l  = number of output v(i);
%
% conditional moments of model-states - q(u)
%--------------------------------------------------------------------------
%   qU.x    = Conditional expectation of hidden states
%   qU.v    = Conditional expectation of causal states
%   qU.w    = Conditional prediction error (states)
%   qU.z    = Conditional prediction error (causes)
%   qU.C    = Conditional covariance: cov(v)
%   qU.S    = Conditional covariance: cov(x)
%
% conditional moments of model-parameters - q(p)
%--------------------------------------------------------------------------
%   qP.P    = Conditional expectation
%   qP.C    = Conditional covariance
%
% conditional moments of hyper-parameters (log-transformed) - q(h)
%--------------------------------------------------------------------------
%   qH.h    = Conditional expectation (cause noise)
%   qH.g    = Conditional expectation (state noise)
%   qH.C    = Conditional covariance
%
% F         = log-evidence = log-marginal likelihood = negative free-energy
%
%__________________________________________________________________________
% Accelerated methods: To accelerate computations one can specify the
% nature of the model equations using:
%
% M(1).E.linear = 0: full        - evaluates 1st and 2nd derivatives
% M(1).E.linear = 1: linear      - equations are linear in x and v
% M(1).E.linear = 2: bilinear    - equations are linear in x, v & x*v
% M(1).E.linear = 3: nonlinear   - equations are linear in x, v, x*v, & x*x
% M(1).E.linear = 4: full linear - evaluates 1st derivatives (for GF)
%
% similarly, for evaluating precisions:
%
% M(1).E.method.h = 0,1  switch for precision parameters (hidden causes)
% M(1).E.method.g = 0,1  switch for precision parameters (hidden states)
% M(1).E.method.v = 0,1  switch for precision (hidden causes)
% M(1).E.method.x = 0,1  switch for precision (hidden states)
%__________________________________________________________________________
%
% spm_LAP implements a variational scheme under the Laplace
% approximation to the conditional joint density q on states u, parameters
% p and hyperparameters (h,g) of an analytic nonlinear hierarchical dynamic
% model, with additive Gaussian innovations.
%
%            q(u,p,h,g) = max <L(t)>q
%
% L is the ln p(y,u,p,h,g|M) under the model M. The conditional covariances
% obtain analytically from the curvature of L with respect to the unknowns.
%__________________________________________________________________________
% Copyright (C) 2010-2013 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_LAP.m 6508 2015-07-25 15:23:25Z karl $

% find or create a DEM figure
%--------------------------------------------------------------------------
try
    DEM.M(1).nograph;
catch
    DEM.M(1).nograph = 0;
end
if ~DEM.M(1).nograph
    Fdem = spm_figure('GetWin','DEM');
end

% check model, data and priors
%==========================================================================
[M,Y,U] = spm_DEM_set(DEM);

% set regularisation
%--------------------------------------------------------------------------
try
    dt = DEM.M(1).E.v;
catch
    dt = 0;
    DEM.M(1).E.v = dt;
end

% number of iterations
%--------------------------------------------------------------------------
try, nD = M(1).E.nD; catch, nD = 1;  end
try, nN = M(1).E.nN; catch, nN = 8;  end

% ensure integration scheme evaluates gradients at each time-step
%--------------------------------------------------------------------------
M(1).E.linear = 4;

% assume precisions are a function of, and only of, hyperparameters
%--------------------------------------------------------------------------
try
    method = M(1).E.method;
catch
    method.h = 1;
    method.g = 1;
    method.x = 0;
    method.v = 0;
end
try method.h; catch, method.h = 0; end
try method.g; catch, method.g = 0; end
try method.x; catch, method.x = 0; end
try method.v; catch, method.v = 0; end

% assume precisions are a function of, and only of, hyperparameters
%--------------------------------------------------------------------------
try
    form = M(1).E.form;
catch
    form = 'Gaussian';
end

% checks for Laplace models (precision functions; ph and pg)
%--------------------------------------------------------------------------
for i  = 1:length(M)
    try
        feval(M(i).ph,M(i).x,M(i + 1).v,M(i).hE,M(i)); method.v = 1;
    catch
        M(i).ph = inline('spm_LAP_ph(x,v,h,M)','x','v','h','M');
    end
    try
        feval(M(i).pg,M(i).x,M(i + 1).v,M(i).gE,M(i)); method.x = 1;
    catch
        M(i).pg = inline('spm_LAP_pg(x,v,h,M)','x','v','h','M');
    end
end

M(1).E.method = method;

% order parameters (d = n = 1 for static models) and checks
%==========================================================================
d   = M(1).E.d + 1;                          % embedding order of q(v)
n   = M(1).E.n + 1;                          % embedding order of q(x)

% number of states and parameters
%--------------------------------------------------------------------------
ns  = size(Y,2);                             % number of samples
nl  = size(M,2);                             % number of levels
nv  = sum(spm_vec(M.m));                     % number of v (casual states)
nx  = sum(spm_vec(M.n));                     % number of x (hidden states)
ny  = M(1).l;                                % number of y (inputs)
nc  = M(end).l;                              % number of c (prior causes)
nu  = nv*d + nx*n;                           % number of generalised states
ne  = nv*n + nx*n + ny*n;                    % number of generalised errors

% precision (R) of generalised errors and null matrices for concatenation
%==========================================================================
s     = M(1).E.s;
Rh    = spm_DEM_R(n,s,form);
Rg    = spm_DEM_R(n,s,form);

if ~nx, Rg = sparse(0,0); end

W     = sparse(nx*n,nx*n);
V     = sparse((ny + nv)*n,(ny + nv)*n);

% fixed priors on states (u)
%--------------------------------------------------------------------------
Px    = kron(spm_DEM_R(n,2),spm_cat(spm_diag({M(1:end).xP})));
Pv    = kron(spm_DEM_R(d,2),spm_cat(spm_diag({M(2:end).vP})));
Pu    = spm_cat(spm_diag({Px Pv}));

% hyperpriors
%--------------------------------------------------------------------------
ph.h  = spm_vec({M.hE M.gE});                % prior expectation of h,g
ph.c  = spm_cat(spm_diag({M.hC M.gC}));      % prior covariances of h,g
Ph    = spm_inv(ph.c);                       % prior precision of h,g

qh.h  = {M.hE};                              % conditional expectation h
qh.g  = {M.gE};                              % conditional expectation g
nh    = length(spm_vec(qh.h));               % number of hyperparameters h
ng    = length(spm_vec(qh.g));               % number of hyperparameters g
nb    = nh + ng;                             % number of hyperparameters

% priors on parameters (in reduced parameter space)
%==========================================================================
pp.c  = cell(nl,nl);
qp.p  = cell(nl,1);
for i = 1:(nl - 1)

    % eigenvector reduction: p <- pE + qp.u*qp.p
    %----------------------------------------------------------------------
    qp.u{i}   = spm_svd(M(i).pC,0);          % basis for parameters
    M(i).p    = size(qp.u{i},2);             % number of qp.p
    qp.p{i}   = sparse(M(i).p,1);            % initial deviates
    pp.c{i,i} = qp.u{i}'*M(i).pC*qp.u{i};    % prior covariance

end
Up    = spm_cat(spm_diag(qp.u));

% priors on parameters
%--------------------------------------------------------------------------
pp.p  = spm_vec(M.pE);
pp.c  = spm_cat(pp.c);
Pp    = spm_inv(pp.c);

% initialise conditional density q(p)
%--------------------------------------------------------------------------
for i = 1:(nl - 1)
    try
        qp.p{i} = qp.p{i} + qp.u{i}'*(spm_vec(M(i).P) - spm_vec(M(i).pE));
    end
end
np    = size(Up,2);

% initialise cell arrays for D-Step; e{i + 1} = (d/dt)^i[e] = e[i]
%==========================================================================
qu.x      = cell(n,1);
qu.v      = cell(n,1);
qu.y      = cell(n,1);
qu.u      = cell(n,1);
[qu.x{:}] = deal(sparse(nx,1));
[qu.v{:}] = deal(sparse(nv,1));
[qu.y{:}] = deal(sparse(ny,1));
[qu.u{:}] = deal(sparse(nc,1));

% initialise cell arrays for hierarchical structure of x[0] and v[0]
%--------------------------------------------------------------------------
x         = {M(1:end - 1).x};
v         = {M(1 + 1:end).v};
qu.x{1}   = spm_vec(x);
qu.v{1}   = spm_vec(v);

% derivatives for Jacobian of D-step
%--------------------------------------------------------------------------
Dx     = kron(spm_speye(n,n,1),spm_speye(nx,nx));
Dv     = kron(spm_speye(d,d,1),spm_speye(nv,nv));
Dy     = kron(spm_speye(n,n,1),spm_speye(ny,ny));
Dc     = kron(spm_speye(d,d,1),spm_speye(nc,nc));
Du     = spm_cat(spm_diag({Dx,Dv}));
Ib     = spm_speye(np + nb,np + nb);
dbdt   = sparse(np + nb,1);

% gradients of generalised weighted errors
%--------------------------------------------------------------------------
dedh   = sparse(nh,ne);
dedg   = sparse(ng,ne);
dedv   = sparse(nv,ne);
dedx   = sparse(nx,ne);
dedhh  = sparse(nh,nh);
dedgg  = sparse(ng,ng);

% curvatures of Gibb's energy w.r.t. hyperparameters
%--------------------------------------------------------------------------
dHdh   = sparse(nh,1);
dHdg   = sparse(ng,1);
dHdp   = sparse(np,1);
dHdu   = sparse(nu,1);

% preclude unnecessary iterations and set switches
%--------------------------------------------------------------------------
if ~np && ~nh && ~ng, nN = 1; end
mnx    = nx*~~method.x;
mnv    = nv*~~method.v;

% preclude very precise states from entering free-energy/action
%--------------------------------------------------------------------------
p      = spm_LAP_eval(M,qu,qh);
ih     = p.h < 8;
ig     = p.g < 8;
ie     = kron(ones(n,1),ih);
ix     = kron(ones(n,1),ig);
iv     = kron(ones(d,1),ih((1:nv) + ny));
je     = find([ie; ix]); ix(1:nx) = 1;
ju     = find([ix; iv]);

% and other useful indices
%--------------------------------------------------------------------------
ix     = (1:nx);
ih     = (1:nb);
iv     = (1:nv) + nx*n;
ip     = (1:np) + nu;

% number of iterations for convergence
%--------------------------------------------------------------------------
convergence = -4;

% Iterate Laplace scheme
%==========================================================================
F      = -Inf;
for iN = 1:nN

    % get time and clear persistent variables in evaluation routines
    %----------------------------------------------------------------------
    tic; clear spm_DEM_eval

    % [re-]set states & their derivatives
    %----------------------------------------------------------------------
    try, qu = Q(1).u; end

    % D-Step: (nD D-Steps for each sample)
    %======================================================================
    for is = 1:ns

        % D-Step: until convergence for static systems
        %==================================================================
        for iD = 1:nD

            % sampling time
            %--------------------------------------------------------------
            ts = is + (iD - 1)/nD;

            % derivatives of responses and inputs
            %--------------------------------------------------------------
            try
                qu.y(1:n) = spm_DEM_embed(Y,n,ts,1,M(1).delays);
                qu.u(1:d) = spm_DEM_embed(U,d,ts);
            catch
                qu.y(1:n) = spm_DEM_embed(Y,n,ts);
                qu.u(1:d) = spm_DEM_embed(U,d,ts);
            end

            % evaluate functions and derivatives
            %==============================================================

            % prediction errors (E) and precision vectors (p)
            %--------------------------------------------------------------
            [E,dE] = spm_DEM_eval(M,qu,qp);
            [p,dp] = spm_LAP_eval(M,qu,qh);

            % gradients of log(det(iS)) dDd...
            %==============================================================

            % get precision matrices
            %--------------------------------------------------------------
            iSh   = diag(exp(p.h));
            iSg   = diag(exp(p.g));
            iS    = blkdiag(kron(Rh,iSh),kron(Rg,iSg));

            % gradients of trace(diag(p)) = sum(p); p = precision vector
            %--------------------------------------------------------------
            dpdx  = n*sum(spm_cat({dp.h.dx; dp.g.dx}));
            dpdv  = n*sum(spm_cat({dp.h.dv; dp.g.dv}));
            dpdh  = n*sum(dp.h.dh);
            dpdg  = n*sum(dp.g.dg);
            dpdx  = kron(sparse(1,1,1,1,n),dpdx);
            dpdv  = kron(sparse(1,1,1,1,d),dpdv);
            dDdu  = [dpdx dpdv]';
            dDdh  = [dpdh dpdg]';

            % gradients precision-weighted generalised error dSd..
            %==============================================================

            % gradients w.r.t. hyperparameters
            %--------------------------------------------------------------
            for i = 1:nh
                diS       = diag(dp.h.dh(:,i).*exp(p.h));
                diSdh{i}  = blkdiag(kron(Rh,diS),W);
                dedh(i,:) = E'*diSdh{i};
            end
            for i = 1:ng
                diS       = diag(dp.g.dg(:,i).*exp(p.g));
                diSdg{i}  = blkdiag(V,kron(Rg,diS));
                dedg(i,:) = E'*diSdg{i};
            end

            % gradients w.r.t. hidden states
            %--------------------------------------------------------------
            for i = 1:mnx
                diV       = diag(dp.h.dx(:,i).*exp(p.h));
                diW       = diag(dp.g.dx(:,i).*exp(p.g));
                diSdx{i}  = blkdiag(kron(Rh,diV),kron(Rg,diW));
                dedx(i,:) = E'*diSdx{i};
            end

            % gradients w.r.t. causal states
            %--------------------------------------------------------------
            for i = 1:mnv
                diV       = diag(dp.h.dv(:,i).*exp(p.h));
                diW       = diag(dp.g.dv(:,i).*exp(p.g));
                diSdv{i}  = blkdiag(kron(Rh,diV),kron(Rg,diW));
                dedv(i,:) = E'*diSdv{i};
            end

            dSdx  = kron(sparse(1,1,1,n,1),dedx);
            dSdv  = kron(sparse(1,1,1,d,1),dedv);
            dSdu  = [dSdx; dSdv];
            dEdh  = [dedh; dedg];
            dEdp  = dE.dp'*iS;
            dEdu  = dE.du'*iS;

            % curvatures w.r.t. hyperparameters
            %--------------------------------------------------------------
            for i = 1:nh
                for j = i:nh
                    diS        = diag(dp.h.dh(:,i).*dp.h.dh(:,j).*exp(p.h));
                    diS        = blkdiag(kron(Rh,diS),W);
                    dedhh(i,j) = E'*diS*E;
                    dedhh(j,i) = dedhh(i,j);
                end
            end
            for i = 1:ng
                for j = i:ng
                    diS        = diag(dp.g.dg(:,i).*dp.g.dg(:,j).*exp(p.g));
                    diS        = blkdiag(V,kron(Rg,diS));
                    dedgg(i,j) = E'*diS*E;
                    dedgg(j,i) = dedgg(i,j);
                end
            end

            % combined curvature
            %--------------------------------------------------------------
            dSdhh = spm_cat({dedhh  [] ;
                             [] dedgg});

            % errors (from prior expectations) (NB pp.p = 0)
            %--------------------------------------------------------------
            Eu    = spm_vec(qu.x(1:n),qu.v(1:d));
            Ep    = spm_vec(qp.p);
            Eh    = spm_vec(qh.h,qh.g) - ph.h;

            % first-order derivatives of Gibb's Energy
            %==============================================================
            dLdu  = dEdu*E + dSdu*E/2 - dDdu/2 + Pu*Eu;
            dLdh  = dEdh*E/2          - dDdh/2 + Ph*Eh;
            dLdp  = dEdp*E                     + Pp*Ep;

            % and second-order derivatives of Gibb's Energy
            %--------------------------------------------------------------
            dLduu = dEdu*dE.du + Pu;
            dLdpp = dEdp*dE.dp + Pp;
            dLdhh = dSdhh/2    + Ph;
            dLduy = dEdu*dE.dy;
            dLduc = dEdu*dE.dc;
            dLdup = dEdu*dE.dp;
            dLdhp = dEdh*dE.dp;
            dLdpu = dLdup';
            dLdph = dLdhp';

            % precision and covariances for entropy
            %--------------------------------------------------------------
            dLdaa = spm_cat({dLduu dLdup  ;
                             dLdpu dLdpp});
            dLdbb = spm_cat({dLdpp dLdph  ;
                             dLdhp dLdhh});

            Cup   = spm_inv(dLdaa);
            Chh   = spm_inv(dLdhh);

            % first-order derivatives of Entropy term
            %==============================================================

            % log-precision
            %--------------------------------------------------------------
            for i = 1:nh
                Luub    = dE.du'*diSdh{i}*dE.du;
                Lpub    = dE.dp'*diSdh{i}*dE.du;
                Lppb    = dE.dp'*diSdh{i}*dE.dp;
                diCdh   = spm_cat({Luub Lpub';
                                   Lpub Lppb});
                dHdh(i) = spm_trace(diCdh,Cup)/2;
            end
            for i = 1:ng
                Luub    = dE.du'*diSdg{i}*dE.du;
                Lpub    = dE.dp'*diSdg{i}*dE.du;
                Lppb    = dE.dp'*diSdg{i}*dE.dp;
                diCdg   = spm_cat({Luub Lpub';
                                   Lpub Lppb});
                dHdg(i) = spm_trace(diCdg,Cup)/2;
            end

            % parameters
            %--------------------------------------------------------------
            for i = 1:np
                Luup    = dE.dup{i}'*dEdu';
                Lpup    = dEdp*dE.dup{i};
                Luup    = Luup + Luup';
                diCdp   = spm_cat({Luup Lpup';
                                   Lpup [] });
                dHdp(i) = spm_trace(diCdp,Cup)/2;
            end

            % hidden states and causes (disabled for stability)
            %--------------------------------------------------------------
            for i = 1:(nu - nu)
                Lppu    = dE.dpu{i}'*dEdp';
                Lupu    = dEdu*dE.dpu{i};
                Lppu    = Lppu + Lppu';
                diCdu   = spm_cat({[]    Lupu;
                                   Lupu' Lppu});
                dHdu(i) = spm_trace(diCdu,Cup)/2;
            end

            % and concatenate
            %--------------------------------------------------------------
            dHdb  = [dHdh; dHdg];
            dHdb  = [dHdp; dHdb];
            dLdb  = [dLdp; dLdh];

            % save conditional moments (and prediction error) at Q{t}
            %==============================================================
            if iD == 1 || ns == 1

                % save means
                %----------------------------------------------------------
                Q(is).e   = E;
                Q(is).E   = diag(diag(iS))*E;
                Q(is).u   = qu;
                Q(is).p   = qp;
                Q(is).h   = qh;

                % and conditional covariances
                %----------------------------------------------------------
                Q(is).u.s = Cup(ix,ix);
                Q(is).u.c = Cup(iv,iv);
                Q(is).p.c = Cup(ip,ip);
                Q(is).h.c = Chh(ih,ih);

                % Free-energy (components)
                %----------------------------------------------------------
                Fc(is,1)  = - E(je)'*iS(je,je)*E(je)/2;
                Fc(is,2)  = - Eu(ju)'*Pu(ju,ju)*Eu(ju)/2;
                Fc(is,3)  = - n*ny*log(2*pi)/2;
                Fc(is,4)  = spm_logdet(iS(je,je))/2;
                Fc(is,5)  = spm_logdet(Pu(ju,ju)*Cup(ju,ju))/2;

                % Free-action (states and parameters)
                %----------------------------------------------------------
                AC(is)    = sum(Fc(is,:))    ...
                          - Ep'*Pp*Ep/2      ...
                          - Eh'*Ph*Eh/2      ...
                          + spm_logdet(Pp)/2 ...
                          + spm_logdet(Ph)/2 ...
                          - spm_logdet(dLdbb)/2;

            end

            % update conditional moments
            %==============================================================

            % prior precision of fluctuations on [hyper] parameters
            %--------------------------------------------------------------
            Kb    = ns*Ib;

            % accumulate curvatures of [hyper] parameters
            %--------------------------------------------------------------
            try
                dLdBB = dLdBB*(1 - 1/ns) + dLdbb/ns;
            catch
                dLdBB = dLdbb + Ib*32;
            end

            % whiten gradient (and curvatures) with regularised precision
            %--------------------------------------------------------------
            Cb    = spm_inv(dLdBB + diag(diag(dLdBB))*exp(dt));
            dLdb  = Cb*dLdb;
            dHdb  = Cb*dHdb;

            % assemble conditional means
            %--------------------------------------------------------------
            q.y  = qu.y(1:n);
            q.x  = qu.x(1:n);
            q.v  = qu.v(1:d);
            q.c  = qu.u(1:d);
            q.p  = qp.p;
            q.h  = qh.h;
            q.g  = qh.g;
            q.d  = dbdt;

            % flow
            %--------------------------------------------------------------
            f.y  =  Dy*spm_vec(q.y)                  ;
            f.u  =  Du*spm_vec(q.x,q.v) - dLdu - dHdu;
            f.c  =  Dc*spm_vec(q.c)                  ;
            f.b  =     spm_vec(q.d)                  ;
            f.d  = -Kb*spm_vec(q.d)     - dLdb - dHdb;

            % and Jacobian
            %--------------------------------------------------------------
            dfdq = {Dy     []        []     []    [] ;
                   -dLduy  Du-dLduu -dLduc  []    [] ;
                    []     []        Dc     []    [] ;
                    []     []        []     []    Ib ;
                    []     []        []    -Ib   -Kb};

            dfdq = spm_cat(dfdq);

            % get eigenvalues of Jacobian (on hidden states)
            %--------------------------------------------------------------
            try DEM.options.eigenvalues;
                DEM.E(:,is) = eig(full(Du - dLduu));
            end

            % update conditional modes of states
            %==============================================================
            dq   = spm_dx(dfdq,spm_vec(f),1/nD);
            q    = spm_unvec(spm_vec(q) + dq,q);

            % unpack conditional means
            %--------------------------------------------------------------
            qu.x(1:n) = q.x;
            qu.v(1:d) = q.v;
            qp.p      = q.p;
            qh.h      = q.h;
            qh.g      = q.g;
            dbdt      = q.d;

        end % D-Step

    end % sequence (ns)

    % Bayesian parameter averaging
    %======================================================================

    % Conditional moments of time-averaged parameters
    %----------------------------------------------------------------------
    Ep    = 0;
    Qp    = 0;
    for i = 1:ns
        P   = spm_inv(Q(i).p.c);
        Ep  = Ep  + P*spm_vec(Q(i).p.p);
        Qp  = Qp + P;
    end
    Ep    = spm_inv(Qp)*Ep;
    Cp    = spm_inv(Qp + (1 - ns)*Pp);

    % conditional moments of hyper-parameters
    %----------------------------------------------------------------------
    Eh    = 0;
    Qh    = 0;
    for i = 1:ns
        P   = spm_inv(Q(i).h.c);
        Eh  = Eh  + P*spm_vec(Q(i).h.h,Q(i).h.g);
        Qh  = Qh + P;
    end
    Eh    = spm_inv(Qh)*Eh - ph.h;
    Ch    = spm_inv(Qh + (1 - ns)*Ph);

    % Free-action of states plus free-energy of parameters
    %======================================================================
    FT    = sum(Fc,2);          %   instantaneous Free energy (of states)
    FC(1) = sum(Fc(:,1));       % - E'*iS*E/2;
    FC(2) = sum(Fc(:,2));       % - Eu'*Pu*Eu/2;
    FC(3) = sum(Fc(:,3));       % - n*ny*log(2*pi)/2;
    FC(4) = sum(Fc(:,4));       %   spm_logdet(iS)/2;
    FC(5) = sum(Fc(:,5));       %   spm_logdet(Pu*Cu)/2;
    FC(6) = -Ep'*Pp*Ep/2;
    FC(7) = -Eh'*Ph*Eh/2;
    FC(8) = spm_logdet(Pp*Cp)/2;
    FC(9) = spm_logdet(Ph*Ch)/2;

    CC(iN,:) = FC;
    S(iN)    = sum(AC);
    Fe       = sum(FC);

    % if F is decreasing, revert [hyper] parameters and slow down
    %----------------------------------------------------------------------
    if Fe < F(iN) && iN > 4

        % save free-energy
        %------------------------------------------------------------------
        F(iN + 1) = F(iN);

        % load current MAP estimates
        %------------------------------------------------------------------
        qp = PQ.qp;
        qh = PQ.qh;

        % decrease update time
        %------------------------------------------------------------------
        dt = max(dt + 2,2);

        % convergence
        %------------------------------------------------------------------
        if dt > 6; convergence = 1; end

    else

        % convergence
        %------------------------------------------------------------------
        if Fe - F(iN) < 1, convergence = convergence + 1; end

        % save free-energy
        %------------------------------------------------------------------
        F(iN)     = Fe;
        F(iN + 1) = Fe;

        % save current MAP estimates
        %------------------------------------------------------------------
        PQ.qp = qp;
        PQ.qh = qh;

        % increase update time
        %------------------------------------------------------------------
        dt    = max(dt - 1,-8);

    end

    % Convergence
    %======================================================================
    if convergence > 0; break, end

    % otherwise save conditional moments (for each time point)
    %======================================================================
    for t = 1:length(Q)

        % states and predictions
        %------------------------------------------------------------------
        v     = spm_unvec(Q(t).u.v{1},v);
        x     = spm_unvec(Q(t).u.x{1},x);
        z     = spm_unvec(Q(t).e(1:(ny + nv)),{M.v});
        Z     = spm_unvec(Q(t).E(1:(ny + nv)),{M.v});
        w     = spm_unvec(Q(t).e((1:nx) + (ny + nv)*n),{M.x});
        X     = spm_unvec(Q(t).E((1:nx) + (ny + nv)*n),{M.x});
        for i = 1:(nl - 1)
            if M(i).m, qU.v{i + 1}(:,t) = spm_vec(v{i});  end
            if M(i).n, qU.x{i}(:,t)     = spm_vec(x{i});  end
            if M(i).n, qU.w{i}(:,t)     = spm_vec(w{i});  end
            if M(i).l, qU.z{i}(:,t)     = spm_vec(z{i});  end
            if M(i).n, qU.W{i}(:,t)     = spm_vec(X{i});  end
            if M(i).l, qU.Z{i}(:,t)     = spm_vec(Z{i});  end
        end
        if    M(nl).l, qU.z{nl}(:,t)    = spm_vec(z{nl}); end
        if    M(nl).l, qU.Z{nl}(:,t)    = spm_vec(Z{nl}); end

        qU.v{1}(:,t)  = spm_vec(Q(t).u.y{1}) - spm_vec(z{1});

        % and conditional covariances
        %------------------------------------------------------------------
        qU.S{t} = Q(t).u.s;
        qU.C{t} = Q(t).u.c;

        % parameters
        %------------------------------------------------------------------
        qP.p{t} = spm_vec(Q(t).p.p);
        qP.c{t} = Q(t).p.c;

        % hyperparameters
        %------------------------------------------------------------------
        qH.p{t} = spm_vec(Q(t).h.h,Q(t).h.g);
        qH.c{t} = Q(t).h.c;

    end

    % graphics (states)
    %----------------------------------------------------------------------
    if exist('Fdem','var')
        spm_figure('Select', Fdem)
        spm_DEM_qU(qU)

        % graphics (parameters and log-precisions)
        %------------------------------------------------------------------
        if np && nb
            subplot(2*nl,2,4*nl - 2)
            plot(1:ns,full(spm_cat(qP.p)))
            set(gca,'XLim',[1 ns])
            title('parameters (modes)','FontSize',16)

            subplot(2*nl,2,4*nl)
            plot(1:ns,full(spm_cat(qH.p)))
            set(gca,'XLim',[1 ns])
            title('log-precision','FontSize',16)

        elseif nb
            subplot(nl,2,2*nl)
            plot(1:ns,full(spm_cat(qH.p)))
            set(gca,'XLim',[1 ns])
            title('log-precision','FontSize',16)

        elseif np
            subplot(nl,2,2*nl)
            plot(1:ns,full(spm_cat(qP.p)))
            set(gca,'XLim',[1 ns])
            title('parameters (modes)','FontSize',16)

        end
        drawnow
    end

    % report (EM-Steps)
    %----------------------------------------------------------------------
    try
        dF = F(iN) - F(iN - 1);
    catch
        dF = 0;
    end
    str{1} = sprintf('LAP: %i (%i)', iN,iD);
    if iN == 1
        str{2} = sprintf('  F0:%.4e', full(F(iN)));
    else
        str{2} = sprintf('F-F0:%.4e', full(F(iN) - F(1)));
    end
    str{3} = sprintf('dF:%.2e',      full(dF));
    str{4} = sprintf('(%.2e sec)',   full(toc));
    fprintf('%-16s%-20s%-14s%-16s\n',str{:})

end
fprintf('%-19sF:%.4e\n', 'LAP: Converged', F(end));
if exist('Fdem','var')
    spm_figure('Focus', Fdem)
end

% Place Bayesian parameter averages in output arguments
%==========================================================================

% Conditional moments of time-averaged parameters
%--------------------------------------------------------------------------
Qp = 0;
Ep = 0;
for i = 1:ns

    % weight in proportion to precisions
    %----------------------------------------------------------------------
    P  = spm_inv(qP.c{i});
    Ep = Ep + P*qP.p{i};
    Qp = Qp + P;

end
Ep     = spm_inv(Qp)*Ep;
Cp     = spm_inv(Qp + (1 - ns)*Pp);
qP.P   = spm_unvec(Up*Ep + pp.p,{M.pE});
qP.C   = Up*Cp*Up';
qP.V   = spm_unvec(diag(qP.C),{M.pE});
qP.U   = Up;

% conditional moments of hyper-parameters
%--------------------------------------------------------------------------
Qh = 0;
Eh = 0;
for i = 1:ns

    % weight in proportion to precisions
    %----------------------------------------------------------------------
    P  = spm_inv(qH.c{i});
    Eh = Eh + P*qH.p{i};
    Qh = Qh + P;

end
Eh     = spm_inv(Qh)*Eh;
Ch     = spm_inv(Qh + (1 - ns)*Ph);
P      = spm_unvec(Eh,{qh.h,qh.g});
qH.h   = P{1};
qH.g   = P{2};
qH.C   = Ch;
P      = spm_unvec(diag(qH.C),P);
qH.V   = P{1};
qH.W   = P{2};

% assign output variables
%--------------------------------------------------------------------------
DEM.M  = M;                   % model
DEM.U  = U;                   % causes

DEM.qU = qU;                  % conditional moments of model-states
DEM.qP = qP;                  % conditional moments of model-parameters
DEM.qH = qH;                  % conditional moments of hyper-parameters

DEM.F  = F(1:iN);             % [-ve] Free-energy
DEM.S  = S(1:iN);             % [-ve] Free-action
DEM.FC = FC;                  % Free-energy components
DEM.CC = CC;                  % over iterations
DEM.FT = FT;                  % over time


text:
function [dx, dy] = GradientSobel(img)

        Sv = [-1 -2 -1;...
               0  0  0;...
               1  2  1];
        Sh = [-1  0  1;
              -2  0  2;
              -1  0  1];
        dx=conv2(img, Sh,'valid');
        dy=conv2(img, Sv,'valid');
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function v3 = r8vec_cross_product_3d ( v1, v2 )

%*****************************************************************************80
%
%% R8VEC_CROSS_PRODUCT_3D computes the cross product of two R8VEC's in 3D.
%
%  Discussion:
%
%    The cross product in 3D can be regarded as the determinant of the
%    symbolic matrix:
%
%          [  i  j  k ]
%      det [ x1 y1 z1 ]
%          [ x2 y2 z2 ]
%
%      = ( y1 * z2 - z1 * y2 ) * i
%      + ( z1 * x2 - x1 * z2 ) * j
%      + ( x1 * y2 - y1 * x2 ) * k
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real V1(3), V2(3), the two vectors.
%
%    Output, real V3(3), the cross product vector.
%
  v3 = zeros ( 3, 1 );

  v3(1) = v1(2) * v2(3) - v1(3) * v2(2);
  v3(2) = v1(3) * v2(1) - v1(1) * v2(3);
  v3(3) = v1(1) * v2(2) - v1(2) * v2(1);

  return
end


text:
function [c, ceq] = lotkaConstraintFCN_models(u,x,N,p,select_model)
%% Constraint function of nonlinear MPC for Lotka-Volterra system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon
%
% Output:
%   c:      inequality constraints applied across prediction horizon
%   ceq:    equality constraints (empty)
%

%% Nonlinear MPC design parameters
% Predator population size always positive: >0, min pop size
zMin = [10];

%% Integrate system
if strcmp(select_model,'DelayDMDc')
     [xk,~] = lsim(p.sys,[p.udelay(1:N);u'],[0:N-1].*p.dt,[p.xdelay(:,1); x]-[p.xmean; p.xmean]);
    xk = xk(:,3:4);
    xk = xk + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'DMDc')
    [xk,~] = lsim(p.sys,[u',0],[0:N].*p.dt,x-p.xmean);
    xk = xk(2:end,:) + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'SINDYc')
    Ns = size(x,1);
    xk = zeros(Ns,N+1); xk(:,1) = x;
    for ct=1:N
        % Obtain plant state at next prediction step.
        xk(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xk(:,ct),u(ct),p.dt,1,[],p);
    end
    xk = xk(:,2:N+1);
elseif strcmp(select_model,'NARX')
    Hu = [u',0];
    Hx = zeros(2,length(Hu)); Hx(:,1) = x;
    [Us,Ui,Si] = preparets(p.net,con2seq(Hu),{},con2seq(Hx));
    xk = p.net(Us,Ui,Si);
    xk = cell2mat(xk);
end

%% Inequality constraints calculation
c = zeros(N,1);
% Apply N population size constraints across prediction horizon, from time
% k+1 to k+N
uk = u(1);
for ct=1:N
    % -z + zMin < 0 % lower bound
    c(ct) = -xk(2,ct)+zMin; %c(2*ct-1)
    % z - zMax < 0 % upper bound
    %c(2*ct) = xk1(1)-zMax;
    % update  input for next step
    if ct<N
        uk = u(ct+1);
    end
end
%% No equality constraints
ceq = [];

author: eurika-kaiser

repo: SINDY-MPC


text:
function map = vega20(N)
% Colormap from MatPlotLib 2.0. For attractive plots using the ColorOrder.
%
%%% Syntax:
%  map = vega20
%  map = vega20(N)
%
% For MatPlotLib 2.0 improved colormaps were created for plot lines of
% categorical data. The new colormaps are introduced here:
% <http://matplotlib.org/2.0.0rc2/users/dflt_style_changes.html>
% <https://github.com/vega/vega/wiki/Scales#scale-range-literals>
% Note that VEGA10 is the new default Line Color Order for MatPlotLib 2.0.
%
% Information on the axes ColorOrder (note that this is NOT the axes COLORMAP):
% <https://www.mathworks.com/help/matlab/creating_plots/defining-the-color-of-lines-for-plotting.html>
% <https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html>
%
%% Examples %%
%
%%% PLOT using matrices:
% N = 20;
% axes('ColorOrder',vega20(N),'NextPlot','replacechildren')
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)n*sin(x+2*n*pi/N), X(:), 1:N);
% plot(X,Y, 'linewidth',4)
%
%%% PLOT in a loop:
% N = 20;
% set(0,'DefaultAxesColorOrder',vega20(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)n*sin(x+2*n*pi/N), X(:), 1:N);
% for n = 1:N
%     plot(X(:),Y(:,n), 'linewidth',4);
%     hold all
% end
%
%%% LINE using matrices:
% N = 20;
% set(0,'DefaultAxesColorOrder',vega20(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)n*cos(x+2*n*pi/N), X(:), 1:N);
% line(X(:),Y)
%
%% Input and Output Arguments %%
%
%%% Inputs (*=default):
% N = NumericScalar, N>=0, an integer to define the colormap length.
%   = *[], use the length of the current figure's colormap (see COLORMAP).
%
%%% Outputs:
% map = NumericMatrix, size Nx3, a colormap of RGB values between 0 and 1.
%
% map = vega20(N)

if nargin<1
	N = size(get(gcf,'colormap'),1);
else
	assert(isscalar(N)&&isreal(N),'First argument must be a real numeric scalar.')
	assert(fix(N)==N&&N>=0,'First argument must be a positive integer.')
end
%
raw = [...
	031,119,180;... #1f77b4
	174,199,232;... #aec7e8
	255,127,014;... #ff7f0e
	255,187,120;... #ffbb78
	044,160,044;... #2ca02c
	152,223,138;... #98df8a
	214,039,040;... #d62728
	255,152,150;... #ff9896
	148,103,189;... #9467bd
	197,176,213;... #c5b0d5
	140,086,075;... #8c564b
	196,156,148;... #c49c94
	227,119,194;... #e377c2
	247,182,210;... #f7b6d2
	127,127,127;... #7f7f7f
	199,199,199;... #c7c7c7
	188,189,034;... #bcbd22
	219,219,141;... #dbdb8d
	023,190,207;... #17becf
	158,218,229;... #9edae5
	];

raw = reshape(raw,2,10,3);
raw = reshape(permute(raw,[2 1 3]),20,3);

%
map = raw(1+mod(0:N-1,size(raw,1)),:) / 255;
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%vega20
% The file <vega_license.txt> contains the colormap license.
% This implmentation: Copyright (c) 2017 Stephen Cobeldick
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%license

author: mtex-toolbox

repo: mtex

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


text:
function [x,minCost]=convexQuadProgIntegerNoConst2D(Q,b,nonNeg)
%%CONVEXQUADPROGINTEGERNOCONST2D Perform quadratic programming on a 2D
%       convex problem with no constraints except that the solutions must
%       be integers that are either unconstrained or are non-negative.
%       This functions solves the optimization problem
%         minimize_x x'*Q*x+2*b'*x
%         subject to x is a 2X1 set of integers (optionally non-negative
%                    integers).
%
%INPUTS: Q An 2X2 real, positive definite symmetric matrix.
%        b An 2X1 real vector.
%   nonNeg This is true if all the x elements should be non-negative and
%          false otherwise. The default if omitted or an empty matrix is
%          passed is false.
%
%OUTPUTS: x The 2X1 optimal solution (or an empty matrix if the optimal
%           cost is infinite).
%   minCost The cost of the optimal constrained solution. If the cost isn't
%           finite, then this is an empty matrix.
%
%One could implement the branch-and bound algorithm for large scale
%problems as in [1]. However, in the 2D case, there are so few branches, it
%isn't worth it. Instead, we simply branch as in [1] (with a
%simple modification for non-negative constraints) and we traverse all 8
%hypotheses and then choose the lowest cost one. In [1], when branching,
%one creates smaller Q and b submatrices to optimize over rather than
%redoing the entire problem with an added eualaity constraint. In the 2D
%case, fixing an index means that the subproblem is just optimizing a scalar
%quadratic equation. Thus, we don't need to explicitely form submatrices;
%we just directly evaluate the scalar solution to the quadratic.
%
%EXAMPLE 1:
%Here is a simple example. We get the unconstrained solution. Then, we get
%the integer constrained solution and see that the cost is higher. Finally,
%we get the integer constrained solution with non-negativity constraints.
%That is NOT just clipping the integer constrained solution to 0. To show
%that, we compute the cost of the clipped solution and we see that it is
%much larger than the solution returned by the function.
% Q=[32,  20;
%    20, 22];
% b=[-178;-1200];
% [xOpt,minCost]=convexQuadProgEqConst(Q,b)
% nonNeg=false;
% [xOptInt,minCostInt]=convexQuadProgIntegerNoConst2D(Q,b,nonNeg)
% nonNeg=true;
% [xNonNegInt,minCostNonNegInt]=convexQuadProgIntegerNoConst2D(Q,b,nonNeg)
% %We show that the zero-clipped solution is much higher than the one with
% %non-negativity constraints.
% xTest=max(xOptInt,0);
% testCost=xTest'*Q*xTest+2*b'*xTest
%
%EXAMPLE 2:
%One might be tempted the think that one can always obtain the globally
%optimal solution by solving the unconstrained problem and then just
%rouning the result or by just trying the 4 points obtained by taking the
%floor and ceiling of each element of the unconstrained solution. However,
%that is not the case. In this example, we show the unconstrained solution
%and then the optimal integer constrained differ by more than 1 index.
% Q=[25, 0.3;
%   0.3, 0.01];
% b=[-156.9;
%    -3.31];
% xUnconst=convexQuadProgEqConst(Q,b)
% [x,minCost]=convexQuadProgIntegerNoConst2D(Q,b)
% %Now, we show that the cost of just rouning up or down the non-integer part
% %of xUnconst does not produce a globally optimal solution. Only the first
% %element of xUnConst is not an integer, so we only need to consider points
% %rounding that one up or down (as opposed to 4 points rounding each
% %element).
% xJustRoundUp=ceil(xUnconst);
% costJustRoundUp=xJustRoundUp'*Q*xJustRoundUp+2*b'*xJustRoundUp
% xJustRoundDown=floor(xUnconst);
% costJustRoundDown=xJustRoundDown'*Q*xJustRoundDown+2*b'*xJustRoundDown
%
%REFERENCES:
%[1] F. Koerner, "An efficient branch and bound algorithm to solve the
%    quadratic integer programming problem," Computing, vol. 30, pp. 253-
%    260, Sep. 1983.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(nonNeg))
    nonNeg=false;
end

xOpt=convexQuadProgEqConst(Q,b);

%Impose the non-negativity constraint. When branching on each value, if the
%non-negativity constraint was enforce,d both branches will be the same.
%However, we aren't going to bother checking to see whether a branch should
%be skipped, because there are so few branches.
if(nonNeg)
    xOpt=max(xOpt,0);
end

x=[];
minCost=Inf;

%Branch on the first index.
for x1=[floor(xOpt(1)),ceil(xOpt(1))]
    %The part of the cost function that depends only on x1 and not x2
    %(which is fixed).
    constCostTerm=2*b(1)*x1+Q(1,1)*x1^2;
    %The optimal solution for x2 given x1 and no constraints.
    x2Node=-(b(2)+Q(1,2)*x1)/(Q(2,2));
    if(nonNeg)
        %Apply a non-negativity constraint.
        x2Node=max(x2Node,0);
    end

    %Now, branch on the second index. With it fixed, we just evaluate the
    %cost (There is nothing else to optimize).
    for x2=[floor(x2Node),ceil(x2Node)]
        costCur=2*b(2)*x2+2*Q(1,2)*x1*x2+Q(2,2)*x2^2+constCostTerm;
        if(costCur<minCost)
            minCost=costCur;
            x=[x1;x2];
        end
    end
end

%Now branch on the second index.
for x2=[floor(xOpt(2)),ceil(xOpt(2))]
    %The part of the cost function that depends only on x2 and not x1
    %(which is fixed).
    constCostTerm=2*b(2)*x2+Q(2,2)*x2^2;
    %The optimal solution for x1 given x2 and no constraints.
    x1Node=-(b(1)+Q(1,2)*x2)/(Q(1,1));

    if(nonNeg)
        %Apply a non-negativity constraint.
        x1Node=max(x1Node,0);
    end

    %Now, branch on the first index.
    for x1=[floor(x1Node),ceil(x1Node)]
        costCur=2*b(1)*x1+Q(1,1)*x1^2+2*Q(1,2)*x1*x2+constCostTerm;
        if(costCur<minCost)
            minCost=costCur;
            x=[x1;x2];
        end
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [x, infos] = anls_nmf(V, rank, in_options)
% Alternative non-negative least squares (ANLS) for non-negative matrix factorization (NMF).
%
% The problem of interest is defined as
%
%       min || V - WH ||_F^2,
%       where
%       {V, W, H} > 0.
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       V           : (m x n) non-negative matrix to factorize
%       rank        : rank
%       in_options  : options
%
%
% References:
%       Jingu Kim, Yunlong He, and Haesun Park,
%       "Algorithms for Nonnegative Matrix and Tensor Factorizations: A Unified View
%       Based on Block Coordinate Descent Framework,"
%       Journal of Global Optimization, 58(2), pp. 285-319, 2014.
%
%       Jingu Kim and Haesun Park.
%       "Fast Nonnegative Matrix Factorization: An Active-set-like Method and Comparisons,"
%       SIAM Journal on Scientific Computing (SISC), 33(6), pp. 3261-3281, 2011.
%
%
% Output:
%       x           : non-negative matrix solution, i.e., x.W: (m x rank), x.H: (rank x n)
%       infos       : log information
%           epoch   : iteration nuber
%           cost    : objective function value
%           optgap  : optimality gap
%           time    : elapsed time
%           grad_calc_count : number of sampled data elements (gradient calculations)
%
%
% This file is part of NMFLibrary
%
% This code calls functions {nnls1_asgivens, nnlsm_activeset, nnlsm_blockpivot}
% written by Jingu Kim. See https://github.com/kimjingu/nonnegfac-matlab.
%
% Ported by H.Kasai on Apr. 04, 2017
%
% Change log:
%
%       Oct. 27, 2017 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jun. 24, 2022 (Hiroyuki Kasai): Added momentum acceleration mode and mofified.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.alg   = 'anls_asgroup';
    local_options.sub_mode = 'std';
    local_options.delta = 0.1;
    local_options.inner_max_epoch = 500;
    local_options.inner_max_epoch_parameter = 0.5;
    local_options.beta0 = 0.5;
    local_options.eta = 1.5;
    local_options.gammabeta = 1.01;
    local_options.gammabetabar = 1.005;
    local_options.momentum_h = 0;
    local_options.momentum_w = 0;
    local_options.scaling = true;
    local_options.warm_restart = false;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % set paramters
    if ~strcmp(options.alg, 'anls_asgroup') && ~strcmp(options.alg, 'anls_asgivens') ...
            && ~strcmp(options.alg, 'anls_bpp')
        fprintf('Invalid algorithm: %s. Therfore, we use anls_asgroup (i.e., ANLS with Active Set Method and Column Grouping).\n', options.alg);
        options.alg = 'anls_asgroup';
    end

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = sprintf('ANLS (%s)', options.alg);
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    if options.scaling
        [W, H] = normalize_WH(V, W, H, rank, 'type1');
    end

    % initialize for ANLS
    [options, beta, betamax] = check_momemtum_setting(options);

    if options.warm_restart
        nV = norm(V, 'fro');
        rel_error = zeros(1, options.max_epoch);
        rel_error(1) = sqrt(nV^2 - 2*sum(sum(V * H' .* W)) + sum(sum( H * H' .* (W'*W)))) / nV;
    end
    W_prev = W;
    H_prev = H;

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        %% update H
        if strcmp(options.alg, 'anls_asgroup')
            ow = 0;

            H = nnlsm_activeset(W'*W, W'*V, ow, 1, H);

        elseif strcmp(options.alg, 'anls_asgivens')
            ow = 0;

            WtV = W' * V;
            for i=1:size(H,2)
                H(:,i) = nnls1_asgivens(W'*W, WtV(:,i), ow, 1, H(:,i));
            end

        elseif strcmp(options.alg, 'anls_bpp')

            H = nnlsm_blockpivot(W'*W, W'*V, 1, H);

        end

        % perform momentum for H
        if strcmp(options.sub_mode, 'momentum')
            [H, H_tmp1, H_tmp2] = do_momentum_h(H, H_prev, beta, epoch, options);
        end

        %% update W
        if strcmp(options.alg, 'anls_asgroup')
            ow = 0;

            W = nnlsm_activeset(H*H', H*V', ow, 1, W');
            W = W';

        elseif strcmp(options.alg, 'anls_asgivens')
            ow = 0;

            HAt = H * V';
            Wt = W';
            for i=1:size(W,1)
                Wt(:,i) = nnls1_asgivens(H*H', HAt(:,i), ow, 1, Wt(:,i));
            end
            W = Wt';

        elseif strcmp(options.alg, 'anls_bpp')

            W = nnlsm_blockpivot(H*H', H*V', 1, W');
            W = W';

        end

        % perform momentum for W
        if strcmp(options.sub_mode, 'momentum')
            [W, H, W_tmp1] = do_momentum_w(W, W_prev, H, H_prev, H_tmp1, beta, epoch, options);
        end

        % perform warm_restart
        if options.warm_restart
            [W, H, W_prev, H_prev, rel_error, beta, betamax, options] = ...
                warm_restart(V, W, H, rank, W_prev, H_prev, W_tmp1, H_tmp1, H_tmp2, rel_error, beta, betamax, epoch, options);
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function expdirs = setUpDir(rootdir,varargin)

[ips,doforce,annotate,features,dosoftlink,frontside,rotatesideview,rotatefrontview] = ...
  myparse(varargin,...
  'ips',[],...
  'doforce',false,...
  'annotate',true,...
  'features',true,...
  'dosoftlink',false,...
  'frontside',false,...
  'rotatesideview',0,...
  'rotatefrontview',0);

if frontside
  createDirFcn = @makeMoviesFrontSide;
else
  createDirFcn = @makeMovies;
end

fprintf('Creating JAABA directories within root directory %s...\n',rootdir);
expdirs = createDirFcn(rootdir,'doforce',doforce,'dosoftlink',dosoftlink,...
  'rotatefrontview',rotatefrontview,'rotatesideview',rotatesideview);

if annotate
  fprintf('Annotate interest points...\n');
  annotateAllDirs(rootdir,'ips',ips,'doforce',doforce,'frontside',frontside);
end

if features
  fprintf('Generating spatiotemporal features...\n');
  genAllFeatures(rootdir,'doforce',doforce,'frontside',frontside);
end

author: kristinbranson

repo: JAABA


text:
function scrollplotdemo

% Created by Steven Lord, slord@mathworks.com
% Uploaded to MATLAB Central
% http://www.mathworks.com/matlabcentral
% 7 May 2002
%
% Permission is granted to adapt this code for your own use.
% However, if it is reposted this message must be intact.

% Generate and plot data
x=0:1e-2:2*pi;
y=sin(x);
dx=2;
%% dx is the width of the axis 'window'
a=gca;
p=plot(x,y);

% Set appropriate axis limits and settings
set(gcf,'doublebuffer','on');
%% This avoids flickering when updating the axis
set(a,'xlim',[0 dx]);
set(a,'ylim',[min(y) max(y)]);

% Generate constants for use in uicontrol initialization
pos=get(a,'position');
Newpos=[pos(1) pos(2)-0.1 pos(3) 0.05];
%% This will create a slider which is just underneath the axis
%% but still leaves room for the axis labels above the slider
xmax=max(x);
S=['set(gca,''xlim'',get(gcbo,''value'')+[0 ' num2str(dx) '])'];
%% Setting up callback string to modify XLim of axis (gca)
%% based on the position of the slider (gcbo)

% Creating Uicontrol
h=uicontrol('style','slider',...
    'units','normalized','position',Newpos,...
    'callback',S,'min',0,'max',xmax-dx);


text:
function mckinnon_test ( )

%*****************************************************************************80
%
%% MCKINNON_TEST works with the McKinnon function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
  global phi
  global tau
  global theta

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MCKINNON_TEST:\n' );
  fprintf ( 1, '  Test COMPASS_SEARCH with the McKinnon function.\n' );
  m = 2;
  delta_tol = 0.00001;
  delta = 0.3;
  k_max = 20000;
%
%  Test 1
%
  a = ( 1.0 + sqrt ( 33.0 ) ) / 8.0;
  b = ( 1.0 - sqrt ( 33.0 ) ) / 8.0;

  phi = 10.0;
  tau = 1.0;
  theta = 15.0;

  x = [ a, b ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '  PHI = %f, TAU = %f, THETA = %f\n', phi, tau, theta );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', mckinnon ( m, x ) );

  [ x, fx, k ] = compass_search ( @mckinnon, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );

  x = [ 0.0, -0.5 ];
  r8vec_print ( m, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', mckinnon ( m, x ) );
%
%  Test 2
%
  a = ( 1.0 + sqrt ( 33.0 ) ) / 8.0;
  b = ( 1.0 - sqrt ( 33.0 ) ) / 8.0;

  phi = 60.0;
  tau = 2.0;
  theta = 6.0;

  x = [ a, b ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '  PHI = %f, TAU = %f, THETA = %f\n', phi, tau, theta );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', mckinnon ( m, x ) );

  [ x, fx, k ] = compass_search ( @mckinnon, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );

  x = [ 0.0, -0.5 ];
  r8vec_print ( m, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', mckinnon ( m, x ) );
%
%  Test 3
%
  a = ( 1.0 + sqrt ( 33.0 ) ) / 8.0;
  b = ( 1.0 - sqrt ( 33.0 ) ) / 8.0;

  phi = 4000.0;
  tau = 3.0;
  theta = 6.0;

  x = [ a, b ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '  PHI = %f, TAU = %f, THETA = %f\n', phi, tau, theta );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', mckinnon ( m, x ) );

  [ x, fx, k ] = compass_search ( @mckinnon, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );

  x = [ 0.0, -0.5 ];
  r8vec_print ( m, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', mckinnon ( m, x ) );

  return
end


text:
% SLAMTB  An EKF-SLAM algorithm with simulator and graphics.
%
%   This script performs multi-robot, multi-sensor, multi-landmark 6DOF
%   EKF-SLAM with simulation and graphics capabilities.
%
%   Please read slamToolbox.pdf in the root directory thoroughly before
%   using this toolbox.
%
%   - Beginners should not modify this file, just edit USERDATA.M and enter
%   and/or modify the data you wish to simulate.
%
%   - More advanced users should be able to create new landmark models, new
%   initialization methods, and possibly extensions to multi-map SLAM. Good
%   luck!
%
%   - Expert users may want to add code for real-data experiments.
%
%   See also USERDATA, USERDATAPNT, USERDATALIN.
%
%   Also consult slamToolbox.pdf in the root directory.

%   Created and maintained by
%   Copyright 2008, 2009, 2010 Joan Sola @ LAAS-CNRS.
%   Copyright 2011, 2012, 2013 Joan Sola.
%   Programmers (for parts of the toolbox):
%   Copyright David Marquez and Jean-Marie Codol @ LAAS-CNRS
%   Copyright Teresa Vidal-Calleja @ ACFR.
%   See COPYING.TXT for full copyright license.

%% OK we start here

% clear workspace and declare globals
clear
global Map

%% I. Specify user-defined options - EDIT USER DATA FILE userData.m

userData;           % user-defined data. SCRIPT.
% userDataPnt;        % user-defined data for points. SCRIPT.
% userDataLin;        % user-defined data for lines. SCRIPT.

%% II. Initialize all data structures from user-defined data in userData.m
% SLAM data
[Rob,Sen,Raw,Lmk,Obs,Tim]     = createSlamStructures(...
    Robot,...
    Sensor,...      % all user data
    Time,...
    Opt);

% Simulation data
[SimRob,SimSen,SimLmk,SimOpt] = createSimStructures(...
    Robot,...
    Sensor,...      % all user data
    World,...
    SimOpt);

% Graphics handles
[MapFig,SenFig]               = createGraphicsStructures(...
    Rob, Sen, Lmk, Obs,...      % SLAM data
    SimRob, SimSen, SimLmk,...  % Simulator data
    FigOpt);                    % User-defined graphic options

%% III. Initialize data logging
% TODO: Create source and/or destination files and paths for data input and
% logs.
% TODO: do something here to collect data for post-processing or
% plotting. Think about collecting data in files using fopen, fwrite,
% etc., instead of creating large Matlab variables for data logging.

% Clear user data - not needed anymore
clear Robot Sensor World Time   % clear all user data

%% IV. Main loop
for currentFrame = Tim.firstFrame : Tim.lastFrame

    % 1. SIMULATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Simulate robots
    for rob = [SimRob.rob]

        % Robot motion
        SimRob(rob) = simMotion(SimRob(rob),Tim);

        % Simulate sensor observations
        for sen = SimRob(rob).sensors

            % Observe simulated landmarks
            Raw(sen) = simObservation(SimRob(rob), SimSen(sen), SimLmk, SimOpt) ;

        end % end process sensors

    end % end process robots

    % 2. ESTIMATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Process robots
    for rob = [Rob.rob]

        % Robot motion
        % NOTE: in a regular, non-simulated SLAM, this line is not here and
        % noise just comes from the real world. Here, the estimated robot
        % is noised so that the simulated trajectory can be made perfect
        % and act as a clear reference. The noise is additive to the
        % control input 'u'.
        Rob(rob).con.u = SimRob(rob).con.u + Rob(rob).con.uStd.*randn(size(Rob(rob).con.uStd));
        Rob(rob) = motion(Rob(rob),Tim);

        Map.t = Map.t + Tim.dt;

        % Process sensor observations
        for sen = Rob(rob).sensors

            % Observe knowm landmarks
            [Rob(rob),Sen(sen),Lmk,Obs(sen,:)] = correctKnownLmks( ...
                Rob(rob),   ...
                Sen(sen),   ...
                Raw(sen),   ...
                Lmk,        ...
                Obs(sen,:), ...
                Opt) ;

            % Initialize new landmarks
            ninits = Opt.init.nbrInits(1 + (currentFrame ~= Tim.firstFrame));
            for i = 1:ninits
                [Lmk,Obs(sen,:)] = initNewLmk(...
                    Rob(rob),   ...
                    Sen(sen),   ...
                    Raw(sen),   ...
                    Lmk,        ...
                    Obs(sen,:), ...
                    Opt) ;
            end

        end % end process sensors

    end % end process robots

    % 3. VISUALIZATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if currentFrame == Tim.firstFrame ...
            || currentFrame == Tim.lastFrame ...
            || mod(currentFrame,FigOpt.rendPeriod) == 0

        % Figure of the Map:
        MapFig = drawMapFig(MapFig,  ...
            Rob, Sen, Lmk,  ...
            SimRob, SimSen, ...
            FigOpt);

        if FigOpt.createVideo
            makeVideoFrame(MapFig, ...
                sprintf('map-%04d.png',currentFrame), ...
                FigOpt, ExpOpt);
        end

        % Figures for all sensors
        for sen = [Sen.sen]
            SenFig(sen) = drawSenFig(SenFig(sen), ...
                Sen(sen), Raw(sen), Obs(sen,:), ...
                FigOpt);

            if FigOpt.createVideo
                makeVideoFrame(SenFig(sen), ...
                    sprintf('sen%02d-%04d.png', sen, currentFrame),...
                    FigOpt, ExpOpt);
            end

        end

        % Do draw all objects
        drawnow;
    end

    % 4. DATA LOGGING
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TODO: do something here to collect data for post-processing or
    % plotting. Think about collecting data in files using fopen, fwrite,
    % etc., instead of creating large Matlab variables for data logging.

end

%% V. Post-processing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter post-processing code here

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function calpak_test48 ( )

%*****************************************************************************80
%
%% CALPAK_TEST48 tests S_TO_YMD_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2013
%
%  Author:
%
%    John Burkardt
%
  n_test = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST48\n' );
  fprintf ( 1, '  S_TO_YMD_COMMON converts a string to a YMD date.\n' );
  fprintf ( 1, '\n' );

  for i_test = 1 : n_test

    if ( i_test == 1 )
      s = '1999-10-31';
      p = 'yyyy-mm-dd';
    elseif ( i_test == 2 )
      s = '01/04/2004';
      p = 'dd/mm/yyyy';
    elseif ( i_test == 3 )
      s = '8/8/88';
      p = 'd/m/yy';
    elseif ( i_test == 4 )
      s = '4 7';
      p = 'd m';
    elseif ( i_test == 5 )
      s = 'On day 1 of month 3 of year 1945';
      p = 'xx xxx d xx xxxxx m xx xxxx yyyy';
    end

    [ y, m, d ] = s_to_ymd_common ( s, p );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %s\n', s );
    fprintf ( 1, '  %s\n', p );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  %4d  %2d  %2d\n', y, m, d );

  end

  return
end


text:
function [supported] = tracker_test(tracker)
% trax_test Test tracker's support for TraX protocol
%
% This function tries to run the tracker and communicate with it using the
% TraX protocol to determine if it supports it. The results are chached
% in the workspace cache directory.
%
% Input:
% - tracker: Tracker structure.
%
% Output:
% - supported (boolean): Is the protocol supported.

% Check if the result of the test is already cached
cache = fullfile(get_global_variable('directory'), 'cache', 'trax');

tracker_hash = md5hash(sprintf('%s-%s-%s', tracker.command, tracker.interpreter, strjoin(tracker.linkpath, '-')));
mkpath(cache);
cache_file = fullfile(cache, sprintf('trax_%s_%s.mat', tracker.identifier, tracker_hash));

supported = [];
if exist(cache_file, 'file')
    load(cache_file);
    if ~isempty(supported)
        return;
    end;
end;

print_text('Testing TraX protocol support for tracker %s.', tracker.identifier);

try

    data.header = false;

    data = tracker_run(tracker, @callback, data);

    supported = data.header;

catch

    supported = false;

end

if supported
    % Only cache if support is detected
    save(cache_file, 'supported');
    print_text('TraX support detected.');
else
	print_text('TraX support not detected.');
end

end

function [image, region, properties, data] = callback(state, data)

    image = [];
    region = [];
    properties = struct();

    data.header = true;
    data.state = state;

end

author: votchallenge

repo: toolkit-legacy


text:
classdef PostProcessor

  methods (Static)

    function posts = PostProcess(curs,postprocessparams,scoreNorm)
      if isempty(postprocessparams);
        posts = curs;
        return;
      end;

      if strcmpi(postprocessparams.method,'hysteresis')
        posts = PostProcessor.ApplyHysteresis(curs,postprocessparams,scoreNorm);
      else
        posts = PostProcessor.ApplyFiltering(curs,postprocessparams);
      end

      posts = PostProcessor.RemoveSmallBouts(posts,postprocessparams);
    end

    % ---------------------------------------------------------------------
    function posts = RemoveSmallBouts(posts,postprocessparams)
      if postprocessparams.blen > 1 && numel(posts)>0,
        if numel(posts)<= postprocessparams.blen
          if nnz(posts>0) > numel(posts)/2,
            posts(:) = 1;
          else
            posts(:) = -1;
          end
          return;
        end

        while true,
          tposts = [posts(:)' 1-posts(end)];
          ends = find(tposts(1:end-1)~=tposts(2:end));
          ends = [1 ends+1];
          blens = ends(2:end)-ends(1:end-1);
          [minblen,smallbout] = min(blens);
          if minblen>=postprocessparams.blen, break; end
          posts(ends(smallbout):ends(smallbout+1)-1) = ...
            1 - posts(ends(smallbout):ends(smallbout+1)-1);
        end
      end
    end

    % ---------------------------------------------------------------------
    function posts = ApplyHysteresis(curs,params,scoreNorm)
      % Use imfill to find the regions.
      if isempty(curs), posts = curs; return; end

      % Select pos bouts that have at least one frame about the high
      % threshold.
      hthresh = curs > params.hystopts(1).value*scoreNorm;
      lthresh = curs > 0;
      if( nnz(hthresh)>0)
        pos = imfill(~lthresh,find(hthresh(:))) & lthresh;
        computeNeg = true;
      else
        pos = false(size(curs));
        computeNeg = false;
      end
      % Select neg bouts that have at least one frame below the low
      % threshold.
      hthresh = curs < params.hystopts(2).value*scoreNorm;
      lthresh = curs < params.hystopts(1).value*scoreNorm;
      if nnz(hthresh)>0 && computeNeg,
        neg = imfill(~lthresh,find(hthresh(:))) & lthresh;
      else
        neg = true(size(curs));
      end

      posts = pos | ~neg;

    end

    % ---------------------------------------------------------------------
    function posts = ApplyFiltering(curs,params)
      % Use filt to find the regions.
      if isempty(curs), posts = curs; return; end
      filts = conv(curs,ones(1,params.filtopts(1).value),'same');
      posts = filts>0;
    end

  end

end


text:
function [fx,dfdx,dfdP] = g_Id_phi(x,P,u,in)
% this function applies the identity mapping
fx = P;
dfdx = [];
dfdP = eye(size(P,1));

author: MBB-team

repo: VBA-toolbox


text:
function x = p18_root ( i )

%*****************************************************************************80
%
%% P18_ROOT returns a root for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the requested root.
%
%    Output, real X, the value of the root.
%
  if ( i == 1 )
    x = 1.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P18_ROOT - Fatal error!\n' );
    fprintf ( 1, '  Illegal root index = %d\n', i );
    error ( 'P18_ROOT - Fatal error!' );
  end

  return
end


text:
% Extract WINDOWSIZE consecutive frames at position SLIDINGWINDOWIDX from a
% given image sequence.
function slidingWindow = assembleSlidingWindows(I, startFrame, windowSize, slidingWindowIdx)

    % Get the index of the reference frame for the current sliding window.
    refFrameIdx = startFrame + slidingWindowIdx - 1;

    % Extract sliding window from the given input sequence.
    slidingWindow.frames = I(:,:,(refFrameIdx - floor(windowSize / 2)):(refFrameIdx + floor(windowSize / 2)));
    % Capture reference frame for the current window.
    slidingWindow.referenceFrame = I(:,:,refFrameIdx);

author: thomas-koehler

repo: SupER


text:
function [img,nr,nc] = getframe_invisible_nocheck(gfdata,sz,removeborder,doresize)

if nargin < 3,
  removeborder = true;
end
if nargin < 4,
  doresize = true;
end

if isfield(gfdata,'grayborder'),
  GRAYBORDER = gfdata.grayborder;
else
  GRAYBORDER = 204;
end

hfig = gfdata.hfig;

% img = hardcopy(gfdata.hardcopy_args{:});
img = print(hfig,'-RGBImage');
set(hfig,'Units',gfdata.units,'Position',gfdata.pos);
if numel(img) == 1,
  fprintf('Could not grab invisible figure. Making visible temporarily.\n');
  set(hfig,'visible','on');
  pause(.1);
  img = hardcopy(gfdata.hardcopy_args{:});
  pause(.1);
  set(hfig,'visible','off');
end

% crop off the extra border
if removeborder,
  isallgray = all(img==GRAYBORDER,3);
end

[nr,nc,tmp] = size(img);
if ~doresize,
  sz(1) = nr; sz(2) = nc;
end

maxfirstcol = nc - sz(2) + 1;
maxfirstrow = nr - sz(1) + 1;
if maxfirstcol <= 0 || maxfirstrow <= 0,
  error('image has size [%d,%d] smaller than input sz = [%d,%d]',nr,nc,sz(1),sz(2));
end
if removeborder,
  firstcol = min(find(~all(isallgray,1),1),maxfirstcol);
  firstrow = min(find(~all(isallgray,2),1),maxfirstrow);
else
  firstcol = 1;
  firstrow = 1;
end
lastcol = firstcol + sz(2) - 1;
lastrow = firstrow + sz(1) - 1;
if removeborder,
  lastcol0 = find(~all(isallgray,1),1,'last');
  lastrow0 = find(~all(isallgray,2),1,'last');
else
  lastcol0 = nc;
  lastrow0 = nr;
end
if lastrow ~= lastrow0 || lastcol ~= lastcol0,
  fprintf('input width = %d, actual width = %d, input height = %d, actual height = %d\n',...
    sz(2),lastcol0-firstcol+1,sz(1),lastrow0-firstrow+1);
end

img = img(firstrow:lastrow,firstcol:lastcol,:);