text:
% RANDN Normally distributed pseudorandom numbers.
%    R = RANDN(N) returns an N-by-N matrix containing pseudorandom values drawn
%    from the standard normal distribution.  RANDN(M,N) or RANDN([M,N]) returns
%    an M-by-N matrix. RANDN(M,N,P,...) or RANDN([M,N,P,...]) returns an
%    M-by-N-by-P-by-... array. RANDN returns a scalar.  RANDN(SIZE(A)) returns
%    an array the same size as A.
%
%    Note: The size inputs M, N, P, ... should be nonnegative integers.
%    Negative integers are treated as 0.
%
%    R = RANDN(..., CLASSNAME) returns an array of normal values of the
%    specified class. CLASSNAME can be 'double' or 'single'.
%
%    R = RANDN(..., 'like', Y) returns an array of normal values of the
%    same class as Y.
%
%    The sequence of numbers produced by RANDN is determined by the settings of
%    the uniform random number generator that underlies RAND, RANDN, and RANDI.
%    RANDN uses one or more uniform random values to create each normal random
%    value.  Control that shared random number generator using RNG.
%
%    Examples:
%
%       Example 1: Generate values from a normal distribution with mean 1
%        and standard deviation 2.
%          r = 1 + 2.*randn(100,1);
%
%       Example 2: Generate values from a bivariate normal distribution with
%       specified mean vector and covariance matrix.
%          mu = [1 2];
%          Sigma = [1 .5; .5 2]; R = chol(Sigma);
%          z = repmat(mu,100,1) + randn(100,2)*R;
%
%       Example 3: Reset the random number generator used by RAND, RANDI, and
%       RANDN to its default startup settings, so that RANDN produces the same
%       random numbers as if you restarted MATLAB.
%          rng('default');
%          randn(1,5)
%
%       Example 4: Save the settings for the random number generator used by
%       RAND, RANDI, and RANDN, generate 5 values from RANDN, restore the
%       settings, and repeat those values.
%          s = rng
%          z1 = randn(1,5)
%          rng(s);
%          z2 = randn(1,5) % z2 contains exactly the same values as z1
%
%       Example 5: Reinitialize the random number generator used by RAND,
%       RANDI, and RANDN with a seed based on the current time.  RANDN will
%       return different values each time you do this.  NOTE: It is usually
%       not necessary to do this more than once per MATLAB session.
%          rng('shuffle');
%          randn(1,5)
%
%    See Replace Discouraged Syntaxes of rand and randn to use RNG to replace
%    RANDN with the 'seed' or 'state' inputs.
%
%    See also RAND, RANDI, RNG, RANDSTREAM, RANDSTREAM/RANDN
%
%    Reference page in Doc Center
%       doc randn
%
%    Other functions named randn
%
%       codistributed/randn        distributed/randn    RandStream/randn
%       codistributor1d/randn      gpuArray/randn       ts/randn
%       codistributor2dbc/randn
%

author: jmaih


text:
function run_demoTracklets(base_dir,calib_dir)
% KITTI RAW DATA DEVELOPMENT KIT
%
% This tool displays the images and the object labels for the benchmark and
% provides an entry point for writing your own interface to the data set.
% Before running this tool, set root_dir to the directory where you have
% downloaded the dataset. 'root_dir' must contain the subdirectory
% 'training', which in turn contains 'image_2', 'label_2' and 'calib'.
% For more information about the data format, please look into readme.txt.
%
% Input arguments:
% base_dir .... absolute path to sequence base directory (ends with _sync)
% calib_dir ... absolute path to directory that contains calibration files
%
% Usage:
%   SPACE: next frame
%   '-':   last frame
%   'x':   +10 frames
%   'y':   -10 frames
%   'q':   quit
%
% Occlusion Coding:
%   green:  not occluded
%   yellow: partly occluded
%   red:    fully occluded
%   white:  unknown

% clear and close everything
clear all; close all; dbstop error; clc;
disp('======= KITTI DevKit Demo =======');

% options (modify this to select your sequence)
% the base_dir must contain:
%   - the data directories (image_00, image_01, ..)
%   - the tracklet file (tracklet_labels.xml)
% the calib directory must contain:
%   - calib_cam_to_cam.txt
%   - calib_velo_to_cam.txt
% cameras:
%   - 0 = left grayscale
%   - 1 = right grayscale
%   - 2 = left color
%   - 3 = right color
if nargin<1
  base_dir  = '/mnt/karlsruhe_dataset/2011_09_26/2011_09_26_drive_0009_sync';
end
if nargin<2
  calib_dir = '/mnt/karlsruhe_dataset/2011_09_26';
end
cam       = 2; % 0-based index

% get image sub-directory
image_dir = fullfile(base_dir, sprintf('/image_%02d/data', cam));

% get number of images for this dataset
nimages = length(dir(fullfile(image_dir, '*.png')));

% set up figure
gh = visualization('init',image_dir);

% read calibration for the day
[veloToCam, K] = loadCalibration(calib_dir);

% read tracklets for the selected sequence
tracklets = readTrackletsMex([base_dir '/tracklet_labels.xml']);

% extract tracklets
% LOCAL OBJECT COORDINATE SYSTEM:
%   x -> facing right
%   y -> facing forward
%   z -> facing up
for it = 1:numel(tracklets)

  % shortcut for tracklet dimensions
  w = tracklets{it}.w;
  h = tracklets{it}.h;
  l = tracklets{it}.l;

  % set bounding box corners
  corners(it).x = [l/2, l/2, -l/2, -l/2, l/2, l/2, -l/2, -l/2]; % front/back
  corners(it).y = [w/2, -w/2, -w/2, w/2, w/2, -w/2, -w/2, w/2]; % left/right
  corners(it).z = [0,0,0,0,h,h,h,h];

  % get translation and orientation
  t{it} = [tracklets{it}.poses(1,:); tracklets{it}.poses(2,:); tracklets{it}.poses(3,:)];
  rz{it} = wrapToPi(tracklets{it}.poses(6,:));
  occlusion{it} = tracklets{it}.poses(8,:);
end

% 3D bounding box faces (indices for corners)
face_idx = [ 1,2,6,5   % front face
             2,3,7,6   % left face
             3,4,8,7   % back face
             4,1,5,8]; % right face

% main loop (start at first image of sequence)
img_idx = 0;
while 1

  % visualization update for next frame
  visualization('update',image_dir,gh,img_idx,nimages);

  % compute bounding boxes for visible tracklets
  for it = 1:numel(tracklets)

    % get relative tracklet frame index (starting at 0 with first appearance;
    % xml data stores poses relative to the first frame where the tracklet appeared)
    pose_idx = img_idx-tracklets{it}.first_frame+1; % 0-based => 1-based MATLAB index

    % only draw tracklets that are visible in current frame
    if pose_idx<1 || pose_idx>(size(tracklets{it}.poses,2))
      continue;
    end

    % compute 3d object rotation in velodyne coordinates
    % VELODYNE COORDINATE SYSTEM:
    %   x -> facing forward
    %   y -> facing left
    %   z -> facing up
    R = [cos(rz{it}(pose_idx)), -sin(rz{it}(pose_idx)), 0;
         sin(rz{it}(pose_idx)),  cos(rz{it}(pose_idx)), 0;
                             0,                      0, 1];

    % rotate and translate 3D bounding box in velodyne coordinate system
    corners_3D      = R*[corners(it).x;corners(it).y;corners(it).z];
    corners_3D(1,:) = corners_3D(1,:) + t{it}(1,pose_idx);
    corners_3D(2,:) = corners_3D(2,:) + t{it}(2,pose_idx);
    corners_3D(3,:) = corners_3D(3,:) + t{it}(3,pose_idx);
    corners_3D      = (veloToCam{cam+1}*[corners_3D; ones(1,size(corners_3D,2))]);

    % generate an orientation vector and compute coordinates in velodyneCS
    orientation_3D      = R*[0.0, 0.7*l; 0.0, 0.0; 0.0, 0.0];
    orientation_3D(1,:) = orientation_3D(1,:) + t{it}(1, pose_idx);
    orientation_3D(2,:) = orientation_3D(2,:) + t{it}(2, pose_idx);
    orientation_3D(3,:) = orientation_3D(3,:) + t{it}(3, pose_idx);
    orientation_3D      = (veloToCam{cam+1}*[orientation_3D; ones(1,size(orientation_3D,2))]);

    % only draw 3D bounding box for objects in front of the image plane
    if any(corners_3D(3,:)<0.5) || any(orientation_3D(3,:)<0.5)
      continue;
    end

    % project the 3D bounding box into the image plane
    corners_2D     = projectToImage(corners_3D, K);
    orientation_2D = projectToImage(orientation_3D, K);
    drawBox3D(gh,occlusion{it}(pose_idx),corners_2D,face_idx,orientation_2D)

    % compute and draw the 2D bounding box from the 3D box projection
    box.x1 = min(corners_2D(1,:));
    box.x2 = max(corners_2D(1,:));
    box.y1 = min(corners_2D(2,:));
    box.y2 = max(corners_2D(2,:));
    drawBox2D(gh,box,occlusion{it}(pose_idx),tracklets{it}.objectType)
  end

  % force drawing and tiny user interface
  waitforbuttonpress;
  key = get(gcf,'CurrentCharacter');
  switch lower(key)
    case 'q',  break;                                 % quit
    case '-',  img_idx = max(img_idx-1,  0);          % previous frame
    case 'x',  img_idx = min(img_idx+100,nimages-1);  % +100 frames
    case 'y',  img_idx = max(img_idx-100,0);          % -100 frames
    otherwise, img_idx = min(img_idx+1,  nimages-1);  % next frame
  end

end

% clean up
close all;

repo: MSCKF


text:
% ----------------------------------------------------------------------------------------------------------------
function net = vid_create_net(varargin)
% Very similar to vanilla AlexNet from MatConvNet examples,
% but with smaller stride at conv1 and no padding
% Used to generate the network described in  the paper
% "Fully-Convolutional Siamese Networks for Object Tracking"
% ----------------------------------------------------------------------------------------------------------------
    opts.exemplarSize = [127 127];
    opts.instanceSize = [255 255];
    opts.scale = 1 ;
    opts.initBias = 0.1 ;
    opts.weightDecay = 1 ;
    %opts.weightInitMethod = 'xavierimproved' ;
    opts.weightInitMethod = 'gaussian';
    opts.batchNormalization = false ;
    opts.networkType = 'simplenn' ;
    opts.strides = [2, 2, 1, 2] ;
    opts.cudnnWorkspaceLimit = 1024*1024*1024 ; % 1GB
    opts = vl_argparse(opts, varargin) ;

    if numel(opts.exemplarSize) == 1
        opts.exemplarSize = [opts.exemplarSize, opts.exemplarSize];
    end
    if numel(opts.instanceSize) == 1
        opts.instanceSize = [opts.instanceSize, opts.instanceSize];
    end

    net = modified_alexnet(struct(), opts) ;

    % Meta parameters
    net.meta.normalization.interpolation = 'bicubic' ;
    net.meta.normalization.averageImage = [] ;
    net.meta.normalization.keepAspect = true ;
    net.meta.augmentation.rgbVariance = zeros(0,3) ;
    net.meta.augmentation.transformation = 'stretch' ;

    % Fill in default values
    net = vl_simplenn_tidy(net) ;

    % Switch to DagNN if requested
    switch lower(opts.networkType)
      case 'simplenn'
        % done
      case 'dagnn'
        net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;
        net.addLayer('top1err', dagnn.Loss('loss', 'classerror'), ...
                     {'prediction','label'}, 'top1err') ;
        net.addLayer('top5err', dagnn.Loss('loss', 'topkerror', ...
                                           'opts', {'topK',5}), ...
                     {'prediction','label'}, 'top5err') ;
      otherwise
        assert(false) ;
    end

end

% --------------------------------------------------------------------
function net = add_block(net, opts, id, h, w, in, out, stride, pad, init_bias)
% --------------------------------------------------------------------
info = vl_simplenn_display(net) ;
fc = (h == info.dataSize(1,end) && w == info.dataSize(2,end)) ;
if fc
  name = 'fc' ;
else
  name = 'conv' ;
end
convOpts = {'CudnnWorkspaceLimit', opts.cudnnWorkspaceLimit} ;
net.layers{end+1} = struct('type', 'conv', 'name', sprintf('%s%s', name, id), ...
                           'weights', {{init_weight(opts, h, w, in, out, 'single'), zeros(out, 1, 'single')}}, ...
                           'stride', stride, ...
                           'pad', pad, ...
                           'learningRate', [1 2], ...
                           'weightDecay', [opts.weightDecay 0], ...
                           'opts', {convOpts}) ;
if opts.batchNormalization
  net.layers{end+1} = struct('type', 'bnorm', 'name', sprintf('bn%s',id), ...
                             'weights', {{ones(out, 1, 'single'), zeros(out, 1, 'single'), zeros(out, 2, 'single')}}, ...
                             'learningRate', [2 1 0.05], ...
                             'weightDecay', [0 0]) ;
end
net.layers{end+1} = struct('type', 'relu', 'name', sprintf('relu%s',id)) ;
end

% --------------------------------------------------------------------
function net = add_block_conv_only(net, opts, id, h, w, in, out, stride, pad, init_bias)
% --------------------------------------------------------------------
info = vl_simplenn_display(net) ;
fc = (h == info.dataSize(1,end) && w == info.dataSize(2,end)) ;
if fc
  name = 'fc' ;
else
  name = 'conv' ;
end
convOpts = {'CudnnWorkspaceLimit', opts.cudnnWorkspaceLimit} ;
net.layers{end+1} = struct('type', 'conv', 'name', sprintf('%s%s', name, id), ...
                           'weights', {{init_weight(opts, h, w, in, out, 'single'), zeros(out, 1, 'single')}}, ...
                           'stride', stride, ...
                           'pad', pad, ...
                           'learningRate', [1 2], ...
                           'weightDecay', [opts.weightDecay 0], ...
                           'opts', {convOpts}) ;
end

% --------------------------------------------------------------------
function net = add_norm(net, opts, id)
% --------------------------------------------------------------------
if ~opts.batchNormalization
  net.layers{end+1} = struct('type', 'normalize', ...
                             'name', sprintf('norm%s', id), ...
                             'param', [5 1 0.0001/5 0.75]) ;
end
end

% --------------------------------------------------------------------
function net = add_dropout(net, opts, id)
% --------------------------------------------------------------------
if ~opts.batchNormalization
  net.layers{end+1} = struct('type', 'dropout', ...
                             'name', sprintf('dropout%s', id), ...
                             'rate', 0.5) ;
end
end

% --------------------------------------------------------------------
function net = modified_alexnet(net, opts)
% --------------------------------------------------------------------

    strides = ones(1, 7);
    strides(1:numel(opts.strides)) = opts.strides(:);

    net.layers = {} ;

    net = add_block(net, opts, '1', 11, 11, 3, 96, strides(1), 0) ;
    net = add_norm(net, opts, '1') ;
    net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                               'method', 'max', ...
                               'pool', [3 3], ...
                               'stride', strides(2), ...
                               'pad', 0) ;

    net = add_block(net, opts, '2', 5, 5, 48, 256, strides(3), 0) ;
    net = add_norm(net, opts, '2') ;
    net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                               'method', 'max', ...
                               'pool', [3 3], ...
                               'stride', strides(4), ...
                               'pad', 0) ;

    net = add_block(net, opts, '3', 3, 3, 256, 384, strides(5), 0) ;
    net = add_block(net, opts, '4', 3, 3, 192, 384, strides(6), 0) ;
    net = add_block_conv_only(net, opts, '5', 3, 3, 192, 256, strides(7), 0) ;

    % Check if the receptive field covers full image

    [ideal_exemplar, ~] = ideal_size(net, opts.exemplarSize);
    [ideal_instance, ~] = ideal_size(net, opts.instanceSize);
    assert(sum(opts.exemplarSize==ideal_exemplar)==2, 'exemplarSize is not ideal.');
    assert(sum(opts.instanceSize==ideal_instance)==2, 'instanceSize is not ideal.');
end

author: bertinetto

repo: siamese-fc


text:
function fx = p54_fun ( n, x )

%*****************************************************************************80
%
%% P54_FUN evaluates the integrand for problem 54.
%
%  Discussion:
%
%    The reference claims that this integrand is more closely approximated
%    by the trapezoid rule than by Gauss-Legendre quadrature.
%
%    Points  Trapezoid  Gauss-Legendre
%     4      1.91667    2.53883
%    12      2.1594     2.25809
%
%    However, the stated results hardly give one confindence in
%    the convergence of the trapezoid results, and I am unable to
%    confirm them, because my results for 4 points give good results
%    (about 1.14) for BOTH Trapezoid and Gauss-Legendre%
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    2 / ( 2 + sin ( 10 * PI * x ) )
%
%  Exact Integral:
%
%    2 / sqrt ( 3 ) = 1.1547...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Prem Kythe, Pratap Puri,
%    Computational Methods for Linear Integral Equations,
%    Birkhaeuser, 2002.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx(1:n) = 2.0 ./ ( 2.0 + sin ( 10.0 * pi * x(1:n) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,y]=selectPolygon()
% function [x,y]=selectPolygon()
% OUTPUT: x,y - xy coordinates of vertices of selected polygon in figure
% NOTE: The coord vectors are set up so that the last and first entries
% are the same
% PURPOSE:
% Function to select a polygonal region in the current figure
% and return the coordinates of the vertices (x,y)
% USAGE:
% - Left click to select point
% - ESC to undo
% - Right click to end and finish loop

% Author: Ranjith Unnikrishnan
%* Carnegie Mellon University, Vision and Mobile Robots Laboratory       *
%* THE MATERIAL EMBODIED IN THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"     *
%* AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,      *
%* INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR      *
%* FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL CARNEGIE MELLON  *
%* UNIVERSITY BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,            *
%* SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY         *
%* KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,        *
%* LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF     *
%* THIRD PARTIES, WHETHER OR NOT CARNEGIE MELLON UNIVERSITY HAS BEEN     *
%* ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON        *
%* ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE     *
%* POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.                      *
%

x=[];
y=[];
done=0;
n_selected=0;
crosses_h=1234.567; % some arbit number
first=1;
xlabel('left mouse=select, ESC=undo last selection, right mouse=finish');

while (done==0)
    [xc,yc,button]=ginput(1);

    switch (button)
        case 1     % left button
            % Select point
            fprintf(1,'Selected point (%.3f,%.3f)\n',xc,yc);
            x=[x;xc];
            y=[y;yc];
            n_selected=n_selected+1;
        case 27      % ESC key
            % Undo last point
            if(n_selected > 0)
                fprintf(1,'Deleted last point\n');
                x(end)=[];
                y(end)=[];
                n_selected=n_selected-1;
            end
        case 3      % Right mouse button
            fprintf(1,'Ending\n');
            done=1;
            if(length(x)>0)
                % Terminate
                x=[x;x(1)];
                y=[y;y(1)];
            end
        otherwise
    end

    if(ishandle(crosses_h))
        %fprintf(1,'crosses_h is a handle\n');
        delete(crosses_h);
    end

    hold on;
    if(length(x)>0)
        crosses_h=plot([x;x(1)],[y;y(1)],'wo-',...
            'markersize',6,'linewidth',2);
    end
    hold off;
end

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
% STD_CHECKSET - check STUDY set consistency
%
% Usage: >> [STUDY, ALLEEG] = std_checkset(STUDY, ALLEEG);
%
% Input:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of EEG datasets included in the STUDY structure
%
% Output:
%   STUDY      - a new STUDY set containing some or all of the datasets in ALLEEG,
%                plus additional information from the optional inputs above.
%   ALLEEG     - an EEGLAB vector of EEG sets included in the STUDY structure
%
% Authors:  Arnaud Delorme & Hilit Serby, SCCN, INC, UCSD, November 2005

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, arno@sccn.ucsd.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 [STUDY, ALLEEG, command] = std_checkset(STUDY, ALLEEG, options)

if nargin < 2
    help std_checkset;
    return;
end
command  = '';

if isempty(STUDY), return; end
studywasempty = 0;

modif = 0;
if ~isfield(STUDY, 'name'),  STUDY.name  = ''; modif = 1; end
if ~isfield(STUDY, 'task'),  STUDY.task  = ''; modif = 1; end
if ~isfield(STUDY, 'notes'), STUDY.notes = ''; modif = 1; end
if ~isfield(STUDY, 'filename'),  STUDY.filename  = ''; modif = 1; end
if ~isfield(STUDY, 'filepath'),  STUDY.filepath  = ''; modif = 1; end
if ~isfield(STUDY, 'history'),   STUDY.history   = ''; modif = 1; end
if ~isfield(STUDY, 'subject'),   STUDY.subject   = {}; modif = 1; end
if ~isfield(STUDY, 'group'),     STUDY.group     = {}; modif = 1; end
if ~isfield(STUDY, 'session'),   STUDY.session   = {}; modif = 1; end
if ~isfield(STUDY, 'run'),       STUDY.run       = {}; modif = 1; end
if ~isfield(STUDY, 'condition'), STUDY.condition = {}; modif = 1; end
if ~isfield(STUDY, 'etc'),       STUDY.etc       = []; modif = 1; end
if ~isfield(STUDY, 'cache'),     STUDY.cache     = []; modif = 1; end
if ~isfield(STUDY, 'etc.warnmemory'), STUDY.etc.warnmemory = 1; modif = 1; end
if ~isfield(STUDY, 'preclust'),  STUDY.preclust  = []; modif = 1; end
if ~isfield(STUDY, 'datasetinfo'), STUDY.datasetinfo = []; modif = 1; end
if ~isfield(STUDY.etc, 'version'), STUDY.etc.version = []; modif = 1; end
if ~isfield(STUDY.preclust, 'erpclusttimes' ),  STUDY.preclust.erpclusttimes = []; modif = 1; end
if ~isfield(STUDY.preclust, 'specclustfreqs' ), STUDY.preclust.specclustfreqs = []; modif = 1; end
if ~isfield(STUDY.preclust, 'erspclustfreqs' ), STUDY.preclust.erspclustfreqs = []; modif = 1; end
if ~isfield(STUDY.preclust, 'erspclusttimes' ), STUDY.preclust.erspclusttimes = []; modif = 1; end
if ~isfield(STUDY.datasetinfo, 'comps') && ~isempty(STUDY.datasetinfo), STUDY.datasetinfo(1).comps = []; modif = 1; end
if ~isfield(STUDY.datasetinfo, 'index') && ~isempty(STUDY.datasetinfo), STUDY.datasetinfo(1).index = []; modif = 1; end
if ~isfield(STUDY.datasetinfo, 'run') && ~isempty(STUDY.datasetinfo), STUDY.datasetinfo(1).run = []; modif = 1; end

% all summary fields
% ------------------
try, subject = unique_bc({ STUDY.datasetinfo.subject });
catch,
    subject = '';
    disp('Important warning: some datasets do not have subject codes; some functions may crash!');
end
try, group = unique_bc({ STUDY.datasetinfo.group });
catch,
    group = {};
    % disp('Important warning: some datasets do not have group codes; some functions may crash!');
end
try, condition = unique_bc({ STUDY.datasetinfo.condition });
catch,
    condition = {};
    disp('Important warning: some datasets do not have condition codes; some functions may crash!');
end
try, session = unique_bc([STUDY.datasetinfo.session]);
catch,
    session = '';
    % disp('Important warning: some datasets do not have session numbers; some functions may crash!');
end
try, run = unique_bc([STUDY.datasetinfo.run]);
catch,
    run = '';
    % disp('Important warning: some datasets do not have session numbers; some functions may crash!');
end
if ~isequal(STUDY.subject,   subject  ), STUDY.subject   = subject;   modif = 1; end
if ~isequal(STUDY.group,     group    ), STUDY.group     = group;     modif = 1; end
if ~isequal(STUDY.condition, condition), STUDY.condition = condition; modif = 1; end
if ~isequal(STUDY.session,   session  ), STUDY.session   = session;   modif = 1; end
if ~isequal(STUDY.run,       run      ), STUDY.run       = run;       modif = 1; end

% check dataset info consistency
% ------------------------------
for k = 1:length(STUDY.datasetinfo)
    if ~strcmpi(STUDY.datasetinfo(k).filename, ALLEEG(k).filename)
        STUDY.datasetinfo(k).filename = ALLEEG(k).filename; modif = 1;
        fprintf('Warning: file name has changed for dataset %d and the study has been updated\n', k);
        fprintf('         to discard this change in the study, reload it from disk\n');
    end
end

% recompute setind array (setind is deprecated but we keep it anyway)
% -------------------------------------------------------------------
setind  = [];
sameica = std_findsameica(ALLEEG);
for index = 1:length(sameica)
    setind(length(sameica{index}):-1:1,index) = sameica{index}';
end
setind(find(setind == 0)) = NaN;
if any(isnan(setind))
    warndlg('Warning: non-uniform set of dataset, some function might not work');
end

% check that dipfit is present in all datasets
% --------------------------------------------
for cc = 1:length(sameica)
    idat = [];
    for tmpi = 1:length(sameica{cc})
        if isfield(ALLEEG(sameica{cc}(tmpi)).dipfit, 'model')
            idat = sameica{cc}(tmpi);
        end
    end
    if ~isempty(idat)
        for tmpi = 1:length(sameica{cc})
            if ~isfield(ALLEEG(sameica{cc}(tmpi)).dipfit, 'model')
                ALLEEG(sameica{cc}(tmpi)).dipfit = ALLEEG(idat).dipfit;
                ALLEEG(sameica{cc}(tmpi)).saved  = 'no';
                fprintf('Warning: no ICA dipoles for dataset %d, using dipoles from dataset %d (same ICA)\n', sameica{cc}(tmpi), idat);
            end
        end
    end
end

% put in fake channels if channel labels are missing
% --------------------------------------------------
chanlabels = { ALLEEG.chanlocs };
if any(cellfun(@isempty, chanlabels))
    if any(~cellfun(@isempty, chanlabels))
        disp('********************************************************************');
        disp(' IMPORTANT WARNING: SOME DATASETS DO NOT HAVE CHANNEL LABELS AND ');
        disp(' SOME OTHERs HAVE CHANNEL LABELS. GENERATING CHANNEL LABELS FOR ');
        disp(' THE FORMER DATASETS (THIS SHOULD PROBABLY BE FIXED BY THE USER).');
        disp('********************************************************************');
    end
    disp('Generating channel labels for all datasets...');
    for currentind = 1:length(ALLEEG)
        for ind = 1:ALLEEG(currentind).nbchan
            ALLEEG(currentind).chanlocs(ind).labels = int2str(ind);
        end
    end
    ALLEEG(currentind).saved = 'no';
end

if length( unique( [ ALLEEG.srate ] )) > 1
    disp('********************************************************************');
    disp(' IMPORTANT WARNING: SOME DATASETS DO NOT HAVE THE SAME SAMPLING ');
    disp(' RATE AND THIS WILL MAKE MOST OF THE STUDY FUNCTIONS CRASH. THIS');
    disp(' SHOULD PROBABLY BE FIXED BY THE USER.');
    disp('********************************************************************');
end

% create STUDY design if it is not present
% ----------------------------------------
rebuild_design = 0;
if ~isfield(STUDY, 'cluster'), STUDY.cluster = []; modif = 1; end
if ~isfield(STUDY, 'changrp'), STUDY.changrp = []; modif = 1; end
if ~studywasempty
    if isfield(STUDY.datasetinfo, 'trialinfo')
        alltrialinfo = { STUDY.datasetinfo.trialinfo };
        if any(cellfun(@isempty, alltrialinfo)) && any(~cellfun(@isempty, alltrialinfo))
            disp('Rebuilding trial information structure for STUDY');
            STUDY  = std_maketrialinfo(STUDY, ALLEEG); % some dataset do not have trialinfo and
            % some other have it, remake it for everybody
        end
    elseif ALLEEG(1).trials > 1
        disp('Rebuilding trial information structure for STUDY');
        STUDY  = std_maketrialinfo(STUDY, ALLEEG); % some dataset do not have trialinfo and
    end
    if ~isfield(STUDY, 'design') || isempty(STUDY.design) || ~isfield(STUDY.design, 'name')
        STUDY  = std_maketrialinfo(STUDY, ALLEEG);
        STUDY  = std_makedesign(STUDY, ALLEEG);
        STUDY  = std_selectdesign(STUDY, ALLEEG,1);
        rebuild_design = 0;
    else
        if isfield(STUDY.design, 'indvar1')
            STUDY  = std_convertdesign(STUDY, ALLEEG);
        end

        % check independent variable
        setFields   = fieldnames(STUDY.datasetinfo);
        if isempty(ALLEEG(1).event)
             trialFields = [];
        else trialFields = fieldnames(ALLEEG(1).event);
        end
        rmInd       = [];
        for indDes = 1:length(STUDY.design)
            for indVar = 1:length(STUDY.design(indDes).variable)
                if isempty(strmatch( STUDY.design(indDes).variable(indVar).label, setFields)) && ...
                    ~isempty(trialFields) && ~isempty(STUDY.design(indDes).variable(indVar).label) && ...
                        isempty(strmatch( STUDY.design(indDes).variable(indVar).label, trialFields))
                        % Missing independent variable

                        if isempty(rmInd)
                            textVal1 = sprintf('The variable ''%s'' defined in the STUDY design is not found in ', STUDY.design(indDes).variable(indVar).label);
                            textVal2 = 'some of the EEG datasets. Some STUDY features will not be functional.';
                            if nargin > 2 && strcmpi(options, 'popup')
                                textVal = [ textVal1 10 textVal2 10 10 'Do you want to remove all designs with missing variables?' ];
                                res = questdlg2(  textVal, 'Issue with studies', 'No', 'Yes', 'Yes');
                                if strcmpi(res, 'Yes')
                                    rmInd = [ rmInd indDes ];
                                end
                            else
                                disp([textVal1 textVal2]);
                            end
                        else
                            rmInd = [ rmInd indDes ];
                        end
                end
            end
        end
        if ~isempty(rmInd)
            STUDY.design(rmInd) = [];
        end

        % convert combined independent variable values
        % between dash to cell array of strings
        % -------------------------------------
        for inddes = 1:length(STUDY.design)
            rmVar = [];
            for indvar = 1:length(STUDY.design(inddes).variable)

                % add pairing info in case it is missing
                if isempty(STUDY.design(inddes).variable(indvar).label)
                    rmVar = [ rmVar indvar ];
                end
                if ~isfield(STUDY.design(inddes).variable, 'pairing') || isempty(STUDY.design(inddes).variable(indvar).pairing)
                    if strcmpi(STUDY.design(inddes).variable(indvar).label, 'group')
                        STUDY.design(inddes).variable(indvar).pairing = 'off';
                    else
                        STUDY.design(inddes).variable(indvar).pairing = 'on';
                    end
                end

                for indval = 1:length(STUDY.design(inddes).variable(indvar).value)
                    STUDY.design(inddes).variable(indvar).value{indval} = convertindvarval(STUDY.design(inddes).variable(indvar).value{indval});
                end
            end
            STUDY.design(inddes).variable(rmVar) = [];
            for indinclude = 1:length(STUDY.design(inddes).include)
                if iscell(STUDY.design(inddes).include{indinclude})
                    for indval = 1:length(STUDY.design(inddes).include{indinclude})
                        STUDY.design(inddes).include{indinclude}{indval} = convertindvarval(STUDY.design(inddes).include{indinclude}{indval});
                    end
                end
            end
        end
    end

    if rebuild_design % in case datasets have been added or removed
        STUDY = std_rebuilddesign(STUDY, ALLEEG);
    end

    % scan design to fix old pairing format
    % ------------------------------------
    for design = 1:length(STUDY.design)
        for var = 1:length(STUDY.design(design).variable)
            if ischar(STUDY.design(design).variable(var).pairing)
                if strcmpi(STUDY.design(design).variable(var).pairing, 'paired')
                    STUDY.design(design).variable(var).pairing = 'on';
                elseif strcmpi(STUDY.design(design).variable(var).pairing, 'unpaired')
                    STUDY.design(design).variable(var).pairing = 'off';
                end
            end
        end
    end

    % add filepath field if absent
    for ind = 1:length(STUDY.design)
        if ~isfield(STUDY.design, 'filepath') || (isnumeric(STUDY.design(ind).filepath) && isempty(STUDY.design(ind).filepath))
            STUDY.design(ind).filepath = '';
            STUDY.saved = 'no';
            modif = 1;
        end
    end

    % check for level field and add it if it not present
    if length(STUDY.design) > 0 && isfield(STUDY.design, 'variable')
        if ~isfield(STUDY.design(1).variable, 'level')
            STUDY = std_addvarlevel(STUDY);
        end
    end

    % check that ICA is present and if it is update STUDY.datasetinfo
    recreateclusters = 0;
    for index = 1:length(STUDY.datasetinfo)
        if max(STUDY.datasetinfo(index).comps) > size(ALLEEG(index).icaweights,1)
            STUDY.datasetinfo(index).comps = [1:size(ALLEEG(index).icaweights,1)];
            recreateclusters = 1;
            modif = 1;
        end
    end
    if recreateclusters
        fprintf('Recreating empty clusters because ICA has changed for some datasets\n');
        STUDY.cluster = [];
    end

    % make channel groups
    % -------------------
    STUDYTMP = std_changroup(STUDY, ALLEEG);
    if length(STUDYTMP.changrp) ~= length(STUDY.changrp)
        STUDY.changrp = STUDYTMP.changrp;
        modif = 1;
    end
    clear STUDYTMP;
end

% check cluster array
% -------------------
if isempty(STUDY.changrp) && isempty(STUDY.cluster)
    rebuild_design = 1;
end
if isfield(STUDY.cluster, 'sets'),
    if max(STUDY.cluster(1).sets(:)) > length(STUDY.datasetinfo)
        disp('Warning: Some datasets had been removed from the STUDY, clusters have been reinitialized');
        STUDY.cluster = []; modif = 1;
    else
    	if isempty(STUDY.cluster(1).sets)
            STUDY.cluster = []; modif = 1;
        end
    end
end
if ~studywasempty
    if isempty(STUDY.cluster)
        modif = 1;
        [STUDY] = std_createclust(STUDY, ALLEEG, 'parentcluster', 'on');
    end
    if length(STUDY.cluster(1).child) == length(STUDY.cluster)-1 && length(STUDY.cluster) > 1
        newchild = { STUDY.cluster(2:end).name };
        if ~isequal(STUDY.cluster(1).child, newchild)
            STUDY.cluster(1).child = newchild;
        end
    end
end

% adapt STUDY design to new framework
% -----------------------------------
if ~isempty(STUDY.design) && ~isfield(STUDY.design(1).variable, 'vartype')
    for inddes = 1:length(STUDY.design)
        for indvar = 1:length(STUDY.design(inddes).variable)
            STUDY.design(inddes).variable(indvar).vartype = 'categorical';
        end
    end
end
if ~isempty(STUDY.design) && isfield(STUDY.design(1), 'cell')
    STUDY.design = rmfield(STUDY.design, 'cell');
end

% determine if there has been any change
% --------------------------------------
if modif
    STUDY.saved = 'no';
    command = '[STUDY ALLEEG] = std_checkset(STUDY, ALLEEG);';
    addToHistory = true;
    % check duplicate
    if length(STUDY.history) >= length(command) && strcmpi(STUDY.history(end-length(command)+1:end), command)
        addToHistory = false;
    end
    if addToHistory
        STUDY.history =  sprintf('%s\n%s',  STUDY.history, command);
    end
end

% convert combined independent variables
% --------------------------------------
function val = convertindvarval(val)
if ischar(val)
    inddash = findstr(' - ', val);
    if isempty(inddash), return; end
    inddash = [ -2 inddash length(val)+1];
    for ind = 1:length(inddash)-1
        newval{ind} = val(inddash(ind)+3:inddash(ind+1)-1);
    end
    val = newval;
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% FILTER_FREQ Calculate center frequencies and bandwidths
%
% Usage
%    [psi_xi, psi_bw, phi_bw] = FILTER_FREQ(filt_opt)
%
% Input
%    filt_opt (struct): The parameters defining the filter bank.
%
% Output
%    psi_xi (numeric): The center frequencies of the wavelet filters.
%    psi_bw (numeric): The bandwidths of the wavelet filters.
%    phi_bw (numeric): The bandwidth of the lowpass filter.
%
% Description
%    Called by WAVELET_1D and WAVELET_LAYER_1D, this function provides the
%    center frequency and bandwidth of a mother wavelet, whose parameters
%    are specified in filt_opt. It operates as a mere disjunction between
%    SPLINE_FREQ_1D, and MORLET_FREQ_1D.
%
% See also
%   MORLET_FREQ_1D, SPLINE_FREQ_1D, WAVELET_1D, WAVELET_LAYER_1D

function [psi_xi, psi_bw, phi_bw] = filter_freq(filter_options)
	if strcmp(filter_options.filter_type,'spline_1d') || ...
        strcmp(filter_options.filter_type,'selesnick_1d')
		[psi_xi,psi_bw,phi_bw] = dyadic_freq_1d(filter_options);
	elseif strcmp(filter_options.filter_type,'morlet_1d') || ...
		strcmp(filter_options.filter_type,'gabor_1d')
		[psi_xi,psi_bw,phi_bw] = morlet_freq_1d(filter_options);
	else
		error('Unknown filter type ''%s''', filter_options.filter_type);
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function v=polyValNewton(z,a,c)
%POLYVALNEWTON Evaluate a set of polynomial given in Newton's form at a
%              desired set of points. This is useful when given a set of
%              Hermite interpolating polynomials for interpolating a
%              vector.
%
%INPUTS: z A numPointsX1 or 1XnumPoints vector of scalar points where the
%          values of the numDim polynomials are desired.
%        a A numDimXn matrix of polynomial coefficients for numDim
%          polynomials, as described below.
%        c A numDim X(n-1) matrix of the control points associated with
%          each of the numDim polynomial in Newton's form, which might be
%          the control points used in the HermiteInterpPoly function.
%
%OUTPUTS: v The numDimXnumPoints matrix of values of the polynomials
%           evaluated at the points in z.
%
%A polynomial function in Newton's form evaluated at point z has the form
%y(z)=a(1)+sum_{k=1}^{n-1}a(k+1)(z-c(1))*(z-c(2))*...*(z-c(k))
%This function just evaluates multiple scalar polynomials.
%
%The algorithm is based on  VALUE from Chapter 19 of [1] to handle vector
%values; a loop is used here.
%
%Newton's polynomial form arises when dealing with Hermite interpolating
%polynomials as is discussed in Section 2.1.3 of [2] and in Chapters 3.3
%and 3.4 of [3].
%
%REFERENCES:
%[1] A. Nijenhuis and H. S. Wilf, Combinatorial Algorithms for Computers
%    and Calculators, 2nd ed. New York: Academic press, 1978.
%[2] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 3rd ed.
%    New York: Springer, 2002.
%[3] R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Boston, MA:
%    Brooks/Cole, 2011.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(a,1);
numPoints=length(z);

v=zeros(numDim,numPoints);
for curDim=1:numDim
   v(curDim,:)=polyValNewtonScalar(z,a(curDim,:),c(curDim,:));
end
end

function v=polyValNewtonScalar(z,a,c)
%Below is VALUE from Chapter 19 of [1]
%
%REFERENCES:
%A. Nijenhuis and H. S. Wilf, Combinatorial Algorithms for Computers
%and Calculators, 2nd ed. New York: Academic press, 1978.

n=length(a);

v=ones(size(z))*a(n);
for k=(n-1):-1:1
    v=(z-c(k)).*v+a(k);
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 Q = demo_rvbpcamv
warning('This demo is deprecated. See folder fa instead')

n = 50;
m = 2;
d = 1;

randn('state', 6);
rand('state', 6);

mu = [0;0];%20 * randn(m,1);
W = [1; 0.5];%orth(randn(m,d)) * diag(10:-1:(10-d+1));
X = orth(randn(d,n)')' * sqrt(n);
s2 = 0.1;

% Generate normal data
Y = W*X + repmat(mu,1,n);

% Add noise
Yn = Y + sqrt(s2)*randn(m,n);

% Generate some outliers
Yno = Yn;
p = (rand(m,n) < 0.05);
%p = zeros(m,n) == 1;
%p(1,1) = true;
%p(2,2) = true;
%p(1,3) = true;
%p(:,4) = true;
Yno(p) = 3 + 2 * rand(sum(p(:)),1);

% Generate some missing values
Ynom = Yno;
pmv = (rand(m,n) < 0.0);
Ynom(pmv) = NaN;

% Use STANDARD PCA for non-outlier data
Q = vbpcamv(Yn, d);
W_pca = Q.W
X_pca = Q.X
mu_pca = Q.mu
s2_pca = 1/Q.tau
%[W_pca, X_pca, mu_pca, s2_pca] = pca_full(Yn, d);
Y_pca = W_pca*X_pca + repmat(mu_pca, 1, n);

% Run different algorithms
Q = vbpcamv(Ynom, d);
W_p = Q.W
X_p = Q.X
mu_p = Q.mu
s2_p = 1/Q.tau
%[W_p, X_p, mu_p, s2_p] = pca_full(Ynom, d);
[W_rp,X_rp,Sv_rp,mu_rp,nu_rp,s2_rp,U_rp] = rppcamv(Ynom, d, 'maxiters', 100);
init.tau = 1e2;
init.nu = 1;
init.mu = 0;%mean(Yno,2)
prior = [];
% $$$ prior.atau = 1e3;
% $$$ prior.atau = 1e0
prior.aw = 1e-10;
prior.bw = 1e-3;
results_rvb = rvbpcamv(Ynom, d, 'maxiters', 200, 'init', init, 'prior', ...
                       prior, 'startupdatehyper', 10, 'startrotate', 1, ...
                       'rotate', true);
W_rvb = results_rvb.W;
X_rvb = results_rvb.X;
S_rvb = results_rvb.Sv;
mu_rvb = results_rvb.mu;
nu_rvb = results_rvb.nu;

% Reconstruct
Y_p = W_p*X_p + repmat(mu_p,1,n);
Y_rp = W_rp*X_rp + repmat(mu_rp,1,n);
Y_rvb = W_rvb*X_rvb + repmat(mu_rvb,1,n);

% Error of the principal subspace
err_W_p = 180 * subspace(W_p, W_pca) / pi
err_W_rp = 180 * subspace(W_rp, W_pca) / pi
err_W_rvb = 180 * subspace(W_rvb, W_pca) / pi

plot_results(Ynom, []);
plot_results(Ynom, Y_p)
plot_results(Ynom, Y_rp)
plot_results(Ynom, Y_rvb)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
function plot_results(Y, Yh)
%subspace2d(Yh, [], Y)
figure
plot(Y(1,:), Y(2,:), 'x', 'markersize', 5, 'color', [.2 .2 .2]);
if ~isempty(Yh)
  hold on
  plot(Yh(1,:), Yh(2,:), 'ko', 'markersize', 3);
  nans = nan*ones(1,size(Y,2));
  plot([Y(1,:);Yh(1,:);nans], [Y(2,:);Yh(2,:);nans], ':', 'color', [.0 .0 .0]);
end
set(gca, 'DataAspectRatioMode', 'manual');
set(gca, 'DataAspectRatio', [1 1 1]);
set(gca, 'XTick', [], 'YTick', []);
set(gcf, 'units', 'centimeters');
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1), pos(2), 5 5])
set(gcf, 'PaperPositionMode', 'auto', 'paperunits', 'centimeters', 'PaperSize', [5 5]);

xmg = 0.05 * (max(Y(1,:)) - min(Y(1,:)));
ymg = 0.05 * (max(Y(2,:)) - min(Y(2,:)));
xl = [min(Y(1,:))-xmg, max(Y(1,:))+xmg];
yl = [min(Y(2,:))-ymg, max(Y(2,:))+ymg];
set(gca, 'xlim', xl, 'ylim', yl);

% $$$ % Mark outliers
% $$$ hold on
% $$$ indeces = sum(p,1)>0;
% $$$ plot(VYno(1,indeces), VYno(2, indeces), 'go', 'MarkerSize', 8);

% $$$ % Mark observations with missing values
% $$$ indeces = sum(pmv,1)>0;
% $$$ scatter(VYno(1,indeces), VYno(2, indeces), 'yo');

return

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function pass = test_isempty( )

% Example 1:
v = ballfunv();
pass(1) = isempty(v);

% Example 2:
f = ballfun(1);
g = ballfun();
v = ballfunv(f,g,f);
pass(2) = ~isempty(v);

% Example 3:
f = ballfun(1);
v = ballfunv(f,f,f);
pass(3) = ~isempty(v);

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef SegmentationLossSemiSupervised < dagnn.Loss
    % SegmentationLossSemiSupervised
    %
    % Similar to dagnn.SegmentationLoss, but also sets pixel weights.
    %
    % Inputs: prediction, labels, labelsImage, [classWeights],
    % isWeaklySupervised, [masksThingsCell]
    % Outputs: loss
    %
    % Note:
    %    - classWeights can be empty, which means they are ignored.
    %    - labels or labelsImage can be empty.
    %
    % Copyright by Holger Caesar, 2016

    properties
        layerFS
        layerWS
    end

    methods
        function outputs = forward(obj, inputs, params) %#ok<INUSD>
            assert(numel(inputs) == 6);
            scoresMap = inputs{1};
            labels = inputs{2};
            labelsImage = inputs{3};
            classWeights = inputs{4};
            isWeaklySupervised = inputs{5};
            masksThingsCell = inputs{6};
            assert(~isempty(isWeaklySupervised));

            if isWeaklySupervised
                outputs = obj.layerWS.forward({scoresMap, labelsImage, classWeights, masksThingsCell}, {});
            else
                outputs = obj.layerFS.forward({scoresMap, labels, classWeights}, {});
            end

            % Combine loss statistics
            imageCount = size(scoresMap, 4);
            n = obj.numAveraged;
            m = n + imageCount;
            obj.average = (n * obj.average + double(gather(outputs{1}))) / m;
            obj.numAveraged = m;
        end

        function [derInputs, derParams] = backward(obj, inputs, ~, derOutputs)
            scoresMap = inputs{1};
            labels = inputs{2};
            labelsImage = inputs{3};
            classWeights = inputs{4};
            isWeaklySupervised = inputs{5};

            if isWeaklySupervised
                derInputs = obj.layerWS.backward({scoresMap, labelsImage, classWeights}, {}, derOutputs);
            else
                derInputs = obj.layerFS.backward({scoresMap, labels, classWeights}, {}, derOutputs);
            end
            derInputs{5} = [];
            derInputs{6} = [];
            derParams = {};
        end

        function obj = SegmentationLossSemiSupervised(varargin)
            obj.load(varargin);
        end

        function forwardAdvanced(obj, layer)
            % Modification: Overrides standard forward pass to avoid giving up when any of
            % the inputs is empty.

            in = layer.inputIndexes;
            out = layer.outputIndexes;
            par = layer.paramIndexes;
            net = obj.net;
            inputs = {net.vars(in).value};

            % clear inputs if not needed anymore
            for v = in
                net.numPendingVarRefs(v) = net.numPendingVarRefs(v) - 1;
                if net.numPendingVarRefs(v) == 0
                    if ~net.vars(v).precious && ~net.computingDerivative && net.conserveMemory
                        net.vars(v).value = [];
                    end
                end
            end

            % call the simplified interface
            outputs = obj.forward(inputs, {net.params(par).value});
            for oi = 1:numel(out)
                net.vars(out(oi)).value = outputs{oi};
            end
        end
    end
end

author: nightrome

repo: matconvnet-calvin


text:
function [fNewLat1, fNewLat2, fNewLon1, fNewLon2] = rotate_xsection(fLat1, fLat2, fLon1, fLon2, fAngle)

    % Compute the center of the given cross-section
    fCenterLat = fLat1 - ((fLat1 - fLat2)/2);
    fCenterLon = fLon1 - ((fLon1 - fLon2)/2);

    % Move the center of the cross-section to the origin
    vPos1 = [fLon1-fCenterLon; fLat1-fCenterLat];
    vPos2 = [fLon2-fCenterLon; fLat2-fCenterLat];

    % Compute angle in radians
    fAngle = fAngle*pi/180;

    % Set up the rotation matrix
    mRotate = [cos(fAngle) -sin(fAngle); sin(fAngle) cos(fAngle)];

    % Rotate the cross-section vectors
    vNewPos1 = mRotate * vPos1;
    vNewPos2 = mRotate * vPos2;

    % Move them back to their previous position
    fNewLon1 = vNewPos1(1) + fCenterLon;
    fNewLat1 = vNewPos1(2) + fCenterLat;
    fNewLon2 = vNewPos2(1) + fCenterLon;
    fNewLat2 = vNewPos2(2) + fCenterLat;

author: CelsoReyes


text:
%% housekeeping
close all
clear
clc
%% Information for the title page

% N.B: Be careful with the "_" it has to be escaped with a "\" in order for
% latex to process it as we wish.

rep_data=struct();
myfields={'name','title','address','author','email','abstract'};

rep_data.name='newkeynesswitch'; % Name under which the report will be saved
rep_data.title='Switches in the US Macroeconomic Data: Policy or Volatility?';
rep_data.address='Your institution';
rep_data.author='first\_name last\_name'; % your name
rep_data.email='first\_name@last\_name.com'; % your email
myabstract={['This report investigates switches in the parameters of a ',...
    ' simple new Keynesian model estimated on US data. 4 variants of the model',...
    'are estimated: (i) the first variant allows for switches in volalitily only; ',...
    '(ii) the second model allows for switches in the policy parameters only; ',...
    '(iii) the third specification allows for switches in both policy parameters ',...
    'and the volatility of shocks, but in a synchronized fashion; (iv) the ',...
    'third variant allows for independent switches in both parameters and the ',...
    'volatility of shocks.']
    ' ' % leave a blank space to start the next sentence in a new line
    'We find ample evidence in favor of switching parameters... '
    };
rep_data.abstract=myabstract;
%% load the estimated models and relabel them if necessary
nk_models={'volatilityOnly','policyOnly','volatilityPolicySame','volatilityPolicyIndependent'};
nk_models_newlabels={'volOnly','polOnly','volPolSame','volPolInd'};

model_objects=rise.empty(0);
for imod=1:numel(nk_models)
    tmp=load([nk_models{imod},filesep,'estimation',filesep,'estimated_model']);
    obj=tmp.obj;
    model_objects(imod,1)=obj;
    model_objects(imod,1).legend=nk_models_newlabels{imod};
end
Tao_report(model_objects,rep_data)


text:
function [v, xml] = loadXML(filename)
%
% Reads a XML file and produces a struct array.
%
% filename = name of the XML file
% v = struct variable containing fields of type char or numeric
%
% filename can also be an URL containing XML

% Read file
if ~strcmp(filename(1:5), 'http:');
    [fid, message] = fopen(filename,'r');
    if fid == -1; error(message); end
    xml = fread(fid, 'uint8=>char');
    fclose(fid);
else
    filename = strrep(filename, '\', '/');
    status =0; trying=0;
    while status == 0
        [xml,status] = urlread(filename,'Timeout',5);
        if status == 0
            trying=1;
            disp(sprintf('Warning: Failure reading %s. Trying again...',filename))
            drawnow
        end
    end
    if trying==1
        disp(sprintf('Ok: Success reading %s',filename))
    end
end

% Remove 'new line' characters from the chain
xml = char(xml(:)');
xml = strrep(xml, char(10), '');
xml = strrep(xml, char(13), '');

%v = xml2struct(xml);
v = xmlparse(xml);

author: CSAILVision

repo: LabelMeToolbox


text:
function gauss_seidel_test01 ( )

%*****************************************************************************80
%
%% GAUSS_SEIDEL_TEST01 tests GAUSS_SEIDEL1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSS_SEIDEL_TEST01:\n' );

  it_num = 400;
  n = 20;

  x_exact = ( 1 : n )';
  a = dif2 ( n );
  b = a * x_exact;

  x = zeros ( n, 1 );
  x_plot(1:n,1) = x;

  step = 1 : it_num + 1;
  e = nan ( it_num+1, 1 );
  xm = nan ( it_num+1, 1 );

  e(1,1) = ( norm ( a * x - b ) ).^2;

  for it = 1 : it_num

    x_new = gauss_seidel1 ( n, a, b, x );

    e(it+1,1) = ( norm ( a * x_new - b ) ).^2;
    x_plot(1:n,it+1) = x_new;
%
%  Display the error.
%
    figure ( 1 )
    plot ( step, log ( e ), 'm-*' )
    title ( 'Log (Error^2)' )
    xlabel ( 'Step' )
    ylabel ( 'Error' )
    grid
%
%  Display the motion.
%
    xm(it,1) = sum ( ( x_new(:) - x(:) ).^2 ) / n;

    figure ( 2 )
    plot ( step, log ( xm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Energy' )
    grid
%
%  Update the solution
%
    x = x_new;

  end
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )

  y = ( 0 : it_num );
  for k = 1 : n
    plot ( x_plot(k,1:it_num+1), y )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution.' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );

  return
end


text:
function H = centeringMatrix(dim)

% CENTERINGMATRIX returns the centering matrix for the given dimensionality.
% FORMAT
% DESC returns the centering matrix for the given dimensionality.
% ARG dim : the dimensionality of the centering matrix.
% RETURN H : the centering matrix for the given dimensionality.
%
% COPYRIGHT : Neil D. Lawrence, 2008
%
% SEEALSO : eye

% NDLUTIL

H = -repmat(1./dim, dim, dim) + speye(dim);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
 function mse = mri_exp_errs(arg, Eh, Llist, ti, zmap)
%function mse = mri_exp_errs(arg, Eh, Llist, ti, zmap)
% evaluate errors of MRI exponential approximations for a range of L values
% see mri_exp_approx()
% in
%	arg	cell	to be passed to mri_exp_approx
%	Eh	[M,K]	K=N or K<<N
%	Llist	[L,1]
%	ti	[M,1]	time samples
%	zmap	[N,1]	rate map
% out
%	mse	[L,K]
%
% Copyright 2004-7-5, Jeff Fessler, The University of Michigan

if nargin < 5, ir_usage, end

nL = length(Llist);
mse = zeros(nL, ncol(Eh));
for ll=1:nL
	L = Llist(ll);
	ticker(arg{1}, ll, nL)
	if size(Eh, 2) == length(zmap) % this one pigs memory
		[B C] = mri_exp_approx(ti, zmap, L, arg{:});
	else % histogram-sized C
		[B C] = mri_exp_approx(ti, zmap, L, 'ctest', 1, arg{:});
	end
	mse(ll,:) = mean(abs(Eh - B * C).^2);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function msm_to_mm_test01 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST01 tests MSM_TO_MM_ARRAY_COMPLEX_GENERAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST01\n' );
  fprintf ( 1, '  Convert an MSM to MM array complex general format.\n' );

  output_filename = 'msm_to_mm_test01.mm';
  a = c8mat_indicator ( 5, 3 );
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'array', 'complex', 'general' );

  return
end


text:
function patch2STL(varargin)

% ------------------------------------------------------------------------
% function patch2STL(fileName,V,F,N,solidName)
%
% This function generates an STL file using the vertices (V) and triangular
% faces (F) provided
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2014/04/08
%
% CHANGE LOG:
% 09/12/2009
% * Created
% 2014/04/08
% * Changed %f to %6.6e
% * Added solidName input option
% * Created varargin based input parsing
% ------------------------------------------------------------------------

%% Parse input

%Check if at least 3 inputs and not more than 6 are provided
if nargin<3 || nargin>6
    error('Wrong number of input arguments');
end

%Get first 3 manditory inputs
fileName=varargin{1};
V=varargin{2};
F=varargin{3};

% Check if faces input is consistent with triangles
if ~size(F,2)==3
    error('Patch data does not represent triangles. Convert to triangles first.');
end

%Check N and solidName entries
switch nargin
    case 3 %N, solidName and writeOpt not provided
        N=[];
        solidName='path2STL';
        writeOpt='w';
    case 4 %solidName and writeOpt not provided
        N=varargin{4};
        solidName='path2STL';
        writeOpt='w';
    case 5
        N=varargin{4};
        solidName=varargin{5};
        writeOpt='w';
    case 6
        N=varargin{4};
        solidName=varargin{5};
        writeOpt=varargin{6};
end

%Check if the solidName is a string
if ~ischar(solidName)
    error('solidName needs to be a string');
end

%% Get face normals
if isempty(N) %if face normals are missing
    [N,~]=patchNormal(F,V);
end

%% Write STL file in ASCII format

fid=fopen(fileName,writeOpt); %Append to existing file

fprintf(fid,['solid ',solidName,' \n']);
hw = waitbar(0,'Writing STL-file, please wait...');

%Start loop for all faces, vertices and normals
for qFace = 1:1:size(F,1)

    %Write face normals
    fprintf(fid,'  facet normal %6.6e %6.6e %6.6e\n',N(qFace,:));
    fprintf(fid,'    outer loop\n');
    fprintf(fid,'      vertex %6.6e %6.6e %6.6e\n',V(F(qFace,1),:));
    fprintf(fid,'      vertex %6.6e %6.6e %6.6e\n',V(F(qFace,2),:));
    fprintf(fid,'      vertex %6.6e %6.6e %6.6e\n',V(F(qFace,3),:));
    fprintf(fid,'    endloop\n');
    fprintf(fid,'  endfacet\n');

    %Increment waitbar
    waitbar(qFace/size(F,1),hw,['Writing STL-file ',num2str(round((qFace/size(F,1)).*100)),'%']);
end
fprintf(fid,['endsolid ',solidName,' \n']);
waitbar(1,hw,'Finished! Closing file.');
drawnow;
fclose(fid); %close file ID
close(hw); %close waitbar

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [E,J]=SynthMeasWatsonSHCylSingleRadIsoV_GPD(x, protocol, fibredir, roots)
% Substrate: Impermeable cylinders with one radius in a homogeneous background.
% Orientation distribution: Watson's distribution with SH approximation
% Signal approximation: Gaussian phase distribution.
% Notes: This version includes an isotropic diffusion compartment with its own
% diffusivity.
%
% [E,J]=SynthMeasWatsonSHCylSingleRadIsoV_GPD(x, protocol, fibredir, roots)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the volume fraction of the intracellular space.
% x(2) is the free diffusivity of the material inside and outside the cylinders.
% x(3) is the hindered diffusivity outside the cylinders in perpendicular directions.
% x(4) is the radius of the cylinders.
% x(5) is the concentration parameter of the Watson's distribution.
% x(6) is the volume fraction of the isotropic compartment.
% x(7) is the diffusivity of the isotropic compartment.
%
% protocol is the object containing the acquisition protocol.
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

fiso = x(6);
dIso = x(7);

% Call the model with no isotropic component to get the anisotropic component.
if(nargout == 1)
    Eaniso=SynthMeasWatsonSHCylSingleRadGPD(x, protocol, fibredir, roots);
    Eiso = SynthMeasIsoGPD(dIso, protocol);
else
    [Eaniso,Janiso]=SynthMeasWatsonSHCylSingleRadGPD(x, protocol, fibredir, roots);
    [Eiso, Jiso] = SynthMeasIsoGPD(dIso, protocol);
end

E = (1-fiso)*Eaniso + fiso*Eiso;

if(nargout>1)

    % Update with anisotropic component.
    J = Janiso*(1-fiso);

    % Add derivatives wrt isotropic fraction.
    J(:,6) = Eiso - Eaniso;

    % Add entry for dIso
    J(:,7) = fiso*Jiso;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function x=adodbinsert(cn,sql)
% [x]=oledbinsert(cn,sql)
%
% adodbinsert   Executes the sql statement against the connection cn
%
% Inputs:
%   cn,     open connection to ADO OLEDB ActiveX Data Source Control
%   sql,    SQL statement to be executed
%
% Output
%   x,      return: 1 - successful or 0 - unsuccessful
%
% Notes: Convert cells to strings using char. Convert cells to numeric
% data using cell2mat() for ints or double(cell2mat()) for floats
%
% This code basis on Tim Myers code (oledb*.m) and use ADO OLE DB instead
% of OWC - Office Web Component
%
% Martin Furlan
% martin.furlan@iskra-ae.com
% January 2007

% Execute insert into database
invoke(cn,'BeginTrans');
try
  invoke(cn,'Execute',sql);
  invoke(cn,'CommitTrans');
  sqlSuccess = 1;
catch
  invoke(cn,'RollbackTrans');
  sqlSuccess = 0;
end

x =  sqlSuccess;

author: Sable

repo: mcbench-benchmarks


text:
function L = le(p1,p2)
%PREAL/LE Overloaded LE (<=) operator for class preal.

global useUnitsFlag

if ~(useUnitsFlag) % If physunits is disabled...
    L=le(double(p1),double(p2)); % ... treat as double.
    return
end

p1=preal(p1);
p2=preal(p2);
tol=0.001;
if isscalar(p1)
    p1=preal(p1.value*ones(size(p2)),p1.units);
elseif isscalar(p2)
    p2=preal(p2.value*ones(size(p1)),p2.units);
elseif ndims(p1)==ndims(p2)&size(p1)==size(p2)
    % do nothing
else
    error('Array dimensions must agree')
end

L=false(size(p1));
for k=1:numel(L)
    if any(abs(p1(k).units-p2(k).units)>tol)
        error('Attempt to compare dimensionally inconsistent preals')
    end
    L(k)=le(p1(k).value,p2(k).value);
end


text:
classdef UnfittedMesh_FlatCurve < UnfittedMesh_AbstractBuilder
    properties (GetAccess = public, SetAccess = private)
        unfittedType = 'BOUNDARY';
        maxSubcells = 2;
        nnodesSubcell = 2;

        subcellsMesher      = SubcellsMesher_Boundary_2D;
        cutPointsCalculator = CutPointsCalculator;
        meshPlotter         = MeshPlotter_Boundary_2D;
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function g = p42_g ( n, x )

%*****************************************************************************80
%
%% P42_G evaluates the gradient for problem 42.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  g(1) = 2.0 * ( x(1) - x(2) ) / ( 1.0 + ( x(1) - x(2) )^2 )^2;

  g(2) = 2.0 * ( x(2) - x(1) ) / ( 1.0 + ( x(1) - x(2) )^2 )^2 ...
    - 0.5 * pi * x(3) * cos ( 0.5 * pi * x(2) * x(3) );

  g(3) = - 0.5 * pi * x(2) * cos ( 0.5 * pi * x(2) * x(3) );

  if ( x(2) ~= 0.0 )

    arg = ( x(1) + 2.0 * x(2) + x(3) ) / x(2);
    term = exp ( - arg^2 );

    g(1) = g(1) + 2.0 * term * ( x(1) + 2.0 * x(2) + x(3) ) / x(2)^2;
    g(2) = g(2) - 2.0 * term * ( x(1) + 2.0 * x(2) + x(3) ) ...
      * ( x(1) + x(3) ) / x(2)^3;
    g(3) = g(3) + 2.0 * term * ( x(1) + 2.0 * x(2) + x(3) ) / x(2)^2;

  end

  return
end


text:
function gotoend(H)
%GOTOEND Skip to end of movie in movie player

feval(H.fcns.goto_end,[],[],H.hfig);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function sqexpKernDisplay(kern, spacing)

% SQEXPKERNDISPLAY Display parameters of the SQEXP kernel.
% FORMAT
% DESC displays the parameters of the pre-built compound squared exponential
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : sqexpKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2004

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('RBF Variance: %2.4f\n', kern.rbfVariance)
fprintf(spacing);
fprintf('RBF inverse width: %2.4f\n', kern.inverseWidth)
fprintf(spacing);
fprintf('White noise Variance: %2.4f\n', kern.whiteVariance)
fprintf(spacing);
fprintf('Bias Kernel Variance: %2.4f\n', kern.biasVariance)


text:
function clearToolPanel(handles)
%
% 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/>.

   t = handles.toolPanel;
    tc = get(t, 'Children');
    delete(tc);

    fnames = fieldnames(handles);
    for i = 1:length(fnames)
        if ~ishandle(handles.(fnames{i})) && ~isstruct(handles.(fnames{i}))
            handles = rmfield(handles, fnames{i});
        end
    end
    guidata(handles.mainframe, handles);

%     for i=1:length(tc)
%         set(tc(i), 'visible', 'off');
%     end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function str = valid_named_set_type(obj, set_type)
%VALID_NAMED_SET_TYPE  Returns a label for the given named set type or empty.
%
%   -----  PRIVATE METHOD  -----
%
%   STR = OBJ.VALID_NAMED_SET_TYPE(SET_TYPE)
%
%   Returns a string corresponding to the type of the named set for
%   valid values of SET_TYPE, otherwise an empty string. Can be used to
%   check whether SET_TYPE is a valid type.
%
%   Valid set types, and the corresponding values are defined as a struct
%   returned by the DEF_SET_TYPES method.
%
%   E.g. The following method ...
%       function obj = def_set_types(obj)
%           obj.set_types = struct(...
%                   'var', 'variable', ...
%                   'lin', 'linear constraint' ...
%               );
%       end
%
%   ... indicates the following valid SET_TYPES and corresponding
%   return values:
%
%       SET_TYPE    STR
%       'var'       'variable'
%       'lin'       'linear constraint'

%   MP-Opt-Model
%   Copyright (c) 2017-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if isfield(obj.set_types, set_type)
    str = obj.set_types.(set_type);
else
    str = '';
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function varargout = process_report_email( varargin )
% PROCESS_REPORT_EMAIL: Send current process report by email.
%
% For calling this function from a script, use directly bst_report.m:
% bst_report('Email', ReportFile, to, subject, isFullReport=1)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2021

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Send report by email';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'File';
    sProcess.Index       = 983;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Scripting#Send_report_by_email';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data', 'results', 'timefreq', 'matrix', 'dipoles', 'pdata', 'presults', 'ptimefreq', 'pmatrix', 'import'};
    sProcess.OutputTypes = {'raw', 'data', 'results', 'timefreq', 'matrix', 'dipoles', 'pdata', 'presults', 'ptimefreq', 'pmatrix', 'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    sProcess.isSeparator = 1;
    % Definition of the options
    % === USERNAME
    sProcess.options.username.Comment = 'Brainstorm username: ';
    sProcess.options.username.Type    = 'text';
    sProcess.options.username.Value   = '';
    % === TO
    sProcess.options.cc.Comment = 'Send copy to (email address): ';
    sProcess.options.cc.Type    = 'text';
    sProcess.options.cc.Value   = '';
    % === SUBJECT
    sProcess.options.subject.Comment = 'Subject: ';
    sProcess.options.subject.Type    = 'text';
    sProcess.options.subject.Value   = 'Process completed';
    % === REPORTFILE
    sProcess.options.reportfile.Comment = 'ReportFile: ';
    sProcess.options.reportfile.Type    = 'text';
    sProcess.options.reportfile.Value   = 'current';
    sProcess.options.reportfile.Hidden  = 1;
    % === FULL REPORT
    sProcess.options.full.Comment = 'Send full report';
    sProcess.options.full.Type    = 'checkbox';
    sProcess.options.full.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Returned files: same as input
    OutputFiles = {sInputs.FileName};
    % Brainstorm username
    username = sProcess.options.username.Value;
    if isempty(username)
        bst_report('error', sProcess, [], 'Invalid Brainstorm username.');
        return;
    end
    % CC address
    cc = sProcess.options.cc.Value;
    if ~isempty(cc) && ~any(cc == '@')
        bst_report('error', sProcess, [], 'Invalid email address.');
        return;
    end
    % Email subject
    subject = sProcess.options.subject.Value;
    if isempty(subject)
        subject = 'Brainstorm report';
    end
    % Report file
    reportfile = sProcess.options.reportfile.Value;
    % Full report
    isFullReport = sProcess.options.full.Value;
    % Send email
    [isOk, resp] = bst_report('Email', reportfile, username, cc, subject, isFullReport);
    % Error handling
    if ~isOk
        bst_report('Error', sProcess, [], ['Email could not be sent: ' 10 resp]);
    end
end

author: brainstorm-tools


text:
function [ lines_ali ] = refitLineSegmentB( lines, vp, vpweight )
%REFITLINESEGMENTA refit direction of line segments
%   lines: original line segments
%   vp: vannishing point
%   vpweight: if set to 0, lines will not change; if set to inf, lines will
%   be forced to pass vp
if ~exist('vpweight', 'var')
    vpweight = 0.1;
end

numSample = 100;
numLine = size(lines,1);
xyz = zeros(numSample+1,3);
wei = ones(numSample+1,1); wei(numSample+1) = vpweight*numSample;
lines_ali = lines;
for i = 1:numLine
    n = lines(i,1:3);
    sid = lines(i,5)*2*pi;
    eid = lines(i,6)*2*pi;
    if eid<sid
        x = linspace(sid,eid+2*pi,numSample);
        x = rem(x,2*pi);
%         x = sid-1:(eid-1+numBins);
%         x = rem(x,numBins) + 1;
    else
        x = linspace(sid,eid,numSample);
    end
%     u = -pi + (x'-1)*uBinSize + uBinSize/2;
    u = -pi + x';
    v = computeUVN(n, u, lines(i,4));
    xyz(1:numSample,:) = uv2xyzN([u v], lines(i,4));
    xyz(numSample+1,:) = vp;
    [ ~, outputNM ] = curveFitting( xyz, wei );
    lines_ali(i,1:3) = outputNM;
end

end

author: zouchuhang

repo: LayoutNet


text:
function [granger, v, n] = ft_connectivity_granger(H, Z, S, varargin)

% FT_CONNECTIVITY_GRANGER computes spectrally resolved granger causality. This
% implementation is loosely based on the code used in Brovelli, et. al., PNAS 101,
% 9849-9854 (2004).
%
% Use as
%   [GRANGER, V, N] = FT_CONNECTIVITY_GRANGER(H, Z, S, ...)
%
% The input data should be
%   H = spectral transfer matrix, Nrpt x Nchan x Nchan x Nfreq (x Ntime),
%       or Nrpt x Nchancmb x Nfreq (x Ntime). Nrpt can be 1.
%   Z = the covariance matrix of the noise, Nrpt x Nchan x Nchan (x Ntime),
%       or Nrpt x Nchancmb (x Ntime).
%   S = the cross-spectral density matrix with the same dimensionality as H.
%
% Additional optional input arguments come as key-value pairs:
%   'dimord'  = required string specifying how to interpret the input data
%               supported values are 'rpt_chan_chan_freq(_time) and
%               'rpt_chan_freq(_time), 'rpt_pos_pos_XXX' and 'rpt_pos_XXX'
%   'method'  = 'granger' (default), or 'instantaneous', or 'total'.
%   'hasjack' = 0 (default) is a boolean specifying whether the input
%               contains leave-one-outs, required for correct variance
%               estimate
%   'powindx' = is a variable determining the exact computation, see below
%
% If the inputdata is such that the channel-pairs are linearly indexed, granger
% causality is computed per quadruplet of consecutive entries, where the convention
% is as follows:
%
%  H(:, (k-1)*4 + 1, :, :, :) -> 'chan1-chan1'
%  H(:, (k-1)*4 + 2, :, :, :) -> 'chan1->chan2'
%  H(:, (k-1)*4 + 3, :, :, :) -> 'chan2->chan1'
%  H(:, (k-1)*4 + 4, :, :, :) -> 'chan2->chan2'
%
% The same holds for the Z and S matrices.
%
% Pairwise block-granger causality can be computed when the inputdata has
% dimensionality Nchan x Nchan. In that case powindx should be specified, as a 1x2
% cell-array indexing the individual channels that go into each 'block'.
%
% See also FT_CONNECTIVITYANALYSIS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Undocumented option: powindx can be a struct. In that case, blockwise
% conditional granger can be computed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2009-2017, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

method  = ft_getopt(varargin, 'method',  'granger');
hasjack = ft_getopt(varargin, 'hasjack', 0);
powindx = ft_getopt(varargin, 'powindx', []);
dimord  = ft_getopt(varargin, 'dimord',  []);

%FIXME speed up code and check
siz = size(H);
if numel(siz)==4
  siz(5) = 1;
end
n   = siz(1);
Nc  = siz(2);

outsum = zeros(siz(2:end));
outssq = zeros(siz(2:end));

% crossterms are described by chan_chan_therest
issquare = length(strfind(dimord, 'chan'))==2 || length(strfind(dimord, 'pos'))==2;

switch method
  case 'granger'

    if issquare && isempty(powindx)
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % data are chan_chan_therest
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%

      for kk = 1:n
        for ii = 1:Nc
          for jj = 1:Nc
            if ii ~=jj
              zc     = reshape(Z(kk,jj,jj,:) - Z(kk,ii,jj,:).^2./Z(kk,ii,ii,:),[1 1 1 1 siz(5)]);
              zc     = repmat(zc,[1 1 1 siz(4) 1]);
              numer  = reshape(abs(S(kk,ii,ii,:,:)),[1 1 siz(4:end)]);
              denom  = reshape(abs(S(kk,ii,ii,:,:)-zc.*abs(H(kk,ii,jj,:,:)).^2),[1 1 siz(4:end)]);
              outsum(jj,ii,:,:) = outsum(jj,ii,:,:) + log(numer./denom);
              outssq(jj,ii,:,:) = outssq(jj,ii,:,:) + (log(numer./denom)).^2;
            end
          end
          outsum(ii,ii,:,:) = 0; %self-granger set to zero
        end
      end

    elseif ~issquare && isempty(powindx)
      %%%%%%%%%%%%%%%%%%%%%%%%%%
      %data are linearly indexed
      %%%%%%%%%%%%%%%%%%%%%%%%%%

      for j = 1:n
        for k = 1:Nc
          %FIXME powindx is not used here anymore
          %iauto1  = sum(powindx==powindx(k,1),2)==2;
          %iauto2  = sum(powindx==powindx(k,2),2)==2;
          %icross1 = k;
          %icross2 = sum(powindx==powindx(ones(Nc,1)*k,[2 1]),2)==2;

          % The following is based on hard-coded assumptions: which is fair
          % to do if the order of the labelcmb is according to the output of
          % ft_connectivity_csd2transfer
          if mod(k-1, 4)==0
            continue; % auto granger set to 0
            %iauto1=k;iauto2=k;icross1=k;icross2=k;
          elseif mod(k-1, 4)==1
            iauto1=k+2;iauto2=k-1;icross1=k;icross2=k+1;
          elseif mod(k-1, 4)==2
            iauto1=k-2;iauto2=k+1;icross1=k;icross2=k-1;
          elseif mod(k-1, 4)==3
            continue; % auto granger set to 0
            %iauto1=k;iauto2=k;icross1=k;icross2=k;
          end

          zc      = Z(j,iauto2,:,:) - Z(j,icross1,:,:).^2./Z(j,iauto1,:,:);
          numer   = abs(S(j,iauto1,:,:));
          denom   = abs(S(j,iauto1,:,:)-zc(:,:,ones(1,size(H,3)),:).*abs(H(j,icross1,:,:)).^2);
          outsum(icross2,:,:) = outsum(icross2,:,:) + reshape(log(numer./denom), [1 siz(3:end)]);
          outssq(icross2,:,:) = outssq(icross2,:,:) + reshape((log(numer./denom)).^2, [1 siz(3:end)]);
        end
      end

    elseif issquare && iscell(powindx)
      %%%%%%%%%%%%%%%%%%%
      % blockwise granger
      %%%%%%%%%%%%%%%%%%%

      % H = transfer function nchan x nchan x nfreq
      % Z = noise covariance  nchan x nchan
      % S = crosspectrum      nchan x nchan x nfreq
      % powindx{k} is a list of indices for block k

      nblock = numel(powindx);
      n      = size(H,1);
      nfreq  = size(H,4);

      outsum = zeros(nblock,nblock,nfreq);
      outssq = zeros(nblock,nblock,nfreq);

      for k = 1:nblock
        for m = (k+1):nblock
          indx  = [powindx{k}(:);powindx{m}(:)];
          indx  = cat(1,indx,setdiff((1:size(Z,2))',indx));
          n1    = numel(powindx{k});
          n2    = numel(powindx{m});
          ntot  = numel(indx);
          indx1 = 1:n1;
          indx1_  = (1:ntot)'; indx1_(indx1) = [];
          indx2   = n1+(1:n2);
          indx12_ = (1:ntot)'; indx12_([indx1(:);indx2(:)]) = [];

          for kk = 1:n
            tmpZ = reshape(Z(kk,indx,indx), [ntot ntot]);

            % projection matrix for therest+block2 -> block1
            P1 = [eye(n1)                                zeros(n1,ntot-n1);
              -tmpZ(indx1_,indx1)/tmpZ(indx1,indx1)      eye(ntot-n1)];

            % projection matrix for therest+block1 -> block2
            P2 = [  eye(n1)    -tmpZ(indx1,indx2)/tmpZ(indx2,indx2) zeros(n1,ntot-n1-n2);
                  zeros(n2,n1)   eye(n2) zeros(n2, ntot-n1-n2);
                  zeros(ntot-n1-n2,n1) -tmpZ(indx12_,indx2)/tmpZ(indx2,indx2) eye(ntot-n1-n2)];

            % invert only once
            for jj = 1:nfreq
              % post multiply transfer matrix with the inverse of the projection matrix
              % this is equivalent to time domain pre multiplication with P
              Sj = reshape(S(kk,indx,indx,jj), [ntot ntot]);
              Zj = tmpZ; %(:,:);
              H1 = reshape(H(kk,indx,indx,jj), [ntot ntot])/P1;
              H2 = reshape(H(kk,indx,indx,jj), [ntot ntot])/P2;
              num1 = abs(det(Sj(indx1,indx1))); % numerical round off leads to tiny imaginary components
              num2 = abs(det(Sj(indx2,indx2))); % numerical round off leads to tiny imaginary components
              denom1 = abs(det(H1(indx1,indx1)*Zj(indx1,indx1)*H1(indx1,indx1)'));
              denom2 = abs(det(H2(indx2,indx2)*Zj(indx2,indx2)*H2(indx2,indx2)'));
              %rH1 = real(H1(indx1,indx1));
              %rH2 = real(H2(indx2,indx2));
              %iH1 = imag(H1(indx1,indx1));
              %iH2 = imag(H2(indx2,indx2));
              %h1 = rH1*Zj(indx1,indx1)*rH1' + iH1*Zj(indx1,indx1)*iH1';
              %h2 = rH2*Zj(indx2,indx2)*rH2' + iH2*Zj(indx2,indx2)*iH2';
              %denom1 = abs(det(h1));
              %denom2 = abs(det(h2));

              outsum(m,k,jj) = log( num1./denom1 )    + outsum(m,k,jj);
              outsum(k,m,jj) = log( num2./denom2 )    + outsum(k,m,jj);
              outssq(m,k,jj) = log( num1./denom1 ).^2 + outssq(m,k,jj);
              outssq(k,m,jj) = log( num2./denom2 ).^2 + outssq(k,m,jj);
            end
          end

        end
      end

    elseif ~issquare && isstruct(powindx) && isfield(powindx, 'n')
      %%%%%%%%%%%%%%%%%%%%%%
      %blockwise conditional
      %%%%%%%%%%%%%%%%%%%%%%

      n     = size(H,1);
      ncmb  = size(H,2);
      nfreq = size(H,3);
      ncnd  = size(powindx.cmbindx,1);

      outsum = zeros(ncnd, nfreq);
      outssq = zeros(ncnd, nfreq);
      for k = 1:n
        tmpS = reshape(S, [ncmb nfreq]);
        tmpH = reshape(H, [ncmb nfreq]);
        tmpZ = reshape(Z, [ncmb 1]);
        tmp  = blockwise_conditionalgranger(tmpS,tmpH,tmpZ,powindx.cmbindx,powindx.n);

        outsum = outsum + tmp;
        outssq = outssq + tmp.^2;
      end

    elseif ~issquare && isstruct(powindx)
      %%%%%%%%%%%%%%%%%%%%%%
      %triplet conditional
      %%%%%%%%%%%%%%%%%%%%%%

      % decode from the powindx struct which rows in the data correspond with
      % the triplets, and which correspond with the duplets
      ublockindx = unique(powindx.blockindx);
      nperblock  = zeros(size(ublockindx));
      for k = 1:numel(ublockindx)
        nperblock(k,1) = sum(powindx.blockindx==ublockindx(k));
      end
      if ~all(ismember(nperblock,[4 9]))
        error('the data should be a mixture of trivariate and bivariate decompositions');
      end
      indx_triplets = ismember(powindx.blockindx, ublockindx(nperblock==9)); ntriplets = sum(indx_triplets)./9;
      indx_duplets  = ismember(powindx.blockindx, ublockindx(nperblock==4)); nduplets  = sum(indx_duplets)./4;

      % this assumes well-behaved powindx.cmbindx
      cmbindx2 = reshape(powindx.cmbindx(indx_duplets, 1), 2, [])';
      cmbindx3 = reshape(powindx.cmbindx(indx_triplets,1), 3, [])';

      cmbindx2 = cmbindx2(1:2:end,:);
      cmbindx3 = cmbindx3(1:3:end,:);

      cmbindx  = powindx.outindx;

      n   = size(H,1);
      siz = size(H);

      outsum = zeros(size(cmbindx,1), size(H,3));
      outssq = outsum;

      % call the low-level function
      for k = 1:n
        H3 = reshape(H(k,indx_triplets,:,:), [3 3 ntriplets siz(3:end)]);
        Z3 = reshape(Z(k,indx_triplets),     [3 3 ntriplets]);

        H2 = reshape(H(k,indx_duplets,:,:), [2 2 nduplets siz(3:end)]);
        Z2 = reshape(Z(k,indx_duplets),     [2 2 nduplets]);

        tmp  = triplet_conditionalgranger(H3,Z3,cmbindx3,H2,Z2,cmbindx2,cmbindx);

        outsum = outsum + tmp;
        outssq = outssq + tmp.^2;
      end

    end

  case 'instantaneous'

    if issquare && isempty(powindx)
      % data are chan_chan_therest
      for kk = 1:n
        for ii = 1:Nc
          for jj = 1:Nc
            if ii ~=jj
              zc1    = reshape(Z(kk,jj,jj,:) - Z(kk,ii,jj,:).^2./Z(kk,ii,ii,:),[1 1 1 1 siz(5)]);
              zc2    = reshape(Z(kk,ii,ii,:) - Z(kk,jj,ii,:).^2./Z(kk,jj,jj,:),[1 1 1 1 siz(5)]);
              zc1    = repmat(zc1,[1 1 1 siz(4) 1]);
              zc2    = repmat(zc2,[1 1 1 siz(4) 1]);
              term1  = abs(S(kk,ii,ii,:,:)) - zc1.*abs(H(kk,ii,jj,:,:)).^2;
              term2  = abs(S(kk,jj,jj,:,:)) - zc2.*abs(H(kk,jj,ii,:,:)).^2;
              numer  = term1.*term2;
              denom  = abs(S(kk,ii,ii,:,:).*S(kk,jj,jj,:,:) - S(kk,ii,jj,:,:).*S(kk,jj,ii,:,:));
              outsum(jj,ii,:,:) = outsum(jj,ii,:,:) + reshape(log(numer./denom),  [1 1 siz(4:end)]);
              outssq(jj,ii,:,:) = outssq(jj,ii,:,:) + reshape((log(numer./denom)).^2, [1 1 siz(4:end)]);
            end
          end
          outsum(ii,ii,:,:) = 0; %self-granger set to zero
        end
      end
    elseif ~issquare && isempty(powindx)
      % data are linearly indexed
      for j = 1:n
        for k = 1:Nc
          %iauto1  = sum(powindx==powindx(k,1),2)==2;
          %iauto2  = sum(powindx==powindx(k,2),2)==2;
          %icross1 = k;
          %icross2 = sum(powindx==powindx(ones(Nc,1)*k,[2 1]),2)==2;
          if mod(k-1, 4)==0
            continue; % auto granger set to 0
            %iauto1=k;iauto2=k;icross1=k;icross2=k;
          elseif mod(k-1, 4)==1
            iauto1=k+2;iauto2=k-1;icross1=k;icross2=k+1;
          elseif mod(k-1, 4)==2
            iauto1=k-2;iauto2=k+1;icross1=k;icross2=k-1;
          elseif mod(k-1, 4)==3
            continue; % auto granger set to 0
            %iauto1=k;iauto2=k;icross1=k;icross2=k;
          end

          zc1     = Z(j,iauto1,:, :) - Z(j,icross2,:, :).^2./Z(j,iauto2,:, :);
          zc2     = Z(j,iauto2,:, :) - Z(j,icross1,:, :).^2./Z(j,iauto1,:, :);
          term1   = abs(S(j,iauto2,:,:)) - zc1(:,:,ones(1,size(H,3)),:).*abs(H(j,icross2,:,:)).^2;
          term2   = abs(S(j,iauto1,:,:)) - zc2(:,:,ones(1,size(H,3)),:).*abs(H(j,icross1,:,:)).^2;
          numer   = term1.*term2;
          denom   = abs(S(j,iauto1,:,:).*S(j,iauto2,:,:) - S(j,icross1,:,:).*S(j,icross2,:,:));

          outsum(icross2,:,:) = outsum(icross2,:,:) + reshape(log(numer./denom), [1 siz(3:end)]);
          outssq(icross2,:,:) = outssq(icross2,:,:) + reshape((log(numer./denom)).^2, [1 siz(3:end)]);
        end
      end
    elseif issquare && iscell(powindx)
      % blockwise granger
      % H = transfer function nchan x nchan x nfreq
      % Z = noise covariance  nchan x nchan
      % S = crosspectrum      nchan x nchan x nfreq
      % powindx{1} is a list of indices for block1
      % powindx{2} is a list of indices for block2
      ft_error('instantaneous causality is not implemented for blockwise factorizations');
    elseif isstruct(powindx)
      %blockwise conditional
      ft_error('blockwise conditional instantaneous causality is not implemented');
    else
      ft_error('not implemented');
    end

  case 'total'

    if issquare && isempty(powindx)
      % data are chan_chan_therest
      for kk = 1:n
        for ii = 1:Nc
          for jj = 1:Nc
            if ii ~=jj
              numer  = abs(S(kk,ii,ii,:,:).*S(kk,jj,jj,:,:));
              denom  = abs(S(kk,ii,ii,:,:).*S(kk,jj,jj,:,:) - S(kk,ii,jj,:,:).*S(kk,jj,ii,:,:));
              outsum(jj,ii,:,:) = outsum(jj,ii,:,:) + reshape(log(numer./denom), [1 1 siz(4:end)]);
              outssq(jj,ii,:,:) = outssq(jj,ii,:,:) + reshape((log(numer./denom)).^2, [1 1 siz(4:end)]);
            end
          end
          outsum(ii,ii,:,:) = 0; %self-granger set to zero
        end
      end
    elseif ~issquare && isempty(powindx)
      % data are linearly indexed
      for j = 1:n
        for k = 1:Nc
          %iauto1  = sum(powindx==powindx(k,1),2)==2;
          %iauto2  = sum(powindx==powindx(k,2),2)==2;
          %icross1 = k;
          %icross2 = sum(powindx==powindx(ones(Nc,1)*k,[2 1]),2)==2;
          if mod(k-1, 4)==0
            continue; % auto granger set to 0
            %iauto1=k;iauto2=k;icross1=k;icross2=k;
          elseif mod(k-1, 4)==1
            iauto1=k+2;iauto2=k-1;icross1=k;icross2=k+1;
          elseif mod(k-1, 4)==2
            iauto1=k-2;iauto2=k+1;icross1=k;icross2=k-1;
          elseif mod(k-1, 4)==3
            continue; % auto granger set to 0
            %iauto1=k;iauto2=k;icross1=k;icross2=k;
          end

          numer   = abs(S(j,iauto1,:,:).*S(j,iauto2,:,:));
          denom   = abs(S(j,iauto1,:,:).*S(j,iauto2,:,:) - S(j,icross1,:,:).*S(j,icross2,:,:));
          outsum(icross2,:,:) = outsum(icross2,:,:) + reshape(log(numer./denom), [1 siz(3:end)]);
          outssq(icross2,:,:) = outssq(icross2,:,:) + reshape((log(numer./denom)).^2, [1 siz(3:end)]);
        end
      end
    elseif issquare && iscell(powindx)
      % blockwise total interdependence

      % S = crosspectrum      nchan x nchan x nfreq
      % powindx{k} is a list of indices for block k

      nblock = numel(powindx);
      n      = size(S,1);
      nfreq  = size(S,4);

      outsum = zeros(nblock,nblock,nfreq);
      outssq = zeros(nblock,nblock,nfreq);

      for k = 1:nblock
        for m = (k+1):nblock
          indx  = [powindx{k}(:);powindx{m}(:)];
          n1    = numel(powindx{k});
          n2    = numel(powindx{m});
          ntot  = n1+n2;
          indx1 = 1:n1;
          indx2 = (n1+1):ntot;

          for kk = 1:n
            for jj = 1:nfreq
              Sj = reshape(S(kk,indx,indx,jj), [ntot ntot]);
              num1 = abs(det(Sj(indx1,indx1))); % numerical round off leads to tiny imaginary components
              num2 = abs(det(Sj(indx2,indx2))); % numerical round off leads to tiny imaginary components
              num  = num1.*num2;
              denom = abs(det(Sj));

              outsum(m,k,jj) = log( num./denom )    + outsum(m,k,jj);
              outsum(k,m,jj) = log( num./denom )    + outsum(k,m,jj);
              outssq(m,k,jj) = log( num./denom ).^2 + outssq(m,k,jj);
              outssq(k,m,jj) = log( num./denom ).^2 + outssq(k,m,jj);
            end
          end

        end
      end

    elseif issquare && isstruct(powindx)
      %blockwise conditional
      ft_error('blockwise conditional total interdependence is not implemented');
    end

  case 'iis'
    ft_warning('THIS IS EXPERIMENTAL CODE, USE AT YOUR OWN RISK!');
    % this is experimental
    if ~issquare && isempty(powindx)
      A = transfer2coeffs(shiftdim(H),(0:size(H,3)-1));
      ncmb = size(A,1)./4;
      iis = coeffs2iis(reshape(A,[2 2 ncmb size(A,2)]),reshape(Z,[2 2 ncmb]));
      iis = repmat(iis(:),[1 4])';
      outsum = iis(:);
      outssq = nan(size(outsum));
    else
      ft_error('iis can only be computed when the input contains sets of bivariate factorizations');
    end

  otherwise
    ft_error('unsupported output requested');
end

granger = outsum./n;
if n>1
  if hasjack
    bias = (n-1).^2;
  else
    bias = 1;
  end
  v = bias*(outssq - (outsum.^2)./n)./(n - 1);
else
  v = [];
end

author: spm


text:
function [system] = LoadTransSystem(model, load_path)
    % Load the hybrid system model

    if nargin < 2
        load_path = [];
    end

    rightStance1 = sys.domains.RightStance(model, load_path);
    rightStance1.UserNlpConstraint = @trans_opt.callback.RightStance1Constraints;
    leftStance1 = sys.domains.LeftStance(model, load_path);
    leftStance1.UserNlpConstraint = @trans_opt.callback.LeftStance1Constraints;

    rightStance2 = copy(rightStance1);
    rightStance2.setName('RightStance2');
    rightStance2.UserNlpConstraint = @trans_opt.callback.RightStance2Constraints;
    leftStance2 = copy(leftStance1);
    leftStance2.setName('LeftStance2');
    leftStance2.UserNlpConstraint = @trans_opt.callback.LeftStance2Constraints;

    leftImpact1 = RigidImpact('LeftImpact', leftStance1, 'leftFootHeight'); % To leftStance
    leftImpact1.addImpactConstraint(struct2array(leftStance1.HolonomicConstraints), load_path);
    leftImpact1.UserNlpConstraint = @trans_opt.callback.LeftImpactConstraints;

    rightImpact1 = RigidImpact('RightImpact', rightStance1, 'rightFootHeight'); % To rightStance
    rightImpact1.addImpactConstraint(struct2array(rightStance1.HolonomicConstraints), load_path);
    rightImpact1.UserNlpConstraint = @trans_opt.callback.RightImpactConstraints;

    leftImpact2 = copy(leftImpact1);
    %     leftImpact2 = RigidImpact('LeftImpact2', leftStance2, 'leftFootHeight'); % To leftStance
    %     leftImpact2.addImpactConstraint(struct2array(leftStance2.HolonomicConstraints), load_path);
    %     leftImpact2.UserNlpConstraint = @trans_opt.callback.LeftImpact2Constraints;

    %     rightImpact2 = RigidImpact('RightImpact2', rightStance2, 'rightFootHeight'); % To rightStance
    %     rightImpact2.addImpactConstraint(struct2array(rightStance2.HolonomicConstraints), load_path);
    %     rightImpact2 = copy(rightImpact1);
    %     rightImpact2.setName('RightImpact2');
    %     rightImpact2.UserNlpConstraint = @trans_opt.callback.RightImpact2Constraints;
    %

    io_control  = IOFeedback('IO');

    system = HybridSystem('atrias');
    system = addVertex(system, 'RightStance1', 'Domain', rightStance1, ...
        'Control', io_control);
    system = addVertex(system, 'LeftStance1', 'Domain', leftStance1, ...
        'Control', io_control);
    system = addVertex(system, 'RightStance2', 'Domain', rightStance2, ...
        'Control', io_control);
    system = addVertex(system, 'LeftStance2', 'Domain', leftStance2, ...
        'Control', io_control);

    srcs = {'RightStance1'
        'LeftStance1'
        'RightStance2'};

    tars = {'LeftStance1'
        'RightStance2'
        'LeftStance2'};

    system = addEdge(system, srcs, tars);
    system = setEdgeProperties(system, srcs, tars, ...
        'Guard', {leftImpact1, rightImpact1, leftImpact2});

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev