text:
function [rcnn_model, rcnn_k_fold_model] = ...
    rcnn_train(imdb, varargin)
% [rcnn_model, rcnn_k_fold_model] = rcnn_train(imdb, varargin)
%   Trains an R-CNN detector for all classes in the imdb.
%
%   Keys that can be passed in:
%
%   svm_C             SVM regularization parameter
%   bias_mult         Bias feature value (for liblinear)
%   pos_loss_weight   Cost factor on hinge loss for positives
%   layer             Feature layer to use (either 5, 6 or 7)
%   k_folds           Train on folds of the imdb
%   checkpoint        Save the rcnn_model every checkpoint images
%   crop_mode         Crop mode (either 'warp' or 'square')
%   crop_padding      Amount of padding in crop
%   net_file          Path to the Caffe CNN to use
%   cache_name        Path to the precomputed feature cache

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

% TODO:
%  - allow training just a subset of classes

ip = inputParser;
ip.addRequired('imdb', @isstruct);
ip.addParamValue('svm_C',           10^-3,  @isscalar);
ip.addParamValue('bias_mult',       10,     @isscalar);
ip.addParamValue('pos_loss_weight', 2,      @isscalar);
ip.addParamValue('layer',           7,      @isscalar);
ip.addParamValue('k_folds',         2,      @isscalar);
ip.addParamValue('checkpoint',      0,      @isscalar);
ip.addParamValue('crop_mode',       'warp', @isstr);
ip.addParamValue('crop_padding',    16,     @isscalar);
ip.addParamValue('net_file', ...
    './data/caffe_nets/finetune_voc_2007_trainval_iter_70k', ...
    @isstr);
ip.addParamValue('cache_name', ...
    'v1_finetune_voc_2007_trainval_iter_70000', @isstr);

ip.parse(imdb, varargin{:});
opts = ip.Results;

opts.net_def_file = './model-defs/rcnn_batch_256_output_fc7.prototxt';

conf = rcnn_config('sub_dir', imdb.name);

% Record a log of the training and test procedure
timestamp = datestr(datevec(now()), 'dd.mmm.yyyy:HH.MM.SS');
diary_file = [conf.cache_dir 'rcnn_train_' timestamp '.txt'];
diary(diary_file);
fprintf('Logging output in %s\n', diary_file);

fprintf('\n\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf('Training options:\n');
disp(opts);
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n');

% ------------------------------------------------------------------------
% Create a new rcnn model
rcnn_model = rcnn_create_model(opts.net_def_file, opts.net_file, opts.cache_name);
rcnn_model = rcnn_load_model(rcnn_model, conf.use_gpu);
rcnn_model.detectors.crop_mode = opts.crop_mode;
rcnn_model.detectors.crop_padding = opts.crop_padding;
rcnn_model.classes = imdb.classes;
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
% Get the average norm of the features
opts.feat_norm_mean = rcnn_feature_stats(imdb, opts.layer, rcnn_model);
fprintf('average norm = %.3f\n', opts.feat_norm_mean);
rcnn_model.training_opts = opts;
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
% Get all positive examples
% We cache only the pool5 features and convert them on-the-fly to
% fc6 or fc7 as required
save_file = sprintf('./feat_cache/%s/%s/gt_pos_layer_5_cache.mat', ...
    rcnn_model.cache_name, imdb.name);
try
  load(save_file);
  fprintf('Loaded saved positives from ground truth boxes\n');
catch
  [X_pos, keys_pos] = get_positive_pool5_features(imdb, opts);
  save(save_file, 'X_pos', 'keys_pos', '-v7.3');
end
% Init training caches
caches = {};
for i = imdb.class_ids
  fprintf('%14s has %6d positive instances\n', ...
      imdb.classes{i}, size(X_pos{i},1));
  X_pos{i} = rcnn_pool5_to_fcX(X_pos{i}, opts.layer, rcnn_model);
  X_pos{i} = rcnn_scale_features(X_pos{i}, opts.feat_norm_mean);
  caches{i} = init_cache(X_pos{i}, keys_pos{i});
end
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
% Train with hard negative mining
first_time = true;
% one pass over the data is enough
max_hard_epochs = 1;

for hard_epoch = 1:max_hard_epochs
  for i = 1:length(imdb.image_ids)
    fprintf('%s: hard neg epoch: %d/%d image: %d/%d\n', ...
            procid(), hard_epoch, max_hard_epochs, i, length(imdb.image_ids));

    % Get hard negatives for all classes at once (avoids loading feature cache
    % more than once)
    [X, keys] = sample_negative_features(first_time, rcnn_model, caches, ...
        imdb, i);

    % Add sampled negatives to each classes training cache, removing
    % duplicates
    for j = imdb.class_ids
      if ~isempty(keys{j})
        if ~isempty(caches{j}.keys_neg)
          [~, ~, dups] = intersect(caches{j}.keys_neg, keys{j}, 'rows');
          assert(isempty(dups));
        end
        caches{j}.X_neg = cat(1, caches{j}.X_neg, X{j});
        caches{j}.keys_neg = cat(1, caches{j}.keys_neg, keys{j});
        caches{j}.num_added = caches{j}.num_added + size(keys{j},1);
      end

      % Update model if
      %  - first time seeing negatives
      %  - more than retrain_limit negatives have been added
      %  - its the final image of the final epoch
      is_last_time = (hard_epoch == max_hard_epochs && i == length(imdb.image_ids));
      hit_retrain_limit = (caches{j}.num_added > caches{j}.retrain_limit);
      if (first_time || hit_retrain_limit || is_last_time) && ...
          ~isempty(caches{j}.X_neg)
        fprintf('>>> Updating %s detector <<<\n', imdb.classes{j});
        fprintf('Cache holds %d pos examples %d neg examples\n', ...
                size(caches{j}.X_pos,1), size(caches{j}.X_neg,1));
        [new_w, new_b] = update_model(caches{j}, opts);
        rcnn_model.detectors.W(:, j) = new_w;
        rcnn_model.detectors.B(j) = new_b;
        caches{j}.num_added = 0;

        z_pos = caches{j}.X_pos * new_w + new_b;
        z_neg = caches{j}.X_neg * new_w + new_b;

        caches{j}.pos_loss(end+1) = opts.svm_C * opts.pos_loss_weight * ...
                                    sum(max(0, 1 - z_pos));
        caches{j}.neg_loss(end+1) = opts.svm_C * sum(max(0, 1 + z_neg));
        caches{j}.reg_loss(end+1) = 0.5 * new_w' * new_w + ...
                                    0.5 * (new_b / opts.bias_mult)^2;
        caches{j}.tot_loss(end+1) = caches{j}.pos_loss(end) + ...
                                    caches{j}.neg_loss(end) + ...
                                    caches{j}.reg_loss(end);

        for t = 1:length(caches{j}.tot_loss)
          fprintf('    %2d: obj val: %.3f = %.3f (pos) + %.3f (neg) + %.3f (reg)\n', ...
                  t, caches{j}.tot_loss(t), caches{j}.pos_loss(t), ...
                  caches{j}.neg_loss(t), caches{j}.reg_loss(t));
        end

        % store negative support vectors for visualizing later
        SVs_neg = find(z_neg > -1 - eps);
        rcnn_model.SVs.keys_neg{j} = caches{j}.keys_neg(SVs_neg, :);
        rcnn_model.SVs.scores_neg{j} = z_neg(SVs_neg);

        % evict easy examples
        easy = find(z_neg < caches{j}.evict_thresh);
        caches{j}.X_neg(easy,:) = [];
        caches{j}.keys_neg(easy,:) = [];
        fprintf('  Pruning easy negatives\n');
        fprintf('  Cache holds %d pos examples %d neg examples\n', ...
                size(caches{j}.X_pos,1), size(caches{j}.X_neg,1));
        fprintf('  %d pos support vectors\n', numel(find(z_pos <  1 + eps)));
        fprintf('  %d neg support vectors\n', numel(find(z_neg > -1 - eps)));
      end
    end
    first_time = false;

    if opts.checkpoint > 0 && mod(i, opts.checkpoint) == 0
      save([conf.cache_dir 'rcnn_model'], 'rcnn_model');
    end
  end
end
% save the final rcnn_model
save([conf.cache_dir 'rcnn_model'], 'rcnn_model');
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
if opts.k_folds > 0
  rcnn_k_fold_model = rcnn_model;
  [W, B, folds] = update_model_k_fold(rcnn_model, caches, imdb);
  rcnn_k_fold_model.folds = folds;
  for f = 1:length(folds)
    rcnn_k_fold_model.detectors(f).W = W{f};
    rcnn_k_fold_model.detectors(f).B = B{f};
  end
  save([conf.cache_dir 'rcnn_k_fold_model'], 'rcnn_k_fold_model');
else
  rcnn_k_fold_model = [];
end
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
function [X_neg, keys] = sample_negative_features(first_time, rcnn_model, ...
                                                  caches, imdb, ind)
% ------------------------------------------------------------------------
opts = rcnn_model.training_opts;

d = rcnn_load_cached_pool5_features(opts.cache_name, ...
    imdb.name, imdb.image_ids{ind});

class_ids = imdb.class_ids;

if isempty(d.feat)
  X_neg = cell(max(class_ids), 1);
  keys = cell(max(class_ids), 1);
  return;
end

d.feat = rcnn_pool5_to_fcX(d.feat, opts.layer, rcnn_model);
d.feat = rcnn_scale_features(d.feat, opts.feat_norm_mean);

neg_ovr_thresh = 0.3;

if first_time
  for cls_id = class_ids
    I = find(d.overlap(:, cls_id) < neg_ovr_thresh);
    X_neg{cls_id} = d.feat(I,:);
    keys{cls_id} = [ind*ones(length(I),1) I];
  end
else
  zs = bsxfun(@plus, d.feat*rcnn_model.detectors.W, rcnn_model.detectors.B);
  for cls_id = class_ids
    z = zs(:, cls_id);
    I = find((z > caches{cls_id}.hard_thresh) & ...
             (d.overlap(:, cls_id) < neg_ovr_thresh));

    % Avoid adding duplicate features
    keys_ = [ind*ones(length(I),1) I];
    if ~isempty(caches{cls_id}.keys_neg) && ~isempty(keys_)
      [~, ~, dups] = intersect(caches{cls_id}.keys_neg, keys_, 'rows');
      keep = setdiff(1:size(keys_,1), dups);
      I = I(keep);
    end

    % Unique hard negatives
    X_neg{cls_id} = d.feat(I,:);
    keys{cls_id} = [ind*ones(length(I),1) I];
  end
end

% ------------------------------------------------------------------------
function [w, b] = update_model(cache, opts, pos_inds, neg_inds)
% ------------------------------------------------------------------------
solver = 'liblinear';
liblinear_type = 3;  % l2 regularized l1 hinge loss
%liblinear_type = 5; % l1 regularized l2 hinge loss

if ~exist('pos_inds', 'var') || isempty(pos_inds)
  num_pos = size(cache.X_pos, 1);
  pos_inds = 1:num_pos;
else
  num_pos = length(pos_inds);
  fprintf('[subset mode] using %d out of %d total positives\n', ...
      num_pos, size(cache.X_pos,1));
end
if ~exist('neg_inds', 'var') || isempty(neg_inds)
  num_neg = size(cache.X_neg, 1);
  neg_inds = 1:num_neg;
else
  num_neg = length(neg_inds);
  fprintf('[subset mode] using %d out of %d total negatives\n', ...
      num_neg, size(cache.X_neg,1));
end

switch solver
  case 'liblinear'
    ll_opts = sprintf('-w1 %.5f -c %.5f -s %d -B %.5f', ...
                      opts.pos_loss_weight, opts.svm_C, ...
                      liblinear_type, opts.bias_mult);
    fprintf('liblinear opts: %s\n', ll_opts);
    X = sparse(size(cache.X_pos,2), num_pos+num_neg);
    X(:,1:num_pos) = cache.X_pos(pos_inds,:)';
    X(:,num_pos+1:end) = cache.X_neg(neg_inds,:)';
    y = cat(1, ones(num_pos,1), -ones(num_neg,1));
    llm = liblinear_train(y, X, ll_opts, 'col');
    w = single(llm.w(1:end-1)');
    b = single(llm.w(end)*opts.bias_mult);

  otherwise
    error('unknown solver: %s', solver);
end

% ------------------------------------------------------------------------
function [W, B, folds] = update_model_k_fold(rcnn_model, caches, imdb)
% ------------------------------------------------------------------------
opts = rcnn_model.training_opts;
num_images = length(imdb.image_ids);
folds = create_folds(num_images, opts.k_folds);
W = cell(opts.k_folds, 1);
B = cell(opts.k_folds, 1);

fprintf('Training k-fold models\n');
for i = imdb.class_ids
  fprintf('\n\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
  fprintf('Training folds for class %s\n', imdb.classes{i});
  fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n');
  for f = 1:length(folds)
    fprintf('Holding out fold %d\n', f);
    [pos_inds, neg_inds] = get_cache_inds_from_fold(caches{i}, folds{f});
    [new_w, new_b] = update_model(caches{i}, opts, ...
        pos_inds, neg_inds);
    W{f}(:,i) = new_w;
    B{f}(i) = new_b;
  end
end

% ------------------------------------------------------------------------
function [pos_inds, neg_inds] = get_cache_inds_from_fold(cache, fold)
% ------------------------------------------------------------------------
pos_inds = find(ismember(cache.keys_pos(:,1), fold) == false);
neg_inds = find(ismember(cache.keys_neg(:,1), fold) == false);

% ------------------------------------------------------------------------
function [X_pos, keys] = get_positive_pool5_features(imdb, opts)
% ------------------------------------------------------------------------
X_pos = cell(max(imdb.class_ids), 1);
keys = cell(max(imdb.class_ids), 1);

for i = 1:length(imdb.image_ids)
  tic_toc_print('%s: pos features %d/%d\n', ...
                procid(), i, length(imdb.image_ids));

  d = rcnn_load_cached_pool5_features(opts.cache_name, ...
      imdb.name, imdb.image_ids{i});

  for j = imdb.class_ids
    if isempty(X_pos{j})
      X_pos{j} = single([]);
      keys{j} = [];
    end
    sel = find(d.class == j);
    if ~isempty(sel)
      X_pos{j} = cat(1, X_pos{j}, d.feat(sel,:));
      keys{j} = cat(1, keys{j}, [i*ones(length(sel),1) sel]);
    end
  end
end

% ------------------------------------------------------------------------
function cache = init_cache(X_pos, keys_pos)
% ------------------------------------------------------------------------
cache.X_pos = X_pos;
cache.X_neg = single([]);
cache.keys_neg = [];
cache.keys_pos = keys_pos;
cache.num_added = 0;
cache.retrain_limit = 2000;
cache.evict_thresh = -1.2;
cache.hard_thresh = -1.0001;
cache.pos_loss = [];
cache.neg_loss = [];
cache.reg_loss = [];
cache.tot_loss = [];

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function J = CalcJacobian_rot(idx)
% Jacobian matrix of current configration in World frame
global uLINK

jsize = length(idx);
target = uLINK(idx(end)).p;   % absolute target position
J = zeros(3,jsize);

for n=1:jsize
    j = idx(n);
    a = uLINK(j).R * uLINK(j).a;  % joint axis vector in world frame
    J(:,n) = cross(a, target - uLINK(j).p);
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function spherical_harmonic_test02 ( )

%*****************************************************************************80
%
%% SPHERICAL_HARMONIC_TEST02 tests SPHERICAL_HARMONIC and SPHERICAL_HARMONIC_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERICAL_HARMONIC_TEST02:\n' );
  fprintf ( 1, '  SPHERICAL_HARMONIC evaluates the\n' );
  fprintf ( 1, '    spherical harmonic function.\n' );
  fprintf ( 1, '  SPHERICAL_HARMONIC_VALUES returns some exact values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      L       M    THETA   PHI   YR   YI\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, l, m, theta, phi, yr, yi ] = spherical_harmonic_values ( n_data );

    if ( n_data == 0 )
      break
    end

    [ c, s ] = spherical_harmonic ( l, m, theta, phi );

    fprintf ( 1, '  %6d  %6d  %6f  %6f  %12f  %12f\n', ...
      l, m, theta, phi, yr, yi );
    fprintf ( 1, '                                      %12f  %12f\n', ...
      c(l+1), s(l+1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fevd_record]=fevd_ols_IV(struct_irf_record,gamma_record,It,Bu,IRFperiods,n)

% function [fevd_record]=bear.fevd(struct_irf_record,gamma_record,It,Bu,IRFperiods,n)
% runs the gibbs sampler to obtain draws from the posterior distribution of FEVD
% inputs:  - cell 'struct_irf_record': record of the gibbs sampler draws for the orthogonalised IRFs
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
% outputs: - cell 'fevd_record': record of the gibbs sampler draws for the FEVD

It=3;
Bu=0;

% this function implements algorithm 3.1.1

% preliminary tasks

% define the time horizon of FEVD as that of IRFs
FEVDperiods=IRFperiods;

% create the first cell
temp=cell(n,n+1);

% now prepare the evaluation of (3.1.13)
% start by filling the first column of every Tij matrix in the cell
% loop over rows of temp
for jj=1:n
   % loop over columns of temp
   for ii=1:n
   % square each IRF element
   temp{jj,ii}(:,1)=struct_irf_record{jj,ii}(:,1).^2;
   end
end

% fill all the other entries of the Tij matrices
% loop over rows of temp
for jj=1:n

   % loop over columns of temp
   for ii=1:n

      % loop over remaining columns
      for kk=2:FEVDperiods

      % define the column as the square of the corresponding column in orthogonalised_irf_record
      % additioned to the value of the preceeding columns, which creates the cumulation
      temp{jj,ii}(:,kk)=struct_irf_record{jj,ii}(:,kk).^2+temp{jj,ii}(:,kk-1);

      end
   end
end

% multiply each matrix in the cell by the variance of the structural shocks
% to do so, loop over simulations (rows of the Tij matrices)
for jj=1:It-Bu
% recover the covariance matrix of structural shocks gamma for this iteration
gamma=reshape(gamma_record,n,n);
% loop over rows of temp
   for ii=1:n
   % loop over columns of temp
      for kk=1:n
      % multiply row jj of the matrix by the variance of the structural shock
      temp{ii,kk}(jj,:)=temp{ii,kk}(jj,:)*gamma(kk,kk);
      end
   end
% then go for next iteration
end

% obtain now the values for Ti, the (n+1)th matrix of each row

% loop over rows of temp
for jj=1:n

% start the summation over Tij matrices
temp{jj,n+1}=temp{jj,1};
   % sum over remaining columns
   for ii=2:n
   temp{jj,n+1}=temp{jj,n+1}+temp{jj,ii};
   end
end

% create the output cell fevd_record
fevd_record=cell(n,n);

% fill the cell
% loop over rows of fevd_record
for jj=1:n
   % loop over columns of fevd_record
   for ii=1:n
   % define the matrix Vfij as the division (pairwise entry) of Tfij by Tfj
   fevd_record{jj,ii}=temp{jj,ii}./temp{jj,n+1};
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
classdef CEC2008_F4 < PROBLEM
% <single> <real> <large/none> <expensive/none>
% Shifted Rastrign's function

%------------------------------- Reference --------------------------------
% K. Tang, X. Yao, P. N. Suganthan, C. MacNish, Y.-P. Chen, C.-M. Chen, and
% Z. Yang, Benchmark functions for the CEC'2008 special session and
% competition on large scale global optimization, Nature Inspired
% Computation and Applications Laboratory, USTC, China, 2007.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2008.mat'),'Data');
            obj.O = Data{4};
            obj.M = 1;
            if isempty(obj.D); obj.D = 100; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 5;
            obj.upper    = zeros(1,obj.D) + 5;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(Z.^2-10*cos(2*pi*Z)+10,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout=ann2rr(varargin)
%
% [RR,tms]=ann2rr(recordName,annotator,N,N0,consecutiveOnly)
%
%    Wrapper to WFDB ANN2RR:
%         http://www.physionet.org/physiotools/wag/ann2rr-1.htm
%
%    Reads a WFDB record and Annotation file to return:
%
%
% RR
%       Nx1 vector of integers representing the duration of the RR
%       interval in samples.
%
% tms
%       Nx1 vector of integers representing the begining of the RR
%       interval in samples.
%
% Required Parameters:
%
% recorName
%       String specifying the name of the record in the WFDB path or
%       in the current directory.
%
% annotator  -
%       String specifying the name of the annotation file in the WFDB path or
%       in the current directory.
%
% Optional Parameters are:
%
% N
%       A 1x1 integer specifying the sample number at which to stop reading the
%       record file (default read all = N).
% N0
%       A 1x1 integer specifying the sample number at which to start reading the
%       annotion file (default 1 = begining of the record).
%
% consecutiveOnly
%       A 1x1 boolean. If true, prints intervals between consecutive valid
%       annotaions only (default =true).
%
%
% Written by Ikaro Silva, 2013
% Last Modified: January, 16, 2013
% Version 1.1
%
% Since 0.0.1
% %Example
%[rr,tm]=ann2rr('challenge/2013/set-a/a01','fqrs');

%endOfHelp

[javaWfdbExec,config]=getWfdbClass('ann2rr');

%Set default pararamter values
inputs={'recordName','annotator','N','N0','consecutiveOnly'};
outputs={'data(:,2)','data(:,1)'};
N=[];
N0=1;
consecutiveOnly=1;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

N0=num2str(N0-1); %-1 is necessary because WFDB is 0 based indexed.
wfdb_argument={'-r',recordName,'-a',annotator,'-f',['s' N0]};

if(~isempty(N))
    wfdb_argument{end+1}='-t';
    %-1 is necessary because WFDB is 0 based indexed.
    wfdb_argument{end+1}=['s' num2str(N-1)];
end

if(consecutiveOnly)
    wfdb_argument{end+1}='-c';
end
wfdb_argument{end+1}='-V';

data=javaWfdbExec.execToDoubleArray(wfdb_argument);
data=wfdbjava2mat(data);
for n=1:nargout
    eval(['varargout{n}=' outputs{n} ';']);
end

author: ikarosilva

repo: wfdb-app-toolbox


text:
clear;
clc
close all;

%% EXAMPLE 2.5(a)
 lat = deg2rad(45);
 lon = deg2rad(45);
 hgt = 0;

 vel_ned = [5 5 5]';
 a_ned = [5 0 0]';
 w_ned = [0.001 0 0]';

[R_meridian, R_transverse, Cecef2ned, Cecef2enu]= ch_earth(lat, lon, hgt);

Cned2ecef = Cecef2ned'

%% Velocity transformation
vel_ecef = Cned2ecef *vel_ned;
vel_ecef
%% Acceleration transformation
a_ecef = Cned2ecef*a_ned;
a_ecef
%% Angular rate transformation
w_ecef = Cned2ecef*w_ned;
w_ecef

%% EXAMPLE 2.5(b)
vel_ecef = [-8.53553391 -1.4644609 4.44089E-16]';
a_ecef = [-2.5 -2.5 3.535534]';
w_ecef = [-0.0005 -0.0005 0.000708]';
lat = deg2rad(45);
lon = deg2rad(45);
h = 0;
[~, ~, Cecef2ned, ~]= ch_earth(lat, lon, hgt);

%% Velocity transformation
v_ned = Cecef2ned*vel_ecef;
v_ned
%% Acceleration transformation
a_ned = Cecef2ned*a_ecef;
a_ned
%% Angular rate transformation
w_ned = Cecef2ned*w_ecef;
w_ned

author: yandld


text:
function nullStats = GiveMeSimpleNullStats(groupLabels,numNulls,cfnParams,meanOverFolds)
% Returns simple permutation-based null samples
%-------------------------------------------------------------------------------

% Check inputs and set defaults:
if nargin < 4
    meanOverFolds = true;
end

%-------------------------------------------------------------------------------
numDataSamples = length(groupLabels);

if meanOverFolds
    nullStats = zeros(numNulls,1);
else
    nullStats = zeros(numNulls,cfnParams.numFolds);
end

for i = 1:numNulls
    shuffledLabels = groupLabels(randperm(numDataSamples));

    % One repeat of 10-fold cross validation using the shuffled labels:
    cvFolds = cvpartition(groupLabels,'KFold',cfnParams.numFolds,'Stratify',true);
    nullAcc_k = zeros(cfnParams.numFolds,1);
    for k = 1:cfnParams.numFolds
        yTrue = groupLabels(cvFolds.test(k));
        yPredict = shuffledLabels(cvFolds.test(k));
        nullAcc_k(k) = BF_LossFunction(yTrue,yPredict,cfnParams.whatLoss,cfnParams.classLabels);
    end

    % Agglomerate over folds:
    if meanOverFolds
        nullStats(i) = mean(nullAcc_k);
    else
        nullStats(i,:) = nullAcc_k;
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [vertices, faces] = steinerPolytope(vectors)
%STEINERPOLYTOPE Create a steiner polytope from a set of vectors
%
%   [VERTICES FACES] = steinerPolygon(VECTORS)
%   Creates the Steiner polytope defined by the set of vectors VECTORS.
%
%   Example
%     % Creates and display a planar Steiner polytope (ie, a polygon)
%     [v f] = steinerPolytope([1 0;0 1;1 1]);
%     fillPolygon(v);
%
%     % Creates and display a 3D Steiner polytope
%     [v f] = steinerPolytope([1 0 0;0 1 0;0 0 1;1 1 1]);
%     drawMesh(v, f);
%     view(3); axis vis3d
%
%   See also
%   meshes3d, drawMesh, steinerPolygon, mergeCoplanarFaces
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2006-04-28
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% History
% 2013-02-22 merge coplanar faces, add management of 2D case, update doc

% compute vectors dimension
nd = size(vectors, 2);

% create candidate vertices
vertices = zeros(1, size(vectors, 2));
for i = 1:length(vectors)
    nv = size(vertices, 1);
    vertices = [vertices; vertices+repmat(vectors(i,:), [nv 1])]; %#ok<AGROW>
end

if nd == 2
    % for planar case, use specific function convhull
    K = convhull(vertices(:,1), vertices(:,2));
    vertices = vertices(K, :);
    faces = 1:length(K);

else
    % Process the general case (tested only for nd==3)

    % compute convex hull
    K = convhulln(vertices);

    % keep only relevant points, and update faces indices
    ind = unique(K);
    for i = 1:length(ind)
        K(K==ind(i)) = i;
    end

    % return results
    vertices = vertices(ind, :);
    faces = K;

    % in case of 3D meshes, merge coplanar faces
    if nd == 3
        faces = mergeCoplanarFaces(vertices, faces);
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [] = draw_crystal_lattice
% File:      draw_crystal_lattice.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2011.01.30 -
% Language:  MATLAB R2012a
% Purpose:   draw crystal lattices
% Copyright: Ioannis Filippidis, 2011-

selection = 'diamond'; % 'diamond', 'copper'

a = 1;
c = a /4;

r_diamond_cubic = c *[0, 0, 0;
        2, 0, 0;
        0, 2, 0;
        2, 2, 0;
        1, 1, 0; %
        1, 0, 1;
        0, 1, 1;
        1, 2, 1;
        2, 1, 1; %
        0, 0, 2;
        2, 0, 2;
        0, 2, 2;
        2, 2, 2;
        1, 1, 2; %
        0.5, 0.5, 0.5;
        1.5, 1.5, 0.5;
        0.5, 1.5, 1.5;
        1.5, 0.5, 1.5]';

r_fcc = c *[0, 0, 0;
        2, 0, 0;
        0, 2, 0;
        2, 2, 0;
        1, 1, 0; %
        1, 0, 1;
        0, 1, 1;
        1, 2, 1;
        2, 1, 1; %
        0, 0, 2;
        2, 0, 2;
        0, 2, 2;
        2, 2, 2;
        1, 1, 2]'; %

i_diamond_cubic =[15, 1;
     15, 5;
     15, 6;
     15, 7; %
     16, 5;
     16, 4;
     16, 8;
     16, 9; %
     17, 7;
     17, 12;
     17, 14;
     17, 8; %
     18, 6;
     18, 9;
     18, 14;
     18, 11];
%

i_fcc = [1, 2;
    2, 4;
    2, 3;
    1, 3;
    1, 4;
    3, 4; %
    1, 10;
    2, 10;
    2, 11;
    1, 11;
    10, 11; %
    3, 12;
    4, 13;
    12, 13;
    4, 12;
    3, 13; %
    10, 12;
    10, 13;
    13, 11;
    12, 11;
    1, 12;
    3, 10; %
    2, 13;
    4, 11];

switch selection
    case 'diamond'
        r = r_diamond_cubic;
        i1 = i_diamond_cubic;
    case 'copper'
        r = r_fcc;
        i1 = i_fcc;
    otherwise
        disp('Unknown material selection.')
end

ax = newax;
hold(ax, 'on')
for i=1:size(r, 2)
    h = drawSphere(r(:, i), 0.05);
    set(h, 'FaceColor', 'r')
end
grid(ax, 'on')
box(ax, 'on')
axis(ax, 'equal')
axis(ax, 'tight')
view(ax, 3)
%maximize(get(ax, 'Parent') )

for i=1:size(i1, 1)
    a = i1(i, 1);
    b = i1(i, 2);
    drawCylinder([r(:, a)', r(:, b)', 0.01] )
end

fname = selection;
fig2u3d(gca, fname, '-pdf');

author: Sable

repo: mcbench-benchmarks


text:
function [kappa,k,theta0,P,b,ebuf,Zf,Zw] = celpana(x,L,M,c,cb,Pidx,bbuf,ebuf,Zf,Zw)
%  celpana --> CELP analyzer (coder).
%
%    [kappa,k,theta0,P,b,ebuf,Zf,Zw] = celpana(x,L,M,c,cb,Pidx,bbuf,ebuf,Zf,Zw)
%
%    The function implements a CELP coder using the following steps:
%
%    (1) Find the reflection coefficients, kappa, using M'th order
%    LP analysis on the signal frame x of length N.
%
%    (2) Find the coefficients of the filter function A(z/c) used in
%    the perceptual weighting filter.
%
%    (3) Determine the excitation parameters k, theta0, P, and b, used
%    to generate the excitation sequence, e(n). The parameters will be
%    estimated in blocks of length L, so N/L values are obtained for
%    the single input frame. Other inputs used here, are the Gaussian
%    codebook given by the L-by-K matrix cb, the pitch search range
%    Pidx(1) < P < Pidx(2), the last estimated b in bbuf, Pidx(2) previous
%    excitation samples buffered in the vector ebuf, and Zf and Zw which
%    are the memory hangover in the filters 1/A(z/c) and W(z) = A(z/c)/A(z),
%    respectively.

N = length(x);                          % Frame length.
J = N/L;                                % Number of sub-blocks.

k      = zeros(J,1);
theta0 = zeros(J,1);
P      = zeros(J,1);
b      = zeros(J,1);

[ar,xi,kappa,ehat] = lpcana(x,M);       % LP analysis of frame.
ac = lpcweight(ar,c);                   % Coefficients of filter A(z/c).

for (j=1:J)                             % Excitation sequence in blocks.
  n = (j-1)*L+1:j*L;
  [k(j),theta0(j),P(j),b(j),ebuf,Zf,Zw] = celpexcit(x(n),cb,Pidx,ar,ac,...
                                                              bbuf,ebuf,Zf,Zw);
  bbuf = b(j);                          % Last estimated b.
end


text:
function x = p19_start ( option, nvar )

%*****************************************************************************80
%
%% P19_START returns a starting point for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  x = [ ...
      0.00218216, ...
      0.03171126, ...
      0.00010562, ...
      0.48301846, ...
      0.48298250, ...
      0.41554567, ...
      0.14949595, ...
      0.43425476, ...
      0.00018983, ...
      0.00051379, ...
    207.02239583, ...
     22.97760417, ...
      1.00000000 ]';

  return
end


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

%%
clear; close all; clc;

%% Syntax
% |dcmFolderAnon(pathName,varargin);|

%% Description
% This function aims to annonymise the dicom data in the folder pathName.
% See also: |dicomanon|

%% Examples
%

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function R = apprRot(Ra)
%R = apprRot(Ra)

i1 = 0.5; i2 = 0.5;
U = Ra(1,:);
V = Ra(2,:);
un = norm(U);
vn = norm(V);
Un = U/un;
Vn = V/vn;

vp = Un*Vn';
up = Vn*Un';

Vc = Vn-vp*Un;  Vc = Vc/norm(Vc);
Uc = Un-up*Vn;  Uc = Uc/norm(Uc);

Ua = i1*Un+i2*Uc; Ua = Ua/norm(Ua);
Va = i1*Vn+i2*Vc; Va = Va/norm(Va);

R = [Ua;Va;cross(Ua,Va)];
if det(R)<0, R(3,:) = -R(3,:); end;

author: TadasBaltrusaitis

repo: OpenFace


text:
function aux_SeisMoment(params, hParentFigure)
% function aux_SeisMoment(params, hParentFigure);
%-------------------------------------------
% Plots cum. number of events  versus time
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% J.Woessner, jowoe@gps.caltech.edu
% last update: 09.11.05

% Get the axes handle of the plotwindow
axes(sv_result('GetAxesHandle', hParentFigure, [], guidata(hParentFigure)));
hold on;
% Select a point in the plot window with the mouse
[fX, fY] = ginput(1);
disp(['X: ' num2str(fX) ' Y: ' num2str(fY)]);
% Plot a small circle at the chosen place
plot(fX,fY,'ok');

% Get closest gridnode for the chosen point on the map
[fXGridNode fYGridNode,  nNodeGridPoint] = calc_ClosestGridNode(params.mPolygon, fX, fY);
plot(fXGridNode, fYGridNode, '*r');
hold off;

% Get the data for the grid node
mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNodeGridPoint}, :);

if (params.nGriddingMode == 1 && params.bMap == 1)
    vDistances_ = sqrt(((mNodeCatalog_(:,1)-fXGridNode)*cos(pi/180*fYGridNode)*111).^2 + ((mNodeCatalog_(:,2)-fYGridNode)*111).^2);
    vSel = (vDistances_ <= params.fRadius);
    fCheckDist = max(vDistances_(vSel, :));
    mNodeCatalog_=mNodeCatalog_(vSel,:);
elseif (params.nGriddingMode == 1 && params.bMap == 0)
    [nRow_, nColumn_] = size(mNodeCatalog_);
    vXSecX_ = mNodeCatalog_(:,nColumn_);  % length along x-section
    vXSecY_ = (-1) * mNodeCatalog_(:,7);
    vDistances_ = sqrt(((vXSecX_ - fXGridNode)).^2 + ((vXSecY_ - fYGridNode)).^2);
    vSel = (vDistances_ <= params.fRadius);
    mNodeCatalog_=mNodeCatalog_(vSel,:);
    fCheckDist = max(vDistances_);
end
% Time period
params.fTimePeriod = params.fTimePeriod/365;

% Sort the catalog according to time
[s,is] = sort(mNodeCatalog_(:,3));
mNodeCatalog_ = mNodeCatalog_(is(:,1),:);

if exist('new_fig','var') &&  ishandle(new_fig)
    set(0,'Currentfigure',new_fig);
    disp('Figure exists');
else
    new_fig=figure_w_normalized_uicontrolunits('tag','bnew','Name','Cumulative FMD and b-value fit','Units','normalized','Nextplot','add','Numbertitle','off');
    new_axs=axes('tag','ax_bnew','Nextplot','add','box','on');
end
set(gca,'tag','ax_bnew','Nextplot','replace','box','on');
axs5=findobj('tag','ax_bnew');
axes(axs5(1));
vTime = mNodeCatalog_(:,3);
vCumNumber = (1:length(vTime));
vCumNumber(length(vTime)) = length(vTime);
plot(vTime,vCumNumber,'k','LineWidth',2);

% Figure refinements
set(gca,'LineWidth',2,'FontSize',12','FontWeight','bold')

xlabel('Time in years ','FontWeight','bold','FontSize',12)
ylabel('Number of events','FontWeight','bold','FontSize',12)

author: CelsoReyes


text:
function coppelia_wait(coppelia, n)
    for i=1:n
        coppelia.sim.simxSynchronousTrigger(coppelia.clientID);
    end
end

repo: ARTE


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_freq(_time)' and
%               'rpt_pos_freq(_time)'
%   'method'  = 'granger' (default), or 'instantaneous', or 'total'
%   'hasjack' = boolean, specifying whether the input contains leave-one-outs,
%               required for correct variance estimate (default = false)
%   '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 CONNECTIVITY, 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$

% check for a supported dimord
ft_checkopt(varargin, 'dimord', 'char');

method  = ft_getopt(varargin, 'method',  'granger');
hasjack = ft_getopt(varargin, 'hasjack', false);
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: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [resources] = getResourcesStruct(resDblArr, numDbl)
%getResourcesStruct Summary of this function goes here
%   Detailed explanation goes here

    r = readPartStrAndDblsFromDblArr(resDblArr, numDbl);

    resources = struct();
    for(i=1:size(r,1))
        resourcePart = r{i,1};
        split = strsplit(resourcePart, '/_/_/');
        resName = lower(split{1});

        if(~isfield(resources,resName))
            resources.(resName) = cell(0);
        end

        rArrOrig = resources.(resName);
        rArr = [split, r{i,2:end}];
        rArrOrig(end+1,:) = rArr;

        resources.(resName) = rArrOrig;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ gd ] = gsp_design_can_dual( g ,tol)
%GSP_DESIGN_CAN_DUAL This function return the canonical dual filters of g
%   Usage:  gd = gsp_design_can_dual( g );
%
%   Inputs parameters:
%       g       : cell array of filters
%       tol     : tolerance for the pseudo-inverse
%
%   Ouputs parameters:
%       g       : cell array of filters
%
%   This function returns the canonical dual filterbank g. Note that it
%   might not be the be the optimal solution in term of computation.
%
%   Example:::
%
%             N = 100;
%             G = gsp_sensor(N);
%             G = gsp_compute_fourier_basis(G);
%             g = gsp_design_abspline(G,8);
%             gd = gsp_design_can_dual(g);
%             paramplot.show_sum = 0;
%             figure(1)
%             gsp_plot_filter(G,g,paramplot);
%             title('Original filters')
%             figure(2)
%             gsp_plot_filter(G,gd,paramplot);
%             title('Canonical dual filters');
%
%             x = rand(N,1);
%             param.method = 'exact';
%             coeff = gsp_filter_analysis(G,g,x,param);
%             xs = gsp_filter_synthesis(G,gd,coeff,param);
%             norm(xs-x)
%
%   See also: gsp_evaluate_can_dual

% Author: Nathanael Perraudin
% Date  : 30 December 2014
% Testing: test_dual

if nargin<2
    tol = 1e-8;
end

Nf = length(g);
gd = cell(Nf,1);

for ii = 1:Nf
    gd{ii} = @(x) can_dual(g,ii,x,tol);
end

end

function ret = can_dual(g,n,x,tol)
    [N1, N2] = size(x);
    x = x(:);
    sol = gsp_evaluate_can_dual( g,x,tol );
    ret = sol(:,n);
    ret = reshape(ret,N1,N2);
end

repo: gspbox


text:
function pass = test_guide(pref)
% Test Chebfun3 guide commands.

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

% commands from guide1:
f = chebfun3(@(x,y,z) 1./(1+x.^2+y.^2+z.^2));
pass(1) = abs(f(0, 0.5, 0.5) - 2/3) < tol;

pass(2) = abs(sum3(f) - 4.28685406230184188268) < tol;

pass(3) = abs(mean3(f) - (1/8)*4.28685406230184188268) < tol;

pass(4) = abs(max3(f) - 1) < tol;

f1 = chebfun(@(x) exp(x));
len1D = length(f1);
f3 = chebfun3(@(x,y,z) exp(x));
[m3, n3, p3] = length(f3);
pass(5) = m3 < 2* len1D && n3 == 1 && p3 == 1;

f3 = chebfun3(@(x,y,z) exp(y));
[m3, n3, p3] = length(f3);
pass(6) = m3 == 1 && n3 < 2* len1D && p3 == 1;

f3 = chebfun3(@(x,y,z) exp(z));
[m3, n3, p3] = length(f3);
pass(7) = m3 == 1 && n3 == 1 && p3 < 2* len1D;

% Is || f * g || <= || f || * || g || ?
f = chebfun3(@(x,y,z) sin(x+y.*z));
g = chebfun3(@(x,y,z) cos(15*exp(z))./(5+x.^3+2*y.^2+z));
pass(8) = max3(f.*g) <= max3(f) * max3(g);

% Test line integration:
curve = chebfun(@(t) [cos(t) sin(t) t/(8*pi)], [0, 8*pi]);
f = chebfun3(@(x,y,z) x+y.*z);
I = integral(f, curve);
exact = -sqrt(1+(8*pi)^2)/(8*pi);
pass(9) = abs(I - exact) < tol;

% Check if the 'trig' flag needs less coefficients for a triply periodic
% function:
ff = @(x,y,z) tanh(3*sin(x))-(sin(y+1/2)).^2+cos(6*z);
dom = [-pi pi -pi pi -pi pi];
fTrig = chebfun3(ff, dom, 'trig');
[m, n, p] = length(fTrig);
fCheb = chebfun3(ff, dom);
[m_fCheb, n_fCheb, p_fCheb] = length(fCheb);
pass(10) = m <= m_fCheb && n <= n_fCheb && p <= p_fCheb;

%% Test Higher-order SVD:
f = chebfun3(@(x,y,z) sin(x+2*y+3*z));
[sv, Score, Scols, Srows, Stubes] = hosvd(f);
sv1 = sv{1};
sv2 = sv{2};
sv3 = sv{3};
% Test decay of singular values:
pass(11) = sv1(2) <= sv1(2);

pass(12) = sv2(2) <= sv2(2);

pass(13) = sv3(2) <= sv3(2);

% Test orthogonality in columns, rows and tubes:
pass(14) = norm(eye(size(Scols,2)) - Scols'*Scols) < tol;

pass(15) = norm(eye(size(Srows,2)) - Srows'*Srows) < tol;

pass(16) = norm(eye(size(Stubes,2)) - Stubes'*Stubes) < tol;

% Test the _all orthogonality_ property of the core tensor:
pass(17) = norm(squeeze(Score(1,:,:) .* Score(2,:,:))) < tol;

pass(18) = norm(squeeze(Score(:,1,:) .* Score(:,2,:))) < tol;

pass(19) = norm(Score(:,:,1) .* Score(:,:,2)) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_failed=test_frft

disp(' ===============  TEST_FRFT ===========');

Lr=[9,10,11,12];

test_failed=0;

% Test the hermite functions and discrete frft
for ii=1:length(Lr)
	L=Lr(ii);
	F=fft(eye(L))/sqrt(L);

	% check if hermite functions are eigenfunctions of F
	V=hermbasis(L,4);
	res=norm(abs(F*V)-abs(V));
	[test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=fprintf('HERMBASIS L:%3i %0.5g %s\n',L,res,fail);

	% Frft of order 1 becomes ordinary DFT
	f1=tester_crand(L,1);
	f2=tester_crand(1,L);

	p=4;
	frf1=dfracft(f1,1,[],p);
	frf2=dfracft(f2,1,2,p);
	res=norm(F*f1-frf1);
	[test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=fprintf('DFRACFT  L:%3i, %0.5g %s\n',L,res,fail);
	res=norm(f2*F-frf2);
	[test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=fprintf('DFRACFT  L:%3i, %0.5g %s\n',L,res,fail);

	frf1=dfracft(f1,1);
	frf2=dfracft(f2,1,2);
	res=norm(F*f1-frf1);
	[test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=fprintf('DFRACFT  L:%3i %0.5g %s\n',L,res,fail);
	res=norm(f2*F-frf2);
	[test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=fprintf('DFRACFT  L:%3i %0.5g %s\n',L,res,fail);

end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun polylog (@var{s}, @var{z})
%% Numerical polylogarithm function
%%
%% Evaluates the polylogarithm of order @var{s} and argument @var{z},
%% in double precision.  Both inputs can be arrays but their sizes
%% must be either the same or scalar.
%%
%% Example:
%% @example
%% @group
%% polylog (2, -4)
%%   @result{} ans = -2.3699
%% @end group
%% @end example
%%
%% @strong{Note} this function may be slow for large numbers of inputs.
%% This is because it is not a native double-precision implementation
%% but rather the numerical evaluation of the Python @code{mpmath} function
%% @code{polylog}.
%%
%% @seealso{@@sym/polylog}
%% @end defun

function y = polylog (s, x)
  if (nargin ~= 2)
    print_usage ();
  end

  if (isequal (size (s), size (x)) || isscalar(s))
    y = zeros (size (x));
  elseif (isscalar (x))
    y = zeros (size( s));
  else
    error ('polylog: inputs S and X must have compatible sizes')
  end

  cmd = { 'Ls = _ins[0]'
          'Lx = _ins[1]'
          'if len(Ls) == 1 and len(Lx) != 1:'
          '    Ls = Ls*len(Lx)'
          'if len(Ls) != 1 and len(Lx) == 1:'
          '    Lx = Lx*len(Ls)'
          'c = [complex(polylog(s, x)) for s,x in zip(Ls, Lx)]'
          'return c,' };
  c = pycall_sympy__ (cmd, num2cell (s(:)), num2cell (x(:)));
  for i = 1:numel (c)
    y(i) = c{i};
  end
end

%!error polylog (1)
%!error polylog (1, 2, 3)

%!error <sizes> polylog ([1 2], [1 2 3])
%!error <sizes> polylog ([1 2], [1; 2])

%!test
%! y = sym(11)/10;
%! t = sym(2);
%! x = 1.1;
%! s = 2;
%! A = polylog (s, x);
%! B = double (polylog (t, y));
%! assert (A, B, -eps);

%!test
%! % maple
%! A = 2.3201804233130983964 - 3.4513922952232026614*1i;
%! B = polylog (2, 3);
%! assert (A, B, -eps)

%!test
%! % maple, complex inputs
%! A = -11.381456201167411758 + 6.2696695219721651947*1i;
%! B = polylog (1+2i, 3+4i);
%! assert (A, B, -eps);

%!test
%! % maple, matrix inputs
%! A1 = 0.47961557317612748431 - 0.52788287823025778869*1i;
%! A2 = -0.0049750526563452645369 - 0.024579343612396884851*1i;
%! B = polylog ([-1-2i -3], [30+40i 40i]);
%! assert ([A1 A2], B, -eps);

%!test
%! % x matrix, s scalar
%! y = [1 2 sym(pi); exp(sym(1)) 5 6];
%! t = sym(2);
%! x = double (y);
%! s = 2;
%! A = polylog (s, x);
%! B = double (polylog (t, y));
%! assert (A, B, -eps);

%!test
%! % s matrix, x scalar
%! t = [1 2 sym(pi); exp(sym(1)) 5 6];
%! y = sym(2);
%! s = double (t);
%! x = 2;
%! A = polylog (s, x);
%! B = double (polylog (t, y));
%! assert (A, B, -eps);

repo: octsympy


text:
function tilefigs(tile,border)
% <cpp> tile figure windows usage: tilefigs ([nrows ncols],border_in pixels)
% Restriction: maximum of 100 figure windows
% Without arguments, tilefigs will determine the closest N x N grid

%Charles Plum                    Nichols Research Corp.
%<cplum@nichols.com>             70 Westview Street
%Tel: (781) 862-9400             Kilnbrook IV
%Fax: (781) 862-9485             Lexington, MA 02173

maxpos  = get (0,'screensize'); % determine terminal size in pixels
maxpos(4) = maxpos(4) - 25;
hands   = get (0,'Children');   % locate fall open figure handles
hands   = sort(hands);          % sort figure handles
numfigs = size(hands,1);        % number of open figures

maxfigs = 100;

if (numfigs>maxfigs)            % figure limit check
        disp([' More than ' num2str(maxfigs) ' figures ... get serious pal'])
        return
end

if nargin == 0
  maxfactor = sqrt(maxfigs);       % max number of figures per row or column
  sq = [2:maxfactor].^2;           % vector of integer squares
  sq = sq(find(sq>=numfigs));      % determine square grid size
  gridsize = sq(1);                % best grid size
  nrows = sqrt(gridsize);          % figure size screen scale factor
  ncols = nrows;                   % figure size screen scale factor
elseif nargin > 0
  nrows = tile(1);
  ncols = tile(2);
  if numfigs > nrows*ncols
    disp ([' requested tile size too small for ' ...
        num2str(numfigs) ' open figures '])
        return
  end
end
if nargin < 2
  border = 0;
else
  maxpos(3) = maxpos(3) - 2*border;
  maxpos(4) = maxpos(4) - 2*border;
end
xlen = fix(maxpos(3)/ncols) - 30; % new tiled figure width
ylen = fix(maxpos(4)/nrows) - 45; % new tiled figure height

% tile figures by postiion
% Location (1,1) is at bottom left corner
pnum=0;
for iy = 1:nrows
  ypos = maxpos(4) - fix((iy)*maxpos(4)/nrows) + border +25; % figure location (row)
  for ix = 1:ncols
        xpos = fix((ix-1)*maxpos(3)/ncols + 1) + border+7;     % figure location (column)
        pnum = pnum+1;
    if (pnum>numfigs)
                break
        else
          figure(hands(pnum))
      set(hands(pnum),'Position',[ xpos ypos xlen ylen ]); % move figure
        end
  end
end
return

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


text:
function [M1,MW,MW1] = perform_wavelet_matching(M1,M,options)

% perform_wavelet_matching - match multiscale histograms
%
% M1 = perform_wavelet_matching(M1,M,options);
%
%   M1 is the image to synthesize.
%   M is the exemplar image.
%
%   This function match the histogram of the image and the histogram
%   of each sub-band of a wavelet-like pyramid.
%
%   To do texture synthesis, one should apply several time this function.
%   You can do it by setting the value of options.niter_synthesis.
%   This leads to the synthesis as described in
%
%       Pyramid-Based Texture Analysis/Synthesis
%       D. Heeger, J. Bergen,
%       Siggraph 1995
%
%   The transform used for synthesis is options.synthesis_method, which can
%   be either 'steerable' 'wavelets-ortho' 'quincunx-ti' 'wavelets-ti'
%       'wavelets-circle'.
%
%   See also perform_wavelet_transform, perform_histogram_equalization
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
niter_synthesis = getoptions(options, 'niter_synthesis', 1);
verb = getoptions(options, 'verb', 0);

if not(isfield(options, 'color_mode'))
    options.color_mode = 'pca';
end
if isfield(options, 'color_mode') && strcmp(options.color_mode, 'pca') && ~isfield(options, 'ColorP') && size(M,3)==3
    [tmp,options.ColorP] = change_color_mode(M,+1,options);
end
rgb_postmatching = getoptions(options, 'rgb_postmatching', 0);

if size(M,3)==3
    options.niter_synthesis = 1;
    options.verb = 0;
    for iter=1:niter_synthesis
        if verb
            progressbar(iter, niter_synthesis);
        end
        % color images
        M  = change_color_mode(M, +1,options);
        M1 = change_color_mode(M1,+1,options);
        for i=1:size(M,3)
            M1(:,:,i) = perform_wavelet_matching(M1(:,:,i),M(:,:,i), options);
        end
        M  = change_color_mode(M, -1,options);
        M1 = change_color_mode(M1,-1,options);
        if rgb_postmatching
            for i=1:size(M,3)
                M1(:,:,i) = perform_histogram_equalization(M1(:,:,i),M(:,:,i));
            end
        end
    end
    return;
end

if size(M,3)>1
    for i=1:size(M,3)
        [M1(:,:,i),MW,MW1] = perform_wavelet_matching(M1(:,:,i),M(:,:,i),options);
    end
    return;
end

n = size(M,1);
n1 = size(M1,1);

synthesis_method = getoptions(options, 'synthesis_method', 'steerable');

m = 2^( ceil(log2(n)) );
m1 = 2^( ceil(log2(n1)) );
M = perform_image_extension(M,m);
M1 = perform_image_extension(M1,m1);

% precompute input
MW = my_transform(M, m, +1, options);

for iter=1:niter_synthesis
    if verb
        progressbar(iter, niter_synthesis);
    end
    % spatial equalization
    M1 = my_equalization(M1,M);
    % forward transforms
    MW1 = my_transform(M1, m1, +1, options);
    % wavelet domain equalization
    MW1 = my_equalization(MW1,MW);
    % backward transform
    M1 = my_transform(MW1, m1, -1, options);
    % spatial equalization
    M1 = my_equalization(M1,M);
end

M1 = M1(1:n1,1:n1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = my_equalization(M,M0)

options.absval = 0;
options.rows = 0;
options.cols = 0;
options.dim3 = 1;

if iscell(M)
    for i=1:min(length(M),length(M0))
        M{i} = my_equalization(M{i},M0{i});
    end
    return;
end
if size(M,3)>1
    for i=1:min(size(M,3),size(M0,3))
        M(:,:,i) = my_equalization(M(:,:,i),M0(:,:,i));
    end
    return;
end
M = perform_histogram_equalization(M,M0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = my_transform(M, n, dir, options)

synthesis_method = getoptions(options, 'synthesis_method', 'steerable');
deltaJ = 5;
if strcmp(synthesis_method, 'steerable')
    deltaJ = 3;
end
Jmax = log2(n)-1;
Jmin = max(Jmax-deltaJ+1,3);

% steerable options
if not(isfield(options, 'nb_orientations'))
    options.nb_orientations = 4;
end
% wave ortho options
options.wavelet_type = 'biorthogonal_swapped';
options.wavelet_vm = 4;

switch synthesis_method
    case 'steerable'
        M = perform_steerable_transform(M, Jmin, options);
    case 'wavelets-ortho'
        if dir==-1
            M = convert_wavelets2list(M, Jmin);
        end
        M = perform_wavelet_transform(M, Jmin, dir, options);
        if dir==1
            M = convert_wavelets2list(M, Jmin);
        end
    case 'quincunx-ti'
        M = perform_quicunx_wavelet_transform_ti(M,Jmin,options);
    case 'wavelets-ti'
        options.wavelet_type = 'biorthogonal';
        options.wavelet_vm = 3;
        M = perform_atrou_transform(M,Jmin,options);
    case 'wavelets-circle'
        if dir==-1
            M = wavecircle2list(M,Jmin);
        end
        wavelet_modulo = getoptions(options, 'wavelet_modulo', 2*pi);
        M = perform_circle_haar_transform(M, Jmin, dir, wavelet_modulo, options);
        if dir==1
            M = wavecircle2list(M,Jmin);
        end
    otherwise
        error('Unknown transform.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function MW = wavecircle2list(M,Jmin)

if not(iscell(M))
    n = size(M,1); Jmax = log2(n)-1;
    MW = {};
    for j=Jmax:-1:Jmin
        MW{end+1} = M(end/2+1:end,:); M(end/2+1:end,:) = [];
        MW{end+1} = M(:,end/2+1:end); M(:,end/2+1:end) = [];
    end
    MW{end+1} = M;
else
    n = size(M{1},1)*2;
    MW = M{end};
    for i=length(M)-1:-1:1
        if mod(i,2)==0
            MW = [MW, M{i}];
        else
            MW = [MW; M{i}];
        end
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
%ISURF SURF feature extractor
%
% SF = ISURF(IM, OPTIONS) returns a vector of SurfPointFeature objects
% representing scale and rotationally invariant interest points in the
% image IM.
%
% The SurfPointFeature object has many properties including:
%  u            horizontal coordinate
%  v            vertical coordinate
%  strength     feature strength
%  descriptor   feature descriptor (64x1 or 128x1)
%  sigma        feature scale
%  theta        feature orientation [rad]
%
% Options::
% 'nfeat',N      set the number of features to return (default Inf)
% 'thresh',T     set Hessian threshold.  Increasing the threshold reduces
%                the number of features computed and reduces computation time.
% 'octaves',N    number of octaves to process (default 5)
% 'extended'     return 128-element descriptor (default 64)
% 'upright'      don't compute rotation invariance
% 'suppress',R   set the suppression radius (default 0).  Features are not
%                returned if they are within R [pixels] of an earlier (stronger)
%                feature.
% Example::
%
% Load the image
%         im = iread('lena.pgm');
% Find the 10 strongest SURF features
%         sf = isurf(im, 'nfeat', 10);
% and overlay them on the original image as blue circles
%         idisp(im);
%         sf.plot_scale()
%
% Notes::
% - Color images, or sequences, are first converted to greyscale.
% - Features are returned in descending strength order
% - If IM is HxWxN it is considered to be an image sequence and F is a cell
%   array with N elements, each of which is the feature vectors for the
%   corresponding image in the sequence.
% - Wraps an M-file implementation of OpenSurf by D. Kroon (U. Twente) or
%   a MEX-file OpenCV wrapper by Petter Strandmark.
% - The sign of the Laplacian is not retained.
% - The SURF algorithm is covered by an extensive suite of international
%   patents including US 8,165,401, EP 1850270 held by Toyota, KU Leuven
%   and ETHZ.  See http://www.kooaba.com/en/plans_and_pricing/ip_licensing
%
% Reference::
% "SURF: Speeded Up Robust Features",
% Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
% Computer Vision and Image Understanding (CVIU),
% Vol. 110, No. 3, pp. 346--359, 2008
%
%
% See also SurfPointFeature, ISIFT, ICORNER.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function features = isurf(im, varargin)

    opt.suppress = 0;
    opt.nfeat = Inf;
    opt.extended = true;
    opt.upright = false;
    opt.octaves = 5;
    opt.thresh = [];
    opt.id = [];

    [opt,arglist] = tb_optparse(opt, varargin);

    if iscell(im)
        % images provided as a cell array, return a cell array
        % of SURF object vectors
        fprintf('extracting SIFT features for %d greyscale images\n', length(im));
        features = {};
        for i=1:length(im)
            sf = isurf(im{i}, 'setopt', opt, 'id', i);

            features{i} = sf;
            fprintf('.');
        end
        fprintf('\n');
        return
    end

    % convert color image to greyscale
    if iscolor(im)
       im = imono(im);
    end

    if ndims(im) > 2

        % TODO sequence of color images..

        % images provided as an array, return a cell array
        % of SURF object vectors
        if opt.verbose
            fprintf('extracting SURF features for %d greyscale images\n', size(im,3));
        end
        features = {};
        for i=1:size(im,3)
            sf = isurf(im(:,:,i), 'setopt', opt);
            for j=1:length(sf)
                sf(j).image_id_ = i;
            end
            features{i} = sf;
            fprintf('.');
        end
        if opt.verbose
            fprintf('\n');
        end
        return
    end

    % do SURF using a static method that wraps the implementation from:
    % OpenSURF for Matlab
    %   written by D.Kroon University of Twente (July 2010)
    %   based on the C++ implementation by Chris Evans
    %try  HACK
        Ipts = SurfPointFeature.surf(im, opt);
%     catch
%         error('Contributed software [OpenSurf] for SURF features is not installed');
%     end

    % Ipts is a structure array with elements x, y, scale, orientation, descriptor

    try
        % sort into descending order of corner strength
        [z,k] = sort([Ipts.strength], 'descend');
        Ipts = Ipts(:,k);
    catch
        % will fail if no features found
        features = [];
        return;
    end

    fprintf('%d corners found (%.1f%%), ', length(Ipts), ...
        length(Ipts)/prod(size(im))*100);

    % allocate storage for the objects
    n = min(opt.nfeat, length(Ipts));

    features = [];
    i = 1;
    while i<=n
        if i > length(Ipts)
            break;
        end

        % enforce separation between corners
        % TODO: strategy of Brown etal. only keep if 10% greater than all within radius
        if (opt.suppress > 0) && (i>1)
            d = sqrt( ([features.v]'-Ipts(i).y).^2 + ([features.u]'-Ipts(i).x).^2 );
            if min(d) < opt.suppress
                continue;
            end
        end
        f = SurfPointFeature(Ipts(i).x, Ipts(i).y, Ipts(i).strength);
        f.scale_ = Ipts(i).scale;
        f.theta_ = Ipts(i).orientation;
        f.descriptor_ = cast(Ipts(i).descriptor, 'single');
        f.image_id_ = opt.id;

        features = [features f];
        i = i+1;
    end
    fprintf(' %d corner features saved\n', i-1);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [p,t,stats,terms,arg] = rmanova(dat,varargin)
% RMANOVA - performs a conventional n-way analysis of variance (ANOVA)
% or a repeated-measures ANOVA. A conventional ANOVA assumes a
% between-subjects design (different groups), a RM ANOVA a within-subjects
% design (all subjects participated in all conditions).
%
%Usage:
% [p,t,stats,terms,arg] = rmanova(dat,<OPT>)
% [p,t,stats,terms,arg] = rmanova(dat,varnames,<OPT>)
%
%Arguments:
% DAT      -  N*F1*F2*F3*...*FN data matrix.
%             First dimension (rows): refers to the subjects, ie each row
%             contains the data of one subject. The second and successive
%             dimensions contain the
%             factors, with the size of the dimension being the number of
%             levels of that factor. Example: for 13 subjects and the
%             factors Speller with 3 levels (Hex,Cake,Center speller) and
%             Attention with 2 levels (overt, covert), the data matrix
%             would have a size of 13*3*2. The entry dat(5,2,1) refers to
%             the fifth subject for Cake Speller (2nd level of Speller) and
%             overt attention (1st level of Attention).
%
% OPT - struct or property/value list of optional properties:
% 'Varnames'   - CELL ARRAY of one or more factors (eg {'Speller'
%              'Attention'}). The order of the factors must correspond
%              to the order of the factors in the DAT matrix.
% 'Design' -  Test design. 'independent' performs a conventional ANOVA.
%             If 'repeated-measures' (default), performs a
%             repeated-measures ANOVA. In the latter case, Subject is
%             included as a random (ie not fixed) effect.
% 'Display'  - if 'on' displays a table with the results (default 'on')
%
% Assumptions:
%   ANOVA   - homogeneity of variance: variances within each group are
%   equal
%   RM ANOVA - sphericity. Variation of population >difference< scores are
%   the same for all differences (for >2 levels of a factor).
%
% All other options are passed to the 'anovan' function.
%
%Returns:
% [p,t,stats,terms]     - 'help anovan' for details
% arg                   - arguments passed to anovan
%
% See also ANOVAN, PLOT_STATS.
%
% Note: If your experiment involves repeated-measures (your subject was run in all
% subconditions) you should use repeated-measures ANOVA (RM-ANOVA), because the
% assumption of independence of samples is violated. Furthermore, RM-ANOVA
% accounts for inter-subject variability and thus has more statistical
% power.

% Author(s): Matthias Treder 2011

varnames = [];
if nargin==2
  varnames = varargin{1};
  varargin={};
elseif  nargin>1 && iscell(varargin{1})
  varnames = varargin{1};
  varargin = varargin(2:end);
end

props = {'Varnames',             varnames,          'CELL{CHAR}|CHAR';
         'Design',          'independent',          '!CHAR(repeated-measures independent)';
         'Display',                'on',            '!CHAR(on off)';
         'Alpha'                    .05,            '!DOUBLE[1]';
         'Model'                    'full',         '!CHAR(linear interaction full)|DOUBLE';
         'SSType'                   3,              '!CHAR(h)|INT';
         'Table'                    1,              '!BOOL';
         'Verbose'                  1,              '!BOOL';
         };

if nargin==0,
  p= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt,isdefault] = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);
misc_checkType(dat,'!DOUBLE');

if ischar(opt.Varnames)
    opt.Varnames = {opt.Varnames};
end

ss = size(dat);
nsbj = ss(1);           % Number of subjects
level = ss(2:end);

%% Factors and levels
% Number of factors
if ~isempty(opt.Varnames)
  nfac = numel(opt.Varnames);
  if nfac ~= ndims(dat)-1
    error('Number of factors %d does not match factors in data %d',nfac,ndims(dat)-1)
  end
else
  nfac = ndims(dat)-1;
end

% If no factor names provided, use default names
if isempty(opt.Varnames)
  vn = num2cell(1:nfac);
  vn = cellfun(@num2str,vn,'UniformOutput',0);
  opt.Varnames = strcat('X',vn);
end

%% Check (RM) ANOVA assumptions // TODO

if strcmp(opt.Design,'independent')
  % Homogenity of variance
  % TODO
elseif strcmp(opt.Design,'repeated-measures')
  % Sphericity
end

%% ANOVA or rmANOVA
if opt.Verbose
  names = cell2mat(strcat(opt.Varnames,',')); names=names(1:end-1);
  fprintf('Performing a %d-way %s ANOVA with factors {%s} and %s = %d levels.\n',...
    nfac,opt.Design,names,str_vec2str(level,'%d','x'),prod(level))
end

if strcmp(opt.Design,'independent')
  % ANOVA model
  design = orthogonal_design(nsbj,level);

elseif strcmp(opt.Design,'repeated-measures')
  % Repeated measures ANOVA
  opt.Varnames = {'Subject' opt.Varnames{:}};  % Add subject as a factor
  if isfield(opt,'random')   % Add subject as random effect
    opt.Random = [1 opt.Random+1];
  else
    opt.Random = 1;
  end
  design = orthogonal_design(1,[nsbj level]);
  % Specify to-be-tested interactions by hand to omit Subject
  model = double(flipud(dec2bin(0:2^(nfac+1)-1))-'0');
  model( model(:,1) & sum(model(:,2:end),2)>0 , :) = []; % Omit all interactions involving Subject
  maineffects = find(sum(model,2)==1);
  interactions = find(sum(model,2)>1);
  opt.Model = model([maineffects; interactions],:); % Bring in right order
else
  error('Unknown design ''%s''',opt.Design)
end

%% Perform ANOVA
arg = opt_structToProplist(rmfield(opt,{'Design','Table','Verbose','Alpha'}));

[p,t,stats,terms] = anovan(dat(:),design.anova,arg{:});

%% Provide output
if opt.Verbose
  fprintf('-------\nResults\n-------\n')
  % Find col indices
  F = find(ismember(t(1,:),'F','legacy'));
  p = find(ismember(t(1,:),'Prob>F','legacy'));
  for ii=2:size(terms,1)+1
    if t{ii,p}<opt.Alpha
      fprintf('''%s'' significant, F = %1.2f, p = %0.4f\n',...
        t{ii,1},t{ii,F},t{ii,p})
    else
      fprintf('''%s'' not significant, F = %1.2f, p = %0.4f\n',...
        t{ii,1},t{ii,F},t{ii,p})
    end
  end
end

function d = orthogonal_design(nRows,nLevels)
% ORTHOGONAL_DESIGN - help function to create matrices with factor values
% assuming an orthogonal design (ie, there is data for each possible
% pairing of subconditions).
% Example: If you have the factors target/nontarget, and electrode
% (Fz, Cz, Pz), and you measure the effect of these variables on 12
% subjects, you would get a 12 x 6 (=2*3) matrix wherein the rows represent
% subjects, and the columns are:
% (1)target-Fz (2)nontarget-Fz (3)target-Cz (4)nontarget-Cz (5)target-Pz (6)nontarget-Pz
% The respective function call would be orthogonal_design(12,[2 3])
%
% Synopsis:
%   D = ORTHOGONAL_DESIGN(NROWS,NLEVELS)
%
% Arguments:
%   NROWS  : number of measurements of each factor (eg number of subjects)
%   NLEVELS: a vector with each field specifying the number of levels of
%     each factor
%
% Returns:
%   A struct with the following fields
%   .mat: a cell array with the corresponding matrices for each variable
%   (except the first)
%   .anova: vector notation of the factor matrices which you can feed
%   directly into anovan()
nFactors = numel(nLevels);
nCols = prod(nLevels);
d = struct('mat',zeros(nRows,nCols,nFactors),'anova',[],'anova_mat',[]);

% Construct first row and then copy n times
for jj=1:nFactors
  repLev = nCols / prod(nLevels(jj:end)); % Number of repeats of a levels

  if jj<nFactors   % Number of repeats of the sequence
    repSeq = prod(nLevels(jj+1:end));
  else repSeq=1;
  end

  row = [];
  for rr=1:repSeq
    for nn=1:nLevels(jj)
      row = [row repmat(nn,[1 repLev])];
    end
  end
  d.mat(1,:,jj)=row;
  %% Extend to all rows
  d.mat(:,:,jj) = repmat(row,[nRows 1]);
end

for jj=1:nFactors
  dd = d.mat(:,:,jj);
  d.anova{jj} = dd(:);
  d.anova_mat = [d.anova_mat dd(:)];
end
end

end

author: bbci


text:
function z = psdinvjmul(xlab,xfrm, y, K) %#ok
% z = psdinvjmul(xlab,xfrm, y, K)
%
% PSDINVJMUL  solves x jmul z = y, with x = XFRM*diag(xlab)*XFRM'
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function nocancel = cinchGenerateData (dt6file, datapath)
nocancel = true;
if ~exist(datapath, 'dir')
        r = questdlg ('CINCH data has not yet been created for this subject. Do this now? (Warning: This could take a minute or two!)');
        if (~strcmp (r,'Yes'))
            nocancel = false;
            return;
        end;
    h = mrvWaitbar(0,'Please wait...');
    dtiConvertDT6ToBinaries (dt6file, datapath);
    close(h);
end;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function h = plotUnitCells(ebsd,d,varargin)
% low level plotting routine for EBSD maps
%

unitCell = ebsd.unitCell;

xy = [ebsd.prop.x(:),ebsd.prop.y(:)];
if check_option(varargin,'region')

  reg = get_option(varargin,'region');

  ind = xy(:,1) > reg(1) & xy(:,1) < reg(2) & xy(:,2) > reg(3) & xy(:,2) < reg(4);

  xy = xy(ind,:);

  if numel(d) == numel(ind) || numel(ind) == size(d,1)
    d = d(ind,:);
  end

end

ax = get_option(varargin,'parent',gca);

if ~isempty(unitCell)

  type = get_flag(varargin,{'unitcell','points','measurements'},'unitcell');

else

  type = 'points';

end

if numel(d) == size(xy,1) || numel(d) == 3*size(xy,1)

  obj.FaceVertexCData = reshape(d,size(xy,1),[]);

  if check_option(varargin,{'transparent','translucent','faceAlpha'})

    s = get_option(varargin,{'transparent','translucent','faceAlpha'},1,'double');
    varargin = delete_option(varargin,'faceAlpha');

    if size(d,2) == 3 % rgb
      obj.FaceVertexAlphaData = s.*(1-min(d,[],2));
    else
      obj.FaceVertexAlphaData = s.*d./max(d);
    end
    obj.AlphaDataMapping = 'none';
    obj.FaceAlpha = 'flat';
  end
  obj.FaceColor = 'flat';
else
  obj.FaceColor = d;
end

obj.EdgeColor = 'none';

switch lower(type)
  case 'unitcell'

    % generate patches
    [obj.Vertices, obj.Faces] = generateUnitCells(xy,unitCell,varargin{:});

  case {'points','measurements'}

    obj.Vertices = xy;
    obj.Faces    = (1:size(xy,1))';

    obj.FaceColor = 'none';
    obj.EdgeColor = 'flat';
    obj.Marker = '.';
    obj.MarkerFaceColor = 'flat';

end

h = optiondraw(patch(obj,'parent',ax),varargin{:});

if ~check_option(varargin,'DisplayName')
  set(get(get(h,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end

if nargout == 0, clear h;end

author: mtex-toolbox

repo: mtex

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


text:
function YESNONA = isconvex(p)
%ISCONVEX Tries to determine if a scalar expression is convex
%
% T = isconvex(p)
%
% p: scalar SDPVAR object
% T: The result (1: convex, 0: concave, NaN: cannot be determined)
%
% NOTE : Under development. Do not trust if you have monomials in your model
%
% Example
% sdpvar x y
% isconvex(x+y) will return 1
% isconvex(x+y^2) will return 1
% isconvex(exp(x+y)) will return 1
% isconvex(-exp(x+y)) will return 0
% isconvex(max(x,exp(x+y))) will return 1
% isconvex(-max(x,exp(x+y))) will return 0
% isconvex(max(x,min(x,-x))) will return NaN

YESNONA = NaN;
[F,failure,cause] = expandmodel([],p,sdpsettings('allownonconvex',0,'allowmilp',0));
if failure == 0
    YESNONA = true;
else
    [F,failure,cause] = expandmodel([],-p,sdpsettings('allownonconvex',0));
    if failure == 0
        % p is nonconvex
        YESNONA = false;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function title = p05_title ( )

%*****************************************************************************80
%
%% P05_TITLE returns a title for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = '#5: The horn.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [rf, vf] = twobody2 (mu, tau, ri, vi)

% solve the two body initial value problem

% Shepperd's method

% input

%  mu  = gravitational constant (km**3/sec**2)
%  tau = propagation time interval (seconds)
%  ri  = initial eci position vector (kilometers)
%  vi  = initial eci velocity vector (kilometers/second)

% output

%  rf = final eci position vector (kilometers)
%  vf = final eci velocity vector (kilometers/second)

% Orbital Mechanics with MATLAB

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

tolerance = 1.0d-12;

u = 0;

imax = 20;

umax = +realmax;

umin = -realmax;

orbits = 0;

tdesired = tau;

threshold = tolerance * abs (tdesired);

r0 = norm(ri);

n0 = dot(ri, vi);

beta = 2 * (mu / r0) - dot(vi, vi);

if (beta ~= 0)
    umax = +1 / sqrt(abs(beta));

    umin = -1 / sqrt(abs(beta));
end

if (beta > 0)
    orbits = beta * tau - 2 * n0;

    orbits = 1 + (orbits * sqrt(beta)) / (pi * mu);

    orbits = floor (orbits / 2);
end

for i = 1:1:imax

    q = beta * u * u;

    q = q / (1 + q);

    n =  0;

    r = 1;

    l = 1;

    s = 1;

    d = 3;

    gcf = 1;

    k = -5;

    gold = 0;

    while (gcf ~= gold)

        k = -k;

        l = l + 2;

        d = d + 4 * l;

        n = n + (1 + k) * l;

        r = d / (d - n * r * q);

        s = (r - 1) * s;

        gold = gcf;

        gcf  = gold + s;

    end

    h0 = 1 - 2 * q;

    h1 = 2 * u * (1 - q);

    u0 = 2 * h0 * h0 - 1;

    u1 = 2 * h0 * h1;

    u2 = 2 * h1 * h1;

    u3 = 2 * h1 * u2 * gcf / 3;

    if (orbits ~= 0)
        u3 = u3 + 2 * pi * orbits / (beta * sqrt(beta));
    end

    r1 = r0 * u0 + n0 * u1 + mu * u2;

    dt = r0 * u1 + n0 * u2 + mu * u3;

    slope = 4 * r1 / (1 + beta * u * u);

    terror = tdesired - dt;

    if (abs (terror) < threshold)
        break;
    end;

    if ((i > 1) && (u  ==  uold) )
        break;
    end

    if ((i > 1) && (dt == dtold) )
        break;
    end

    uold  = u;

    dtold = dt;

    ustep = terror / slope;

    if (ustep > 0)
        umin = u;

        u = u + ustep;

        if (u > umax)
            u = (umin + umax) / 2;
        end
    else
        umax = u;

        u = u + ustep;

        if (u < umin)
            u = (umin + umax) / 2;
        end
    end

    if (i == imax)
        fprintf('\n\nmax iterations in twobody2 function');
    end

end

usaved = u;

f = 1.0 - (mu / r0) * u2;

gg = 1.0 - (mu / r1) * u2;

g  =  r0 * u1 + n0 * u2;

ff = -mu * u1 / (r0 * r1);

% final position and velocity vectors

for i = 1:1:3
    rf(i) = f  * ri(i) + g  * vi(i);

    vf(i) = ff * ri(i) + gg * vi(i);
end

repo: mcbench-benchmarks