text:
function stats = glmfit_multilevel(Y, X1, X2, varargin)
% :Usage:
% ::
%
%     stats = glmfit_multilevel(Y, X1, X2, varargin)
%
% Mixed-effects models differ in their assumptions and implementation
% details. glmfit_multilevel is a fast and simple option for running a
% two-level mixed effects model with participant as a random effect.
% It implements a random-intercept, random-slope model across 2nd-level units
% (e.g., participants). it fits regressions for individual 2nd-level
% units (e.g., participants), and then (optionally) uses a precision-weighted least
% squares approach to model group effects. It thus treats participants as a
% random effect. This is appropriate when 2nd-level units are participants
% and 1st-level units are observations (e.g., trials) within participants.
% glmfit_multilevel was designed with this use case in mind.
%
% Options:
% glmfit_multilevel includes some options that are not included in many
% mixed effects models, including:
% - bootstrapping or sign permutation for inference
% - robust regression (*needs code update*)
% - AR(p) autoregressive model (*needs code update*)
%
% Requirements: glmfit_multilevel requires enough 1st-level units to fit a
% separate model for each  2nd-level unit (participant). If this is not the
% case, other models (igls.m, LMER, etc.) are preferred.
%
% Degrees of freedom: glmfit_multilevel is conservative in the sense that the degrees of
% freedom in the group statistical test is always based on the number of
% subjects - 2nd-level parameters.
% The df is never higher than the sample size, which you would have with
% mixed effects models that estimate the df from the data. This causes
% problems in many other packages, particularly when there are many 1st-level
% observations and they are uncorrelated, resulting in large and undesirable
% estimated df.
%
% Correlated effects: The correlations across 1st-level observations are not measured, and
% 1st-level obs are assumed to be IID. This is valid when generalizing across
% 2nd-level units, but may not be fully efficient (powerful) if 1st-level
% units are correlated.
%
% :Inputs:
%
%   **Y:**
%        Is data in either:
%           -cell array, one cell per subject.
%            Column vector of subject outcome data in each cell.
%           -Matrix
%            One column per subject, with vector of subject outcome
%            data in that column
%
%   **X1 and X2:**
%        are first and 2nd level design matrices
%          - X1 in cell array, one cell per subject
%            design matrix for each subject in each cell.
%            *columns must code for the same variable for all subjects*
%
%         - X2 in rect. matrix
%         - can be empty (intercept only)
%
% E.g., with one 2nd-level predictor:
% stats = glmfit_multilevel(Y, X1, X2, ...
% 'names', {'Int' 'Temp'}, 'beta_names', {'2nd-level Intercept (overall group effect)' '2nd-lev predictor: Group membership'});
%
% :Output:
%
%   **stats:**
%        is structure with results.
%        Intercept is always added as first column!
%        (do not add intercept to input predictors)
%
% See glmfit_general.m for varargin variable input options.
%
% :Examples:
% ::
%
%    len = 200; sub = 20;
%    x = zeros(len,sub);
%    x(11:20,:) = 2;                   % create signal
%    x(111:120,:) = 2;
%    c = normrnd(0.5,0.1,sub,1);       % slope between-subjects variations
%    d = normrnd(3,0.2,sub,1);         % intercept between-subjects variations
%    % Create y: Add between-subjects error (random effects) and measurement noise
%    % (within-subjects error)
%    for i=1:sub, y(:,i) = d(i) + c(i).*x(:,i) + normrnd(0,0.5,len,1);
%    end;
%
%    for i = 1:size(y, 2), YY{i} = y(:, i); end
%    for i = 1:size(y, 2), XX{i} = x(:, i); end
%
%    % one-sample t-test, weighted by inv of btwn + within vars
%    stats = glmfit_multilevel(YY, XX, [], 'verbose', 'weighted');
%
%    statsg = glmfit_multilevel(y, x, covti, 'names', {'L1 Intercept' 'L1 Slope'},...
%    'beta_names', {'Group Average', 'L2_Covt'});
%
% :Input Options:
%
% General Defaults
%  - case 'names', Names of first-level predictors, starting
%    with 'Intercept', in cell array
%
%  - case 'analysisname', analysisname = varargin{i+1}; varargin{i+1} = [];
%  - case 'beta_names', beta_names = Names of 2nd-level predictors, starting
%    with 'Intercept', in cell array
%
% Estimation Defaults
%  - case 'robust', robust_option = 'yes';
%  - case {'weight', 'weighted', 'var', 's2'}, weight_option = 'unweighted';
%  - case {'nocenter'}, do not force centering of 2nd-level predictors
%
% Inference defaults
%  - case {'boot1', 'boot', 'bootstrap'}, inference_option = 'bootstrap';
%  - case {'sign perm', 'signperm', 'sign'}, inference_option = 'signperm';
%  - case {'t-test', 'ttest'}, inference_option = 't-test';
%
% Display control defaults
%  - case 'plots', doplots = 1; plotstr = 'plots';
%  - case 'noplots', doplots = 0; plotstr = 'noplots';
%
%  - case {'dosave', 'save', 'saveplots'}, dosave = 1; savestr = 'save';
%  - case 'verbose', verbose = 1; verbstr = 'verbose';
%  - case 'noverbose', verbose = 0; verbstr = 'noverbose';
%
%  - case {'savefile', 'savefilename'}, savefilename = varargin{i + 1}; varargin{i+1} = [];
%
% Bootstrap defaults
%  - case 'nresample', nresample = varargin{i+1};
%  - case {'pvals', 'whpvals_for_boot'}, whpvals_for_boot = varargin{i+1};
%
% Sign perm defaults
%  - case {'permsign'}, permsign = varargin{i+1};
%

%    Programmer's notes:
%    9/2/09: Tor and Lauren: Edited to drop NaNs within-subject, and drop
%    subject only if there are too few observations to estimate.
% ..

 % Convert from matrix form to cells
 % Matrix: Y is a col vector, X is one predictor column per subject
 % -------------------------------------------------------------------

  if ~iscell(Y)
    N = size(Y, 2);
    for i = 1:N
      YY{i} = Y(:, i);
    end
    Y = YY;
    clear YY;
  end

  if ~iscell(X1)
    N2 = size(X1, 2);
    if N ~= N2, error('Sizes of X and Y do not match'); end
    for i = 1:N
      XX{i} = X1(:, i);
    end
    X1 = XX;
    clear XX
  end

  N = length(Y);
  if N ~= length(X1)
    error('Enter one cell per subject for each of X and Y');
  end

 % first level: SETUP
 % -------------------------------------------------------------------

 % Check variances and exclude subjects with no variance in any Y or X
 [Y, X1, X2] = check_variances_and_exclude(Y, X1, X2);
 N = length(Y);

 % set up first-level X matrix (sample)
  if any(strcmp(varargin, 'noint')) % no-intercept version
    X1tmp = X1{1};
  else
    X1tmp = setup_X_matrix(X1{1}); % intercept first
  end

  k = size(X1tmp, 2);             % num predictors; assumed to be the same!!

 % Second level: SETUP
 % Need to remove 2nd-level units with NaN data at first level
 % -------------------------------------------------------------------
  wh_omit = false(1, N);
  for i = 1:N
     %if any(isnan(Y{i})) || any(isnan(X1{i}(:))), wh_omit(i) = 1; end

    can_be_nans = length(Y{i}) - k - 1;  % up to this many can be NaN, still leaving 1 degree of freedom
    if can_be_nans < 0, warning('Warning: you might be overparameterized!  Seems like you have more predictors than observations'); end
    if sum(isnan(Y{i}) | any(isnan(X1{i}), 2)) > can_be_nans
      wh_omit(i) = 1;
    end
  end

  if any(wh_omit)
    if isempty(X2), X2 = ones(N, 1); end

    Y(wh_omit) = [];
    X1(wh_omit) = [];
    X2(wh_omit, :) = [];
    N = length(Y);
  end

  beta = zeros(k, N);
  sterr = zeros(k, N);
  t = zeros(k, N);
  p = zeros(k, N);
  dfe = zeros(1, N);
%phi = zeros(arorder, N);

 % first level: ESTIMATE
 % -------------------------------------------------------------------
  for i = 1:N
    [beta(:, i), sterr(:, i), t(:, i), p(:, i), dfe(:, i), phi(:,i), V{i}] = first_level_model(Y{i}, X1{i}, varargin{:});
       % V{i} is var/cov matrix (xtxi)*sigmasq
  end

  varnames = {'beta' 't' 'p' 'dfe' 'phi'};
  first_level = create_struct(varnames);
  first_level.ste = sterr;

 % second level: Finish SETUP and ESTIMATE
 % -------------------------------------------------------------------
 % set up second-level X matrix: intercept first
  X2 = setup_X_matrix(X2, beta(1,:)');

% set up second-level options
% names of outcomes become beta_names here b/c second level test on 1st
% level betas
  [beta_names1, analysisname, beta_names2, robust_option, weight_option, inference_option, ...
   verbose, dosave, doplots, ...
   verbstr, savestr, plotstr, ...
   targetu, nresample, whpvals_for_boot, ...
   permsign] = ...
  setup_inputs(beta', X2, varargin{:});

  switch weight_option
    case 'weighted'
  % Note: R & B-style : replaced sterr' with V

      stats = glmfit_general( ...
	      beta', X2, ...
	      'analysisname', analysisname, 'names', beta_names1, 'beta_names', beta_names2, ...
	      verbstr, savestr, plotstr, ...
	      weight_option, V, inference_option, 'dfwithin', dfe', ...
	      'targetu', targetu, 'nresample', nresample, ...
	      'whpvals_for_boot', whpvals_for_boot, 'permsign', permsign);

    case 'unweighted'

      stats = glmfit_general( ...
	      beta', X2, ...
	      'analysisname', analysisname, 'names', beta_names1, 'beta_names', beta_names2, ...
	      verbstr, savestr, plotstr, ...
	      weight_option, inference_option, ...
	      'targetu', targetu, 'nresample', nresample, ...
	      'whpvals_for_boot', whpvals_for_boot, 'permsign', permsign);

    otherwise
      error('Problem with weight_option. Please select either weighted or unweighted.')
  end

  stats.first_level = first_level;

  if doplots
    scn_stats_helper_functions('xyplot', X1, Y, weight_option, 'names', beta_names1, 'nostats');
    xlabel('X'); ylabel('Y');
  end

% _________________________________________________________________________
%
%
%
% * Inline functions
%
%
%
%__________________________________________________________________________

  function newstruct = create_struct(varnames)
    newstruct = struct();
    for i = 1:length(varnames)
      eval(['newstruct.' varnames{i} ' = ' varnames{i} ';']);
    end
  end

end %End of glmfit_multilevel

function [b, sterr, t, p, dfe, phi, V] = first_level_model(y, X, varargin)

% defaults
% -------------------------------------------------------------------
verbose = 0;
verbstr = 'noverbose';
arorder = 0;                    % or Zero for no AR
interceptstr = 'intercept';

% optional inputs
% -------------------------------------------------------------------
for varg = 1:length(varargin)
    if ischar(varargin{varg})
        switch varargin{varg}

            % reserved keywords
            case 'verbose all', verbose = 1; verbstr = 'verbose';
            case 'verbose', % do nothing
            case {'ar', 'arorder'} , arorder = varargin{varg+1};
                %otherwise, disp(['Unknown input string option: ' varargin{varg}]);

            case 'noint', interceptstr = 'noint';
        end
    end
end

k = size(X, 2) + 1;

[whnan X y] = nanremove(X, y);

if isempty(X)
    % no data

    [b, t, p, sterr] = deal(NaN * ones(k, 1));
    dfe = deal(NaN);
    if arorder, phi = NaN * ones(arorder, 1); else phi = NaN; end
    V = [];

    return

end

% set up X matrix: intercept first
if ~strcmp(interceptstr, 'noint')
    X = setup_X_matrix(X, y);
end

if arorder
    % if we have missing observations or redundant columns, let's
    % regularize a bit so we can still estimate this, using a ridge prior
    % The degree of regularization is arbitrary.
    % V is used to estimate variance components and re-weight.
    if rank(X) < size(X, 2)
        disp('WARNING! RANK DEFICIENT.  THIS FUNCTION WILL RETURN AN ERROR.')
        X = [X; eye(size(X, 2)) ./ size(X, 1)];
        if ~strcmp(interceptstr, 'noint'), X(:, 1) = 1; end
        y = [y; ones(size(X, 2), 1) .* nanmean(y)];
    end

    [t, dfe, b, phi, sigma, sterr] = fit_gls(y, X, [], arorder);
    p = 2 * (1 - tcdf(abs(t), dfe));  % two-tailed

    V = inv(X' * X) * sigma .^ 2; % Var/Cov mtx, Precision^-1, used in weighted est. and empirical bayes

else
    % if we have missing observations or redundant columns, let's
    % regularize a bit so we can still estimate this, using a ridge prior
    % The degree of regularization is arbitrary.
    if rank(X) < size(X, 2)
        disp('WARNING! RANK DEFICIENT.  THIS FUNCTION WILL RETURN AN ERROR.')
        X = [X; eye(size(X, 2)) ./ size(X, 1)];
        if ~strcmp(interceptstr, 'noint'), X(:, 1) = 1; end
        y = [y; ones(size(X, 2), 1) .* nanmean(y)];
    end

    stats = glmfit_general(y, X, verbstr, interceptstr);
    t = stats.t; dfe = stats.dfe; b = stats.beta; phi = NaN; sterr = stats.ste; p = stats.p;

    V = inv(X' * X) * stats.var; % Var/Cov mtx, Precision^-1, used in weighted est. and empirical bayes
end

end

function X = setup_X_matrix(X, y)
  % set up X matrix: intercept first
  [n, k] = size(X);
  if n == 0
    n = size(y, 1);
  end
  equal_x = false(1, k);
  for i = 1:k
    if all(X(:, i) == X(1, i))
    % weights are equal
      equal_x(i) = 1;
    end
  end
  if any(equal_x)
    error('Warning: some columns of X have no variance.  Do not enter intercept in X; it will be added automatically as the first predictor.');
    X(:, equal_x) = [];
  end
  X = [ones(n, 1) X];
end

% -------------------------------------------------------------------------
% Setup inputs, print info to screen if verbose
% THIS IS FOR SECOND-LEVEL MODEL -- NOT FIRST
% -------------------------------------------------------------------------
function [names, analysisname, beta_names, robust_option, weight_option, inference_option, ...
    verbose, dosave, doplots, ...
    verbstr, savestr, plotstr, ...
    targetu, nresample, whpvals_for_boot, ...
    permsign] = ...
    setup_inputs(Y, X, varargin)

% Initial compliance checks
% ----------------------------------------------------------------
[n, k] = size(X);
nvars = size(Y, 2);
nobs_tmp = size(Y, 1);

% Check sizes
if nobs_tmp ~= n, error('Y and X must have same number of rows!'); end

% Check intercept
if ~(all(X(:, 1) == X(1)))
    error('First column of X must be an intercept column. (e.g., all ones)');
end

beta_names = cell(1, k);
for i = 1:k
    beta_names{i} = sprintf('2nd-level B%02d', i);
end

% Defaults
% ----------------------------------------------------------------

% General Defaults
names = cell(1, nvars);     % variable names, columns of Y
for i = 1:nvars, names{i} = ['1st-level B' num2str(i)]; end

analysisname = 'Second Level of Multilevel Model ';

% Estimation Defaults
robust_option = 'no';               % robust IRLS; 'no' 'yes'
weight_option = 'unweighted';       % 'weighted' 'unweighted'
force_centering = true;             % force 2nd-level predictor centering

% Inference defaults
inference_option = 't-test';        % 't-test' 'bootstrap' 'signperm'

% Display control defaults
verbose = 1;                % verbose output
dosave = 0;                 % save figures at end
doplots = 0;                % make plots
plotstr = 'noplots';
savestr = 'nosave';
verbstr = 'verbose';

savefilename = 'glmfit_general_output.txt';

% Bootstrap defaults
targetu = .20;              % proportion contribution of boot procedure to p-value
nresample = 1000;         % initial bootstrap samples
whpvals_for_boot = 1:size(Y,2);   % indices of p-values, the min of which is used to determine boot samples needed
% lower p-values require more boot samples
% for the p-vals to be meaningful.

% Sign perm defaults
permsign = [];              % empty: setup new sign permutation indices
% if entered: keep same permutation matrix
% across repeated calls (much faster! but
% reduces accuracy in simulations!)

% Inputs
% ----------------------------------------------------------------
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % General Defaults
            case 'names', names = varargin{i+1};  varargin{i+1} = [];
            case 'analysisname', analysisname = varargin{i+1}; varargin{i+1} = [];
            case 'beta_names', beta_names = varargin{i+1}; varargin{i+1} = [];

                % Estimation Defaults
            case 'robust', robust_option = 'yes';
            case {'weight', 'weighted', 'var', 's2'}, weight_option = 'weighted';
            case {'nocenter', 'nocentering'}, force_centering = false;

                % Inference defaults
            case {'boot1', 'boot', 'bootstrap'}, inference_option = 'bootstrap';
            case {'sign perm', 'signperm', 'sign'}, inference_option = 'signperm';
            case {'t-test', 'ttest'}, inference_option = 't-test';

                % Display control defaults
            case {'plot', 'plots'}, doplots = 1; plotstr = 'plots';
            case {'noplots', 'noplot'}, doplots = 0; plotstr = 'noplots';

            case {'dosave', 'save', 'saveplots'}, dosave = 1; savestr = 'save';
            case 'verbose', verbose = 1; verbstr = 'verbose';
            case 'noverbose', verbose = 0; verbstr = 'noverbose';

            case {'savefile', 'savefilename'}, savefilename = varargin{i + 1}; varargin{i+1} = [];

                % Bootstrap defaults
            case 'nresample', nresample = varargin{i+1};
            case {'pvals', 'whpvals_for_boot'}, whpvals_for_boot = varargin{i+1};

                % Sign perm defaults
            case {'permsign'}, permsign = varargin{i+1};

            case 'intercept'

            otherwise
                fprintf('Warning! Unknown input string option: %s', varargin{i});

        end
    end
end

% all the other checking, etc. is done in glmfit_general

is_intercept = all(abs(diff(X, 1, 1)) < 100 * eps);

if force_centering
    X(:, ~is_intercept) = X(:, ~is_intercept) - mean(X(:, ~is_intercept));
end

is_centered = abs(mean(X)) < 100*eps;

any_noncentered = any(~is_intercept & ~is_centered);
intercept_beta_name = 'Within-person average effects';

if any_noncentered
    disp('WARNING!!! Some 2nd-level predictors are not mean-centered.')
    disp('Within-person effects and P-values are NOT interpretable as the')
    disp('average within-person effect');

    intercept_beta_name = 'Within-person effects when all 2nd-level predictors are zero';

end

% fix names by adding intercept if needed
% This is over-rided by user input
if length(beta_names) == k - 1
    if ~isrow(beta_names), beta_names = beta_names'; end
    beta_names = {intercept_beta_name beta_names{:}};
end

end

% -------------------------------------------------------------------------
% Check variances
% -------------------------------------------------------------------------

function [Y, X1, X2] = check_variances_and_exclude(Y, X1, X2)
xvar = cellfun(@var, X1, 'UniformOutput', false);
xvar = cat(1, xvar{:});
wh = any(xvar < 100 * eps, 2);
if any(wh)
    disp('Warning! Some participants have no variance in X column(s) and will be excluded')
    fprintf('Participant numbers:');
    fprintf('%d ', find(wh));
    fprintf('\n');
end

wh_out = wh;

yvar = cellfun(@var, Y, 'UniformOutput', false);
yvar = cat(1, yvar{:});
wh = any(yvar < 100 * eps, 2);
if any(wh)
    disp('Warning! Some participants have no variance in Y and will be excluded')
    fprintf('Participant numbers:');
    fprintf('%d ', find(wh));
    fprintf('\n');
end

wh_out = wh_out | wh;

X1(wh_out) = [];
Y(wh_out) = [];
if ~isempty(X2), X2(wh_out,  :) = []; end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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

%%
clear; close all; clc;

%% Syntax
% |[C]=grayColor(varargin);|

%% Description
% UNDOCUMENTED
%% 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 Obs = matchFeature(Sen,Raw,Obs)

% MATCHFEATURE  Match feature.
% 	Obs = MATCHFEATURE(Sen,Raw,Obs) matches one feature in Raw to the predicted
% 	feature in Obs.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

switch Obs.ltype(4:6)
    case 'Pnt'
        rawDataLmks = Raw.data.points;
        R = Sen.par.cov;
    case 'Lin'
        rawDataLmks = Raw.data.segments;
        R = blkdiag(Sen.par.cov,Sen.par.cov);
    otherwise
        error('??? Unknown landmark type ''%s''.',Obs.ltype);
end

switch Raw.type

    case {'simu','dump'}

        id  = Obs.lid;
        idx = find(rawDataLmks.app==id);

        if ~isempty(idx)
            Obs.meas.y   = rawDataLmks.coord(:,idx);
            Obs.meas.R   = R;
            Obs.measured = true;
            Obs.matched  = true;
        else
            Obs.meas.y   = zeros(size(Obs.meas.y));
            Obs.meas.R   = R;
            Obs.measured = false;
            Obs.matched  = false;
        end

    case 'image'
        error('??? Feature matching for Raw data type ''%s'' not implemented yet.', Raw.type)

    otherwise

        error('??? Unknown Raw data type ''%s''.',Raw.type)

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function points3d(varargin)
%POINTS3D Description of functions operating on 3D points.
%
%   Points are represented by their 3 Cartesian coordinates:
%   P = [X Y Z];
%
%   Arrays of points consist in N*3 arrays, each row being a point.
%
%   See also
%   isCoplanar, distancePoints, boundingBox3d
%   anglePoints3d, angleSort3d, sphericalAngle
%   sph2cart2, cart2sph2, cart2cyl, cyl2cart
%   transformPoint3d, clipPoints3d
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function x = daub14_transform_inverse ( n, y )

%*****************************************************************************80
%
%% DAUB14_TRANSFORM_INVERSE inverts the DAUB14 transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%    N must be a power of 2 and at least 4.
%
%    Input, real Y(N), the transformed vector.
%
%    Output, real X(N), the original vector.
%
  c = [ ...
     7.785205408500917e-02; ...
     3.965393194819173e-01; ...
     7.291320908462351e-01; ...
     4.697822874051931e-01; ...
    -1.439060039285649e-01; ...
    -2.240361849938749e-01; ...
     7.130921926683026e-02; ...
     8.061260915108307e-02; ...
    -3.802993693501441e-02; ...
    -1.657454163066688e-02; ...
     1.255099855609984e-02; ...
     4.295779729213665e-04; ...
    -1.801640704047490e-03; ...
     3.537137999745202e-04 ];
  p = 13;
  x(1:n,1) = y(1:n);
  m = 4;
  q = floor ( ( p - 1 ) / 2 );

  while ( m <= n )

    z(1:m,1) = 0.0;

    j = 1;

    mh = floor ( m / 2 );

    for i = - q + 1 : mh - q

      for k = 0 : 2 : p - 1
        i0 = i4_wrap ( i      + k / 2,     1,      mh );
        i1 = i4_wrap ( i + mh + k / 2,     mh + 1, m  );
        z(j,1)   = z(j,1)   + c(p-k)   * x(i0) + c(k+2) * x(i1);
        z(j+1,1) = z(j+1,1) + c(p-k+1) * x(i0) - c(k+1) * x(i1);
      end

      j = j + 2;

    end

    x(1:m,1) = z(1:m);

    m = m * 2;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function SO3VF = power(SO3VF,a)
% overloads |SO3VF.^2|
%
% Syntax
%   SO3VF = SO3VF.^2
%
% Input
%  SO3VF - @SO3VectorField
%  a - double
%
% Output
%  SO3VF - @SO3VectorField
%

if isa(a,'vector3d')
  a = a.xyz;
end

if ~isnumeric(a)
  error('The exponent has to be numeric.')
end

SO3VF = SO3VectorFieldHandle(@(rot) SO3VF.eval(rot).^(a(:)),SO3VF.SRight,SO3VF.SLeft);

end

author: mtex-toolbox

repo: mtex

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


text:
function q = conjexp(p)

% conjexp
%
%   1/p+1/q=1

q = p/(p-1);
if p==1
    q = Inf;
end
if p==Inf
    q = p;
end

end

author: optimaltransport


text:
function get_seed_test ( )

%*****************************************************************************80
%
%% GET_SEED_TEST tests GET_SEED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GET_SEED_TEST\n' );
  fprintf ( 1, '  GET_SEED picks an initial seed value for R8_UNIFORM_01.\n' );
  fprintf ( 1, '  The value chosen should vary over time, because\n' );
  fprintf ( 1, '  the seed is based on reading the clock.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This is just the "calendar" clock, which does\n' );
  fprintf ( 1, '  not change very fast, so calling GET_SEED several\n' );
  fprintf ( 1, '  times in a row may result in the same value.\n' );

  seed = 12345678;
  seed_old = seed;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial seed is %12d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Next 3 values of R8_UNIFORM_01:\n' );
  fprintf ( 1, '\n' );

  for j = 1 : 3
    [ r, seed ] = r8_uniform_01 ( seed );
    fprintf ( 1, '  %f\n', r );
  end

  for i = 1 : 4

    while ( 1 )

      seed = get_seed ( );

      if ( seed ~= seed_old )
        seed_old = seed;
        break;
      end

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  New seed from GET_SEED is %12d\n', seed );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Next 3 values of R8_UNIFORM_01:\n' );
    fprintf ( 1, '\n' );

    for j = 1 : 3
      [ r, seed ] = r8_uniform_01 ( seed );
      fprintf ( 1, '  %f\n', r );
    end

  end

  return
end


text:
% Stability domain of second order extrapolated BDF scheme

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory is given in Section 5.8 of
% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program generates Fig. 5.7  in the book

figure(1), clf, hold on

		% Plot of stability domain
theta = 0:0.01:0.5; theta = theta*pi;
a = (6*cos(theta)-9)./(2*cos(2*theta) - 4*cos(theta)) - 3/2;
b = (sin(theta).*(2-cos(theta)))./(cos(2*theta) - 2*cos(theta));
z = - a - i*b; plot(z)

		% Plot of oval
a = 1; b =(2*a/4)^(1/4); theta = 0:0.002:1; theta = theta*pi;
z =  - a*(1-cos(theta)) + i*b*(sin(theta)).^0.5; plot(z,'--')

		% Plot of parabola
yy = 0:0.02:1.8; b = 1; par = -(yy/b).^2; plot(par,yy,'.')

repo: Aero-matlab


text:
%% OMAS GRAPHICS (OMAS_graphics.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This utility provides tools for interacting with 3D models within OMAS
% based on the 'stlTools' library. This function does in no way take credit
% for the design of the functions inclosed, only for their arrangment as
% part of the OMAS package.

% Author: James A. Douthwaite 08/03/18

classdef OMAS_graphics
    %% GENERAL GEOMETRY MANIPULATION
    methods (Static, Access = public)
        % SORT POINT CLOUD IN A CLOCK WISE ABOUT A VECTOR
        function [verticesCW] = sortVerticesCW(planarNormal,centroid,vertexData)
            % This function takes a list of vertex data and sorts them into
            % angularly in relation to the planar normal vector.

            % INPUT HANDLING
            if numel(centroid) == 2
                centroid = [centroid;0];
            end
            % THE REFERENCE
            unit_reference = [1;0;0];

            % DEFINE THE VECTORS RELATIVE TO THE CENTROID
            vertexVectors = vertexData - centroid';

            tempMatrix = zeros(size(vertexVectors,1),2);
            for i = 1:size(vertexVectors,1)
                unit_v = vertexVectors(i,:)/norm(vertexVectors(i,:));
                vertexNormal = cross(unit_reference,unit_v);
                % THE ROTATION DIRECTION
                tempMatrix(i,1) = -sign(dot(vertexNormal,planarNormal));
                % THE ROTATION MAGNITUDE
                tempMatrix(i,2) = acos(dot(unit_reference,unit_v));
                % MAP TO CCW ROTATION
                if tempMatrix(i,1) < 0
                    tempMatrix(i,2) = 2*pi - tempMatrix(i,2) ;
                end
            end
            % SORT THE VERTICE BY ANGLE RELATIVE
            [~,sortedIndices] = sortrows(tempMatrix(:,2),1);
            % GET THE SORTED VERTEX ARRAY
            verticesCW = vertexData(sortedIndices,:);
        end
        % SORT POINT CLOUD IN A COUNTER-CLOCK WISE ABOUT A VECTOR
        function [verticesCCW] = sortVerticesCCW(planarNormal,centroid,vertexData)
            % This function takes a list of vertex data and sorts them into
            % angularly in relation to the planar normal vector.

            % INPUT HANDLING
            if numel(centroid) < 3
                centroid = [centroid;0];
            end
            % THE REFERENCE
            unit_reference = [1;0;0];

            % DEFINE THE VECTORS RELATIVE TO THE CENTROID
            vertexVectors = vertexData - centroid';

            tempMatrix = zeros(size(vertexVectors,1),2);
            for i = 1:size(vertexVectors,1)
                unit_v = vertexVectors(i,:)/norm(vertexVectors(i,:));
                vertexNormal = cross(unit_reference,unit_v);
                % THE ROTATION DIRECTION
                tempMatrix(i,1) = sign(dot(vertexNormal,planarNormal));
                % THE ROTATION MAGNITUDE
                tempMatrix(i,2) = acos(dot(unit_reference,unit_v));
                % MAP TO CCW ROTATION
                if tempMatrix(i,1) < 0
                    tempMatrix(i,2) = 2*pi - tempMatrix(i,2) ;
                end
            end
            % SORT THE VERTICE BY ANGLE RELATIVE
            [~,sortedIndices] = sortrows(tempMatrix(:,2),1);
            % GET THE SORTED VERTEX ARRAY
            verticesCCW = vertexData(sortedIndices,:);

%             % DEBUG PLOTS
%             figure(3)
%             hold on; axis equal; grid on;
%             scale = 10;
%             q = quiver3(centroid(1),centroid(2),centroid(3),...
%                         scale*unit_reference(1),scale*unit_reference(2),scale*unit_reference(3),'r');
%             q = quiver3(centroid(1),centroid(2),centroid(3),...
%                         scale*planarNormal(1),scale*planarNormal(2),scale*planarNormal(3),'g');
%             for i = 1:size(verticesCCW,1)
%                 plot3(verticesCCW(i,1),verticesCCW(i,2),verticesCCW(i,3),'ro');
%             end
        end
        % PLANAR PROHECTION OF A 3D GEOMETRY
        function [planarVertices] = geometryPlanarProjection(planarNormal,centroid,vertexData)
            % This function returns a vertex projections on a plane defined
            % by the normal vector provided

            % INPUT HANDLING
            assert(size(planarNormal,1) == size(centroid,1),'The planar normal and centroid must be [3x1]');
            assert(size(planarNormal,2) == 1,'The planar normal and centroid must be [3x1]');
            assert(size(vertexData,2) == 3,'The geometry data must be [:x3].');

            % DEFINE THE VECTORS BETWEEN THE CENTROID AND VERTICES
            vertexData = vertexData - centroid';
            planarVertices = zeros(size(vertexData,1),3);
            for i = 1:size(vertexData,1)
                % The vertex planar projections
                [projection] = OMAS_geometry.vectorPlanarProjection(...
                                    planarNormal,...
                                    vertexData(i,:)');
                % REFORMAT THE PROJECTIONS
                planarVertices(i,:) = projection' + centroid';
            end
            % REMOVE DUPLICATES
            [planarVertices,~] =  unique(planarVertices, 'rows');
            % We cannot gaurantee the geometry will be the same one the
            % unique command is ran. The face data may reference vertices
            % that have been moved or deleted.
        end
        % DEFINE SCALE OF GEOMETRY STRUCTURE
        function [geometry] = scale(geometry,scale)
            % This function scales a geometry in accordance to a provided
            % scalar or vector of dimensional scalars. A geometry is
            % defined as a vertices, faces structure.

            assert(numel(scale) == 1 || numel(scale) == 3,'Please provide a valid scale value.');
            % Apply the scaling either by dimension or unilaterally.
            if numel(scale) == 3
                geometry.vertices = geometry.vertices*diag(scale);
            else
                geometry.vertices = geometry.vertices*(eye(3)*scale);
            end
        end
        % NORMALISE THE STL FILE
        function [geometry] = normalise(geometry)
            % This function normalises the STL geometry to allow it to be
            % scaled appropriately.
            % GEOMETRY CHECK
            assert(size(geometry.vertices,2) == 3 && size(geometry.faces,2) == 3,...
                'The provided patch object has invalid vertex and faces assignments.');

            % GET THE VERTEX MAGNITUDES RELATIVE TO THE ORIGIN
            verticesNorms = abs(geometry.vertices);
            % MAXIMUMS IN THE FIRST AND SECOND DIMENSIONS
            dimMaximal = max(max(abs(verticesNorms),[],1),[],2);
            % SCALE GEOMETRY SPECIFIED BY VIRTUAL.radius
            geometry.vertices = (geometry.vertices/dimMaximal);
        end
        % REMOVE FACE/VERTEX DUPLICATES
        function [geometry] = removeDuplicateVertices(geometry)
            % PATCHSLIM removes duplicate vertices in surface meshes.
            %
            % This function finds and removes duplicate vertices.
            %
            % USAGE: [v, f]=patchslim(v, f)
            %
            % Where v is the vertex list and f is the face list specifying vertex
            % connectivity.
            %
            % v contains the vertices for all triangles [3*n x 3].
            % f contains the vertex lists defining each triangle face [n x 3].
            %
            % This will reduce the size of typical v matrix by about a factor of 6.
            %
            % For more information see:
            %  http://www.esmonde-white.com/home/diversions/matlab-program-for-loading-stl-files
            %
            % Francis Esmonde-White, May 2010

            if ~isfield(geometry,'vertices')
                error('The geometry does not have a specified "vertex" list.');
            end
            if ~isfield(geometry,'faces')
                error('The geometry does not have a specified "face" triangulation list.');
            end

            % REMOVE DUPLICATES
            [geometry.vertices,~,indexn] =  unique(geometry.vertices, 'rows');
            geometry.faces = indexn(geometry.faces);
        end
        % CALCULATE GEOMETRY SURFACE NORMALS
        function [normals]  = normals(geometry)
            % This function computes the surface normals of a defined
            % geometry structure
            normals = zeros(size(geometry.faces));
            for face = 1:size(geometry.faces,1)
                % MEMBERS IF THE PLANE
                memberID_A = geometry.faces(face,1);
                memberID_B = geometry.faces(face,2);
                memberID_C = geometry.faces(face,3);
                % SURFACE DEFINING VECTORS
                BA = geometry.vertices(memberID_B,:) - geometry.vertices(memberID_A,:);
                CA = geometry.vertices(memberID_C,:) - geometry.vertices(memberID_A,:);
                % THE NORMAL
                normals(face,:) = cross(BA,CA);
                normals(face,:) = normals(face,:)/norm(normals(face,:));
            end
        end
    end

    %% UNIVERSAL DRAWING MECHANISMS
    methods (Static)
        % GET HIT-BOX GEOMETRY
        function [hitBoxGeometry] = getHitBoxGeometry(VIRTUAL,geometry)
            % Input check
            assert(isa(VIRTUAL.hitBoxType,'uint8'),'Object hit-box type must be of type "uint8".');
            % Derive hit-box geometry
            switch VIRTUAL.hitBoxType
                case OMAS_hitBoxType.none
                    % Assemble the spherical constraint volume
                    hitBoxGeometry = [];
                case OMAS_hitBoxType.spherical
                    % Assemble the spherical constraint volume
                    hitBoxGeometry = OMAS_graphics.defineSphere(zeros(3,1),VIRTUAL.radius,10);
                case OMAS_hitBoxType.AABB
                    % Assemble the AABB constraint volume
                    R = OMAS_geometry.quaternionToRotationMatrix(VIRTUAL.quaternion);
                    % Rotate the body before defining constraint volume
                    geometry.vertices = R*geometry.vertices;
                    % DEFINE AN ALIGNED AABB CUBOID
                    hitBoxGeometry = OMAS_graphics.defineCuboid(min(geometry.vertices),...
                                                                  max(geometry.vertices));
                case OMAS_hitBoxType.OBB
                    % DEFINE AN ALIGNED OBB CUBOID
                    hitBoxGeometry = OMAS_graphics.defineCuboid(min(geometry.vertices),...
                                                                  max(geometry.vertices));
                    % Assemble the OBB constraint volume
                    R = OMAS_geometry.quaternionToRotationMatrix(VIRTUAL.quaternion);
                    % Rotate the hit-box to be aligned with the geometry
                    hitBoxGeometry.vertices = R*hitBoxGeometry.vertices;
                otherwise
                    error('Hit-box type not recognised');
            end
        end
        % DEFINE CUBOID (MINOR) FROM RADIUS
        function [geometry,minExtents,maxExtents] = defineCuboidFromRadius(center,radius)
            % This function creates a set of vertices for a cube
            % encapsuated in a sphere of a given radius.
            assert(numel(center) == 3,'The cuboid center must be a cartesian vector.');
            assert(numel(radius) == 1,'The radius of the sphere must be a scalar and non-zero.');
            % RATE OF DEMENSIONAL EXPANSION
            h = radius/1.7321;
            % DEFINE THE CUBOID EXTENTS
            minExtents = center - h;
            maxExtents = center + h;
            % DEFINE THE CUBOID VERTICES
            [geometry] = OMAS_graphics.defineCuboid(minExtents,maxExtents);
        end
        % DEFINE CUBOID
        function [geometry,minExtents,maxExtents] = defineCuboid(minExtents,maxExtents)
            % Return a matrix of point defining a cuboid scaled to that of
            % a dimensions provided.

            % Define vertex data from limits
            vertices = zeros(8,3);
            vertices(1,:) = [maxExtents(1),maxExtents(2),maxExtents(3)];
            vertices(2,:) = [maxExtents(1),maxExtents(2),minExtents(3)];
            vertices(3,:) = [maxExtents(1),minExtents(2),minExtents(3)];
            vertices(4,:) = [maxExtents(1),minExtents(2),maxExtents(3)];
            vertices(5,:) = [minExtents(1),minExtents(2),minExtents(3)];
            vertices(6,:) = [minExtents(1),minExtents(2),maxExtents(3)];
            vertices(7,:) = [minExtents(1),maxExtents(2),maxExtents(3)];
            vertices(8,:) = [minExtents(1),maxExtents(2),minExtents(3)];
            geometry.vertices = vertices;
            % Define face connectivity matrix
            geometry.faces =  [  1     2     7
                                 1     4     2
                                 1     7     4
                                 2     3     8
                                 2     4     3
                                 2     8     7
                                 3     4     6
                                 3     5     8
                                 3     6     5
                                 4     7     6
                                 5     6     8
                                 6     7     8];
        end
        % DRAW SPHERE
        function [geometry] = defineSphere(position,radius,faces)
            % This function returns a sphere coordinate cloud at a given
            % position in space, of a given radius. 'figureHandle' is used
            % to provide the function context.

            % INPUT CHECKING
            if nargin < 3
                faces = 10;
            end
            % DEFINE SPHERE TRIANGULATION
            [X,Y,Z] = sphere(faces);
            X = X.*radius + position(1);
            Y = Y.*radius + position(2);
            Z = Z.*radius + position(3);
            % CONVERT TO PATCH OBJECT
            [geometry.faces,geometry.vertices,~] = surf2patch(X,Y,Z,'triangles');
            % DEFINE CENTROID
            geometry.centroid = position'; %reducepatch(P, R)
        end
        % DRAW UNIT TRIAD
        function [triadHandle] = drawTriad(figureHandle,position,R,scale)
            % Draw a unit triad at a cartesian position, rotated by R.
            % FigureHandl Here is used to bring the current figure to the
            % function context.

            % INPUT CHECL
            if nargin < 4
                scale = 1;
            end

            colourVector = 'rgb';
            triadVectors = scale*eye(3);
            for axis = 1:size(triadVectors,2)
                triadVectors(:,axis) = R*triadVectors(:,axis);
                % Draw vectors
                triadHandle(axis) = quiver3(gca,position(1),position(2),position(3),triadVectors(1,axis),triadVectors(2,axis),triadVectors(3,axis),colourVector(axis));
            end
        end
    end

    %% STL IMPORTATION
    methods (Static, Access = public)
        % PREPARE THE PATCH FILES FOR VISUALS
        function [geometry,successFlag] = importStlFromFile(filename)
            % This function prepares object STL files for presentation. The assumption
            % is that object stl file has the same name as the object being simulated.

            try
                % GET THE STL FILE
                geometry = OMAS_graphics.stlRead(filename);
                successFlag = 1;
            catch
                % NO STL WAS FOUND BY THAT FILE NAME (fail quietly)
                geometry = [];
                successFlag = 0;
            end
            % ENSURE UNIQUE VERTICES
            if successFlag
                geometry = OMAS_graphics.removeDuplicateVertices(geometry);
            end
        end
        % READ AND STL INTO VERTICIES, FACES
        function [geometry, name] = stlRead(fileName)
            %STLREAD reads any STL file not depending on its format
            %V are the vertices
            %F are the faces
            %N are the normals
            %NAME is the name of the STL object (NOT the name of the STL file)

            [format,isSuccessful] = OMAS_graphics.stlGetFormat(fileName);

            if strcmp(format,'ascii')
                [geometry,name] = OMAS_graphics.stlReadAscii(fileName);
            elseif strcmp(format,'binary')
                [geometry,name] = OMAS_graphics.stlReadBinary(fileName);
            end
        end
        % IDENTIFY STL TYPE
        function [format,isSuccessful] = stlGetFormat(fileName)
            %STLGETFORMAT identifies the format of the STL file and returns 'binary' or
            %'ascii'

            fid = fopen(fileName);
            if fid == -1        % Unsuccessful file load
                isSuccessful = 0;
            else
                isSuccessful = 1;
            end

            % Check the file size first, since binary files MUST have a size of 84+(50*n)
            fseek(fid,0,1);         % Go to the end of the file
            fidSIZE = ftell(fid);   % Check the size of the file
            if rem(fidSIZE-84,50) > 0
                format = 'ascii';
            else
                % Files with a size of 84+(50*n), might be either ascii or binary...
                % Read first 80 characters of the file.
                % For an ASCII file, the data should begin immediately (give or take a few
                % blank lines or spaces) and the first word must be 'solid'.
                % For a binary file, the first 80 characters contains the header.
                % It is bad practice to begin the header of a binary file with the word
                % 'solid', so it can be used to identify whether the file is ASCII or
                % binary.
                fseek(fid,0,-1);                                           % go to the beginning of the file
                header = strtrim(char(fread(fid,80,'uchar')'));            % trim leading and trailing spaces
                isSolid = strcmp(header(1:min(5,length(header))),'solid'); % take first 5 char
                fseek(fid,-80,1);                                          % go to the end of the file minus 80 characters
                tail = char(fread(fid,80,'uchar')');
                isEndSolid = findstr(tail,'endsolid');

                % Double check by reading the last 80 characters of the file.
                % For an ASCII file, the data should end (give or take a few
                % blank lines or spaces) with 'endsolid <object_name>'.
                % If the last 80 characters contains the word 'endsolid' then this
                % confirms that the file is indeed ASCII.
                if isSolid && isEndSolid
                    format = 'ascii';
                else
                    format = 'binary';
                end
            end
            fclose(fid);
        end
        % INTERPRET A BINARY STL FILE
        function [geometry, name] = stlReadBinary(fileName)
            %STLREADBINARY reads a STL file written in BINARY format
            %V are the vertices
            %F are the faces
            %N are the normals
            %NAME is the name of the STL object (NOT the name of the STL file)

            %=======================
            % STL binary file format
            %=======================
            % Binary STL files have an 84 byte header followed by 50-byte records, each
            % describing a single facet of the mesh.  Technically each facet could be
            % any 2D shape, but that would screw up the 50-byte-per-facet structure, so
            % in practice only triangular facets are used.  The present code ONLY works
            % for meshes composed of triangular facets.
            %
            % HEADER:
            % 80 bytes:  Header text
            % 4 bytes:   (int) The number of facets in the STL mesh
            %
            % DATA:
            % 4 bytes:  (float) normal x
            % 4 bytes:  (float) normal y
            % 4 bytes:  (float) normal z
            % 4 bytes:  (float) vertex1 x
            % 4 bytes:  (float) vertex1 y
            % 4 bytes:  (float) vertex1 z
            % 4 bytes:  (float) vertex2 x
            % 4 bytes:  (float) vertex2 y
            % 4 bytes:  (float) vertex2 z
            % 4 bytes:  (float) vertex3 x
            % 4 bytes:  (float) vertex3 y
            % 4 bytes:  (float) vertex3 z
            % 2 bytes:  Padding to make the data for each facet 50-bytes in length
            %   ...and repeat for next facet...

            fid = fopen(fileName);
            header = fread(fid,80,'int8'); % reading header's 80 bytes
            name = deblank(native2unicode(header,'ascii')');
            if isempty(name)
                name = 'Unnamed Object'; % no object name in binary files!
            end
            nfaces = fread(fid,1,'int32');  % reading the number of facets in the stl file (next 4 byters)
            nvert = 3*nfaces; % number of vertices
            % reserve memory for vectors (increase the processing speed)
            n = zeros(nfaces,3);
            v = zeros(nvert,3);
            f = zeros(nfaces,3);
            for i = 1 : nfaces % read the data for each facet
                tmp = fread(fid,3*4,'float'); % read coordinates
                n(i,:) = tmp(1:3); % x,y,z components of the facet's normal vector
                v(3*i-2,:) = tmp(4:6); % x,y,z coordinates of vertex 1
                v(3*i-1,:) = tmp(7:9); % x,y,z coordinates of vertex 2
                v(3*i,:) = tmp(10:12); % x,y,z coordinates of vertex 3
                f(i,:) = [3*i-2 3*i-1 3*i]; % face
                fread(fid,1,'int16'); % Move to the start of the next facet (2 bytes of padding)
            end
            fclose(fid);
            % Geometric structure
            geometry = struct('vertices',v,'faces',f,'normals',n);
            [geometry] = OMAS_graphics.removeDuplicateVertices(geometry);
        end
        % INTERPRET AN ASCII STL FILE
        function [geometry, n, name] = stlReadAscii(fileName)
            %STLREADASCII reads a STL file written in ASCII format
            %V are the vertices
            %F are the faces
            %N are the normals
            %NAME is the name of the STL object (NOT the name of the STL file)

            %======================
            % STL ascii file format
            %======================
            % ASCII STL files have the following structure.  Technically each facet
            % could be any 2D shape, but in practice only triangular facets tend to be
            % used.  The present code ONLY works for meshes composed of triangular
            % facets.
            %
            % solid object_name
            % facet normal x y z
            %   outer loop
            %     vertex x y z
            %     vertex x y z
            %     vertex x y z
            %   endloop
            % endfacet
            %
            % <Repeat for all facets...>
            %
            % endsolid object_name

            fid = fopen(fileName);
            cellcontent = textscan(fid,'%s','delimiter','\n');             % read all the file and put content in cells
            content = cellcontent{:}(logical(~strcmp(cellcontent{:},''))); % remove all blank lines
            fclose(fid);

            % read the STL name
            line1 = char(content(1));
            if (size(line1,2) >= 7)
                name = line1(7:end);
            else
                name = 'Unnamed Object';
            end

            % read the vector normals
            normals = char(content(logical(strncmp(content,'facet normal',12))));
            n = str2double(normals(:,13:end));

            % read the vertex coordinates (vertices)
            vertices = char(content(logical(strncmp(content,'vertex',6))));
            v = str2double(vertices(:,7:end));
            nvert = size(vertices,1);                       % number of vertices
            nfaces = sum(strcmp(content,'endfacet'));       % number of faces
            if (nvert == 3*nfaces)
                f = reshape(1:nvert,[3 nfaces])';           % create faces
            end
            % Geometric structure
            geometry = struct('vertices',v,'faces',f,'normals',n);
            % slim the file (delete duplicated vertices)
            [geometry] = OMAS_graphics.removeDuplicateVertices(geometry);
        end
        % PLOT THE STL
        function stlPlot(geometry, name)
            %STLPLOT is an easy way to plot an STL object
            %V is the Nx3 array of vertices
            %F is the Mx3 array of faces
            %NAME is the name of the object, that will be displayed as a title

            figure;
            patch(geometry,...
                'FaceColor',[0.8 0.8 1.0], ...
                'EdgeColor','none',...
                'FaceLighting','gouraud',...
                'AmbientStrength', 0.15);

            % Add a camera light, and tone down the specular highlighting
            camlight('headlight');
            material('dull');

            % Fix the axes scaling, and set a nice view angle
            axis('image');
            view([-135 35]);
            grid on;
            title(name);
        end
    end
end

repo: OpenMAS


text:
function c = u2nonucfmt(cu, p)
%U2NONUCFMT Uniform to non-uniform filterbank coefficient format
%   Usage:  c=u2nonucfmt(cu,pk)
%
%   Input parameters:
%         cu   : Uniform filterbank coefficients.
%
%   Output parameters:
%         c    : Non-uniform filterbank coefficients.
%         p    : Numbers of copies of each filter.
%
%   `c = u2nonucfmt(cu,pk)` changes the coefficient format from
%   uniform filterbank coefficients *cu* (`M=sum(p)` channels) to
%   non-uniform coefficients *c* (`numel(p)` channels)  such that each
%   channel of *c* consinst of `p(m)` interleaved channels of *cu*.
%
%   The output *c* is a cell-array in any case.
%
%   See also: nonu2ufilterbank
%
%   References: akkva2003

complainif_notenoughargs(nargin,2,mfilename);

if isempty(cu)
   error('%s: cu must be non-empty.',upper(mfilename));
end

if iscell(cu)
    if any(cellfun(@isempty,cu));
      error('%s: Elements of cu must be non-empty.',upper(mfilename));
    end

    M = numel(cu);
    W = size(cu{1},2);

    Lc = size(cu{1},1);
    if any(Lc~=cellfun(@(cEl)size(cEl,1),cu))
        error('%s: Coefficient subbands do not have an uniform length',...
              upper(mfilename));
    end
elseif isnumeric(cu)
    M = size(cu,2);
    W = size(cu,3);
    Lc = size(cu,1);
else
    error('%s: cu must be a cell array or numeric.',upper(mfilename));
end

if isempty(p) || ~isvector(p)
   error('%s: p must be a non-empty vector.',upper(mfilename));
end

if sum(p) ~= M
    error(['%s: Total number of filters in p does not comply with ',...
           'number of channels'],upper(mfilename));
end

Mnonu = numel(p);
c = cell(Mnonu,1);
p = p(:);
pkcumsum = cumsum([1;p]);
crange = arrayfun(@(pEl,pcEl)pcEl:pcEl+pEl-1,p,pkcumsum(1:end-1),'UniformOutput',0);

if iscell(cu)
    for m=1:Mnonu
        ctmp = [cu{crange{m}}].';
        c{m} = reshape(ctmp(:),W,[]).';
    end
else
    for m=1:Mnonu
        c{m} = zeros(p(m)*Lc,W,assert_classname(cu));
        for w=1:W
            c{m}(:,w) = reshape(cu(:,crange{m},w).',1,[]);
        end
    end
end

% Post check whether there is the same number of coefficients
if sum(cellfun(@(cEl) size(cEl,1),c)) ~= M*Lc
    error(['%s: Invalid number of coefficients in subbands.'],upper(mfilename));
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [AR, lambdamax]=sim_genRndVARcoeffs(varargin)
% Generate random coefficients for a VAR model of given order and with a
% specified degree of interaction sparsity.
% Coefficients are drawn from a normal distribution with specified standard
% deviation (sigma).
%
% The function will repeatedly generate random models until it obtains a
% stable VAR process (all eigenvalues of the system matrix are less than or
% equal to 1 in magnitude).
%
%
% See Also: sim_genTVARcoeffs()
%
% References:
% [1] This function is adapted from Stefan Haufe's gen_ar_sech.m
%
% Author: Tim Mullen 2011-2013, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

g = arg_define(varargin,...
        arg({'nchs','NumChans'},10,[],'Number of channels (variables)'), ...
        arg({'morder','ModelOrder'},10,[],'Model order'), ...
        arg({'sigma','Sigma'},0.2,[0 Inf],'Scale (std) of random VAR parameters. This determines the variance of the determinstic component of the generated process.'), ...
        arg({'maxinter','MaxNumInter'},0.5,[],'Maximum number of interactions. If >= 1, this is the expected number of interacting variables (chosen randomly). If < 1, this represents the proportion of interacting variable (out of nchs^2-nchs max possible interactions). This can also be a vector specifying the linear (column-major) index of each nonzero (non-diagonal) element of the interaction (connectivity) matrix'), ...
        arg({'probinter','PerLagInteractProb'},1,[0 1],'Percent of nonzero AR coefficients for a given interaction. This applies to auto-interactions as well as cross-interactions. For instance, a value of 1 will generate nonzero gaussian iid AR filter coefficients. A value of 0.5 will generate a mixed sparse AR filter with 50% zero coefficients and remaining nonzero iid gaussian coefficients'), ...
        arg({'maxtrys','MaxAttempts','maxattempts'},Inf,[],'Max number of retries for stable model gen.'), ...
        arg({'maxlambda','MaxLambdaMag'},1,[0 1],'Largest allowable eigenvalue (magnitude) of system matrix. This determines the stability of the process. The smaller this value the more stable the process. All eigenvalues must be <= 1 in magnitude for a stable process.') ...
        );

% RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));

if g.maxinter < 1
    g.maxinter = ceil(g.maxinter*(g.nchs^2-g.nchs));
end

% generate canonical interaction indices
if length(g.maxinter) == 1
    inddiag     = linspace(1, g.nchs^2, g.nchs);
    indndiag    = setdiff(1:g.nchs^2, inddiag);
    per         = randperm(length(indndiag));
    indndiag    = indndiag(per);
    indndiag    = indndiag(1:g.maxinter);
    ind         = [inddiag indndiag];
else
    inddiag     = linspace(1, g.nchs^2, g.nchs);
    ind         = unique([inddiag g.maxinter]);
end
ind = sort(ind);
% generate interaction indices for all model orders
ind = tensorsum(0:g.nchs^2:g.nchs^2*(g.morder-1),ind);

% generate random coefficients for interacting variables
lambdamax = Inf; ntrys = 0;
while lambdamax > g.maxlambda && ntrys < g.maxtrys
    AR      = zeros(g.nchs,g.nchs*g.morder);
    AR(ind) = double(rand(length(ind), 1) <= g.probinter).*randn(length(ind), 1)*g.sigma;

%     for k=1:g.morder
%       aloc = zeros(g.nchs);
%       aloc(ind) = double(rand(length(ind), 1) < g.probinter).*randn(length(ind), 1)*g.sigma;
%       AR(:,:,(k-1)*g.nchs+1:k*g.nchs) = aloc;
%     end

    % check model stability
    E         = eye(g.nchs*g.morder);
    AA        = [AR;E(1:end-g.nchs,:)];
    lambda    = eig(AA);
    lambdamax = max(abs(lambda));

    ntrys = ntrys + 1;
end
if ntrys>=g.maxtrys
    fprintf('Unable to find a stable model. Maximum number of attempts exceeded\n');
    AR = [];
    return;
end

repo: REST


text:
classdef VideoWriter < handle
    %VIDEOWRITER  Video Writer class
    %
    % The class provides API for writing video files or image sequences.
    %
    % Note: When querying a property that is not supported by the backend
    % used by the cv.VideoWriter instance, value 0 is returned.
    %
    % ## Example
    % Here is how to write to a video file:
    %
    %     vid = cv.VideoWriter('myvideo.mpg', [640,480]);
    %     vid.write(im);         % add a frame
    %     clear vid;             % finish
    %
    % See also: cv.VideoWriter.VideoWriter, cv.VideoWriter.write,
    %  cv.VideoCapture, VideoWriter, avifile, movie2avi,
    %  vision.VideoFileWriter
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Current quality (0..100%) of the encoded videostream. Can be
        % adjusted dynamically in some codecs.
        Quality

        % Number of stripes for parallel encoding. -1 for auto detection
        NStripes
    end

    properties (Dependent, SetAccess = private)
        % Size of just encoded video frame; note that the encoding order may
        % be different from representation order.
        FrameBytes
    end

    methods
        function this = VideoWriter(varargin)
            %VIDEOWRITER  VideoWriter constructor
            %
            %     vid = cv.VideoWriter()
            %     vid = cv.VideoWriter(filename, frameSize)
            %     vid = cv.VideoWriter(filename, frameSize, 'OptionName',optionValue, ...)
            %
            % ## Output
            % * __vid__ New instance of the VideoWriter
            %
            % The constructor initializes a video writer. When arguments are
            % supplied, the constructor calls the cv.VideoWriter.open method
            % right after creating the object.
            %
            % Accepts the same inputs and options as the cv.VideoWriter.open
            % method.
            %
            % See also: cv.VideoWriter.open, cv.VideoWriter.write
            %
            this.id = VideoWriter_(0, 'new');
            if nargin > 0
                this.open(varargin{:});
            end
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     vid.delete()
            %
            % The method first calls cv.VideoWriter.release to close the
            % already opened file.
            %
            % See also: cv.VideoWriter.release
            %
            if isempty(this.id), return; end
            VideoWriter_(this.id, 'delete');
        end

        function successFlag = open(this, filename, frameSize, varargin)
            %OPEN  Initializes or reinitializes video writer
            %
            %     successFlag = vid.open(filename, frameSize)
            %     successFlag = vid.open(filename, frameSize, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __filename__ Name of the output video file
            % * __frameSize__ Size of the video frames `[width, height]`.
            %
            % ## Output
            % * __successFlag__ bool, true if video writer has been
            %   successfully initialized.
            %
            % ## Options
            % * __API__ The parameter allows to specify API backends to use.
            %   Can be used to enforce a specific writer implementation if
            %   multiple are available: e.g. 'FFMPEG' or 'GStreamer'
            %   * __Any__ (default) Auto detect
            %   * __VfW__ Video For Windows
            %   * __QuickTime__ QuickTime
            %   * __AVFoundation__ AVFoundation framework for iOS
            %   * __MediaFoundation__ Microsoft Media Foundation
            %   * __GStreamer__ GStreamer
            %   * __FFMPEG__ FFMPEG library
            %   * __Images__ OpenCV Image Sequence (e.g. `img_%02d.jpg`)
            %   * __MotionJPEG__ Built-in OpenCV MotionJPEG codec
            %   * __MediaSDK__ Intel Media SDK
            % * __FourCC__ 4-character code of codec used to compress the
            %   frames. List of codes can be obtained at [FOURCC]. FFMPEG
            %   backend with MP4 container natively uses other values as
            %   FourCC code: see [ObjectType], so you may receive a warning
            %   message from OpenCV about fourcc code conversion. Examples are:
            %   * __PIM1__ MPEG-1 codec
            %   * __MJPG__ Motion-JPEG codec (default)
            %   * __MP42__ MPEG-4 (Microsoft)
            %   * __DIV3__ DivX MPEG-4 Part 2
            %   * __DIVX__ DivX codec
            %   * __XVID__ XVID MPEG-4 Part 2
            %   * __U263__ H263
            %   * __I263__ ITU H.263
            %   * __FLV1__ Sorenson Spark (Flash Video)
            %   * __X264__ H.264
            %   * __AVC1__ MPEG-4 Part 10/H.264 (Apple)
            %   * __WMV1__ Windows Media Video 7 (Microsoft)
            %   * __WMV2__ Windows Media Video 8 (Microsoft)
            %   * __-1__ prompts with codec selection dialog (Windows only)
            % * __FPS__ Framerate of the created video stream. default 25.
            % * __Color__ If true, the encoder will expect and encode color
            %   frames, otherwise it will work with grayscale frames (the flag
            %   is currently supported on Windows only). default true
            %
            % The method first calls cv.VideoWriter.release to close the
            % already opened file.
            %
            % The method opens a video writer:
            %
            % * On Linux FFMPEG is used to write videos;
            % * On Windows FFMPEG or VFW is used;
            % * On MacOSX QTKit is used.
            %
            % ### Tips
            % * With some backends `FourCC=-1` pops up the codec selection
            %   dialog from the system.
            % * To save image sequence use a proper filename (eg.
            %   `img_%02d.jpg`) and `FourCC=0` OR `FPS=0`. Use uncompressed
            %   image format (eg. `img_%02d.BMP`) to save raw frames.
            % * Most codecs are lossy. If you want lossless video file you
            %   need to use a lossless codecs (eg. FFMPEG FFV1, Huffman HFYU,
            %   Lagarith LAGS, etc...)
            % * If FFMPEG is enabled, using `codec=0; fps=0;` you can create
            %   an uncompressed (raw) video file.
            %
            % ## References
            % [FOURCC]:
            % > Video Codecs by [FOURCC](http://www.fourcc.org/codecs.php)
            %
            % [ObjectType]:
            % > [Codecs](http://www.mp4ra.org/codecs.html)
            %
            % See also: cv.VideoWriter.VideoWriter, cv.VideoWriter.isOpened
            %
            successFlag = VideoWriter_(this.id, 'open', filename, frameSize, varargin{:});
        end

        function retval = isOpened(this)
            %ISOPENED  Returns true if video writer has been successfully initialized
            %
            %     retval = vid.isOpened()
            %
            % ## Output
            % * __val__ bool, return value
            %
            % Returns true if video writer has been successfully initialized.
            %
            % See also: cv.VideoWriter.open
            %
            retval = VideoWriter_(this.id, 'isOpened');
        end

        function release(this)
            %RELEASE  Closes the video writer
            %
            %     vid.release()
            %
            % The method is automatically called by subsequent
            % cv.VideoWriter.open and by the cv.VideoWriter destructor.
            %
            % See also: cv.VideoWriter.open
            %
            VideoWriter_(this.id, 'release');
        end

        function write(this, frame, varargin)
            %WRITE  Writes the next video frame
            %
            %     vid.write(frame)
            %     vid.write(frame, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __frame__ The written frame
            %
            % ## Options
            % * __FlipChannels__ in case the input is color image, flips the
            %   color order from MATLAB's RGB to OpenCV's BGR order.
            %   default true
            %
            % The method writes the specified image to video file. It must
            % have the same size as has been specified when opening the video
            % writer.
            %
            % See also: cv.VideoWriter.open
            %
            VideoWriter_(this.id, 'write', frame, varargin{:});
        end

        function value = get(this, prop)
            %GET  Returns the specified VideoWriter property
            %
            %     value = vid.get(prop)
            %
            % ## Input
            % * __prop__ Property identifier. It can be specified as a string
            %   (one of the recognized properties), or directly as its
            %   corresponding integer code.
            %
            % ## Output
            % * __value__ Value for the specified property (as a `double`).
            %   Value 0 is returned when querying a property that is not
            %   supported by the backend used by the VideoWriter instance.
            %
            % See also: cv.VideoWriter.set
            %
            value = VideoWriter_(this.id, 'get', prop);
        end

        function set(this, prop, value)
            %SET  Sets a property in the VideoWriter
            %
            %     vid.set(prop, value)
            %
            % ## Input
            % * __prop__ Property identifier. It can be specified as a string
            %   (one of the recognized properties), or directly as its
            %   corresponding integer code.
            % * __value__ Value of the property (as a `double`).
            %
            % On failure (unsupported property by backend), the function
            % issues a warning.
            %
            % See also: cv.VideoWriter.get
            %
            VideoWriter_(this.id, 'set', prop, value);
        end
    end

    %% Getters/Setters
    methods
        function value = get.Quality(this)
            value = VideoWriter_(this.id, 'get', 'Quality');
        end
        function set.Quality(this, value)
            VideoWriter_(this.id, 'set', 'Quality', value);
        end

        function value = get.NStripes(this)
            value = VideoWriter_(this.id, 'get', 'NStripes');
        end
        function set.NStripes(this, value)
            VideoWriter_(this.id, 'set', 'NStripes', value);
        end

        function value = get.FrameBytes(this)
            value = VideoWriter_(this.id, 'get', 'FrameBytes');
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function copyDataToCluster(localLoc, remoteLoc, clusterHost)
%COPYDATATOCLUSTER Copies files or directories from the local machine
% to a location on a remote host.

% Copyright 2006 The MathWorks, Inc.

% Use scp to copy files.
copyCmd = sprintf('scp -r %s "%s:%s"', localLoc, clusterHost, remoteLoc);

[s, r] = system(copyCmd);
if s ~= 0
    error('distcomp:scheduler:FailedRemoteOperation', ...
        ['Failed to copy files from "%s"\n' ...
        'to "%s" on the host "%s".\n' ...
        'Command Output:\n' ...
        '"%s"\n' ...
        ], localLoc, remoteLoc, clusterHost, r);
end


text:
classdef ReferenceCellOfPathComputer < handle

    properties (Access = private)
        isInCell
    end

    properties (Access = private)
        pathVertexes
        mesh
    end

    methods (Access = public)

        function obj = ReferenceCellOfPathComputer(cParams)
            obj.init(cParams)

        end

        function cellsR = compute(obj)
            cellsR = obj.computeFirstReferenceCell();
            vertex = obj.pathVertexes;
            for ivertex = 2:length(vertex)
                vI   = vertex(ivertex);
                vOld = vertex(ivertex-1);
                cells    = obj.computeCellsOfVertex(vI);
                cellsOld = obj.computeCellsOfVertex(vOld);
                commonCells = intersect(cells,cellsOld);
                cellsR(ivertex) = commonCells(1);
            end

        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.pathVertexes = cParams.pathVertexes;
            obj.mesh         = cParams.mesh;
        end

        function cell = computeFirstReferenceCell(obj)
            firstVertex  = obj.pathVertexes(1);
            cells =  obj.computeCellsOfVertex(firstVertex);
            cell  = cells(1);
        end

        function cells = computeCellsOfVertex(obj,vertex)
            vertexInCell  = obj.mesh.connec;
            itIs = (vertexInCell == vertex);
            cells = any(itIs,2);
            cells = find(cells);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ n, x, indx ] = i4vec_index_delete_dupes ( n, x, indx )

%*****************************************************************************80
%
%% I4VEC_INDEX_DELETE_DUPES deletes duplicates from an indexed sorted I4VEC.
%
%  Discussion:
%
%    The output quantities N2, X2, and INDX2 are computed from the
%    input quantities by sorting, and eliminating duplicates.
%
%    The output arrays should be dimensioned of size N, unless the user
%    knows in advance what the value of N2 will be.
%
%    The output arrays may be identified with the input arrays.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the input list.
%
%    Input, integer X(N), the list.
%
%    Input, integer INDX(N), the sort index of the list.
%
%    Output, integer N, the number of unique entries in X.
%
%    Output, integer X(N), a copy of the list which has
%    been sorted, and made unique.
%
%    Output, integer INDX(N), the sort index of the new list.
%
  i = 0;
  n3 = 0;

  while ( 1 )

    i = i + 1;

    if ( n < i )
      break
    end

    if ( 1 < i )
      if ( x(indx(i)) == x3(n3) )
        continue
      end
    end

    n3 = n3 + 1;
    x3(n3) = x(indx(i));

  end
%
%  Copy data into output arrays.
%
  n = n3;
  x(1:n) = x3(1:n3);
  indx = i4vec_indicator1 ( n );

  return
end


text:
function [X, info] = dominant_invariant_subspace_complex(A, p)
% Returns a unitary basis of the dominant invariant p-subspace of A.
%
% function X = dominant_invariant_subspace_complex(A, p)
%
% Input: A complex, Hermitian matrix A of size nxn and an integer p < n.
% Output: A complex, unitary matrix X of size nxp such that trace(X'*A*X)
%         is maximized. That is, the columns of X form a unitary basis
%         of a dominant subspace of dimension p of A.
%
% The optimization is performed on the complex Grassmann manifold, since
% only the space spanned by the columns of X matters.
%
% See dominant_invariant_subspace for more details in the real case.
%
% See also: dominant_invariant_subspace grassmanncomplexfactory

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, June 30, 2015
% Contributors:
%
% Change log:
%
%   Xiaowen Jiang Aug. 31, 2021
%       Added AD to compute the egrad and the ehess

    % Generate some random data to test the function
    if ~exist('A', 'var') || isempty(A)
        A = randn(128) + 1i*randn(128);
        A = (A+A')/2;
    end
    if ~exist('p', 'var') || isempty(p)
        p = 3;
    end

    % Make sure the input matrix is Hermitian
    n = size(A, 1);
    assert(size(A, 2) == n, 'A must be square.');
    assert(norm(A-A', 'fro') < n*eps, 'A must be Hermitian.');
	assert(p<=n, 'p must be smaller than n.');

    % Define the cost and its derivatives on the complex Grassmann manifold
    Gr = grassmanncomplexfactory(n, p);
    problem.M = Gr;
    problem.cost  = @(X)    -real(trace(X'*A*X));
    problem.egrad = @(X)    -2*A*X;
    problem.ehess = @(X, H) -2*A*H;

    % An alternative way to compute the egrad and the ehess is to use
    % automatic differentiation provided in the deep learning toolbox
    % (slower). AD does not support complex numbers if the Matlab version
    % is R2021a or earlier. The cost function should be defined differently
    % In this case. See complex_example_AD.m and manoptADhelp.m for more
    % information.
    % problem.cost = @cost_complex;
    %    function f = cost_complex(X)
    %        AX = cprod(A,X);
    %        Xtransp = ctransp(X);
    %        product = cprod(Xtransp,AX);
    %        f = -creal(ctrace(product));
    %    end
    % call manoptAD to automatically obtain the egrad and the ehess
    % problem = manoptAD(problem);

    % If the version of Matlab installed is R2021b or later, specify the
    % cost function in the normal way and call manoptAD. Notice that
    % the function trace is not supported for AD so far. Replace it with
    % ctrace described in the file manoptADhelp.m
    % problem.cost  = @(X)    -real(ctrace(X'*A*X));
    % problem = manoptAD(problem);

    % Execute some checks on the derivatives for early debugging.
    % These can be commented out.
    % checkgradient(problem);
    % pause;
    % checkhessian(problem);
    % pause;

    % Issue a call to a solver. A random initial guess will be chosen and
    % default options are selected except for the ones we specify here.
    options.Delta_bar = 8*sqrt(p);
    [X, costX, info, options] = trustregions(problem, [], options); %#ok<ASGLU>

    fprintf('Options used:\n');
    disp(options);

    % For our information, Manopt can also compute the spectrum of the
    % Riemannian Hessian on the tangent space at (any) X. Computing the
    % spectrum at the solution gives us some idea of the conditioning of
    % the problem. If we were to implement a preconditioner for the
    % Hessian, this would also inform us on its performance.
    %
    % Notice that (typically) all eigenvalues of the Hessian at the
    % solution are positive, i.e., we find an isolated minimizer. If we
    % replace the Grassmann manifold by the Stiefel manifold, hence still
    % optimizing over orthonormal matrices but ignoring the invariance
    % cost(XQ) = cost(X) for all Q orthogonal, then we see
    % dim O(p) = p(p-1)/2 zero eigenvalues in the Hessian spectrum, making
    % the optimizer not isolated anymore.
    if Gr.dim() < 512
        evs = hessianspectrum(problem, X);
        stairs(sort(evs));
        title(['Eigenvalues of the Hessian of the cost function ' ...
               'at the solution']);
        xlabel('Eigenvalue number (sorted)');
        ylabel('Value of the eigenvalue');
    end

end

author: NicolasBoumal

repo: manopt


text:
% FT_POSTAMBLE_SAVEVAR is a helper script that optionally saves the output
% FieldTrip data structures to a *.mat file on disk. This is useful for
% batching and for distributed processing. This makes use of the
% cfg.outputfile variable.
%
% Use as
%   ft_postamble savevar data
%   ft_postamble savevar source mri
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_POSTAMBLE_SAVEFIG

% Copyright (C) 2011-2019, Robert Oostenveld, DCCN
%
% 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$

% the output data should be saved to a MATLAB file
if (isfield(cfg, 'outputfile') && ~isempty(cfg.outputfile)) || exist('Fief7bee_reproducescript', 'var')

  if isfield(cfg, 'outputlock') && ~isempty(cfg.outputlock)
    mutexlock(cfg.outputlock);
  end

  if exist('Fief7bee_reproducescript', 'var')
    % write the output variable(s) to a MATLAB file
    iW1aenge_now = datestr(now, 30);
    cfg.outputfile = {};
    for i=1:numel(iW1aenge_postamble)
      cfg.outputfile{i} = fullfile(Fief7bee_reproducescript, sprintf('%s_%s_output_%s.mat', ...
        iW1aenge_now, FjmoT6aA_highest_ft, iW1aenge_postamble{i}));
    end

    % write the large configuration fields to a MATLAB file
    % this applies to layout, event, sourcemodel, headmodel, grad, etc.
    % note that this is here, rather than in the (seemingly more logical)
    % ft_preamble_loadvar, because this code depends on cfg.callinfo (which
    % is only present at postamble stage)
    cfg = save_large_cfg_fields(cfg, FjmoT6aA_highest_ft, Fief7bee_reproducescript, iW1aenge_now);

    % write a snippet of MATLAB code with the user-specified configuration and function call
    reproducescript(Fief7bee_reproducescript, fullfile(Fief7bee_reproducescript, 'script.m'),...
      FjmoT6aA_highest_ft, cfg, isempty(iW1aenge_postamble));

    % instruct savevar() to also write hashes of output files to the
    % hashfile
    Fief7bee_hashfile = fullfile(Fief7bee_reproducescript, 'hashes.mat');

  elseif (isfield(cfg, 'outputfile') && ~isempty(cfg.outputfile))
    % keep the output file as it is
    Fief7bee_hashfile = [];
  else
    % don't write to an output file
    cfg.outputfile = {};
    Fief7bee_hashfile = [];
  end

  if isequal(iW1aenge_postamble, {'varargout'}) && ~iscell(cfg.outputfile)
    % this should be a cell-array, oterwise it cannot be matched with varargout
    cfg.outputfile = {cfg.outputfile};
  end

  % save the output data structures to a MATLAB .mat file
  if iscell(cfg.outputfile)
    % iW1aenge_postamble is a cell-array containing the variable names
    if isequal(iW1aenge_postamble, {'varargout'})
      % the output is in varargout
      for tmpindx=1:length(cfg.outputfile)
        savevar(cfg.outputfile{tmpindx}, 'data', varargout{tmpindx}, Fief7bee_hashfile);
      end % for
      clear tmpindx
    else
      % the output is in explicitly named variables
      for tmpindx=1:length(cfg.outputfile)
        savevar(cfg.outputfile{tmpindx}, iW1aenge_postamble{tmpindx}, eval(iW1aenge_postamble{tmpindx}), Fief7bee_hashfile);
      end % for
      clear tmpindx
    end
  else
    % iW1aenge_postamble{1} contains the name of the only variable
    savevar(cfg.outputfile, iW1aenge_postamble{1}, eval(iW1aenge_postamble{1}), Fief7bee_hashfile);
  end

  if isfield(cfg, 'outputlock') && ~isempty(cfg.outputlock)
    mutexunlock(cfg.outputlock);
  end

  if ~ft_nargout
    % do not return the output variable "ans"
    clear(iW1aenge_postamble{1});
  end
end

author: spm


text:
function c = imcmass(im)
% IMCMASS  Centre of mass of an image.
%
% IMCMASS computes the centre of mass of a grayscale image weighting the
% pixel coordinates by the pixel intensities.
%
% C = IMCMASS(IM)
%
%   IM is a grayscale 2D image.
%
%   C is the center of mass of the image in [row, col] format.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2015 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 1);
nargoutchk(0, 1);

if (size(im, 3) > 1)
    error('IM must be a grayscale image')
end

% grid for image coordinates
[grow, gcol] = ndgrid(1:size(im, 1), 1:size(im, 2));

% centre of mass
c = single(im(:));
c = sum([c .* grow(:), c .* gcol(:)]) / sum(c);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [data,ctf] = readepochs(folder,varargin);
%
%
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%      <                                                       >
%      <                      DISCLAIMER:                      >
%      <                                                       >
%      <  THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. >
%      <  THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   >
%      <                     OFFICIAL USE.                     >
%      <                                                       >
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>
%
% function to read set time windows (epochs) around event markers.
%
[marker_info] = readmarkerfile(folder);
ctf = ctf_read_res4(folder,1);
CHAN = ctf.sensor.index.megsens;
trials = [1:ctf.setup.number_trials];
win = 'alltimes';
for i = 1:size(varargin,2)

  if ~iscell(varargin{i})

    switch num2str(varargin{i})

      case 'megsens',
        CHAN = ctf.sensor.index.megsens(varargin{i+1});

      case 'refsens',
        CHAN = ctf.sensor.index.refsens(varargin{i+1});

      case 'eegsens',
        CHAN = ctf.sensor.index.eegsens(varargin{i+1});

      case 'othersens',
        CHAN = ctf.sensor.index.othersens(varargin{i+1});

      case 'vc',
        CHAN = ctf.sensor.index.vc(varargin{i+1});

      case 'allchans',
        CHAN = [1:ctf.setup.number_channels];

      case 'trials',
        trials = varargin{i+1};

      case 'markers',
        marks = varargin{i+1};

      case 'window',
        win = varargin{i+1};

    end
  end
end

if ~exist('marks','var'),

  ctf = ctf_read_meg4(folder,ctf,CHAN,win,trials);

  epochs = cell(1,1);

  epochs{1} = zeros(size(ctf.data{1},1),size(ctf.data{1},2),size(ctf.data,1));

  for i = 1:size(ctf.data,1),

    epochs{1}(:,:,i) = ctf.data{i};

  end

else

  nm = size(marks,2);

  epochs = cell(1,nm);

  for mkr = 1:nm,

    mk=ismember(marker_info.marker_names,marks(mkr));

    marker_info.marker_names(mk);

    nsamp=marker_info.number_samples(mk)

    nss=0;

    for ns=1:nsamp

      tr = marker_info.trial_times{mk}(ns,1);

      if ismember(tr,trials)

        nss=nss+1;

        tim=marker_info.trial_times{mk}(ns,2);

        times=win+tim;

        ctf = ctf_read_ds(folder,ctf,CHAN,times,tr);

        temp=ctf.data{1};

        if nss==1,
          epochs{mkr}=zeros(size(temp,1),size(temp,2),1);
        end

        epochs{mkr}(:,:,nss)=temp;

      end
    end
  end
end
data = struct('epochs',{epochs},'setup',{ctf.setup},'sensor',{ctf.sensor});
return

author: PatternRecognition

repo: OpenBMI


text:
function output = callpenlab(model)

model = yalmip2nonlinearsolver(model);

if ~model.derivative_available
    disp('Derivate-free call to penlab not supported')
    error('Derivate-free call to penlab not implemented')
end
if model.options.savedebug
    save penlabdebug model
end

penm = [];
penm.Nx=length(model.linearindicies);

if length(nnz(model.K.s)>0)
    top = 1 + model.K.f + model.K.l + sum(model.K.q);
    for i = 1:length(model.K.s)
        model.vecF{i} = model.F_struc(top:top+model.K.s(i)^2-1,:);
        top = top + model.K.s(i)^2;
    end
end

showprogress('Calling PENLAB',model.options.showprogress);

% These are needed to avoid recomputation due to penlab double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G = [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

penm.userdata=model;

penm.objfun =  @(x,Y,userdata) penlab_callback_f(x,userdata);
penm.objgrad = @(x,Y,userdata) penlab_callback_df(x,userdata);
penm.objhess = @(x,Y,userdata) penlab_callback_df2(x,userdata);

if length(model.b)>0 | length(model.beq)>0
    penm.lbg = [model.beq;model.b-inf];
    penm.ubg = [model.beq;model.b];
    penm.NgLIN = length(model.beq) + length(model.b);
    penm.confun =  @(x,Y,userdata) penlab_callback_con(x,userdata);
    penm.congrad = @(x,Y,userdata) penlab_callback_dcon(x,userdata);
    penm.conhess = @(x,Y,k,userdata) penlab_callback_dcon2(x,k,userdata);
end

if model.K.s(1)>0
    penm.NALIN=length(model.K.s);
    penm.lbA=zeros(find(model.K.s)>0,1);
    penm.mconfun  = @(x,Y,k,userdata) penlab_callback_matrixG(x, k,userdata);
    penm.mcongrad = @(x,Y,k,i,userdata) penlab_callback_matrixdG(x,k,i,userdata);
end

prob=penlab(penm);

solvertime = tic;
prob.solve();
solvertime = toc(solvertime);
x = prob.x;

% Duals currently not supported
lambda = [];

switch 0
    case {0,1}
        problem = 0;
    case {2}
        problem = 1;
    case {-1}
        problem = 3;
    case {3,4,-2,-3}
        problem = 4;
    case {-11,-12,-13}
        problem = 7;
    case {-10,-100,-101,-102,-199}
        problem = 11;
    otherwise
        problem = -1;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(x,D_struc,[],problem,'PENLAB',solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function dFdz = deriv(F,z0,Fz0,type)
%DERIV Approximate gradient and Jacobian.
%   For real x0, deriv(F,x0) computes the gradient dF(x0)/dx or Jacobian
%   dF(x0)/dx^T of a real-valued function F(x) in the real variables x,
%   assuming F is analytic in x when x is complex. If this assumption is
%   not satisfied (e.g., when f is function of x' instead of x.'), use the
%   finite difference methods below. The input variables x may be a scalar,
%   vector, matrix, tensor or even a (nested) cell array of tensors.
%
%   For complex z0, deriv(F,z0) and deriv(F,z0,Fz0) are equivalent to
%   deriv(F,z0,Fz0,'gradient') if F(z0) is a real scalar, and
%   deriv(F,z0,Fz0,'Jacobian') otherwise. The input variables z may be a
%   scalar, vector, matrix, tensor or even a cell array of tensors.
%
%   deriv(f,z0,fz0,'gradient') approximates the real gradient df(z0)/dx or
%   scaled conjugate cogradient 2*conj(df(z0)/dz) of a real-valued function
%   f in the real variables x or complex variables z, respectively. The
%   scaled conjugate cogradient is defined as twice the complex conjugate
%   of the partial derivative of f w.r.t. the complex variables z, while
%   treating conj(z) as constant. The real gradient is returned if z0 is
%   real, else the scaled conjugate cogradient is returned. Note that the
%   former is equal to the latter if both z0 and f(z) are real. The input
%   fz0 should supply f(z0). If fz0 is [], it is computed as f(z0).
%
%   deriv(F,z0,Fz0,'Jacobian') approximates the Jacobian dF(z0)/dz^T in the
%   real or complex variables z. If F(z) is nonanalytic in z and depends on
%   conj(z), use 'Jacobian-C' instead. The input Fz0 should supply F(z0).
%   If Fz0 is [], it is computed as F(z0).
%
%   deriv(F,z0,Fz0,'Jacobian-C') approximates the complex Jacobian
%   [dF(z0)/dz^T dF(z0)/dconj(z)^T] of a function F in the complex
%   variables z. The partial derivatives in the complex Jacobian are w.r.t.
%   the complex variables z and conj(z), while treating conj(z) and z as
%   constant, respectively. The input Fz0 should supply F(z0). If Fz0 is
%   [], it is computed as F(z0).
%
%   See also gradient, diff.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables", SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.

% If possible, use the i-trick.
dim = structure(z0);
z0 = serialize(z0);
F = @(z)serialize(F(deserialize(z,dim)));
if nargin < 3 && all(isreal(z0))
    % Assume F is analytic and both F and z0 are real-valued; return dF/dx.
    e = 2^floor(log2(nthroot(realmin(class(z0)),3)));
    p = zeros(size(z0));
    p(1) = e;
    dFdz = imag(F(z0+p*1i))/e;
    dFdz = [dFdz zeros(numel(dFdz),length(z0)-1)];
    for n = 2:length(z0)
        p(n-1) = 0; p(n) = e;
        dFdz(:,n) = imag(F(z0+p*1i))/e;
    end
    if size(dFdz,1) == 1
        dFdz = deserialize(dFdz(:),dim);
    end
    return;
elseif nargin < 3 || isempty(Fz0)
    Fz0 = F(z0);
end
Fz0 = serialize(Fz0);

% Define type of derivative.
if nargin < 4
    if numel(Fz0) == 1 && isreal(Fz0), type = 'gradient';
    else type = 'Jacobian'; end
end
dFdzIsGrad = any(strfind(type,'gradient'));
dFdzIsRealGrad = dFdzIsGrad && all(isreal(z0));
dFdzIsAnaJac = ~dFdzIsGrad && ~any(strfind(type,'-C'));

% Else use finite (complex) derivatives.
e = sqrt(eps(class(z0)));
p = zeros(size(z0));
p(1) = e*max(1,abs(real(z0(1))));
dFdx = (F(z0+p)-Fz0)/p(1);
if ~dFdzIsRealGrad || ~dFdzIsAnaJac
    p(1) = e*max(1,abs(imag(z0(1))))*1i;
    dFdy = (F(z0+p)-Fz0)/imag(p(1));
end
switch strtok(type,'-')
    case 'gradient'
        if dFdzIsRealGrad, dFdz = dFdx;
        else dFdz = real(dFdx)+real(dFdy)*1i; end
        dFdz = [dFdz zeros(numel(dFdz),length(z0)-1)];
    case 'Jacobian'
        if dFdzIsAnaJac
            dFdz = dFdx;
            dFdz = [dFdz zeros(numel(dFdz),length(z0)-1)];
        else
            dFdz = 0.5*cat(3,dFdx-dFdy*1i,dFdx+dFdy*1i);
            dFdz = cat(2,dFdz,zeros(size(dFdz,1),length(z0)-1,2));
        end
end
for n = 2:length(z0)
    p(n-1) = 0; p(n) = e*max(1,abs(real(z0(n))));
    dFdx = (F(z0+p)-Fz0)/p(n);
    if ~dFdzIsRealGrad || ~dFdzIsAnaJac
        p(n) = e*max(1,abs(imag(z0(n))))*1i;
        dFdy = (F(z0+p)-Fz0)/imag(p(n));
    end
    switch strtok(type,'-')
        case 'gradient'
            if dFdzIsRealGrad, dFdz(:,n) = dFdx;
            else dFdz(:,n) = real(dFdx)+real(dFdy)*1i; end
        case 'Jacobian'
            if dFdzIsAnaJac, dFdz(:,n) = dFdx;
            else dFdz(:,n,:) = 0.5*cat(3,dFdx-dFdy*1i,dFdx+dFdy*1i); end
    end
end
dFdz = dFdz(:,:);
if size(dFdz,1) == 1 && dFdzIsGrad
    dFdz = deserialize(dFdz(:),dim);
end

end

function [z,offset] = deserialize(z,dim,offset)
    if iscell(dim)
        v = z;
        z = cell(size(dim));
        if nargin < 3, offset = 0; end
        for i = 1:numel(z)
            if iscell(dim{i})
                [z{i},offset] = deserialize(v,dim{i},offset);
            else
                n = prod(dim{i}(:));
                z{i} = reshape(v(offset+(1:n)),dim{i});
                offset = offset+n;
            end
        end
    elseif ~isempty(dim)
        z = reshape(z,dim);
    end
end

function z = serialize(z)
    if iscell(z)
        for i = find(cellfun(@iscell,z(:).'))
            z{i} = serialize(z{i});
        end
        s = cellfun(@numel,z(:)); o = [0; cumsum(s)];
        c = z; z = zeros(o(end),1);
        for i = 1:length(s), z(o(i)+(1:s(i))) = c{i}(:); end
    else
        z = z(:);
    end
end

function dim = structure(z)
    if iscell(z)
        dim = cellfun(@size,z,'UniformOutput',false);
        for i = find(cellfun(@iscell,z(:).'))
            dim{i} = structure(z{i});
        end
    else
        dim = size(z);
        if numel(z) == dim(1), dim = []; end
    end
end

author: andrewssobral

repo: lrslibrary