text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function dwindowt
%DWINDOWT Unit test for the function DWINDOW.

%	O. Lemoine - March 1996.

N=5000;

% Rectangular window
h=tftb_window(N,'rect');
Dh1=[1;zeros(N-2,1);-1];
Dh2=dwindow(h);
if any(abs(Dh1-Dh2)>sqrt(eps)),
 error('dwindow test 1 failed');
end;

% Hanning window
h=tftb_window(N,'hanning');
Dh1=pi*sin(2*pi*(1:N)'/(N+1))/(N+1);
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)),
 error('dwindow test 2 failed');
end;

% Bartlett window
h=tftb_window(N,'bartlett');
Dh1=[2*ones(N/2-1,1);1;-1;-2*ones(N/2-1,1)]/N;
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)*10),
 error('dwindow test 3 failed');
end;

% Papoulis window
h=tftb_window(N,'papoulis');
Dh1=pi*cos(pi*(1:N)'/(N+1))/(N+1);
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)*10),
 error('dwindow test 4 failed');
end;

% Gaussian window
K=0.02; h=tftb_window(N,'gauss',K);
t=linspace(-1,1,N)';
Dh1=4*log(K)*t.*exp(log(K)*t.^2)/N;
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)*10),
 error('dwindow test 5 failed');
end;

N=4977;

% Rectangular window
h=tftb_window(N,'rect');
Dh1=[1;zeros(N-2,1);-1];
Dh2=dwindow(h);
if any(abs(Dh1-Dh2)>sqrt(eps)),
 error('dwindow test 6 failed');
end;

% Hanning window
h=tftb_window(N,'hanning');
Dh1=pi*sin(2*pi*(1:N)'/(N+1))/(N+1);
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)),
 error('dwindow test 7 failed');
end;

% Papoulis window
h=tftb_window(N,'papoulis');
Dh1=pi*cos(pi*(1:N)'/(N+1))/(N+1);
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)*10),
 error('dwindow test 8 failed');
end;

% Gaussian window
K=0.02; h=tftb_window(N,'gauss',K);
t=linspace(-1,1,N)';
Dh1=4*log(K)*t.*exp(log(K)*t.^2)/N;
Dh2=dwindow(h);
if any(abs(Dh1(2:N-1)-Dh2(2:N-1))>sqrt(eps)*10),
 error('dwindow test 9 failed');
end;

author: HeLiangHIT


text:
setenv('MKL_NUM_THREADS','1')
setenv('MKL_SERIAL','YES')
setenv('MKL_DYNAMIC','NO')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function vectarrow(p0,p1)
%Arrowline 3-D vector plot.
%   vectarrow(p0,p1) plots a line vector with arrow pointing from point p0
%   to point p1. The function can plot both 2D and 3D vector with arrow
%   depending on the dimension of the input
%
%   Example:
%       3D vector
%       p0 = [1 2 3];   % Coordinate of the first point p0
%       p1 = [4 5 6];   % Coordinate of the second point p1
%       vectarrow(p0,p1)
%
%       2D vector
%       p0 = [1 2];     % Coordinate of the first point p0
%       p1 = [4 5];     % Coordinate of the second point p1
%       vectarrow(p0,p1)
%
%   See also Vectline

%   Rentian Xiong 4-18-05
%   $Revision: 1.0

  if max(size(p0))==3
      if max(size(p1))==3
          x0 = p0(1);
          y0 = p0(2);
          z0 = p0(3);
          x1 = p1(1);
          y1 = p1(2);
          z1 = p1(3);
          plot3([x0;x1],[y0;y1],[z0;z1]);   % Draw a line between p0 and p1

          p = p1-p0;
          alpha = 0.1;  % Size of arrow head relative to the length of the vector
          beta = 0.1;  % Width of the base of the arrow head relative to the length

          hu = [x1-alpha*(p(1)+beta*(p(2)+eps)); x1; x1-alpha*(p(1)-beta*(p(2)+eps))];
          hv = [y1-alpha*(p(2)-beta*(p(1)+eps)); y1; y1-alpha*(p(2)+beta*(p(1)+eps))];
          hw = [z1-alpha*p(3);z1;z1-alpha*p(3)];

          hold on
          plot3(hu(:),hv(:),hw(:))  % Plot arrow head
          grid on
          xlabel('x')
          ylabel('y')
          zlabel('z')
          hold off
      else
          error('p0 and p1 must have the same dimension')
      end
  elseif max(size(p0))==2
      if max(size(p1))==2
          x0 = p0(1);
          y0 = p0(2);
          x1 = p1(1);
          y1 = p1(2);
          plot([x0;x1],[y0;y1]);   % Draw a line between p0 and p1

          p = p1-p0;
          alpha = 0.1;  % Size of arrow head relative to the length of the vector
          beta = 0.1;  % Width of the base of the arrow head relative to the length

          hu = [x1-alpha*(p(1)+beta*(p(2)+eps)); x1; x1-alpha*(p(1)-beta*(p(2)+eps))];
          hv = [y1-alpha*(p(2)-beta*(p(1)+eps)); y1; y1-alpha*(p(2)+beta*(p(1)+eps))];

          hold on
          plot(hu(:),hv(:))  % Plot arrow head
          grid on
          xlabel('x')
          ylabel('y')
          hold off
      else
          error('p0 and p1 must have the same dimension')
      end
  else
      error('this function only accepts 2D or 3D vector')
  end

author: Sable

repo: mcbench-benchmarks


text:
function [ret]=realTime_startVelControlJoints( t_Kuka )
%% Applicable to KST 1.6

%% Syntax:
% realTime_startVelControlJoints( t_Kuka )

%% About:
% This function is used to turn on the soft real-time control at joints
% velocities level

%% Check also
% ((sendJointsVel)), refer also to the function
% ((realTime_stopVelControlJoints)).

%% Arreguments
% t_Kuka: is the TCP/IP connection object

%% Return value:
% ret: 'true' if the joints positions message has been received and processed
% successfully by the server, otherwise 'false' is returned.

% Copy right, Mohammad SAFEEA, 1st of April 2018

theCommand='stVelDcJoint_';
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);

[ret]=checkAcknowledgment(message);
if ret==true
	disp('Acknowledge, realtime control turned on');
end
delay(0.5); % introduce some time delay, so the robot turns on the direct servo before starting the motion
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function result = is_orthonormal_set(A)
%
% Orthonormal Sets
%
% is_orthonormal_set(A) determines if a set of vectors in Euclidean n-space
% is orthonormal. The matrix A is formed from these vectors as its columns.
% That is, the subspace spanned by the set of vectors is the column space
% of A. The value 1 is returned if the set is orthonormal. The value 0 is
% returned if the set is not orthonormal. An error is returned if a set
% containing only zero vectors is attempted to be determined for
% orthonormality.
%
% For example, if the set of row vectors (u1,u2,...} is to be determined
% for orthonormality, set A to be equal to [u1' u2' ...].
%
%  Function written by Anthony Russo, downloaded from MatlabCentral.
%

matrix_size = size(A);

m = matrix_size(1,1);
n = matrix_size(1,2);

tolerance = 10^-10;

if A == zeros(m,n)
    error('The set that contains just zero vectors cannot be orthonormal.');
elseif n == 1
    if norm(A(1:m,1)) == 1
        result = 1;
    else
        result = 0;
    end
else
    if is_orthogonal_set(A) == 1
        length_counter = 0;

        for i = 1:n
            if abs(norm(A(1:m,i)) - 1) <= tolerance
                length_counter = length_counter + 1;
            end
        end

        if length_counter == n
            result = 1;
        else
            result = 0;
        end
    else
        result = 0;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% EXAMPLE 8: Band-pass filtering.
%
%   Reference:
%     [1] A. Nekkanti, O. T. Schmidt, Frequency–time analysis, low-rank reconstruction and denoising of turbulent flows using SPOD,
%         Journal of Fluid Mechanics 926, A26, 2021
%
% O. T. Schmidt (oschmidt@ucsd.edu)
% Last revision: 5-Sep-2022

clc, clear variables
addpath('utils')
disp('Loading the entire test database might take a second...')
load(fullfile('jet_data','jetLES.mat'),'p','p_mean','x','r','dt');

%   trapezoidal quadrature weights for cylindrical coordinates
intWeights      = trapzWeightsPolar(r(:,1),x(1,:));

%% SPOD
%   We will use standard parameters: a Hann window of length 256 and 50%
%   overlap.
nDFT            = 256;
nOvlp           = nDFT/2;
[L,P,f,~,A]     = spod(p,nDFT,intWeights,nOvlp,dt);

figure
loglog(f,L)
title('SPOD of full data')
xlabel('frequency'), ylabel('SPOD mode energy')
ylims           = ylim;

%% Filtering
%   An band-pass filter that focusses on the 'low-rank' portion of the
%   spectrum is implemented by setting to zero the SPOD expansion
%   coefficients in the range f<=0.08 and f>=1.0.
f_lowpass       = 1;
f_highpass      = 0.08;
A(f>=f_lowpass|f<=f_highpass,:,:) ...
                = 0;

%   The inverse SPOD using the modified SPOD expansion coefficients yields
%   the band-pass filtered data.
nt              = size(p,1);
p_rec           = invspod(P,A,nDFT,nOvlp);

%% Animate
%   Animate the original, filtered, and removed data.
figure
for t_i=1:1:30
    subplot(3,1,1)
    pcolor(x,r,squeeze(p(t_i,:,:))-p_mean); shading interp, axis equal tight
    title('Original data')
    if t_i==1; pmax = max(abs(caxis)); end, caxis(0.5*pmax*[-1 1]), colorbar
    subplot(3,1,2)
    pcolor(x,r,squeeze(p_rec(t_i,:,:))); shading interp, axis equal tight
    title('Reconstructed data')
    caxis(0.5*pmax*[-1 1]), colorbar
    subplot(3,1,3)
    pcolor(x,r,(squeeze(p(t_i,:,:))-p_mean)-squeeze(p_rec(t_i,:,:))); shading interp, axis equal tight
    title('Filtered/removed component')
    caxis(0.5*pmax*[-1 1]), colorbar
    drawnow
end

%% SPOD of filtered data
%   The effect of the band-pass filter should to a large degree remove
%   high-frequency components above the filter cut-on frequency.
[L,P,f,~,A]     = spod(p_rec,nDFT,intWeights,nOvlp,dt);

figure
loglog([f_lowpass f_lowpass],ylims,'k--'); hold on
loglog([f_highpass f_highpass],ylims,'k:');
loglog(f,L)
title('SPOD of filtered data')
xlabel('frequency'), ylabel('SPOD mode energy')
ylim(ylims);
legend('f_{low-pass}','f_{high-pass}')

author: SpectralPOD


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

p = sigmoid(X * theta);

idx_1 = find(p >= 0.5);
idx_0 = find(p < 0.5);

p(idx_1) = ones(size(idx_1));
p(idx_0) = zeros(size(idx_0));

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

end

repo: machine-learning-coursera


text:
function l = l1dn_inverse ( n, h )

%*****************************************************************************80
%
%% L1DN_INVERSE stores the inverse of the 1D DN Laplacian.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real L(N,N), the inverse of the Laplacian matrix.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1DN_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1DN_INVERSE - Fatal error!' );
  end

  l = zeros ( n, n );

  for j = 1 : n
    for i = 1 : n
      l(i,j) = min ( i, j ) * h * h;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ gmm_c ] = compress_gmm(gmm,Xtrain,threashold,options)
%COMPRESS_GMM Compresses a GMM in which there are possibly duplicate
%parameters.
%
%   input ------------------------------------------------------
%
%       o gmm: strcut, gmm.Priors
%                      gmm.Mu
%                      gmm.Sigma
%
%       o Xtrain: (N x D), training data used to retrain the parameters of
%       the GMM after a few have been removed
%
%   output -----------------------------------------------------
%
%       o gmm_c: struct, GMM with less than or equal number of parameters
%       as gmm.
%

GMModel.ComponentProportion = gmm.Priors;
GMModel.mu                  = gmm.Mu;
GMModel.Sigma               = gmm.Sigma;

gmm_c = merge_gmm_components2(GMModel,Xtrain,threashold,options);

end

author: epfl-lasa


text:
function pz=lpcpp2pz(pp)
%LPCPP2PZ LPC: Convert power spectrum polynomial in cos(w) to power spectrum zeros PZ=(RP)
% pp is a polynomial such that |polyval(ra,e^jw)| = polyval(pp,cos(w))

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

pz=roots(pp);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function varargout = dirwalk(topPath, visitor, recursionLimit, varargin)
%DIRWALK Generate the file names in a directory tree by walking the tree
%
% Description:
%   Function DIRWALK generates the file names in a directory tree by walking the tree
%   top-down. For each directory in the tree rooted at directory topPath.
%   For each directory of tree you can call "Visitor Function" for files processing.
%
%
% Using:
%   [pathNames, dirNames, fileNames] = dirwalk(topPath)
%   dirwalk(topPath, visitor)
%   dirwalk(topPath, visitor, varargin)
%   varargout = dirwalk(topPath, visitor, varargin)
%   [visitorOutput1, visitorOutput2, ..., visitorOutputN] = dirwalk(topPath, visitor)
%   [...] = dirwalk(topPath, visitor, visitorInput1, visitorInput2, ..., visitorInputN)
%   [...] = dirwalk(topPath, visitor, varargin)
%c
%
% Input:
%   topPath -- (Required) Top path name (Root path name)
%
%   visitor -- (Optional) Function handle. The function will be called
%                         when visiting each directory of tree.
%
%   Signatures of visitor function:
%   	visitor(pathName, DirListing)
%       visitor(pathName, DirListing, in1, in2, ..., inN)
%       visitor(pathName, DirListing, varargin)
%       varargout = visitor(...)
%      	[out1, out2, ..., outN] = visitor(...)
%
%   Input arguments:
%       pathName   -- (Required) Path to visited directory. String.
%                                (passed within DIRWALK)
%
%       DirListing -- (Required) Visited directory listing.
%                                Array of structs output of function DIR.
%                                (passed within DIRWALK)
%
%       varargin   -- (Optional) Other input arguments. (Passed outside)
%
%
% Output:
%   dirPaths  -- (Default visitor output) Visited path names. Cell array of strings.
%   dirNames  -- (Default visitor output) Directory names in visited paths (without '.' and '..'). Cell array of cell arrays.
%   fileNames -- (Default visitor output) File names in visited paths. Cell array of cell arrays.
%   visitorOutputs -- Visitor function outputs
%
%
% Examples:
%   topPath = fullfile(matlabroot, 'toolbox', 'matlab', 'demos');
%
%   [pathNames, dirNames, fileNames] = dirwalk(topPath);
%
%   dirwalk(topPath, @(y,x) disp(strcat(y, filesep, {x.name}')))
%
%
% See Also DIR, LS
%

% -------------------------------------------------------------------------
%   Version   : 1.1
%   Author    : Evgeny Pr aka iroln <esp.home@gmail.com>
%   Created   : 10.10.10
%   Updated   : 03.07.11
%
%   Copyright : Evgeny Prilepin (c) 2010-2011
% -------------------------------------------------------------------------

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

if (nargin < 2)
    error(nargoutchk(0, 3, nargout));
    recursionLimit = Inf;
    visitor = @default_visitor;
    visitorNumOutputs = 3;
    isUseDefaultVisitor = true;
elseif nargin < 3
    recursionLimit = Inf;
else
    visitorNumOutputs = nargout;
    isUseDefaultVisitor = false;
end

validateattributes(topPath, {'char'}, {'row'}, ...
    mfilename('fullpath'), '"Top Path Name"', 1)

if ~isempty(visitor)
    validateattributes(visitor, {'function_handle'}, {'scalar'}, ...
        mfilename('fullpath'), '"Visitor Function Handle"', 2)
end

varargout = dir_tree_helper(topPath, ...
    visitor, visitorNumOutputs, isUseDefaultVisitor, recursionLimit, varargin{:});
%--------------------------------------------------------------------------

%==========================================================================
function visitorOutputs = dir_tree_helper(topPath, visitorHandle, ...
    visitorNumOutputs, isUseDefaultVisitor, recursionLimit, varargin)
%DIR_TREE_HELPER Helper function for get directory tree

% Allocate memory for tree listing
preallocDirItems = 100000;
outputs = cell(preallocDirItems, visitorNumOutputs);
counter = 1;
depth = 0;

% Get tree listing
[counter, outputs] = dir_tree_listing(topPath, counter, depth, outputs, ...
    visitorHandle, visitorNumOutputs, isUseDefaultVisitor, recursionLimit, varargin{:});

% Remove extra
if (counter < preallocDirItems)
    outputs(counter+1:end, :) = [];
end

% Construct visitor function outputs
visitorOutputs = cell(1, visitorNumOutputs);
for i = 1:visitorNumOutputs
    visitorOutputs(i) = {outputs(:,i)};
end
%--------------------------------------------------------------------------

%==========================================================================
function [counter, outputs] = dir_tree_listing(topPath, counter, depth, outputs, ...
    visitorHandle, visitorNumOutputs, isUseDefaultVisitor, recursionLimit, varargin)
%DIR_TREE_LISTING Generate dir tree listing

% Get listing of current directory
Listing = dir(topPath);

% remove .* hidden directories here
for i =1:length(Listing)
    if strcmp(Listing(i).name(1),'.') & length(Listing(i).name)>2
       bad(i) = 1;
    else
        bad(i) = 0;
    end
end
if ~isempty(Listing)
Listing(find(bad)) = [];
end
% Call Visitor function
cwd = pwd;
cd(topPath)
disp(varargin)
visitorOutputs = visitor_call_helper(topPath, Listing, ...
    visitorHandle, visitorNumOutputs, varargin{:});
cd(cwd);

if isUseDefaultVisitor
    dirNames = visitorOutputs{2};
else
    dirNames = get_dir_file_names(Listing);
end

if (visitorNumOutputs > 0)
    outputs(counter,:) = visitorOutputs;
end

% Recursive walking directories in current directory
for i = 1:length(dirNames)
    nextRootPath = fullfile(topPath, dirNames{i});
    if depth + 1 <= recursionLimit
        [counter, outputs] = dir_tree_listing(nextRootPath, counter+1, depth+1, outputs, ...
            visitorHandle, visitorNumOutputs, isUseDefaultVisitor, recursionLimit, varargin{:});
    end
end

%--------------------------------------------------------------------------

%==========================================================================
function visitorOutputs = visitor_call_helper(topPath, Listing, ...
    visitor, numOutputs, varargin)
%VISITOR_CALL_HELPER Helper function for call visitor function

if (numOutputs == 0)
    % Visitor function without output arguments
    visitorOutputs = {};
%     visitor(topPath, Listing, varargin{:});

    visitor();
%     visitorCalling = sprintf('[] = visitor(topPath, Listing, varargin{:});');
%     eval(visitorCalling);
else
    % Construct eval command for call visitor function with any number output arguments
    outputNumbers = num2cell(1:numOutputs);

    %FIXME: It is not recommended, but otherwise doesn't work
    outputArgs = deblank(sprintf('outputs{%d} ', outputNumbers{:}));
    visitorCalling = sprintf('[%s] = visitor(topPath, Listing, varargin{:});', outputArgs);
    eval(visitorCalling);

    visitorOutputs = outputs;
end
%--------------------------------------------------------------------------

%==========================================================================
function [pathName, dirNames, fileNames] = default_visitor(rootPath, Listing)
%DEFAULT_VISITOR Default Visitor function
%
% Default visitor function return 3 output arguments:
%   pathNames -- visited paths names
%   dirNames  -- Directories names in visited directories
%   fileNames -- Files names in visited directories
%

pathName = rootPath;
[dirNames, fileNames] = get_dir_file_names(Listing);
%--------------------------------------------------------------------------

%==========================================================================
function [dirNames, fileNames] = get_dir_file_names(Listing)
%GET_DIR_FILE_NAMES

names = {Listing.name}';
isDirs = [Listing.isdir];
dirNames = names(isDirs);
fileNames = names(~isDirs);

% Exclude special directories '.' and '..'
inds = ~strcmp(dirNames, '.') & ~strcmp(dirNames, '..');
dirNames = dirNames(inds);
%--------------------------------------------------------------------------

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function linpack_d_test20 ( )

%*****************************************************************************80
%
%% TEST20 tests DPPFA and DPPDI.
%
%  Discussion:
%
%    DPPFA factors a packed positive definite symmetric matrix,
%    and DPPDI can compute the determinant or the inverse.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST20\n' );
  fprintf ( 1, '  For a positive definite symmetric packed matrix,\n' );
  fprintf ( 1, '  DPPFA factors the matrix.\n' );
  fprintf ( 1, '  DPPDI computes the inverse or determinant.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the matrix A.
%
  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      if ( i == j - 1 )
        a(k) = -1.0;
      elseif ( i == j )
        a(k) = 2.0;
      else
        a(k) = 0.0;
      end
    end
  end
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, info ] = dppfa ( a, n );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Error, DPPFA returns INFO = %d\n', info );
    return
  end
%
%  Invert the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get the determinant and inverse.\n' );

  job = 11;
  [ a, det ] = dppdi ( a, n, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1,'  Determinant  = %f * 10 ^ %f\n', det(1), det(2) );
%
%  DPPDI produces only the 'upper half triangle' of the inverse,
%  which is actually symmetric.  Thus, the lower half could be
%  produced by copying from the upper half.  However, the first row
%  of A, as returned, is exactly the first row of the inverse.
%
  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      b(i,j) = a(k);
      b(j,i) = a(k);
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Inverse:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', b(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function gi = slgraphinfo(G, conds)
%SLGRAPHINFO Extracts basic information of a given graph representation
%
% $ Syntax $
%   - gi = slgraphinfo(G)
%   - gi = slgraphinfo(G, conds)
%
% $ Arguments $
%   - G:        The input graph
%   - conds:    The cell array of conditions to be checked
%   - gi:       the information of the graph
%               a struct with the following fields:
%               - type: the type of the graph:
%                       - 'ge': a general graph
%                       - 'bi': a bigraph
%               - form: the form of the graph representation:
%                       - 'edgeset': edget set representation
%                       - 'adjlist': adjacency list representation
%                       - 'adjmat':  adjacency matrix representation
%               - n:    the number of (source) nodes
%               - nt:   the number of (target) nodes
%               - valued: whether the graph has values on edges
%
% $ Description $
%   - gi = slgraphinfo(G) extracts the basic information of the graph
%     representation and examines the integrity of it.
%
%   - gi = slgraphinfo(G, conds) verifies whether the graph conforms
%     to special conditions. The acceptable conditions:
%       - 'square': a square graph, it can be a general graph or a bigraph
%                   with n == nt. All graph that is square can be
%                   considered as a general graph
%       - 'edgeset': the representation is an edge set
%       - 'adjlist': the representation is an adjacency list
%       - 'adjmat':  the representation is an adjacency matrix
%       - [n]:       the n is equal to specified value
%       - [n, nt]:   the n and nt are equal to specified value
%       - 'numeric': the value type is numeric (only take effect for adjmat)
%       - 'logical': the value type is logical (only take effect for adjmat)
%
% $ Remarks $
%   - For a general graph, n and nt are equal, being the number of nodes.
%
%   - When multiple conditions on representation form are specified, only
%     the last one takes effect.
%
% $ History $
%   - Created by Dahua Lin, on Sep 9, 2006
%

%% parse conditions

if nargin < 2 || isempty(conds)
    cc = {};
    fc = [];
else
    nconds = length(conds);
    tf = false(1, nconds);
    pf = 0;
    for i = 1 : nconds
        curcond = conds{i};
        if ischar(curcond)
            switch curcond
                case {'adjmat', 'edgeset', 'adjlist'}
                    tf(i) = 1;
                    pf = i;
            end
        end
    end
    if pf > 0
        fc = conds{pf};
    else
        fc = [];
    end
    cc = conds(~tf);
end

%% parse representation form and delegate

if isnumeric(G) || islogical(G)

    if ~isempty(fc)
        chk_form('adjmat', fc);
    end
    gi = ginfo_adjmat(G);

elseif isstruct(G)
    if isfield(G, 'edges') && ~isfield(G, 'targets')

        if ~isempty(fc)
            chk_form('edgeset', fc);
        end
        gi = ginfo_edgeset(G);

    elseif isfield(G, 'targets') && ~isfield(G, 'edges')

        if ~isempty(fc)
            chk_form('adjlist', fc);
        end
        gi = ginfo_adjlist(G);

    else
        report_unreg();
    end
else
    report_unreg();
end

%% check conditions

if ~isempty(cc)
    ncc = numel(cc);
    for i = 1 : ncc
        chk_cond(G, gi, cc{i});
    end
end

%% Information function for different forms

function gi = ginfo_edgeset(G)

gi.form = 'edgeset';

if ~isfield(G, 'n') || ~isfield(G, 'edges')
    chkerr('The edgeset representation of a graph should have the fields n and edges');
end

ncols = size(G.edges, 2);
if ~isempty(G.edges) && ncols ~= 2 && ncols ~= 3
    chkerr('The non-empty edges should have 2 or 3 columns');
end

if isfield(G, 'nt')
    gi.type = 'bi';
    gi.n = G.n;
    gi.nt = G.nt;

else
    gi.type = 'ge';
    gi.n = G.n;
    gi.nt = G.n;

end

gi.valued = (ncols == 3);

function gi = ginfo_adjlist(G)

gi.form = 'adjlist';

if ~isfield(G, 'n') || ~isfield(G, 'targets')
    chkerr('The adjlist representation of a graph should have the fields n and targets');
end

tars = G.targets;
if ~iscell(tars)
    chkerr('The targets should be a cell array');
end
if ~isequal(size(tars), [1, G.n]) && ~isequal(size(tars), [G.n, 1])
    chkerr('The targets should be a cell array of size 1 x n or n x 1');
end

if isfield(G, 'nt')
    gi.type = 'bi';
    gi.n = G.n;
    gi.nt = G.nt;

else
    gi.type = 'ge';
    gi.n = G.n;
    gi.nt = G.n;

end

n = G.n;
gi.valued = false;
for i = 1 : n
    if ~isempty(tars{i})
        ncols = size(tars{i}, 2);
        gi.valued = (ncols > 1);
        break;
    end
end

function gi = ginfo_adjmat(G)

gi.form = 'adjmat';

gi.n = size(G, 1);
gi.nt = size(G, 2);

if gi.n == gi.nt
    gi.type = 'ge';
else
    gi.type = 'bi';
end

gi.valued = isnumeric(G);

%% Condition checking function

function chk_cond(G, gi, cond)

if ischar(cond)
    switch cond
        case 'square'
            if gi.n ~= gi.nt
                conderr('The graph is square');
            end
        case 'logical'
            if isnumeric(G)
                conderr('The value type is logical');
            end
        case 'numeric'
            if islogical(G)
                conderr('The value type is numeric');
            end
        otherwise
            error('sltoolbox:invalidarg', ...
                'Unknown condition: %s', cond);
    end
else
    if isscalar(cond)
        if gi.n ~= cond
            conderr(sprintf('n = %d', cond));
        end
    elseif isvector(cond) && length(cond) == 2
        if gi.n ~= cond(1) || gi.nt ~= cond(2)
            conderr(sprintf('n = %d, nt = %d', ...
                cond(1), cond(2)));
        end
    else
        error('sltoolbox:invalidarg', ...
            'Invalid condition is specified');
    end
end

%% Auxiliary function

function chk_form(df, cf)

if ~strcmp(df, cf)
    error('sltoolbox:invalidarg', ...
        'The given graph is like the form %s, which is not the required form %s', ...
        df, cf);
end

function report_unreg()

error('sltoolbox:invalidarg', ...
    'The graph form is unrecognizable');

function chkerr(msg)

error('sltoolbox:invalidarg', ...
    'Invalid graph representation: %s', msg);

function conderr(msg)

error('sltoolbox:invalidarg', ...
    'The given graph does not meet the requirement: %s', msg);

author: lmthang


text:
%%
 %  Copyright (c) 2014, Facebook, Inc.
 %  All rights reserved.
 %
 %  This source code is licensed under the BSD-style license found in the
 %  LICENSE file in the root directory of this source tree. An additional grant
 %  of patent rights can be found in the PATENTS file in the same directory.
 %
 %%
function level2_model=train_level2_cnn(train_labels, train_features, ...
    train_gfeatures, config)

level2_features_val = generate_level2_features_cnn(train_features,train_gfeatures,config);
clear train_features

tic;
fprintf('Training attribute classifier');
for attr_id = 1:size(train_labels,2)
    fprintf('%s ',config.ATTR_NAME{attr_id});
    %ignore those with labels = 0
    idx  = find(train_labels(:,attr_id) ~= 0);
    labels = (train_labels(idx,attr_id)>0)*2-1;
    level2_model{attr_id}.svm_weights = liblinear_do_train(labels,level2_features_val(idx,:));
end
fprintf('\n');
toc;
end

author: facebookarchive

repo: pose-aligned-deep-networks


text:
classdef dme_load < mp.dm_element
%MP.DME_LOAD  MATPOWER data model class for load data

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

    properties
        bus     %% bus index vector (all loads)
        pd      %% active power demand (p.u.) for constant power loads that are on
        qd      %% reactive power demand (p.u.) for constant power loads that are on
        pd_i    %% active power demand (p.u.) for constant current loads that are on
        qd_i    %% reactive power demand (p.u.) for constant current loads that are on
        pd_z    %% active power demand (p.u.) for constant impedance loads that are on
        qd_z    %% reactive power demand (p.u.) for constant impedance loads that are on
    end     %% properties

    methods
        function name = name(obj)
            name = 'load';
        end

        function label = label(obj)
            label = 'Load';
        end

        function label = labels(obj)
            label = 'Loads';
        end

        function name = cxn_type(obj)
            name = 'bus';
        end

        function name = cxn_idx_prop(obj)
            name = 'bus';
        end

        function names = main_table_var_names(obj)
            names = horzcat( main_table_var_names@mp.dm_element(obj), ...
                {'bus', 'pd', 'qd', 'pd_i', 'qd_i', 'pd_z', 'qd_z', ...
                'p', 'q'});
        end

        function vars = export_vars(obj)
            vars = {};
        end

        function nr = count(obj, dm)
            nr = count@mp.dm_element(obj, dm);
            if nr
                obj.bus = obj.tab.source_uid;
            end
        end

        function obj = update_status(obj, dm)
            %% get bus status info
            bs = dm.elements.bus.tab.status;    %% bus status

            %% update status of loads at isolated/offline buses
            obj.tab.status = obj.tab.status & bs(obj.bus);

            %% call parent to fill in on/off
            update_status@mp.dm_element(obj, dm);
        end

        function obj = build_params(obj, dm)
            obj.pd   = obj.tab.pd(obj.on) / dm.base_mva;
            obj.qd   = obj.tab.qd(obj.on) / dm.base_mva;
            obj.pd_i = obj.tab.pd_i(obj.on) / dm.base_mva;
            obj.qd_i = obj.tab.qd_i(obj.on) / dm.base_mva;
            obj.pd_z = obj.tab.pd_z(obj.on) / dm.base_mva;
            obj.qd_z = obj.tab.qd_z(obj.on) / dm.base_mva;
        end

        function TorF = pp_have_section_sum(obj, mpopt, pp_args)
            TorF = true;
        end

        function obj = pp_data_sum(obj, dm, rows, out_e, mpopt, fd, pp_args)
            %% call parent
            pp_data_sum@mp.dm_element(obj, dm, rows, out_e, mpopt, fd, pp_args);

            %% print load summary
            fprintf(fd, '  %-29s %12.1f MW', 'Total load', ...
                                            sum(obj.tab.p));
            if mpopt.model(1) ~= 'D'    %% AC model
                fprintf(fd, ' %12.1f MVAr', sum(obj.tab.q));
            end
            fprintf(fd, '\n');
            if obj.n ~= obj.nr
                fprintf(fd, '  %-29s %12.1f MW', '  online', ...
                                                sum(obj.tab.p(obj.on)));
                if mpopt.model(1) ~= 'D'    %% AC model
                    fprintf(fd, ' %12.1f MVAr', sum(obj.tab.q(obj.on)));
                end
                fprintf(fd, '\n');
            end
        end

        function TorF = pp_have_section_det(obj, mpopt, pp_args)
            TorF = true;
        end

        function h = pp_get_headers_det(obj, dm, out_e, mpopt, pp_args)
            h = [ pp_get_headers_det@mp.dm_element(obj, dm, out_e, mpopt, pp_args) ...
                {   '                             Power Consumption', ...
                    'Load ID    Bus ID   Status   P (MW)   Q (MVAr)', ...
                    '--------  --------  ------  --------  --------' } ];
            %%       1234567 123456789 -----1 12345678.0 1234567.9
        end

        function f = pp_get_footers_det(obj, dm, out_e, mpopt, pp_args)
            f = {'                            --------  --------',
                sprintf('%18s Total:%10.1f %9.1f', ...
                    '', sum(obj.tab.p(obj.on)), sum(obj.tab.q(obj.on)))};
        end

        function str = pp_data_row_det(obj, dm, k, out_e, mpopt, fd, pp_args)
            if obj.tab.status(k) && abs(obj.tab.p(k)) > 1e-5
                p = sprintf('%10.1f', obj.tab.p(k));
            else
                p = '       -  ';
            end
            if obj.tab.status(k) && abs(obj.tab.q(k)) > 1e-5
                q = sprintf('%9.1f', obj.tab.q(k));
            else
                q = '      -  ';
            end
            str = sprintf('%7d %9d %6d %10s %9s', ...
                obj.tab.uid(k), obj.tab.bus(k), obj.tab.status(k), p, q);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef LaunchVehicleGroundObjectSet < matlab.mixin.SetGet
    %LaunchVehicleGroundObjectSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        groundObjs(1,:) LaunchVehicleGroundObject = LaunchVehicleGroundObject.empty(1,0);

        lvdData LvdData
    end

    methods
        function obj = LaunchVehicleGroundObjectSet(lvdData)
            obj.lvdData = lvdData;
        end

        function addGroundObj(obj, groundObj)
            obj.groundObjs(end+1) = groundObj;
            groundObj.groundObjs = obj;
        end

        function removeGroundObj(obj, groundObj)
            obj.groundObjs([obj.groundObjs] == groundObj) = [];
        end

        function [listBoxStr, groundObjs] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.groundObjs))
                listBoxStr{end+1} = obj.groundObjs(i).name; %#ok<AGROW>
            end

            groundObjs = obj.groundObjs;
        end

        function groundObjects = getGroundObjsForInds(obj, inds)
            groundObjects = obj.groundObjs(inds);
        end

        function inds = getIndsForGroundObjs(obj, grdObjs)
            inds = find(ismember(obj.groundObjs, grdObjs));
        end

        function grndObj = getGroundObjAtInd(obj, ind)
            if(ind > 0 && ind <= length(obj.groundObjs))
                grndObj = obj.groundObjs(ind);
            else
                grndObj = LaunchVehicleGroundObject.empty(1,0);
            end
        end

        function numGrndObjs = getNumGroundObj(obj)
            numGrndObjs = length(obj.groundObjs);
        end

        function [grdObjAzGAStr, grdObjs] = getGrdObjAzGraphAnalysisTaskStrs(obj)
            grdObjs = obj.groundObjs;

            grdObjAzGAStr = cell(1,length(grdObjs));
            A = length(grdObjs);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(grdObjs))
                grdObjAzGAStr{i} = sprintf(sprintf('Ground Object %s Azimuth to S/C - "%s"',formSpec, grdObjs(i).name), i);
            end
        end

        function [grdObjElevGAStr, grdObjs] = getGrdObjElevGraphAnalysisTaskStrs(obj)
            grdObjs = obj.groundObjs;

            grdObjElevGAStr = cell(1,length(grdObjs));
            A = length(grdObjs);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(grdObjs))
                grdObjElevGAStr{i} = sprintf(sprintf('Ground Object %s Elevation to S/C - "%s"',formSpec, grdObjs(i).name), i);
            end
        end

        function [grdObjRngGAStr, grdObjs] = getGrdObjRangeGraphAnalysisTaskStrs(obj)
            grdObjs = obj.groundObjs;

            grdObjRngGAStr = cell(1,length(grdObjs));
            A = length(grdObjs);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(grdObjs))
                grdObjRngGAStr{i} = sprintf(sprintf('Ground Object %s Range to S/C - "%s"',formSpec, grdObjs(i).name), i);
            end
        end

        function [grdObjLoSGAStr, grdObjs] = getGrdObjLoSGraphAnalysisTaskStrs(obj)
            grdObjs = obj.groundObjs;

            grdObjLoSGAStr = cell(1,length(grdObjs));
            A = length(grdObjs);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(grdObjs))
                grdObjLoSGAStr{i} = sprintf(sprintf('Ground Object %s Line of Sight to S/C - "%s"',formSpec, grdObjs(i).name), i);
            end
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            for(i=1:length(obj.groundObjs))
                grdObj = obj.groundObjs(i);

                if(isempty(grdObj.groundObjs))
                    grdObj.groundObjs = obj;
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [x_min,x_max,y_min,y_max] = handle_cross_boundary(x_min,x_max,y_min,y_max,img_size)
    row = img_size(1);
    col = img_size(2);
    if x_min < 1
        x_min = 1;
    end
    if y_min < 1
        y_min = 1;
    end
    if x_max > row
        x_max = row;
    end
    if y_max > col
        y_max = col;
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [S,K,s,w,t,dfdx] = spm_dcm_mtf(P,M,U)
% computes transfer functions using the system's eigenspectrum
% FORMAT [S,K,s,w,t,dfdx] = spm_dcm_mtf(P,M,[U])
%
% P - model parameters
% M - model (with flow M.f and expansion point M.x and M.u)
% U - induces expansion around steady state (from spm_dcm_neural_x(P,M))
%
% S    - modulation transfer functions (complex)
% K    - Volterra kernels (real)
% s    - eigenspectrum (complex)
% w    - frequencies (Hz) = M.Hz
% t    - time (seconds)   = M.pst
% dfdx - Jacobian
%
% This routine uses the eigensolution of a dynamical systems Jacobian to
% complete the first-order Volterra terminals and transfer functions  in
% peristimulus and frequency space respectively.  The advantage of using
% the-solution is that unstable modes (eigenvectors of the Jacobian) can be
% conditioned (suppressed). Furthermore, this provides for a
% computationally efficient and transparent evaluation of the transfer
% functions that draws on linear signal processing theory in frequency
% space.
%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_mtf.m 6856 2016-08-10 17:55:05Z karl $

% get local linear approximation
%==========================================================================

% solve for steady-state - if required
%--------------------------------------------------------------------------
if nargin > 2
    M.x   = spm_dcm_neural_x(P,M);
end

% check expansion points
%--------------------------------------------------------------------------
try, M.x; catch, M.x = spm_dcm_x_neural(P,M.dipfit.model); end
try, M.u; catch, M.u = sparse(M.m,1);                      end

% frequencies and peristimulus time of interest
%--------------------------------------------------------------------------
w      = (1:64)';
t      = (0:64 - 1)'/64;
try, w = M.Hz(:);       end
try, t = M.pst(:);      end
try, t = M.dt*(1:M.N)'; end

% delay operator - if not specified already
%--------------------------------------------------------------------------
if isfield(M,'D')

    dfdx = spm_diff(M.f,M.x,M.u,P,M,1);
    dfdu = spm_diff(M.f,M.x,M.u,P,M,2);
    D    = M.D;

else

    if nargout(M.f) == 4
        [f,dfdx,D,dfdu] = feval(M.f,M.x,M.u,P,M);

    elseif nargout(M.f) == 3
        [f,dfdx,D]      = feval(M.f,M.x,M.u,P,M);
        dfdu            = spm_diff(M.f,M.x,M.u,P,M,2);

    elseif nargout(M.f) == 2
        [f,dfdx]        = feval(M.f,M.x,M.u,P,M);
        dfdu            = spm_diff(M.f,M.x,M.u,P,M,2);
        D               = 1;
    else
        dfdx            = spm_diff(M.f,M.x,M.u,P,M,1);
        dfdu            = spm_diff(M.f,M.x,M.u,P,M,2);
        D               = 1;
    end
end

% Jacobian and eigenspectrum
%==========================================================================
if nargout(M.g) == 2
    [g,dgdx] = feval(M.g,M.x,M.u,P,M);
else
    dgdx     = spm_diff(M.g,M.x,M.u,P,M,1);
end
dfdx  = D*dfdx;
dfdu  = D*dfdu;
try
   [v,s] = eig(full(dfdx),'nobalance');
catch
   v  = eye(size(dfdx));
   s  = NaN(size(dfdx));
end
s     = diag(s);

% condition unstable eigenmodes
%--------------------------------------------------------------------------
if max(w) > 1
    s = 1j*imag(s) + real(s) - exp(real(s));
else
    s = 1j*imag(s) + min(real(s),-1/32);
end

% Transfer functions
%==========================================================================

% transfer functions (FFT of kernel)
%--------------------------------------------------------------------------
nw    = size(w,1);            % number of frequencies
nt    = size(t,1);            % number of time bins
ng    = size(dgdx,1);         % number of outputs
nu    = size(dfdu,2);         % number of inputs
nk    = size(v,2);            % number of modes
S     = zeros(nw,ng,nu);
K     = zeros(nt,ng,nu);

% derivatives over modes
%--------------------------------------------------------------------------
dgdv  = dgdx*v;
dvdu  = pinv(v)*dfdu;
for j = 1:nu
    for i = 1:ng
        for k = 1:nk

            % transfer functions (FFT of kernel)
            %--------------------------------------------------------------
            Sk       = 1./(1j*2*pi*w - s(k));
            S(:,i,j) = S(:,i,j) + dgdv(i,k)*dvdu(k,j)*Sk;

            % kernels
            %--------------------------------------------------------------
            if nargout > 1
                Kk       = exp(s(k)*t);
                K(:,i,j) = K(:,i,j) + real(dgdv(i,k)*dvdu(k,j)*Kk);
            end

        end
    end
end

return

% NOTES: internal consistency with explicit Fourier transform of kernels
%==========================================================================

% augment and bi-linearise (with intrinsic delays)
%--------------------------------------------------------------------------
M.D       = D;
[M0,M1,L] = spm_bireduce(M,P);

% project onto spatial modes
%--------------------------------------------------------------------------
try,    L = M.U'*L; end

% kernels
%--------------------------------------------------------------------------
N         = length(t);
dt        = (t(2) - t(1));
[K0,K1]   = spm_kernels(M0,M1,L,N,dt);

% Transfer functions (FFT of kernel)
%--------------------------------------------------------------------------
S1        = fft(K1)*dt;
w1        = ((1:N) - 1)/(N*dt);
j         = w1 < max(w);
S1        = S1(j,1,1);
w1        = w1(j);

subplot(2,2,1), plot(t,K(:,1,1),t,K1(:,1,1));
title('kernels','fontsize',16)
xlabel('peristimulus time')

subplot(2,2,2), plot(w,real(S(:,1,1)),    w1,real(S1(:,1,1))); hold on
subplot(2,2,2), plot(w,imag(S(:,1,1)),':',w1,imag(S1(:,1,1)),':'); hold off
title('transfer functions','fontsize',16)
xlabel('frequency');

author: spm


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This is a testing environment for the files in the folder kernel/data
% 1. Based on data/contents, a list of required files is generated and it
%    is verified, that all files are present; additional files are listed.
% 2. All c-files are compiled.
% 3. All files are executed.
% 4. Check administrative modul distance.m
% 5. Test explicitly several specific distances
%==============================================================================

FAIRcheckFiles(mfilename);

%% Test module distance, start with syntax
distance('reset','distance','MIcc',...
  'tol',1e-7,'minT',0,'maxT',256,'nT',60,'minR',0,'maxR',256,'nR',60);
distance('disp');
[scheme,parameter] = distance
distance('clear');
distance('disp');

%% 5. setup test cases and run these
fprintf(' ... %s\n','generate hand data and setup transformation model:')
setup2DhandData
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e0);
level = 4; omega = ML{level}.omega; m = ML{level}.m;
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega(1,:),'out',0);
xc  = getCellCenteredGrid(omega(1,:),m);
Rc  = imgModel(R,omega(end,:),xc);

FAIRfigure(2);
subplot(1,2,1); viewImage2Dsc(T,omega,m,'title','template','colormap','gray')
subplot(1,2,2); viewImage2Dsc(R,omega,m,'title','reference','colormap','gray')

%% test several distances
distances = {'SSD','NCC','MIspline','MImex','NGFdot'};
clear fig; l = 0;

for k=1:length(distances);

  distance('reset','distance',distances{k});
  switch distance,
    case {'SSD','NCC'},
    case 'MIcc',
      %setup default parameter
      distance('set','tol',1e-7,...
        'minT',0,'maxT',256,'nT',60,'minR',0,'maxR',256,'nR',40);
      distance('disp');
    case 'NGFdot',
      %setup default parameter
      distance('set','edge',100);
      distance('disp');
  end;

  type = {'function','derivative','Gauss-Newton'};
  for k = 1:length(type),
    fctn = @(x) parseDistance(type{k},T,Rc,omega(end,:),m,x);
    l = l+1;
    fig(l)=checkDerivative(fctn,xc);
    ylabel([distance,'-',type{k}])
    builtin('pause',1);
  end;
end;

testEnd;
%==============================================================================


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deffn  Command assume {@var{x} @var{cond}}
%% @deffnx Command assume {@var{x} @var{cond} @var{cond2} @dots{}}
%% @deffnx Command assume {@var{x} @var{y} @dots{} @var{cond} @var{cond2} @dots{}}
%% Specify assumptions for a symbolic variable (replace existing).
%%
%% Example:
%% @example
%% @group
%% syms n x y
%% assume n integer
%% assume x y real
%% assumptions
%%   @result{} ans =
%%     @{
%%       [1,1] = n: integer
%%       [1,2] = x: real
%%       [1,3] = y: real
%%     @}
%% @end group
%% @end example
%%
%% To clear assumptions on a variable, use @code{assume x clear}, for example:
%% @example
%% @group
%% assume x y clear
%% assumptions
%%   @result{} ans =
%%     @{
%%       [1,1] = n: integer
%%     @}
%% @end group
%% @end example
%%
%% For more precise control over assumptions, @pxref{@@sym/assume}
%% and @pxref{@@sym/assumeAlso}.
%%
%% @seealso{@@sym/assume, @@sym/assumeAlso, assumptions, sym, syms}
%% @end deffn

function varargout = assume(varargin)

  assert (nargout == 0, 'assume: use functional form if you want output');

  assert (nargin > 1, 'assume: general algebraic assumptions are not supported');

  %% Find symbol/assumptions boundary and verify input
  valid_asm = assumptions ('possible');
  lastvar = -1;
  for n = 1:nargin
    assert (ischar (varargin{n}), 'assume: command form expects string inputs only')
    if (ismember (varargin{n}, valid_asm))
      if (lastvar < 0)
        lastvar = n - 1;
      end
    elseif (strcmp (varargin{n}, 'clear'))
      assert (n == nargin, 'assume: "clear" should be the final argument')
      assert (lastvar < 0, 'assume: should not combine "clear" with other assumptions')
      lastvar = n - 1;
    elseif (lastvar > 0)
      error('assume: cannot have symbols after assumptions')
    else
      assert (isvarname (varargin{n}), 'assume: only symbols can have assumptions')
    end
  end

  if (lastvar < 0)
    error ('assume: no assumptions were given')
  end
  if (lastvar == 0)
    error ('assume: cannot have only assumptions w/o symbols')
  end

  asm = varargin((lastvar+1):end);
  vars = varargin(1:lastvar);

  %% loop over each variable
  for n = 1:length (vars)
    vals = evalin ('caller', vars{n});
    newvals = cell(1, numel (vals));
    [newvals{:}] = assume (vals, asm{:});

    for i = 1:length (newvals)
      newx = newvals{i};
      xstr = newx.flat;

      % ---------------------------------------------
      % Muck around in the caller's namespace, replacing syms
      % that match 'xstr' (a string) with the 'newx' sym.
      context = 'caller';
      % ---------------------------------------------
      S = evalin(context, 'whos');
      evalin(context, '[];');  % clear 'ans'
      for i = 1:numel(S)
        obj = evalin(context, S(i).name);
        [newobj, flag] = symreplace(obj, xstr, newx);
        if flag, assignin(context, S(i).name, newobj); end
      end
      % ---------------------------------------------

    end
  end
end

%!error <if you want output>
%! a = assume('a', 'real')

%!error <cannot have only assumptions>
%! assume positive integer

%!error <no assumptions>
%! assume x y

%!error <should be the final>
%! assume x clear real

%!error <algebraic assumptions are not supported>
%! assume a>0

%!error <only symbols>
%! assume 'x/pi' integer

%!test
%! syms x
%! assume x positive
%! a = assumptions(x);
%! assert(strcmp(a, 'x: positive'))
%! assume x even
%! a = assumptions(x);
%! assert(strcmp(a, 'x: even'))

%!test
%! % multiple assumptions
%! syms x
%! assume x positive integer
%! [tilde, a] = assumptions(x, 'dict');
%! assert(a{1}.integer)
%! assert(a{1}.positive)

%!test
%! % does workspace
%! syms x positive
%! x2 = x;
%! f = sin(x);
%! assume x negative
%! a = assumptions(x);
%! assert(strcmp(a, 'x: negative'))
%! a = assumptions(x2);
%! assert(strcmp(a, 'x: negative'))
%! a = assumptions(f);
%! assert(strcmp(a, 'x: negative'))

%!error <undefined>
%! % does not create new variable x
%! clear x
%! assume x real

%!error <undefined>
%! % no explicit variable named x
%! clear x
%! f = 2*sym('x');
%! assume x real

%!test
%! % clear does workspace
%! syms x positive
%! f = 2*x;
%! assume x clear
%! assert (isempty (assumptions (f)));
%! assert (isempty (assumptions ()));

%!test
%! syms x y
%! f = sin (2*x);
%! assume x y real
%! assert (strcmp (assumptions (x), 'x: real'))
%! assert (strcmp (assumptions (y), 'y: real'))
%! assert (strcmp (assumptions (f), 'x: real'))

%!test
%! syms x y
%! f = sin (2*x);
%! assume x y positive even
%! assert (strcmp (assumptions (x), 'x: positive, even') || strcmp (assumptions (x), 'x: even, positive'))
%! assert (strcmp (assumptions (y), 'y: positive, even') || strcmp (assumptions (y), 'y: even, positive'))
%! assert (strcmp (assumptions (f), 'x: positive, even') || strcmp (assumptions (f), 'x: even, positive'))

%!test
%! % works from variable names not symbols
%! syms x y
%! a = [x y];
%! assume a real
%! assert (strcmp (assumptions (x), 'x: real'))
%! assert (strcmp (assumptions (y), 'y: real'))

%!test
%! % works from variable names not symbols
%! y = sym('x');
%! f = 2*y;
%! assume y real
%! assert (strcmp (assumptions (f), 'x: real'))

%!test
%! % matrix of symbols
%! syms a b c d
%! A = [a b; c d];
%! assume A real
%! assert (strcmp (assumptions (a), 'a: real'))
%! assert (strcmp (assumptions (b), 'b: real'))
%! assert (strcmp (assumptions (c), 'c: real'))
%! assert (strcmp (assumptions (d), 'd: real'))

%!test
%! % assume after symfun
%! clear x
%! syms f(x)
%! assume x real
%! assert (~ isempty (assumptions (formula (f))))
%! assert (~ isempty (assumptions (argnames (f))))

repo: octsympy


text:
function [rotations,flipLR] = getFlatRotations(view)
%
%  [rotations,flipLR] = getFlatRotations(view)
%
% Gets rotation field from FLAT
% Returns a [2x1]  pair of rotations if sucessful.
% Returns [0 0] if no such field present.
% Also returns a pair of flags indicating whether a L/R flip happens.

if (~strcmp(view.viewType,'Flat'))
    error('getFlatRotation called for non-flat view');
end

if (~isfield(view,'rotateImageDegrees'))
    rotations=[0 0];
else
    rotations=view.rotateImageDegrees;
end

if (~isfield(view,'flipLR'))
    flipLR = [0 0]; % If this flag is set, coordinates are flipped L/R : The rotation matrix becomes -c s; s c
    return;
else
    flipLR = view.flipLR;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug1425

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivityanalysis ft_connectivity_corr

% the bug pertains to a non-specific error when trying to do coherence computation on single trial data

% reproduce the issue
data = [];
data.trial{1} = randn(3,100);
data.time{1}  = (0:99)./100;
data.label    = {'chan01';'chan02';'chan03'};

cfg = [];
cfg.method = 'mtmfft';
cfg.output = 'fourier';
cfg.taper  = 'hanning';
freq = ft_freqanalysis(cfg, data);

cfg = [];
cfg.method = 'coh';
cfg.channelcmb = {'chan01' 'chan02'};
coh = ft_connectivityanalysis(cfg, freq);

% issue has been reproduced. It can be tracked down to ft_checkdata (fixcsd) where the singleton first dimension is not removed, leading to an incorrect dimensionality of the numeric data with respect to the dimord


text:
clear all; close all;
I=imread('circbw.tif');
J=watershed(I, 8);
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function out = ICC(cse,typ,dat)
%function to work out ICCs according to shrout & fleiss' schema (Shrout PE,
%Fleiss JL. Intraclass correlations: uses in assessing rater reliability.
%Psychol Bull. 1979;86:420-428).
%
% 'dat' is data whose rows represent different ratings or raters & whose
% columns represent different cases or targets being measured. Each target
% is assumed to be a random sample from a population of targets.
%
% 'cse' is either 1,2,3. 'cse' is: 1 if each target is measured by a
% different set of raters from a population of raters, 2 if each target is
% measured by the same raters, but that these raters are sampled from a
% population of raters, 3 if each target is measured by the same raters and
% these raters are the only raters of interest.
%
% 'typ' is either 'single' or 'k' & denotes whether the ICC is based on a
% single measurement or on an average of k measurements, where k = the
% number of ratings/raters.
%
% This has been tested using the example data in the paper by shrout & fleiss.
%
% Example: out = ICC(3,'k',S_Fdata)
% returns ICC(3,k) of data 'S_Fdata' to double 'out'.
%
% Kevin Brownhill, Imaging Sciences, KCL, London kevin.brownhill@kcl.ac.uk
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%number of raters/ratings
k = size(dat,2);
%number of targets
n = size(dat,1);
%mean per target
mpt = mean(dat,2);
%mean per rater/rating
mpr = mean(dat);
%get total mean
tm = mean(mpt);
%within target sum sqrs
WSS = sum(sum(bsxfun(@minus,dat,mpt).^2));
%within target mean sqrs
WMS = WSS / (n * (k - 1));
%between rater sum sqrs
RSS = sum((mpr - tm).^2) * n;
%between rater mean sqrs
RMS = RSS / (k - 1);
% %get total sum sqrs
% TSS = sum(sum((dat - tm).^2));
%between target sum sqrs
BSS = sum((mpt - tm).^2) * k;
%between targets mean squares
BMS = BSS / (n - 1);
%residual sum of squares
ESS = WSS - RSS;
%residual mean sqrs
EMS = ESS / ((k - 1) * (n - 1));
switch cse
    case 1
        switch typ
            case 'single'
                out = (BMS - WMS) / (BMS + (k - 1) * WMS);
            case 'k'
                out = (BMS - WMS) / BMS;
            otherwise
               error('Wrong value for input typ')
        end
    case 2
        switch typ
            case 'single'
                out = (BMS - EMS) / (BMS + (k - 1) * EMS + k * (RMS - EMS) / n);
            case 'k'
                out = (BMS - EMS) / (BMS + (RMS - EMS) / n);
            otherwise
               error('Wrong value for input typ')
        end
    case 3
        switch typ
            case 'single'
                out = (BMS - EMS) / (BMS + (k - 1) * EMS);
            case 'k'
                out = (BMS - EMS) / BMS;
            otherwise
               error('Wrong value for input typ')
        end
    otherwise
        error('Wrong value for input cse')
end


text:
% POP_REJMENU - Main menu for rejecting trials in an EEG dataset
%
% Usage: >> pop_rejmenu(INEEG, typerej);
%
% Inputs:
%   INEEG      - input dataset
%   typerej    - data to reject on (0 = component activations;
%                1 = raw electrode data). {Default: 1 = reject on raw data}
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGLAB, POP_EEGPLOT, pop_eegthresh, POP_REJTREND
% POP_REJKURT, POP_JOINTPROB, POP_REJSPEC

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

function pop_rejmenu( EEG, icacomp );

if icacomp == 0
	if isempty( EEG.icasphere )
		disp('Error: you must first run ICA on the data'); return;
	end
end

if icacomp == 1 	rejtitle = 'Reject trials using data statistics - pop_rejmenu()'; tagmenu = 'rejtrialraw';
else            	rejtitle = 'Reject trials using component activity statistics - pop_rejmenu()'; tagmenu = 'rejtrialica';
end;

if ~isempty( findobj('tag', tagmenu))
	error('cannot open two identical windows; close the first one first');
end

figure('visible', 'off', 'numbertitle', 'off', 'name', rejtitle, 'tag', tagmenu);

% definition of callbacks
% -----------------------
checkstatus = [ 'rejstatus = get( findobj(''parent'', gcbf, ''tag'', ''rejstatus''), ''value'');' ...
                'if rejstatus == 3,' ...
                '    EEG.reject.disprej = {};' ...
				'    if get( findobj(''parent'', gcbf, ''tag'', ''IManual''), ''value''), EEG.reject.disprej{1} = ''manual''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IThresh''), ''value''), EEG.reject.disprej{2} = ''thresh''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IConst''), ''value''), EEG.reject.disprej{3} = ''const''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IEnt''), ''value''), EEG.reject.disprej{4} = ''jp''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IKurt''), ''value''), EEG.reject.disprej{5} = ''kurt''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IFreq''), ''value''), EEG.reject.disprej{6} = ''freq''; end;' ...
                'end;' ...
                'rejstatus = rejstatus-1;' ]; % from 1-3 range, go to 0-2 range

if icacomp, ICAPREFIX = '';
else        ICAPREFIX = 'ica';
end
% tmp_comall is used when returning from eegplot
tmp_comall =      [ 'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''mantrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejmanual)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''threshtrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejthresh)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''freqtrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejfreq)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''consttrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejconst)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''enttrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejjp)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''kurttrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejkurt)));' ];

cb_manual = [ 	checkstatus ...
				'pop_eegplot( EEG, ' int2str( icacomp ) ', rejstatus, 0,''' tmp_comall ''');' ...
				'clear rejstatus;' ];

% -----------------------------------------------------

cb_compthresh = [ '  posthresh = get( findobj(''parent'', gcbf, ''tag'', ''threshpos''), ''string'' );' ...
                  '  negthresh = get( findobj(''parent'', gcbf, ''tag'', ''threshneg''), ''string'' );' ...
                  '  startime  = get( findobj(''parent'', gcbf, ''tag'', ''threshstart''), ''string'' );' ...
                  '  endtime   = get( findobj(''parent'', gcbf, ''tag'', ''threshend''), ''string'' );' ...
                  '  elecrange = get( findobj(''parent'', gcbf, ''tag'', ''threshelec''), ''string'' );' ...
                  checkstatus ...
				  '[EEG Itmp LASTCOM] = pop_eegthresh( EEG,' int2str(icacomp) ...
				  ', elecrange, negthresh, posthresh, str2num(startime)/1000, str2num(endtime)/1000, rejstatus, 0,''' tmp_comall ''');' ...
				  'EEG = eegh(LASTCOM, EEG);' ...
				  'clear com Itmp elecrange posthresh negthresh startime endtime rejstatus;' ];

%                  '  set(findobj(''parent'', gcbf, ''tag'', ''threshtrial''), ''string'', num2str(EEG.trials - length(Itmp)));' ...
%				  'eegh(LASTCOM);' ...
%				  'clear Itmp elecrange posthresh negthresh startime endtime rejstatus;' ];

% -----------------------------------------------------
cb_compfreq   = [ '  posthresh = get( findobj(''parent'', gcbf, ''tag'', ''freqpos''), ''string'' );' ...
                  '  negthresh = get( findobj(''parent'', gcbf, ''tag'', ''freqneg''), ''string'' );' ...
                  '  startfreq = get( findobj(''parent'', gcbf, ''tag'', ''freqstart''), ''string'' );' ...
                  '  endfreq   = get( findobj(''parent'', gcbf, ''tag'', ''freqend''), ''string'' );' ...
                  '  elecrange = get( findobj(''parent'', gcbf, ''tag'', ''freqelec''), ''string'' );' ...
                  checkstatus ...
				  '[EEG Itmp LASTCOM] = pop_rejspec( EEG,' int2str(icacomp) ...
				  ', elecrange, negthresh, posthresh, startfreq, endfreq, rejstatus, 0,''' tmp_comall ''');' ...
				  'EEG = eegh(LASTCOM, EEG);' ...
				  'clear Itmp elecrange posthresh negthresh startfreq endfreq rejstatus;' ];

% -----------------------------------------------------
cb_compconstrej = [ '  minslope   = get( findobj(''parent'', gcbf, ''tag'', ''constpnts''), ''string'' );' ...
					'  minstd    = get( findobj(''parent'', gcbf, ''tag'', ''conststd''), ''string'' );' ...
					'  elecrange = get( findobj(''parent'', gcbf, ''tag'', ''constelec''), ''string'' );' ...
                    checkstatus ...
				    '[rej LASTCOM] = pop_rejtrend(EEG,' int2str(icacomp) ', elecrange, ''' ...
					int2str(EEG.pnts) ''', minslope, minstd, rejstatus, 0,''' tmp_comall ''');' ...
				    'EEG = eegh(LASTCOM, EEG);' ...
				    'clear rej elecrange minslope minstd rejstatus;' ];

% -----------------------------------------------------
cb_compenthead = [ '  locthresh  = get( findobj(''parent'', gcbf, ''tag'', ''entloc''), ''string'' );', ...
				   '  globthresh = get( findobj(''parent'', gcbf, ''tag'', ''entglob''), ''string'' );', ...
				   '  elecrange  = get( findobj(''parent'', gcbf, ''tag'', ''entelec''), ''string'' );' ];
cb_compenttail = [ '  set( findobj(''parent'', gcbf, ''tag'', ''entloc''), ''string'', num2str(locthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''entglob''), ''string'', num2str(globthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''enttrial''), ''string'', num2str(nrej) );' ...
				   'EEG = eegh(LASTCOM, EEG);' ...
				   'clear nrej elecrange locthresh globthresh rejstatus;' ];
cb_compentplot =  [ cb_compenthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_jointprob( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, 0, 0,''' tmp_comall ''');', ...
				    cb_compenttail ];
cb_compentcalc =  [ cb_compenthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_jointprob( EEG, ' int2str(icacomp) ...
					', [ str2num(elecrange) ], str2num(locthresh), str2num(globthresh), 0, 0);', ...
 				    cb_compenttail ];
cb_compenteeg  =  [ cb_compenthead ...
                    checkstatus ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_jointprob( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, rejstatus, 0, 1,''' tmp_comall ''');', ...
 				    cb_compenttail ];

% -----------------------------------------------------
cb_compkurthead =[ '  locthresh  = get( findobj(''parent'', gcbf, ''tag'', ''kurtloc''), ''string'' );', ...
				   '  globthresh = get( findobj(''parent'', gcbf, ''tag'', ''kurtglob''), ''string'' );', ...
				   '  elecrange  = get( findobj(''parent'', gcbf, ''tag'', ''kurtelec''), ''string'' );' ];
cb_compkurttail =[ '  set( findobj(''parent'', gcbf, ''tag'', ''kurtloc''), ''string'', num2str(locthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''kurtglob''), ''string'', num2str(globthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''kurttrial''), ''string'', num2str(nrej) );' ...
				   'EEG = eegh(LASTCOM, EEG);' ...
				   'clear nrej elecrange locthresh globthresh rejstatus;' ];
cb_compkurtplot = [ cb_compkurthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_rejkurt( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, 0, 0,''' tmp_comall ''');', ...
				    cb_compkurttail ];
cb_compkurtcalc = [ cb_compkurthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_rejkurt( EEG, ' int2str(icacomp) ...
					', [ str2num(elecrange) ], str2num(locthresh), str2num(globthresh), 0, 0);', ...
 				    cb_compkurttail ];
cb_compkurteeg  = [ cb_compkurthead ...
                    checkstatus ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_rejkurt( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, rejstatus, 0, 1,''' tmp_comall ''');', ...
 				    cb_compkurttail ];

% -----------------------------------------------------
cb_reject =      [ 'set( findobj(''parent'', gcbf, ''tag'', ''rejstatus''), ''value'', 3);' ... % force status to 3
				    checkstatus ...
				   '[EEG LASTCOM] = eeg_rejsuperpose(EEG,' int2str(icacomp) ',1,1,1,1,1,1,1); EEG = eegh(LASTCOM, EEG);' ...
                   'if isempty(find(EEG.reject.rejglobal)),' ...
                   '   warndlg2(strvcat(''No epoch selected...'', ''When using thresholding, click update'',''marks in the EEG plotting window''));' ...
				   'else,' ...
                   '   [EEG LASTCOM] = pop_rejepoch( EEG, EEG.reject.rejglobal, 1);' ...
				   '   if ~isempty(LASTCOM), ' ...
				   '      EEG = eegh(LASTCOM, EEG); [ALLEEG EEG CURRENTSET LASTCOM] = pop_newset(ALLEEG, EEG, CURRENTSET); eegh(LASTCOM);' ...
				   '   end; eeglab redraw; close(gcbf);' ...
                   'end;' ];

cb_clear =       [ 'close gcbf; EEG = rmfield( EEG, ''reject''); EEG.reject.rejmanual = [];' ...
				   'EEG=eeg_checkset(EEG); pop_rejmenu(EEG,' int2str(icacomp) ');' ];

cb_close =       [ 'close gcbf;' ...
				   'disp(''Marks stored in dataset'');' ...
				   '[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);' ...
				   'eegh(''[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);'');'];

lisboxoptions = { 'string', [ 'Show only the new trials marked for rejection by the measure selected above|' ...
                  'Show previous and new trials marked for rejection by the measure selected above|' ...
                  'Show all trials marked for rejection by the measure selected above or checked below'], 'tag', 'rejstatus', 'value', 1, 'callback', ...
                  [ 'if get(gcbo, ''value'') == 3,' ...
                    '   set(findobj(''parent'', gcbf, ''style'', ''checkbox''), ''enable'', ''on'');' ...
                    'else' ...
                    '   set(findobj(''parent'', gcbf, ''style'', ''checkbox''), ''enable'', ''off'');' ...
                    'end;' ] };

	chanliststr = [fastif(icacomp,'Electrode(s)','Component(s)') ];
	chanlistval = fastif(icacomp, [ '1:' int2str(EEG.nbchan) ], [ '1:' int2str(size(EEG.icaweights,1)) ]);

	% assess previous rejections
	% --------------------------
	sizeman = 0;
	sizethresh = 0;
	sizetrend = 0;
	sizejp = 0;
	sizekurt = 0;
	sizespec = 0;
	if icacomp == 1
		if ~isempty(EEG.reject.rejmanual), sizeman   = length(find(EEG.reject.rejmanual)); end
		if ~isempty(EEG.reject.rejconst),  sizetrend = length(find(EEG.reject.rejconst)); end
		if ~isempty(EEG.reject.rejjp),     sizejp    = length(find(EEG.reject.rejjp)); end
		if ~isempty(EEG.reject.rejkurt),   sizekurt  = length(find(EEG.reject.rejkurt)); end
		if ~isempty(EEG.reject.rejfreq),   sizespec  = length(find(EEG.reject.rejfreq)); end
	else
		if ~isempty(EEG.reject.icarejmanual), sizeman   = length(find(EEG.reject.icarejmanual)); end
		if ~isempty(EEG.reject.icarejconst),  sizetrend = length(find(EEG.reject.icarejconst)); end
		if ~isempty(EEG.reject.icarejjp),     sizejp    = length(find(EEG.reject.icarejjp)); end
		if ~isempty(EEG.reject.icarejkurt),   sizekurt  = length(find(EEG.reject.icarejkurt)); end
		if ~isempty(EEG.reject.icarejfreq),   sizespec  = length(find(EEG.reject.icarejfreq)); end;
	end

    stdl = [0.25 1.2 0.8 1.2 0.8]; % standard line
	titl = [0.9 0.18 1.55]; % title line
    geometry = { 	[0.883 0.195 .2 .45 .4 .4] ...
					[1] titl stdl stdl stdl stdl ...
					[1] titl stdl stdl stdl  ...
					[1] titl stdl stdl stdl  ...
					[1] titl stdl stdl stdl  ...
					[1] titl stdl stdl stdl stdl ...
					[1] [1] [1] [1 1 1] [1 1 1]  ...
					[1] [1 1 1]};

	listui = {{'Style', 'text', 'string', 'Mark trials by appearance', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butmanual', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejmanualcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejmanualcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] },...
	{ } { 'Style', 'pushbutton', 'string', fastif(icacomp,'Scroll Data','Scroll Acts.'), 'callback', cb_manual }, ...
	{ 'Style', 'text', 'string', 'Marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizeman), 'tag', 'mantrial' }, ...
	{ }, ...
	... % ---------------------------------------------------------------------------
	{ 'Style', 'text', 'string', 'Find abnormal values', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butthresh',  ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejthreshcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejthreshcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', ['Upper limit(s) ' fastif(icacomp,'(uV)', '(std. dev.)')] }, ...
	{ 'Style', 'edit', 'string', '25', 'tag', 'threshpos' }, ...
	{ 'Style', 'text', 'string', ['Lower limit(s) ' fastif(icacomp,'(uV)', '(std. dev.)')] }, ...
	{ 'Style', 'edit', 'string', '-25', 'tag', 'threshneg' }, ...
	...
	{ }, { 'Style', 'text', 'string', 'Start time(s) (ms)' }, ...
	{ 'Style', 'edit', 'string', int2str(EEG.xmin*1000), 'tag', 'threshstart' }, ...
	{ 'Style', 'text', 'string', 'Ending time(s) (ms)' }, ...
	{ 'Style', 'edit', 'string', int2str(EEG.xmax*1000), 'tag', 'threshend' }, ...
	...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'threshelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizethresh), 'tag', 'threshtrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calc / Plot', 'callback', cb_compthresh }, ...
	{ }, { },{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_eegthresh'');' }, ...
	... % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find abnormal trends', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'buttrend', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejconstcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejconstcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', ['Max slope ' fastif(icacomp, ...
           '(uV/epoch)', ...
           '(std. dev./epoch)') ] }, ...
	{ 'Style', 'edit', 'string', '50', 'tag', 'constpnts' }, ...
	{ 'Style', 'text', 'string', 'R-squared limit (0 to 1)' }, ...
	{ 'Style', 'edit', 'string', '0.3', 'tag', 'conststd'  }, ...
    ...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'constelec'  }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizetrend), 'tag', 'consttrial'  }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calc / Plot', 'callback', cb_compconstrej }, ...
	{ }, { }, { 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_rejtrend'');' }, ...
	...  % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find improbable data', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '',  'tag', 'butjp', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejjpcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejjpcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', fastif(icacomp, 'Single-channel limit (std. dev.)', 'Single-comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', fastif(icacomp, '5', '20'), 'tag', 'entloc' }, ...
	{ 'Style', 'text', 'string', fastif(icacomp, 'All channels limit (std. dev.)', 'All comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', '5', 'tag', 'entglob' }, ...
    ...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'entelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizejp), 'tag', 'enttrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calculate', 'callback', cb_compentcalc }, ...
	{ 'Style', 'pushbutton', 'string', 'Scroll Data', 'callback', cb_compenteeg }, ...
	{ 'Style', 'pushbutton', 'string', 'PLOT', 'callback', cb_compentplot }, ...
	{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_jointprob'');' }, ...
	...  % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find abnormal distributions', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butkurt', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejkurtcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejkurtcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', fastif(icacomp, 'Single-channel limit (std. dev.)', 'Single-comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', fastif(icacomp, '5', '20'), 'tag', 'kurtloc' }, ...
	{ 'Style', 'text', 'string', fastif(icacomp, 'All channels limit (std. dev.)', 'All comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', '5', 'tag', 'kurtglob' }, ...
    ...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'kurtelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizekurt), 'tag', 'kurttrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calculate', 'callback', cb_compkurtcalc }, ...
	{ 'Style', 'pushbutton', 'string', 'Scroll Data', 'callback', cb_compkurteeg }, ...
	{ 'Style', 'pushbutton', 'string', 'PLOT', 'callback', cb_compkurtplot }, ...
	{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_rejkurt'');' }, ...
	... % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find abnormal spectra (slow)', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butspec', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejfreqcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejfreqcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', 'Upper limit(s) (dB)' }, ...
	{ 'Style', 'edit', 'string', '25', 'tag', 'freqpos' }, ...
	{ 'Style', 'text', 'string', 'Lower limit(s) (dB)' }, ...
	{ 'Style', 'edit', 'string', '-25', 'tag', 'freqneg' }, ...
	...
	{ }, { 'Style', 'text', 'string', 'Low frequency(s) (Hz)' }, ...
	{ 'Style', 'edit', 'string', '0', 'tag', 'freqstart' }, ...
	{ 'Style', 'text', 'string', 'High frequency(s) (Hz)' }, ...
	{ 'Style', 'edit', 'string', '50', 'tag', 'freqend' }, ...
	...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'freqelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizespec), 'tag', 'freqtrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calc / Plot', 'callback', cb_compfreq }, ...
	{ }, { },{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_rejspec'');' }, ...
    ...
    {}, ...  % ---------------------------------------------------------------------------
	...
	{ 'Style', 'text', 'string', 'Plotting options', 'fontweight', 'bold' }, ...
	{ 'style', 'listbox', lisboxoptions{:}, 'value', 3 }, ...
	{ 'style', 'checkbox', 'String', 'Abnormal appearance', 'tag', 'IManual', 'value', 1, ...
      'callback', 'set(gcbo, ''value'', 1); disp(''This checkbox must be set so that you may mark/unmark data epochs'');'}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal values', 'tag', 'IThresh', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal trends', 'tag', 'IConst', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Improbable epochs', 'tag', 'IEnt', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal distributions', 'tag', 'IKurt', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal spectra', 'tag', 'IFreq', 'value', 1}, ...
	...
	{ }, ...
 	{ 'Style', 'pushbutton', 'string', 'CLOSE (KEEP MARKS)', 'callback', cb_close }, ...
 	{ 'Style', 'pushbutton', 'string', 'CLEAR ALL MARKS', 'callback', cb_clear  }, ...
	{ 'Style', 'pushbutton', 'string', 'REJECT MARKED TRIALS', 'callback', cb_reject }};

	allh = supergui( gcf, geometry, [], listui{:});

%	{ 'style', 'checkbox', 'String', ['Include ' fastif(icacomp, 'ica data', 'raw data')], 'tag', 'IOthertype', 'value', 1}, { }, ...
set(gcf, 'userdata', { allh });
set(findobj('parent', gcf', 'tag', 'butmanual'), 'backgroundcolor', EEG.reject.rejmanualcol);
set(findobj('parent', gcf', 'tag', 'butthresh'), 'backgroundcolor', EEG.reject.rejthreshcol);
set(findobj('parent', gcf', 'tag', 'buttrend'),  'backgroundcolor', EEG.reject.rejconstcol);
set(findobj('parent', gcf', 'tag', 'butjp'),     'backgroundcolor', EEG.reject.rejjpcol);
set(findobj('parent', gcf', 'tag', 'butkurt'),   'backgroundcolor', EEG.reject.rejkurtcol);
set(findobj('parent', gcf', 'tag', 'butspec'),   'backgroundcolor', EEG.reject.rejfreqcol);
set( findobj('parent', gcf, 'tag', 'rejstatus'), 'style', 'popup');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function T = convert_to_table(CPD, domain, evidence)
% CONVERT_TO_TABLE Convert a discrete CPD to a table
% T = convert_to_table(CPD, domain, evidence)
%
% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.
% The resulting table can easily be converted to a potential.

CPT = CPD.table;
odom = domain(~isemptycell(evidence(domain)));
vals = cat(1, evidence{odom});
map = find_equiv_posns(odom, domain);
index = mk_multi_index(length(domain), map, vals);
T = CPT(index{:});

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt