text:
function [ fitted_shape, fitted_app ] = fit_2d(AAM, init_shape, image_data, max_it)
% [ fitted_shape fitted_app ] = fit_2d(AAM, init_shape, image_data, max_it)
%
% Fit an Active Appearance Model to an image, starting from the provided initial shape.
% Based on the 'Active Appearance Models Revisited' paper by Iain Matthews and Simon Baker.
%
%                         PARAMETERS
%
% AAM Active Appearance Models structure as returned by build_model_2d
%
% init_shape Shape to be used to start the fitting process.
%            It is A Nx2 matrix containing N landmarks over image_data,
%            with image coordinates stored on rows: [i j]
%
% image_data The image to be fitted.
%            It is A HxWxC floating point image, with pixels arranged such
%            that app_data(i,j,c) actually returns the pixel value
%            associated with position (i,j). Must be in floating point
%            format, ranging between 0 and 1.
%
% max_it Maximum number of iterations to be performed.
%
%
%                        RETURN VALUES
%
% fitted_shape The shape produced by the fitting process. It has the
%              same format as init_shape.
%              Due to the local nature of the optimization process, there is no
%              guarantee of a meaningful result. Hopefully, the closer
%              the initial guess is, the better the fit.
%
%
% fitted_app (Optional) The fitted appearance. It is the best approximation of
%            the appearance that the model can do.
%            It is a MHxMWxC image matrix, where MH = size(AAM.warp_map, 1) and
%            MW = size(AAM.warp_map, 2) and C is the number of colors and
%            is the same as image_data.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

	if isa(image_data, 'float') == 0
		warning 'Image data is not in floating point format, assuming unsigned integer data.'
		image_data = double(image_data) / 255;
	end

	% Number of points in the shape
	np = size(init_shape, 1);

	% Number of colors
	nc = size(image_data, 3);

	% Size of the model
	modelh = size(AAM.warp_map, 1);
	modelw = size(AAM.warp_map, 2);

	% Check that the shape has the right size
	if np ~= (size(AAM.s0, 1) / 2) || nc ~= (size(AAM.A0, 1) / (modelw * modelh))
		msg = strcat('Number of landmarks in init_shape and/or number of ', ...
		             ' colors in image_data are not consistent with the model.');
		error(msg)
	end

	% Initialize current shape
	cur_shape = init_shape;

	% Used to estimate convergence, it's the root mean
	% square change on the shapes
	%rms_change = inf;

	iter = 1;

	%while (rms_change > 1 || iter < 7) && iter <= max_it
	while iter <= max_it
		% Warp the appearance to compute the error image, as a column vector
		%warped_app = reshape(pa_warp(AAM, cur_shape, image_data), [], 1);

		% Error image, simply the difference between the image and A0
		error_img = reshape(pa_warp(AAM, cur_shape, image_data), [], 1) - AAM.A0;

		% In the first 5 iterations, try to align the shape to the image using only
		% the affine transform. This will maximize convergence by using smaller error
		% images in determining the incremental shape parameters, giving smaller but
		% more meaningful updates.
		if iter > 5 || max_it < 10
			% Incremental parameter updates which should minimize the error.
			% The first four will be the delta_q and the latter ones are the
			% delta_p
			delta_qp = AAM.R * error_img;

			% Damping factor, used to scale parameters updates, thus avoiding
			% large updates to parameters associated with less important
			% modes.
			% FIXME: this improves fitting stability, but at the expense of precision.
			% damping = min(AAM.shape_eiv / AAM.shape_eiv(min(round((iter-6)/4+1),numel(AAM.shape_eiv))), ones(size(AAM.shape_eiv)));
			% damping = AAM.shape_eiv / AAM.shape_eiv(1);
			%
			% Compute the incremental warp using the *negated* incremental
			% shape parameters and reshape to prepare for the affine transform
			% d_s0 = reshape(AAM.s0 - sum(AAM.s * (damping .* delta_qp(5:end)), 2), np, 2);

			d_s0 = reshape(AAM.s0 - sum(AAM.s * delta_qp(5:end), 2), np, 2);

			% Compute the incremental affine transformation using the *negated*
			% parameters
			[ A trans ] = to_affine(AAM, -delta_qp(1:4)');

			d_s0 = d_s0 * A + repmat(trans, np, 1);

			comp_warp = pa_warp_composition(AAM, cur_shape, d_s0);
		else
			delta_q = AAM.R(1:4,:) * error_img;

			% Compute the incremental affine transformation using the *negated*
			% parameters
			[ A trans ] = to_affine(AAM, -delta_q');
			d_s0 = reshape(AAM.s0, np, 2);
			d_s0 = (d_s0 * A + repmat(trans, np, 1)) - d_s0;
			comp_warp = cur_shape + d_s0;
		end

%figure(1);
%hold off;
%imshow(image_data);
%hold on;
%triplot(AAM.shape_mesh, cur_shape(:,2), cur_shape(:,1), 'b');
%triplot(AAM.shape_mesh, comp_warp(:,2), comp_warp(:,1), 'r');
%
%pause;

		% Convergence estimate is simply given by relative average vertex rms_change
		%rms_change = sqrt(sum(sum((cur_shape - comp_warp).^2)) / np);
		%pause

		iter = iter + 1;
		% Update current shape
		cur_shape = comp_warp;
	end

	% Set return value
	fitted_shape = cur_shape;

	if nargout > 1
		% Error image, simply the difference between the warped image and A0
		error_img = reshape(pa_warp(AAM, cur_shape, image_data), [], 1) - AAM.A0;

		% Project the error image onto the appearance modes and use the resulting parameters
		% to build the reconstructed appearance
		fitted_app = reshape(AAM.A0 + AAM.A * (AAM.A' * error_img), [modelh modelw nc]);
	end


text:
function num = licenseuser(username,verbose)
%LICENSEUSER Displays the names of users who have a license checked out.
% NUM = LICENSEUSER(USERNAME) displays the feature toolboxes that are
%  currently checked out by USERNAME from a shared network license.
%  LICENSEUSER uses the Flexlm (R) license manager utility to get
%  information about the license status.
%
% USERNAME is the username on the host system. Default, if USERNAME is not
%  given, is to try and get the username of the current user from the host.
% VERBOSE (0|1|2) is the level of messages printed on the command line.
%
% NUM returns the total number of licenses checked out by USERNAME.
%
% Example:
%  To see a list of the toolboxes that user wsmith is using, use:
%
%  >> licenseuser('wsmith')
%
%
% Note that the only way to release a license is to exit Matlab.
%
%   See also LICENSEUSE, LICENSE, VER.
%
%
% Matt Hopcroft
%  mhopeng@gmail.com
%

% MH Nov2012
% v2.01 handle Win/Mac/Linux environment username
% v2.0  update for current Matlab releases (& sync with licenseuse)
% MH Jun2011
% v1.1  return total number of licenses checked out
% MH Apr2011
% v1.0
%

%% Welcome
if nargin < 2
    verbose=1;
end
if verbose >= 1, fprintf(1,'licenseuser v2.01: %s\n',datestr(clock)); end

if nargin < 1 || isempty(username)
    % use this ridiculous logic to handle Win/Mac/Linux (v2.01)
    username=[];
    try %#ok<TRYNC>
        username = getenv('USER');
    end
    if isempty(username)
        try %#ok<TRYNC>
            username = getenv('USERNAME');
        end
    end
end

if isempty(username)
    error('licenseuser: Could not determine username');
end

% get the list of all active licenses ("-a")
fprintf(1,'licenseuser: Finding licenses checked out to "%s"...\n',username);

cdir=cd;

%% Determine path to lmutil
if isunix % v2.11
    lmExe='lmutil';
else
    lmExe='lmutil.exe';
end

lmPath=[matlabroot filesep 'bin' filesep computer('arch') filesep]; % prior to R2010b
if ~exist([lmPath lmExe],'file')
    lmPath=[matlabroot filesep 'etc' filesep computer('arch') filesep]; % R2010b+
    if ~exist([lmPath lmExe],'file')
        fprintf(1,'ERROR: "%s" not found in bin/ or etc/ locations!\n\n',lmExe);
        return
    end
end
if verbose >= 3, fprintf(1,'License manager utilities found at:\n  %s\n',lmPath); end

cd(lmPath);

%% Execute license manager command
[ltstat, lt]=system(['.' filesep 'lmutil lmstat -c "' matlabroot filesep 'licenses' filesep 'network.lic" -a']);

cd(cdir);

if ltstat ~= 0
    fprintf(1,'ERROR: "%s"\n',strtrim(lt));
    if strfind(lt, 'No such file')
        fprintf(1,'On Linux you may need to upgrade the package "lsb"\n  (see %s)\n',...
            'https://www.mathworks.com/support/solutions/en/data/1-GLXUHV/index.html?product=ML&solution=1-GLXUHV');
    end
    error('OS returned error %d when executing lmutil',ltstat);
else
    if verbose >= 2, disp(lt); end % display the license manager status
end

%% Process the list

% tline = 1; tlineMatch={};
% k=0; k2=0;

tlines=textscan(lt,'%s','delimiter','\n'); tlines=tlines{1};
tIndC=strfind(tlines,'Users of');
%length(tIndC)
uIndC=strfind(tlines,username);
% create numeric binary arrays
for i=1:length(tIndC)
    if isempty(tIndC{i})
        tInd(i)=0;                                                              %#ok<*AGROW>
    else
        tInd(i)=1;
    end
    if isempty(uIndC{i})
        uInd(i)=0;
    else
        uInd(i)=1;
    end
end

tlbInd=find(tInd);
usrInd=find(uInd);

if ~isempty(usrInd)
    num = length(usrInd);
    fprintf(1,'User %s has %d licenses checked out:\n',username,num);

    for i=usrInd
        tlbx=tlbInd(tlbInd<i);
        %tlines{tlbx(end)}
        toolbox=textscan(tlines{tlbx(end)},'Users of %s','Delimiter',':'); toolbox=toolbox{1};
        fprintf(1,' %s',toolbox{1});
        hostname=textscan(tlines{i},[username '%s']); hostname=hostname{1};
        fprintf(1,' (on host %s)\n',hostname{1});
    end
    fprintf(1,'\n');
else
    fprintf(1,'No licenses found for username "%s".\n',username);
    num = 0;
end

author: Sable

repo: mcbench-benchmarks


text:
%% Testing Will's adiff
clc
clear
%Array of functions to check
funs = {'abs','acos','asin','atan','cos','cosh','exp','log','log10','sin','sinh','sqrt','tan','tanh'};

%Check against both numerical and symbolic gradients
n = length(funs);
ajac = zeros(n,1);
njac = zeros(n,1);
sjac = zeros(n,1);

x0 = abs(0.5+0.2*randn);

%Test each one
for i = 1:n
    fun = str2func(['@(x)' funs{i} '(x)'])
    ajac(i) = autoJac(fun,x0);
    njac(i) = mklJac(fun,x0);
    s = symJac(fun);
    sjac(i) = s(x0);
end

%Display Results
raw = [ajac njac sjac]
acc = [abs(sjac-ajac) abs(sjac-njac)]

%% Index Checking

%% Vection addition1
%double + adiff
clc
clear
a = 2;
b = [1 3 5; 2 4 6];
x = adiff(2);

%2+x1
y2 = a+x(1)

%2+[x..]
y = a+x

%[1..]+x1
y = b+x(1)

%[1..]+[x..]
z = b+x

%% Vector addition2
%adiff + double
clc
clear
a = 2;
b = [1 3 5; 2 4 6];
x = adiff(2);

%x1+2
x(1)+a

%[x..]+2
y = x+a

%x1+[1..]
x(1)+b

%[x..]+[1..]
z = x+b

%% Vector addition3
%adiff + adiff
clc
clear all
x = adiff(2);

%x1+x1
y = x(1)+x(1)

%x1+[x..]
y2 = x(1)+x

%[x..]+x1
z = x+x(1)

%[x..]+[x..]
z2 = x+x

%% Loop Checking
clc
% clear all
b = 1:3;
x = adiff(1:3);
j = x(1)*x(2);
for i = 1:3
    j = j + b(i)*x(i);
end
j

%% Test1 2D
clear all
clc
x0 = [-1;2];
x = adiff(x0);

fun = @(x) 1/(27*sqrt(3)) * ((x(1) - 3)^2 - 9)*x(2)^3;

b = fun(x)

fun(x0)
[f,j] = adiffget(b)

s = symJac(fun);
j2 = s(x0)

%% Test2 3D
clear all
clc
x0 = [-1;2;4];
x = adiff(x0);

fun = @(x) 9 - 8*x(1) - 6*x(2) - 4*x(3) + 2*x(1)^2 + 2*x(2)^2 + x(3)^2 + 2*x(1)*x(2) + 2*x(1)*x(3);

b = fun(x)

fun(x0)
[f,j] = adiffget(b)

s = symJac(fun);
j2 = s(x0)

%% Test3 4D
clear all
clc
x0 = [-1;2;4;-5];
x = adiff(x0);

fun = @(x) 100*(x(2)-x(1)^2)^2 + (1-x(1))^2 + 90*(x(4)-x(3)^2)^2 + (1-x(3))^2 + 10.1*((x(2)-1)^2 + (x(4)-1)^2) + 19.8*(x(2)-1)*(x(4)-1);

b = fun(x)

fun(x0)
[f,j] = adiffget(b)

s = symJac(fun);
j2 = s(x0)

%% Test4 5D
clear all
clc
a = 100;
x0 = [-1;2;4;-5;7];
x = adiff(x0);

fun = @(x) (x(1)-x(2))^2 + (x(2)-x(3))^3 + (x(3)-x(4))^4 + (x(4)-x(5))^4;

fun(x0)
b = fun(x)

[f,j] = adiffget(b)

s = symJac(fun);
j2 = s(x0)

%% Subsasgn, Empty LHS
clear all
clc
x = adiff(1:4);
X = adiff(2);
a = [1;2;3;4];

%Simple indexing
y = x(1)*x(2)

%Simple assign
y2(2) = x(3)*2*x(4)
%Simple 2D assign
y22(2) = x(3)*2*x(4)

% %Multi index, indexing
z = x(1:3)

%Multi index, assign
y3(1:2) = [x(1);x(3)]

%Multi index, assign, longer vec
y4(2:3) = [x(3);x(4)]

%Logical Index
y5(a>2) = x(a>2)

%% Subsasgn, Existing LHS
clear
clc
x = adiff(1:4);
a = [1;2;3;4];

%Add element to vector
y = x(1)
y(2) = x(2)
y(3) = x(3)
y(5) = x(4)

%Add element to indexed vector
y2(1) = x(1)

%Add element to existing vector
% y2(:,2) = x(1:2)
% y2(:,3) = x;

%% Subsagn, Element Numbers
clear
clc
x = adiff(2);

% X([8 9]) = [x(1); x(2)]

% X(:,1:2) = barvec(2)

%% Multiline Equation + multi index
clear
clc
x = adiff(1:4);

f = 3*x.^2 - x(1)
g = sum(f)
z1(1:4) = f
% z1(5) = g
% y = z1(5)*2
%
% %Multi Assign
% c(1,1) = 2*x(1);
% c(2:3,1) = 3*x(2:3).^2;
% c(4,1) = 0.1*x(4)

%% Multi index subs
clear
clc
x = adiff(1:3);

y = sum(x(1:2).*x(2:3))

% y.vecstr

%% Multi2
clear all
clc

b = [1:4]'
x = adiff(1:4);

f = @(x) b(1:2).*x(1:2) + b(3:end).*x(3:end).^2
y = f(x)

x0 = ones(4,1);
f(x0)

%% Multi assign
clear
clc

x = adiff(1:4);

f([2 3]) = 2*x([3 2])
f([1 4]) = 5*x([1 4])

%% Constraint Test
clc
clear
x = adiff(1:5);

nlcon = @(x) [x(1) + 2*x(2) + 3*x(3) - 6;
              x(2) + 2*x(3) + 3*x(4) - 6;
              x(3) + 2*x(4) + 3*x(5) - 6];

b = nlcon(x)

[f,j] = adiffget(b)

s = symJac(nlcon);
j2 = s(ones(5,1))

%% Constraint Test 2
clc
clear
a = 100;
x = adiff(1:5);

c(3) = x(1) + 2*x(2) + 3*x(3) - 6;
c(2) = x(2) + 2*x(3) + 3*x(4) - 6;
c(1) = x(3) + 2*x(4) + 3*x(5) - 6;
c

% y = c*x(1:3)

[f,j] = adiffget(c)

%% DIW

clc
clear all

x = adiff([1 1])

length(x)
size(x)
ndims(x)

n = length(x);
j = @(x) sum((1-x).^2) + sum(100*(x(2:n) - x(1:n-1).^2).^2)

[f,jac] = adiffget(j(x))

mklJac(j,ones(2,1))

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function cleanup(obj)
% remove temporary files created by prepare
%
% cleanup(obj)
%
% Inputs:
%   obj                 MOcovMFileCollection instance
%
% Notes:
% - this function should be called after coverage has been determined.
%
% See also: rewrite_mfiles, prepare

    notify(obj.monitor,'Cleanup');
    if ~isempty(obj.orig_path)
        notify(obj.monitor,'','Resetting path');
        path(obj.orig_path);
    end

    if ~isempty(obj.temp_dir)
        msg=sprintf('Removing temporary files in %s',obj.temp_dir);
        notify(obj.monitor,'',msg);

        if mocov_util_platform_is_octave()
            % GNU Octave requires, by default, confirmation when using
            % rmdir. The state of confirm_recursive_rmdir is stored,
            % and set back to its original value when leaving this
            % function.
            confirm_val=confirm_recursive_rmdir(false);
            cleaner=onCleanup(@()confirm_recursive_rmdir(confirm_val));
        end

        rmdir(obj.temp_dir,'s');
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function M = femMultiplyCell(varargin)
%+========================================================================+
%|                                                                        |
%|              OPENFEM - LIBRARY FOR FINITE ELEMENT METHOD               |
%|           openFem is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018.          |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : femMultiplyCell.m                             |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & François Alouges            |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Multiply finite elements cells                |
%|  `---'  |                                                              |
%+========================================================================+

if (nargin == 1)
    M = varargin{1};

elseif (nargin == 2)
    A = varargin{1};
    B = varargin{2};

    if iscell(A) && iscell(B)
        M = A{1} * B{1} + ...
            A{2} * B{2} + ...
            A{3} * B{3} ;

    elseif iscell(A) && ~iscell(B)
        M{1} = A{1} * B ;
        M{2} = A{2} * B ;
        M{3} = A{3} * B ;

    elseif ~iscell(A) && iscell(B)
        M{1} = A * B{1} ;
        M{2} = A * B{2} ;
        M{3} = A * B{3} ;

    elseif ~iscell(A) && ~iscell(B)
        M = A * B;

    else
        error('femMultiplyCell.m : unavailable case')
    end

elseif (nargin == 3)
    A = varargin{1};
    B = varargin{2};
    C = varargin{3};

    if iscell(A) && iscell(B) && iscell(C)
        M = A{1} * B{2} * C{3} - ...
            A{1} * B{3} * C{2} + ...
            A{2} * B{3} * C{1} - ...
            A{2} * B{1} * C{3} + ...
            A{3} * B{1} * C{2} - ...
            A{3} * B{2} * C{1} ;

    elseif iscell(A) && iscell(B) && ~iscell(C)
        M = ( A{1} * B{1} + ...
            A{2} * B{2} + ...
            A{3} * B{3} ) * C;

    elseif iscell(A) && ~iscell(B) && iscell(C)
        M = A{1} * B * C{1} + ...
            A{2} * B * C{2} + ...
            A{3} * B * C{3} ;

    elseif ~iscell(A) && iscell(B) && iscell(C)
        M = A * ( B{1} * C{1} + ...
            B{2} * C{2} + ...
            B{3} * C{3} );

    elseif ~iscell(A) && ~iscell(B) && iscell(C)
        M{1} = A * B * C{1};
        M{2} = A * B * C{2};
        M{3} = A * B * C{3};

    elseif ~iscell(A) && iscell(B) && ~iscell(C)
        M{1} = A * B{1} * C;
        M{2} = A * B{2} * C;
        M{3} = A * B{3} * C;

    elseif iscell(A) && ~iscell(B) && ~iscell(C)
        M{1} = A{1} * B * C;
        M{2} = A{2} * B * C;
        M{3} = A{3} * B * C;

    elseif ~iscell(A) && ~iscell(B) && ~iscell(C)
        M = A * B * C;
    else
        error('femMultiplyCell.m : unavailable case')
    end

else
    error('femMultiplyCell.m : unavailable case')
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function s = addtls(s)
%ADDTLS Add the t location-scale distribution.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.7 $  $Date: 2004/01/24 09:35:08 $

j = length(s) + 1;
s(j).name = 't location-scale';
s(j).code = 'tlocationscale';
s(j).pnames = {'mu' 'sigma' 'nu'};
s(j).pdescription = {'location' 'scale' 'shape'};
s(j).prequired = [false false false];
s(j).fitfunc = @tlsfit;
s(j).likefunc = @tlslike;
s(j).cdffunc = @tlscdf;
s(j).pdffunc = @tlspdf;
s(j).invfunc = @tlsinv;
s(j).statfunc = @tlsstat;
s(j).loginvfunc = [];
s(j).logcdffunc = [];
s(j).hasconfbounds = false;
s(j).censoring = true;
s(j).paramvec = true;
s(j).support = [-Inf Inf];
s(j).closedbound = [false false];
s(j).iscontinuous = true;
s(j).islocscale = false;
s(j).uselogpp = false;

% ==== t Location-Scale distribution functions ====

% these distribution functions do not yet handle arrays of parameters

function y = tlspdf(x, mu, sigma, nu)
%TLSPDF T location-scale probability density function (pdf).
sigma(sigma <= 0) = NaN;
nu(nu <= 0) = NaN;

y = tpdf((x - mu)./sigma,nu)./sigma;

function p = tlscdf(x ,mu, sigma, nu)
%TLSCDF T location-scale cumulative distribution function (cdf).
sigma(sigma <= 0) = NaN;
nu(nu <= 0) = NaN;

p = tcdf((x - mu)./sigma,nu);

function x = tlsinv(p, mu, sigma, nu)
%TLSINV Inverse of the t location-scale cumulative distribution function (cdf).
sigma(sigma <= 0) = NaN;
nu(nu <= 0) = NaN;

x = tinv(p,nu).*sigma + mu;

function r = tlsrnd(mu, sigma, nu, varargin)
%TLSRND Random arrays from the t location-scale distribution.
sigma(sigma <= 0) = NaN;
nu(nu <= 0) = NaN;

[err, sizeOut] = statsizechk(3,mu,sigma,nu,varargin{:});
if err > 0
    error('stats:tlsrnd:InconsistentSizes','Size information is inconsistent.');
end

r = mu + sigma.*trnd(nu,sizeOut);

function [m,v] = tlsstat(mu, sigma, nu)
%TLSSTAT Mean and variance for the t location-scale distribution.
sigma(sigma <= 0) = NaN;
nu(nu <= 0) = NaN;

if nu <= 1
    m = NaN;
else
    m = mu;
end
if nu <= 2
    v = Inf;
else
    v = sigma.^2 .* nu ./ (nu - 2);
end

function [nlogL,acov] = tlslike(params,data,cens,freq)
%TLSLIKE Negative log-likelihood for the t location-scale distribution.
if nargin < 4 || isempty(freq), freq = ones(size(data)); end
if nargin < 3 || isempty(cens), cens = zeros(size(data)); end

nlogL = tls_nloglf(params, data, cens, freq);
if nargout > 1
    acov = mlecov(params, data, 'nloglf',@tls_nloglf, 'cens',cens, 'freq',freq);
end

% ==== t location-scale fitting functions ====

function [phat,pci] = tlsfit(x,alpha,cens,freq,opts)

if nargin < 2 || isempty(alpha), alpha = .05; end
if nargin < 3 || isempty(cens), cens = zeros(size(x)); end
if nargin < 4 || isempty(freq), freq = ones(size(x)); end
if nargin < 5, opts = []; end

% Robust estimators for the mean and std dev of a normal, and method
% of moments on t-kurtosis for nu
xunc = x(cens == 0);
k = max(kurtosis(xunc), 4);
start = [median(xunc), 1.253.*mad(xunc), 2.*(2.*k-3)./(k-3)];

% The default options include turning fminsearch's display off.  This
% function gives its own warning/error messages, and the caller can turn
% display on to get the text output from fminsearch if desired.
options = statset(statset('tlsfit'), opts);
tolBnd = options.TolBnd;
options = optimset(options);

% Maximize the log-likelihood with respect to mu, sigma, and nu.
[phat,nll,err,output] = ...
    fminsearch(@tls_nloglf, start, options, x, cens, freq, tolBnd);
if (err == 0)
    % fminsearch may print its own output text; in any case give something
    % more statistical here, controllable via warning IDs.
    if output.funcCount >= options.MaxFunEvals
        wmsg = 'Maximum likelihood estimation did not converge.  Function evaluation limit exceeded.';
    else
        wmsg = 'Maximum likelihood estimation did not converge.  Iteration limit exceeded.';
    end
    if phat(3) > 100 % degrees of freedom became very large
       wmsg = sprintf('%s\n%s', wmsg, ...
                      'The normal distribution might provide a better fit.');
    end
    warning('stats:tlsfit:IterOrEvalLimit',wmsg);
elseif (err < 0)
    error('stats:tlsfit:NoSolution',...
          'Unable to reach a maximum likelihood solution.');
end

if nargout > 1
    acov = mlecov(phat, x, 'nloglf',@tls_nloglf, 'cens',cens, 'freq',freq);
    probs = [alpha/2; 1-alpha/2];
    se = sqrt(diag(acov))';

    % Compute the CI for mu using a normal approximation for muhat.
    pci(:,1) = norminv(probs, phat(1), se(1));

    % Compute the CI for sigma using a normal approximation for
    % log(sigmahat), and transform back to the original scale.
    % se(log(sigmahat)) is se(sigmahat) / sigmahat.
    logsigci = norminv(probs, log(phat(2)), se(2)./phat(2));
    pci(:,2) = exp(logsigci);

    % Compute the CI for nu using a normal distribution for nuhat.
    pci(:,3) = norminv(probs, phat(3), se(3));
end

function nll = tls_nloglf(parms, x, cens, freq, tolBnd)
%TLS_NLOGLF Objective function for t location-scale maximum likelihood.
mu = parms(1);
sigma = parms(2);
nu = parms(3);

% Restrict sigma and nu to the open interval (0, Inf).
if nargin > 4
    if sigma < tolBnd || nu < tolBnd
        nll = Inf;
        return
    end
end

t = (x - mu) ./ sigma;
w = nu + (t.^2);
logw = log(w);

L = -.5.*(nu+1).*logw + gammaln(.5.*(nu+1)) - gammaln(.5.*nu) + 0.5.*nu.*log(nu) - log(sigma) - .5.*log(pi);
ncen = sum(freq.*cens);
if ncen > 0
    cen = (cens == 1);
    if nu < 1e7  % Use the standard formula
        Scen = betainc(nu ./ w(cen), .5.*nu, 0.5) ./ 2;

        % Reflect for negative t.
        reflect = (t(cen) < 0);
        Scen(reflect) = 1 - Scen(reflect);

    else  % Use a normal approximation.
        Scen = log(0.5 * erfc(t(cen) ./ sqrt(2)));
    end
    L(cen) = log(Scen);
end
nll = -sum(freq .* L);

% Don't yet have dbetainc, so can't compute an analytic gradient with censoring.
%
% if nargout > 1
%     dL1 = (nu+1).*t./(w.*sigma);
%     dL2 = t.*dL1 - 1./sigma;
%     dL3 = .5.*(-logw - (nu+1)./w + psi(.5.*(nu+1)) - psi(.5.*nu) + log(nu) + 1);
%     if ncen > 0
% %         dL1(cen) = ;
% %         dL2(cen) = ;
% %         dL3(cen) = ;
%     end
%     ngrad = -[sum(freq .* dL1) sum(freq .* dL2) sum(freq .* dL3)];
% end

author: zouchuhang

repo: LayoutNet


text:
%% Machine Learning Online Class - Exercise 3 | Part 1: One-vs-all

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     lrCostFunction.m (logistic regression cost function)
%     oneVsAll.m
%     predictOneVsAll.m
%     predict.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Setup the parameters you will use for this part of the exercise
input_layer_size  = 400;  % 20x20 Input Images of Digits
num_labels = 10;          % 10 labels, from 1 to 10
                          % (note that we have mapped "0" to label 10)

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  You will be working with a dataset that contains handwritten digits.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

load('ex3data1.mat'); % training data stored in arrays X, y
m = size(X, 1);

% Randomly select 100 data points to display
rand_indices = randperm(m);
sel = X(rand_indices(1:100), :);

displayData(sel);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============ Part 2a: Vectorize Logistic Regression ============
%  In this part of the exercise, you will reuse your logistic regression
%  code from the last exercise. You task here is to make sure that your
%  regularized logistic regression implementation is vectorized. After
%  that, you will implement one-vs-all classification for the handwritten
%  digit dataset.
%

% Test case for lrCostFunction
fprintf('\nTesting lrCostFunction() with regularization');

theta_t = [-2; -1; 1; 2];
X_t = [ones(5,1) reshape(1:15,5,3)/10];
y_t = ([1;0;1;0;1] >= 0.5);
lambda_t = 3;
[J grad] = lrCostFunction(theta_t, X_t, y_t, lambda_t);

fprintf('\nCost: %f\n', J);
fprintf('Expected cost: 2.534819\n');
fprintf('Gradients:\n');
fprintf(' %f \n', grad);
fprintf('Expected gradients:\n');
fprintf(' 0.146561\n -0.548558\n 0.724722\n 1.398003\n');

fprintf('Program paused. Press enter to continue.\n');
pause;
%% ============ Part 2b: One-vs-All Training ============
fprintf('\nTraining One-vs-All Logistic Regression...\n')

lambda = 0.1;
[all_theta] = oneVsAll(X, y, num_labels, lambda);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 3: Predict for One-Vs-All ================

pred = predictOneVsAll(all_theta, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);

author: Ayatans

repo: Machine-Learning-homework


text:
function headmodel = ft_headmodel_duneuro(mesh, varargin)

% FT_HEADMODEL_DUNEURO creates a volume conduction model of the head
% using the finite element method (FEM) for EEG and MEG. Different source models
% are implemented, including the St. Venant, the subtraction and partial
% integration model. This function takes as input a mesh with tetrahedral
% or hexahedral elements and corresponding conductivities and returns
% as output a volume conduction model which can be used to compute EEG/MEG
% leadfields.
%
% Use as
%   headmodel = ft_headmodel_duneuro(mesh,'conductivity', conductivities, ...)
%   headmodel = ft_headmodel_duneuro(mesh,'grid_filename', grid_filename, 'tensors_filename', tensors_filename, ...)
%
% Required input arguments should be specified in key-value pairs and have
% to include either
%   grid_filename   = string, filename for grid in "msh" fileformat (see here: https://gmsh.info/doc/texinfo/gmsh.html#File-formats)
%   tensors_filename= string, filename for conductivities, txt file with conductivity values
%
% or
%   conductivity    = vector, conductivity values for tissues
%
% Optional input arguments are passed with
%   duneuro_settings = (optional) struct, which can contain the following fields
%
%   type            = string, 'fitted' (default)
%   solver_type     = string, 'cg' (default)
%   electrodes      = string, 'closest_subentity_center' (default)
%   subentities     = string, e.g. '1 2 3' (default) or '3'
%   forward         = string, 'venant' (default), 'partial_integration'
%   intorderadd     = string, e.g. '2' (default)
%   intorderadd_lb  = string, e.g. '2' (default)
%   initialization  = string, e.g. 'closest_vertex' (default)
%   numberOfMoments = string, e.g. '3' (default)
%   referenceLength = string, e.g. '20' (default)
%   relaxationFactor= string, e.g. '1e-6' (default)
%   restrict        = string, e.g. 'true' (default)
%   weightingExponent= string, e.g. '1' (default)
%   post_process    = string, e.g. 'true' (default)
%   subtract_mean   = string, e.g. 'true' (default)
%   reduction       = string, e.g. '1e-10' (default)
%   intorderadd_meg = integer, e.g.'0' (default)
%   mixedMoments    = logical, e.g. 'true' (default)
%   meg_type        = string, e.g. 'physical' (default)
%   meg_eneablecache= logical, e.g. 'false (default)

ft_hastoolbox('duneuro', 1);

grid_filename   = ft_getopt(varargin, 'grid_filename');
tensors_filename= ft_getopt(varargin, 'tensors_filename');
conductivity    = ft_getopt(varargin, 'conductivity');
duneuro_settings = ft_getopt(varargin, 'duneuro_settings');

% get the optional arguments and set defaults
if(isfield(duneuro_settings, 'type')) type = duneuro_settings.type;
else type = 'fitted';
end

if(isfield(duneuro_settings, 'solver_type')) solver_type = duneuro_settings.solver_type;
else solver_type = 'cg';
end

if(isfield(duneuro_settings, 'electrodes')) electrodes = duneuro_settings.electrodes;
else electrodes = 'closest_subentity_center';
end

if(isfield(duneuro_settings, 'subentities')) subentities = duneuro_settings.subentities;
else subentities = '1 2 3';
end

if(isfield(duneuro_settings, 'forward')) forward = duneuro_settings.forward;
else forward = 'venant';
end

if(isfield(duneuro_settings, 'intorderadd')) intorderadd = duneuro_settings.intorderadd;
else intorderadd = '2';
end

if(isfield(duneuro_settings, 'intorderadd_lb')) intorderadd_lb = duneuro_settings.intorderadd_lb;
else intorderadd_lb = '2';
end

if(isfield(duneuro_settings, 'initialization')) initialization = duneuro_settings.initialization;
else initialization = 'closest_vertex';
end

if(isfield(duneuro_settings, 'numberOfMoments')) numberOfMoments = duneuro_settings.numberOfMoments;
else numberOfMoments = '3';
end

if(isfield(duneuro_settings, 'referenceLength')) referenceLength = duneuro_settings.referenceLength;
else referenceLength = '20';
end

if(isfield(duneuro_settings, 'relaxationFactor')) relaxationFactor = duneuro_settings.relaxationFactor;
else relaxationFactor = '1e-6';
end

if(isfield(duneuro_settings, 'restrict')) restrict = duneuro_settings.restrict;
else restrict = 'true';
end

if(isfield(duneuro_settings, 'weightingExponent')) weightingExponent = duneuro_settings.weightingExponent;
else weightingExponent = '1';
end

if(isfield(duneuro_settings, 'post_process')) post_process = duneuro_settings.post_process;
else post_process = 'true';
end

if(isfield(duneuro_settings, 'subtract_mean')) subtract_mean = duneuro_settings.subtract_mean;
else subtract_mean = 'true';
end

if(isfield(duneuro_settings, 'reduction')) reduction = duneuro_settings.reduction;
else reduction = '1e-15';
end

if(isfield(duneuro_settings, 'intorderadd_meg')) intorderadd_meg = duneuro_settings.intorderadd_meg;
else intorderadd_meg = '2';
end

if(isfield(duneuro_settings, 'mixedMoments')) mixedMoments = duneuro_settings.mixedMoments;
else mixedMoments = 'true';
end

if(isfield(duneuro_settings, 'meg_type')) reduction = duneuro_settings.meg_type;
else meg_type = 'physical';
end

if(isfield(duneuro_settings, 'meg_eneablecache')) meg_eneablecache = duneuro_settings.meg_eneablecache;
else meg_eneablecache = 'false';
end

% start with an empty volume conductor
mesh = ft_datatype_parcellation(mesh);
headmodel = [];
if isfield(mesh,'pos')
  headmodel.pos = mesh.pos;
else
  error('Vertex field is required!')
end

if isfield(mesh,'tet')
  headmodel.tet = mesh.tet;
elseif isfield(mesh,'hex')
  headmodel.hex = mesh.hex;
  %%assume fieldtrip ordering, reorder to conform to DUNE standard
  headmodel.hex(:,3) = mesh.hex(:,4);
  headmodel.hex(:,4) = mesh.hex(:,3);
  headmodel.hex(:,7) = mesh.hex(:,8);
  headmodel.hex(:,8) = mesh.hex(:,7);
else
  error('Element field with tetrahedral or hexahedral information is required!')
end

if isfield(mesh,'tissue')
  headmodel.tissue = mesh.tissue;
else
  error('No tissue labels declared!')
end

if (isempty(grid_filename) || isempty(tensors_filename)) && (isempty(conductivity))
  error('Either a filename for the grid and a filename for the conductivity or a conductivity array must be supplied')
end

if ~isempty(conductivity)
  if (size(conductivity,1) ~= 1 && size(conductivity,1) ~= 9)  || size(conductivity,2) ~= length(unique(headmodel.tissue))
    error('Conductivity has wrong dimension!')
  end
end

if ~isfield(mesh,'tissuelabel')
  numlabels = size(unique(mesh.tissue),1);
  headmodel.tissuelabel = {};
  ulabel = unique(labels);
  for i = 1:numlabels
    headmodel.tissuelabel{i} = num2str(ulabel(i));
  end
else
  headmodel.tissuelabel = mesh.tissuelabel;
end

% create driver object
cfg                 = [];
cfg.type            = type;
cfg.solver_type     = solver_type;
cfg.meg.intorderadd = intorderadd_meg;
cfg.meg.type        = meg_type;
cfg.meg.cache.enable = meg_eneablecache;

if isfield(mesh,'tet')
  cfg.element_type = 'tetrahedron';
elseif isfield(mesh,'hex')
  cfg.element_type = 'hexahedron';
end

if(~isempty(grid_filename) && (~isempty(tensors_filename)))
  cfg.volume_conductor.grid.filename = grid_filename;
  cfg.volume_conductor.tensors.filename = tensors_filename;
else
  if isfield(headmodel,'tet')
    cfg.volume_conductor.grid.elements = (uint64(headmodel.tet))' - 1;
  elseif isfield(headmodel,'hex')
    cfg.volume_conductor.grid.elements = (uint64(headmodel.hex))' - 1;
  end
  cfg.volume_conductor.grid.nodes = headmodel.pos';

  %make sure labels start at 0 %TODO: faster for 6C anisotropy
  utissue = sort(unique(headmodel.tissue));
  for i = 1:size(headmodel.tissue,1)
    headmodel.tissue(i,1) =  find(utissue == headmodel.tissue(i,1));
  end
  cfg.volume_conductor.tensors.labels = uint64(headmodel.tissue -1);
  if(size(conductivity,1) == 1)
    cfg.volume_conductor.tensors.conductivities = conductivity;
  elseif (size(conductivity,1) == 9)
    cfg.volume_conductor.tensors.tensors = conductivity;
  end
end

try
    headmodel.driver = duneuro_meeg(cfg);
catch
    warning('To generate a valid Mex File, please compile duneuro-matlab on your machine following the instructions on the DUNEuro Wiki page.');
    return
end

headmodel.type              = 'duneuro';
headmodel.forward           = forward;
headmodel.electrodes        = electrodes;
headmodel.subentities       = subentities;
headmodel.intorderadd       = intorderadd;
headmodel.intorderadd_lb    = intorderadd_lb;
headmodel.initialization    = initialization;
headmodel.numberOfMoments   = numberOfMoments;
headmodel.referenceLength   = referenceLength;
headmodel.relaxationFactor  = relaxationFactor;
headmodel.restrict          = restrict;
headmodel.weightingExponent = weightingExponent;
headmodel.mixedMoments      = mixedMoments;
headmodel.post_process      = post_process;
headmodel.subtract_mean     = subtract_mean;
headmodel.reduction         = reduction;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function results = PTKMultipleRegionAnalysis(slice_bins, roi, context_mask, context_no_airways, distance_label, reporting)
    % PTKMultipleRegionAnalysis. Calculates metrics based on bins
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.

    results = PTKMetrics.empty;

    if isempty(context_mask) || ~context_mask.ImageExists
        return;
    end

    bin_image = slice_bins.BinImage;
    bin_regions = slice_bins.BinRegions;

    % Reduce all images to a consistent size
    bin_image.CropToFit;
    roi.ResizeToMatch(bin_image);
    context_mask.ResizeToMatch(bin_image);
    context_no_airways.ResizeToMatch(bin_image);

    reporting.UpdateProgressAndMessage(0, 'Calculating metrics for each bin');

    % Iterate over each bin
    for bin_index = 1 : numel(bin_regions)
        reporting.UpdateProgressStage(bin_index, numel(bin_regions));
        region = bin_regions(bin_index);
        bin_number = region.RegionNumber;
        bin_colour_index = region.ColormapIndex;
        distance_from_origin = region.DistanceFromOrigin;
        coordinates = region.Coordinates;

        % Create a mask for this bin
        mask = bin_image.BlankCopy;
        mask.ChangeRawImage(bin_image.RawImage == bin_colour_index & context_mask.RawImage);
        mask.CropToFit;

        % Create a mask for this bin excluding the airways
        no_airways_mask = bin_image.BlankCopy;
        no_airways_mask.ChangeRawImage(bin_image.RawImage == bin_colour_index & context_no_airways.RawImage);
        no_airways_mask.ResizeToMatch(mask);

        roi_reduced = roi.Copy;
        roi_reduced.ResizeToMatch(mask);
        bin_results = PTKComputeAirTissueFraction(roi_reduced, mask, reporting);
        [emphysema_results, ~] = PTKComputeEmphysemaFromMask(roi_reduced, no_airways_mask);
        bin_results.Merge(emphysema_results);
        bin_results.AddMetric('RegionNumber', bin_number, 'Region Number');
        bin_results.AddMetric('DistanceFromLungBaseMm', distance_from_origin, distance_label);
        bin_results.AddMetric('RegionCoordinateX', coordinates.CoordX, 'Region Coordinates X (mm)');
        bin_results.AddMetric('RegionCoordinateY', coordinates.CoordY, 'Region Coordinates Y (mm)');
        bin_results.AddMetric('RegionCoordinateZ', coordinates.CoordZ, 'Region Coordinates Z (mm)');

        results(end + 1) = bin_results;
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%%%%%%%%%%% author: Szymon Leja 2008-01-11
function [y]=dec2q15(x,form)

option = {'mode' , 'roundmode', 'overflowmode', 'format'};
value   = {'fixed', 'ceil'     , 'saturate'    ,  [16 15]};
q = quantizer( option, value ) ;

if form=='bin'
    y=num2bin(q,x);
end;

if form=='hex'
    y=num2hex(q,x);
end;

end


text:
function f = cumsum3(f)
%CUMSUM3   Triple indefinite integral of a CHEBFUN3.
%   F = CUMSUM3(F) returns the triple indefinite integral of a CHEBFUN3.
%   That is
%                  z  y  x
%                 /  /  /
%   CUMSUM3(F) = |  |  |   F(x,y,z) dx dy dz
%               /  /  /
%              e  c  a
%
%   where [a,b] x [c,d] x [e,g] is the domain of F.
%
% See also CHEBFUN3/CUMSUM, CHEBFUN3/CUMSUM2, CHEBFUN3/SUM, CHEBFUN3/SUM2
% and CHEBFUN3/SUM3.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Check for empty:
if ( isempty(f) )
    f = [];
    return
end

f.cols = cumsum(f.cols);     % CUMSUM along the cols.
f.rows = cumsum(f.rows);     % CUMSUM along the rows.
f.tubes = cumsum(f.tubes);   % CUMSUM along the tubes.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear all
m=40; n = 400; nrhs=1; A=randn(m,n); b=randn(m,nrhs);

fprintf('\nUnderdetermined system of size %d-by-%d, %d right-hand side(s) \n',m,n,nrhs);

fprintf('     \n');
fprintf('     m=40; n = 400; nrhs=1;\n');
fprintf('     A=randn(m,n); b=randn(m,nrhs);\n');
fprintf('     \n');
fprintf('     [x]=gels(A,b) \n');
fprintf('     [x]=gels(A,b,''svd'') \n');
fprintf('     [x]=gels(A,b,''dc'') \n');
fprintf('     [x]=gels(A,b,''y'') \n');
fprintf('     \n');

tic; [x]=gels(A,b); time=toc;
fprintf('     GELS         :: time =%4.2fs\tnorm(x-A''*(A''\\x))=%4.2e\tnorm(A*x-b)=%4.2e\n',time,norm(x-A'*(A'\x)),norm(A*x-b));

tic; [x]=gels(A,b,'svd'); time=toc;
fprintf('     GELS SVD     :: time =%4.2fs\tnorm(x-A''*(A''\\x))=%4.2e\tnorm(A*x-b)=%4.2e\n',time,norm(x-A'*(A'\x)),norm(A*x-b));

tic; [x]=gels(A,b,'dc'); time=toc;
fprintf('     GELS DC      :: time =%4.2fs\tnorm(x-A''*(A''\\x))=%4.2e\tnorm(A*x-b)=%4.2e\n',time,norm(x-A'*(A'\x)),norm(A*x-b));

tic; [x]=gels(A,b,'y'); time=toc;
fprintf('     GELS Y       :: time =%4.2fs\tnorm(x-A''*(A''\\x))=%4.2e\tnorm(A*x-b)=%4.2e\n',time,norm(x-A'*(A'\x)),norm(A*x-b));

tic; x=A\b; time=toc;
fprintf('     Matlab \\     :: time =%4.2fs\t                          \tnorm(A*x-b)=%4.2e\n',time,norm(A*x-b));

clear all
m=400; n = 40; nrhs=1; A=randn(m,n); b=randn(m,nrhs);

fprintf('\nOverdetermined system of size %d-by-%d, %d right-hand side(s) \n',m,n,nrhs);

fprintf('     \n');
fprintf('     m=%d; n=%d; nrhs=%d;\n',m,n,nrhs);
fprintf('     A=randn(m,n); b=randn(m,nrhs);\n');
fprintf('     \n');
fprintf('     [x]=gels(A,b) \n');
fprintf('     [x]=gels(A,b,''svd'') \n');
fprintf('     [x]=gels(A,b,''dc'') \n');
fprintf('     [x]=gels(A,b,''y'') \n');
fprintf('     \n');

tic; [x]=gels(A,b); time=toc;
fprintf('     GELS         :: time =%f\tnorm(x-A''*(A''\\x))=%e\tnorm(A''*(A*x-b))=%e\n',time,norm(x-A'*(A'\x)),norm(A'*(A*x-b)));

tic; [x]=gels(A,b,'svd'); time=toc;
fprintf('     GELS SVD     :: time =%f\tnorm(x-A''*(A''\\x))=%e\tnorm(A''*(A*x-b))=%e\n',time,norm(x-A'*(A'\x)),norm(A'*(A*x-b)));

tic; [x]=gels(A,b,'dc'); time=toc;
fprintf('     GELS DC      :: time =%f\tnorm(x-A''*(A''\\x))=%e\tnorm(A''*(A*x-b))=%e\n',time,norm(x-A'*(A'\x)),norm(A'*(A*x-b)));

tic; [x]=gels(A,b,'y'); time=toc;
fprintf('     GELS Y       :: time =%f\tnorm(x-A''*(A''\\x))=%e\tnorm(A''*(A*x-b))=%e\n',time,norm(x-A'*(A'\x)),norm(A'*(A*x-b)));

tic; x=A\b; time=toc;
fprintf('     Matlab \\     :: time =%f\tnorm(x-A''*(A''\\x))=%e\tnorm(A''*(A*x-b))=%e\n',time,norm(x-A'*(A'\x)),norm(A'*(A*x-b)));
fprintf('     \n');

repo: Aero-matlab


text:
function machar_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests R8_MACHAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  R8_MACHAR computes double\n' );
  fprintf ( 1, '  precision machine constants.\n' );

  [ ibeta, it, irnd, ngrd, machep, negep, iexp, ...
    minexp, maxexp, eps, epsneg, xmin, xmax ] = r8_machar ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  IBETA is the internal base for machine arithmetic.\n' );
  fprintf ( 1, '    IBETA =  %d\n', ibeta );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  IT is the number of digits, base IBETA, in the\n' );
  fprintf ( 1, '  floating point significand.\n' );
  fprintf ( 1, '    IT =     %d\n', it );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  IRND reports on floating point addition rounding:\n' );
  fprintf ( 1, '  0, for chopping;\n' );
  fprintf ( 1, '  1, for non-IEEE rounding;\n' );
  fprintf ( 1, '  2, for IEEE rounding;\n' );
  fprintf ( 1, '  3, for chopping with partial underflow;\n' );
  fprintf ( 1, '  4, for non-IEEE rounding with partial underflow.\n' );
  fprintf ( 1, '  5, for IEEE rounding with partial underflow.\n' );
  fprintf ( 1, '    IRND =   %d\n', irnd );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NGRD is the number of guard digits for floating point\n' );
  fprintf ( 1, '  multiplication with truncating arithmetic.\n' );
  fprintf ( 1, '    NGRD =   %d\n', ngrd );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  MACHEP is the largest negative integer such that\n' );
  fprintf ( 1, '  1.0 < 1.0 + BETA^MACHEP.\n' );
  fprintf ( 1, '    MACHEP = %d\n', machep );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NEGEPS is the largest negative integer such that\n' );
  fprintf ( 1, '  1.0 - BETA^NEGEPS < 1.0:\n' );
  fprintf ( 1, '    NEGEP =  %d\n', negep );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  IEXP is the number of bits reserved for the exponent\n' );
  fprintf ( 1, '  of a floating point number:\n' );
  fprintf ( 1, '    IEXP =   %d\n', iexp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  MINEXP is the most negative power of BETA such that\n' );
  fprintf ( 1, '  BETA^MINEXP is positive and normalized.\n' );
  fprintf ( 1, '    MINEXP = %d\n', minexp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  MAXEXP is the smallest positive power of BETA that\n' );
  fprintf ( 1, '  overflows:\n' );
  fprintf ( 1, '    MAXEXP = %d\n', maxexp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  EPS is a small positive floating point number\n' );
  fprintf ( 1, '  such that 1.0 < 1.0 + EPS.\n' );
  fprintf ( 1, '    EPS    = %26.16e\n', eps );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  EPSNEG is a small positive floating point number\n' );
  fprintf ( 1, '  such that 1.0 - EPSNEG < 1.0.\n' );
  fprintf ( 1, '    EPSNEG = %26.16e\n', epsneg );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  XMIN is the smallest positive normalized floating\n' );
  fprintf ( 1, '  point power of the radix:\n' );
  fprintf ( 1, '    XMIN =   %26.16e\n', xmin );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  XMAX is the largest finite floating point number:\n' );
  fprintf ( 1, '    XMAX   = %26.16e\n', xmax );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [DAT newFileName] = PTimport(csvFname, BBLFileName)
%% [DAT newFileName] = PTimport(filename, BBLFileName)
%  Imports log data in multiple formats (fw=selected firmware/logfile type). Default is .csv, but if using .bbl or .bfl files,

% ----------------------------------------------------------------------------------
% "THE BEER-WARE LICENSE" (Revision 42):
% <brian.white@queensu.ca> wrote this file. As long as you retain this notice you
% can do whatever you want with this stuff. If we meet some day, and you think
% this stuff is worth it, you can buy me a beer in return. -Brian White
% ----------------------------------------------------------------------------------

newFileName=csvFname;
try % speed things up a bit with this method
    hdr=detectImportOptions(csvFname);
    T=readtable(csvFname, 'HeaderLines', 2);
    T.Properties.VariableNames=hdr.VariableNames;
catch
    T=readtable(csvFname);
end

delete(csvFname);% we r done with this-dont wanna leave junk on main directory

% get header info from .bbl/.bfl file directly
a=strfind(BBLFileName,'.');
fid = fopen(BBLFileName, 'r');
c = textscan(fid, '%s', 'Delimiter','\n','CollectOutput',1); %c2 = native2unicode(c)
fclose(fid);

logStartPoints=[];
logEndPoints=[];

a=strfind(c{1},'Firmware version');
if isempty(find(cellfun(@(a)~isempty(a)&&a>0,a)))
    a=strfind(c{1},'Firmware revision');
    logStartPoints=find(cellfun(@(a)~isempty(a)&&a>0,a));
else
    a=strfind(c{1},'Firmware version');
    logStartPoints=find(cellfun(@(a)~isempty(a)&&a>0,a));
end

endStr={'throttle_boost_cutoff'; 'debug_mode' ; 'motor_pwm_rate' ; 'use_unsynced_pwm' ; 'gyro_32khz_hardware_lpf' ; 'gyro_hardware_lpf' ; 'yaw_deadband' ; 'deadband' ; 'pidsum_limit_yaw' ; 'pidsum_limit' ; 'energyCumulative (mAh)'; 'looptime'};% strings not in all revisions, try a few
m=1; a=strfind(c{1},endStr{m});
while isempty(find(cellfun(@(a)~isempty(a)&&a>0,a)))
    a=strfind(c{1},endStr{m});
    m=m+1;
end
logEndPoints = find(cellfun(@(a)~isempty(a)&&a>0,a));

relevantLogNum=str2num(csvFname(end-5:end-4));
s=c{1}(logStartPoints(relevantLogNum):logEndPoints(relevantLogNum));
n=1;
for m=1:size(s,1)
    if size(strsplit(char(s(m)),':'),2)==2
        SetupInfo(n,:)=(strsplit(char(s(m)),':'));
        a=char(SetupInfo(n,1));
        if strcmp(a(1:2),'H ')
            SetupInfo(n,1)={a(3:end)};
        end
        n=n+1;
    end
end

%
DAT.T=T;
DAT.mainFname=BBLFileName;
DAT.csvFname=csvFname;
DAT.SetupInfo=SetupInfo;
try
DAT.blackbox_decode_results=result;
catch
end

end

repo: PIDtoolbox


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM%
% 2D Multilevel LDDMM Example using stationary velocity field and
% diffusion regularizer.The example is described in detail in
% Section 4.2 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
% =========================================================================

close all; clear all; clc;
setup2Ddisc2CData

%% run affine pre-registration
imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',.1);

alpha = [4e2 0];
parametric = 0;
pad  = .5;
N    = 3;
mV     = @(m) ceil(1*m);
minLevel = 5;
maxLevel = 7;

%% run multilevel LDDMM registration

% 1) setup grid for velocities (padded)
omegaV = omega; omegaV(1:2:end) = omegaV(1:2:end)-pad;  omegaV(2:2:end) = omega(2:2:end)+pad;

% 2) setup regularizer
% regularizer('reset','regularizer','mfDiffusionCC','nt',0,'alpha',alpha,'HessianShift',1e-2);
regularizer('reset','regularizer','mfDiffusionCC','nt',0,'alpha',alpha,'HessianShift',1e-2);

NPIRpara         = optPara('NPIR-GN');
NPIRpara.maxIter = 40;
NPIRpara.scheme  = @GaussNewtonLDDMM;
[vc,~,wc,his] = MLLDDMM(ML,'minLevel',minLevel,'maxLevel',maxLevel,...
    'omegaV',omegaV,'mV',mV,'N',N,'parametric',parametric,'NPIRpara',NPIRpara,'plots',1);

%% show results
yc = getTrafoFromVelocityRK4(vc,getNodalGrid(omega,m),'omega',omegaV,'m',m,'tspan',[1,0],'N',N);
Topt = linearInterMex(dataT,omega,center(yc,m));
Jac = geometry(yc,m,'Jac','omega',omega);
D0  = distance(dataT(:),dataR(:),omega,m);
DOpt = distance(Topt(:),dataR(:),omega,m);

fig = figure(); clf;
fig.Name = sprintf('LDDMM Results: %s',mfilename);

subplot(2,3,1);
viewImage(dataR,omega,m);
title('reference');

subplot(2,3,4);
viewImage(dataT,omega,m);
hold on;
plotGrid(yc,omega,m,'spacing',4)
title('template');

subplot(2,3,2);
viewImage(Topt,omega,m);
title('T(yc)')

subplot(2,3,3);
viewImage(dataT(:)-dataR(:),omega,m);
title('init. residual, SSD=100%');

subplot(2,3,5);
viewImage2Dsc(Jac,omega,m);
title(sprintf('Jac, min=%1.2f max=%1.2f',min(Jac(:)),max(Jac(:))));

subplot(2,3,6);
viewImage(Topt(:)-dataR(:),omega,m);
title(sprintf('opt residual, SSD=%1.2f%%',100*DOpt/D0));


text:
function dx = tora_dx(u, x)

tora_lpv

p = [x(3) x(4)];

AB = querylpv(lpv, p);

% dx = A x + B u
dx = AB * [x(:); u];


text:
function test_suite=test_assert_less_than()
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_assert_less_than_exceptions
    assertExceptionThrown(@()assertLessThan([1 2],[1,2,3]),...
                          'assertLessThan:sizeNotEqual');
    assertExceptionThrown(@()assertLessThan([0,0],[1;1]),...
                          'assertLessThan:sizeNotEqual');
    assertExceptionThrown(@()assertLessThan([1],'a'), ...
                          'assertLessThan:classNotEqual');
    assertExceptionThrown(@()assertLessThan(1,1), ...
                          'assertLessThan:notLessThan');
    assertExceptionThrown(@()assertLessThan([1,2,3],[2,3,3]), ...
                          'assertLessThan:notLessThan');
    assertExceptionThrown(@()assertLessThan(nan,nan), ...
                          'assertLessThan:notLessThan');
    assertExceptionThrown(@()assertLessThan(nan,1), ...
                          'assertLessThan:notLessThan');
    assertExceptionThrown(@()assertLessThan(1,nan), ...
                          'assertLessThan:notLessThan');
    assertExceptionThrown(@()assertLessThan(ones(2,2,2,2),ones(2,2,2,2)*nan), ...
                          'assertLessThan:notLessThan');

    a = zeros(2,2,2,2);
    assertExceptionThrown(@()assertLessThan(ones(2,2,2,2),a), ...
                          'assertLessThan:notLessThan');
    a(1,1,1,1) = 1;
    assertExceptionThrown(@()assertLessThan(ones(2,2,2,2),a), ...
                          'assertLessThan:notLessThan');

function test_assert_less_than_passes
    assertLessThan(0,1);
    assertLessThan([0,0],[1,1]);
    assertLessThan(zeros(2,2,2,2),ones(2,2,2,2));


text:
% absolute change in body orientation
function [data,units] = compute_absdtheta(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  data{i} = abs(trx(fly).dtheta);
end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
function test_bug2197

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_selectdata

freq = [];
freq.labelcmb = {
  '1' '1'
  '1' '2'
  '1' '3'
  '1' '4'
  '1' '5'
  '1' '6'
  %  '2' '1'
  %  '2' '2'
  '2' '3'
  '2' '4'
  '2' '5'
  '2' '6'
  };
freq.cumtapcnt = ones(15,1);
freq.freq      = 1:7;
freq.crsspctrm = randn(15,10,7);
% freq.crsspctrmdimord = 'rpt_chancmb_freq'; % this is not needed

cfg = [];
cfg.avgoverrpt = 'yes';
output = ft_selectdata(cfg, freq);
assert(size(output.crsspctrm,1)==10)

cfg = [];
cfg.avgoverchancmb = 'yes';
output = ft_selectdata(cfg, freq);
assert(size(output.crsspctrm,2)==1)

cfg = [];
cfg.avgoverrpt = 'yes';
cfg.avgoverchancmb = 'yes';
output = ft_selectdata(cfg, freq);
assert(size(output.crsspctrm,1)==1)
assert(size(output.crsspctrm,2)==7) % frequencies

cfg = [];
cfg.avgoverrpt = 'yes';
cfg.avgoverchancmb = 'yes';
cfg.foilim = [3 5];
output = ft_selectdata(cfg, freq);
assert(size(output.crsspctrm,1)==1)
assert(size(output.crsspctrm,2)==3) % frequencies


text:
function g = whitehKernGradient(kern, x, varargin)

% WHITEHKERNGRADIENT Gradient of WHITEH kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% whiteh noise
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO whitehKernParamInit, kernGradient, whitehKernDiagGradient, kernGradX
%
% COPYRIGHT : Mauricio A. Alvarez, Neil D. Lawrence, 2009

% KERN

if nargin < 4
  g(1, 1) = sum(diag(varargin{end}).*(1./x(:,end)));
else
  g = 0;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function cfg = topoplot_common(cfg, varargin)

% TOPOPLOT_COMMON is shared by FT_TOPOPLOTTFR, FT_TOPOPLOTER and FT_TOPOPLOTIC, which
% serve as placeholder for the documentation and for the pre/postamble.

% Copyright (C) 2005-2011, F.C. Donders Centre
%
% 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$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPERS NOTE: This code is organized in a similar fashion for multiplot/singleplot/topoplot
% and for ER/TFR and should remain consistent over those 6 functions.
% Section 1: general cfg handling that is independent from the data
% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data
% Section 3: select the data to be plotted and determine min/max range
% Section 4: do the actual plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Section 1: general cfg handling that is independent from the data

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos
cfg = ft_checkconfig(cfg, 'unused',     {'cohtargetchannel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'cohrefchannel' 'refchannel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'zparam', 'parameter'});

% check for option-values to be renamed
cfg = ft_checkconfig(cfg, 'renamedval', {'electrodes',     'dotnum',      'numbers'});
cfg = ft_checkconfig(cfg, 'renamedval', {'zlim',           'absmax',      'maxabs'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedforward', 'outflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedback',    'inflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'highlight',      'yes',         'on'});

% check for renamed options
cfg = ft_checkconfig(cfg, 'renamed',     {'matrixside',    'directionality'});
cfg = ft_checkconfig(cfg, 'renamed',     {'electrodes',    'marker'});
cfg = ft_checkconfig(cfg, 'renamed',     {'emarker',       'markersymbol'});
cfg = ft_checkconfig(cfg, 'renamed',     {'ecolor',        'markercolor'});
cfg = ft_checkconfig(cfg, 'renamed',     {'emarkersize',   'markersize'});
cfg = ft_checkconfig(cfg, 'renamed',     {'efontsize',     'markerfontsize'});
cfg = ft_checkconfig(cfg, 'renamed',     {'hlmarker',      'highlightsymbol'});
cfg = ft_checkconfig(cfg, 'renamed',     {'hlcolor',       'highlightcolor'});
cfg = ft_checkconfig(cfg, 'renamed',     {'hlmarkersize',  'highlightsize'});
cfg = ft_checkconfig(cfg, 'renamed',     {'maplimits',     'zlim'});

% old ft_checkconfig adapted partially from topoplot.m (backwards backwards compatability)
cfg = ft_checkconfig(cfg, 'renamed',     {'grid_scale',    'gridscale'});
cfg = ft_checkconfig(cfg, 'renamed',     {'interpolate',   'interpolation'});
cfg = ft_checkconfig(cfg, 'renamed',     {'numcontour',    'contournum'});
cfg = ft_checkconfig(cfg, 'renamed',     {'electrod',      'marker'});
cfg = ft_checkconfig(cfg, 'renamed',     {'electcolor',    'markercolor'});
cfg = ft_checkconfig(cfg, 'renamed',     {'emsize',        'markersize'});
cfg = ft_checkconfig(cfg, 'renamed',     {'efsize',        'markerfontsize'});
cfg = ft_checkconfig(cfg, 'renamed',     {'headlimits',    'interplimits'});

% check for forbidden options
cfg = ft_checkconfig(cfg, 'forbidden',  {'hllinewidth', ...
  'headcolor', ...
  'hcolor', ...
  'hlinewidth', ...
  'contcolor', ...
  'outline', ...
  'highlightfacecolor', ...
  'showlabels'});

if ft_platform_supports('griddata-v4')
  default_interpmethod = 'v4';
else
  % Octave does not support 'v4', and 'cubic' is not yet implemented
  default_interpmethod = 'linear';
end

% set other config defaults
cfg.xlim              = ft_getopt(cfg, 'xlim',             'maxmin');
cfg.ylim              = ft_getopt(cfg, 'ylim',             'maxmin');
cfg.zlim              = ft_getopt(cfg, 'zlim',             'maxmin');
cfg.style             = ft_getopt(cfg, 'style',            'both');
cfg.gridscale         = ft_getopt(cfg, 'gridscale',         67);
cfg.interplimits      = ft_getopt(cfg, 'interplimits',     'head');
cfg.interpolation     = ft_getopt(cfg, 'interpolation',     default_interpmethod);
cfg.contournum        = ft_getopt(cfg, 'contournum',        6);
cfg.colormap          = ft_getopt(cfg, 'colormap',         'default');
cfg.colorbar          = ft_getopt(cfg, 'colorbar',         'no');
cfg.colorbartext      = ft_getopt(cfg, 'colorbartext',    '');
cfg.shading           = ft_getopt(cfg, 'shading',          'flat');
cfg.comment           = ft_getopt(cfg, 'comment',          'auto');
cfg.fontsize          = ft_getopt(cfg, 'fontsize',          8);
cfg.fontweight        = ft_getopt(cfg, 'fontweight',       'normal');
cfg.baseline          = ft_getopt(cfg, 'baseline',         'no'); % to avoid warning in timelock/freqbaseline
cfg.trials            = ft_getopt(cfg, 'trials',           'all', 1);
cfg.interactive       = ft_getopt(cfg, 'interactive',      'yes');
cfg.hotkeys           = ft_getopt(cfg, 'hotkeys',          'yes');
cfg.renderer          = ft_getopt(cfg, 'renderer',          []); % let MATLAB decide on the default
cfg.marker            = ft_getopt(cfg, 'marker',           'on');
cfg.markersymbol      = ft_getopt(cfg, 'markersymbol',     'o');
cfg.markercolor       = ft_getopt(cfg, 'markercolor',       [0 0 0]);
cfg.markersize        = ft_getopt(cfg, 'markersize',        2);
cfg.markerfontsize    = ft_getopt(cfg, 'markerfontsize',    8);
cfg.highlight         = ft_getopt(cfg, 'highlight',        'off');
cfg.highlightchannel  = ft_getopt(cfg, 'highlightchannel', 'all', 1); % highlight may be 'on', making highlightchannel {} meaningful
cfg.highlightsymbol   = ft_getopt(cfg, 'highlightsymbol',  '*');
cfg.highlightcolor    = ft_getopt(cfg, 'highlightcolor',    [0 0 0]);
cfg.highlightsize     = ft_getopt(cfg, 'highlightsize',     6);
cfg.highlightfontsize = ft_getopt(cfg, 'highlightfontsize', 8);
cfg.labeloffset       = ft_getopt(cfg, 'labeloffset',       0.005);
cfg.maskparameter     = ft_getopt(cfg, 'maskparameter',     []);
cfg.component         = ft_getopt(cfg, 'component',         []);
cfg.directionality    = ft_getopt(cfg, 'directionality',    []);
cfg.channel           = ft_getopt(cfg, 'channel',          'all');
cfg.refchannel        = ft_getopt(cfg, 'refchannel',        []);
cfg.figurename        = ft_getopt(cfg, 'figurename',        []);
cfg.interpolatenan    = ft_getopt(cfg, 'interpolatenan',   'yes');
cfg.commentpos        = ft_getopt(cfg, 'commentpos',       'layout');
cfg.scalepos          = ft_getopt(cfg, 'scalepos',         'layout');
cfg.figure            = ft_getopt(cfg, 'figure',           'yes');

% the user can either specify a single group of channels for highlighting
% which are all to be plotted in the same style, or multiple groups with a
% different style for each group. The latter is used by ft_clusterplot.
if iscell(cfg.highlightchannel) && ~isempty(cfg.highlightchannel) && ischar(cfg.highlightchannel{1})
  % it is a single cell-array with channels names, e.g. {'C1', 'Cz', 'C2'}
  cfg.highlightchannel = {cfg.highlightchannel};
elseif isnumeric(cfg.highlightchannel)
  % it is a numeric selection of channels, e.g. [1 2 3 4]
  cfg.highlightchannel = {cfg.highlightchannel};
elseif ischar(cfg.highlightchannel)
  % it is a single channel or single channel group, e.g. 'all'
  cfg.highlightchannel = {{cfg.highlightchannel}};
end
if ~iscell(cfg.highlight)
  cfg.highlight = {cfg.highlight};
end
if ~iscell(cfg.highlightsymbol)
  cfg.highlightsymbol = {cfg.highlightsymbol};
end
if ~iscell(cfg.highlightcolor)
  cfg.highlightcolor = {cfg.highlightcolor};
end
if ~iscell(cfg.highlightsize)
  cfg.highlightsize = {cfg.highlightsize};
end
if ~iscell(cfg.highlightfontsize)
  cfg.highlightfontsize = {cfg.highlightfontsize};
end
% make sure all cell-arrays for options are sufficiently long
ncellhigh = length(cfg.highlightchannel);
if length(cfg.highlight)          < ncellhigh,   cfg.highlight{ncellhigh}          = [];  end
if length(cfg.highlightsymbol)    < ncellhigh,   cfg.highlightsymbol{ncellhigh}    = [];  end
if length(cfg.highlightcolor)     < ncellhigh,   cfg.highlightcolor{ncellhigh}     = [];  end
if length(cfg.highlightsize)      < ncellhigh,   cfg.highlightsize{ncellhigh}      = [];  end
if length(cfg.highlightfontsize)  < ncellhigh,   cfg.highlightfontsize{ncellhigh}  = [];  end
% make sure all cell-arrays for options are not too long
cfg.highlight         (ncellhigh+1:end) = [];
cfg.highlightsymbol   (ncellhigh+1:end) = [];
cfg.highlightcolor    (ncellhigh+1:end) = [];
cfg.highlightsize     (ncellhigh+1:end) = [];
cfg.highlightfontsize (ncellhigh+1:end) = [];
% then default all empty cells
for icell = 1:ncellhigh
  if isempty(cfg.highlight{icell}),          cfg.highlight{icell} = 'on';          end
  if isempty(cfg.highlightsymbol{icell}),    cfg.highlightsymbol{icell} = 'o';     end
  if isempty(cfg.highlightcolor{icell}),     cfg.highlightcolor{icell} = [0 0 0];  end
  if isempty(cfg.highlightsize{icell}),      cfg.highlightsize{icell} = 6;         end
  if isempty(cfg.highlightfontsize{icell}),  cfg.highlightfontsize{icell} = 8;     end
end

% for backwards compatability
if strcmp(cfg.marker, 'highlights')
  ft_warning('cfg.marker=''highlights'' is no longer supported, please use cfg.highlight')
  cfg.marker = 'off';
end

% check if the colormap is in the proper format
if ~isequal(cfg.colormap, 'default')
  if ischar(cfg.colormap)
    cfg.colormap = ft_colormap(cfg.colormap);
  elseif iscell(cfg.colormap)
    cfg.colormap = ft_colormap(cfg.colormap{:});
  elseif isnumeric(cfg.colormap) && size(cfg.colormap,2)~=3
    ft_error('cfg.colormap must be Nx3');
  end
  % the actual colormap will be set below
end

Ndata = numel(varargin);
for indx=1:Ndata

  % open a new figure, or add it to the existing one
  open_figure(keepfields(cfg, {'figure', 'position', 'visible', 'renderer', 'figurename', 'title'}));

  % apply the same colormap to all figures
  if ~isempty(cfg.colormap)
    set(gcf,  'colormap', cfg.colormap);
  end

  if iscell(cfg.dataname)
    dataname = cfg.dataname{indx};
  else
    dataname = cfg.dataname;
  end

  data = varargin{indx};
  data = ft_checkdata(data, 'datatype', {'comp', 'timelock', 'freq'});

  %% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data

  dtype  = ft_datatype(data);
  hastime = isfield(data, 'time');

  % Set x/y/parameter defaults according to datatype and dimord
  switch dtype
    case 'timelock'
      xparam = ft_getopt(cfg, 'xparam', 'time');
      yparam = ft_getopt(cfg, 'yparam', '');
      if isfield(data, 'trial')
        cfg.parameter = ft_getopt(cfg, 'parameter', 'trial');
      elseif isfield(data, 'individual')
        cfg.parameter = ft_getopt(cfg, 'parameter', 'individual');
      elseif isfield(data, 'avg')
        cfg.parameter = ft_getopt(cfg, 'parameter', 'avg');
      end
    case 'freq'
      if hastime
        xparam = ft_getopt(cfg, 'xparam', 'time');
        yparam = ft_getopt(cfg, 'yparam', 'freq');
        cfg.parameter = ft_getopt(cfg, 'parameter', 'powspctrm');
      else
        xparam = 'freq';
        yparam = '';
        cfg.parameter = ft_getopt(cfg, 'parameter', 'powspctrm');
      end
    case 'comp'
      % Add a pseudo-axis with the component numbers
      data.comp = 1:size(data.topo,2);
      if ~isempty(cfg.component)
        % make a selection of components
        data.comp  = data.comp(cfg.component);
        data.topo  = data.topo(:,cfg.component);
        try, data.label     = data.label(cfg.component); end
        try, data.unmixing  = data.unmixing(cfg.component,:); end
      end
      % Rename the field with topographic label information
      data.label      = data.topolabel;
      data.topodimord = 'chan_comp';
      data = removefields(data, {'topolabel', 'unmixing', 'unmixingdimord'}); % not needed any more
      xparam = 'comp';
      yparam = '';
      cfg.parameter = ft_getopt(cfg, 'parameter', 'topo');
    otherwise
      % if the input data is not one of the standard data types, or if the functional
      % data is just one value per channel: in this case xparam, yparam are not defined
      % and the user should define the parameter
      if ~isfield(data, 'label'),     ft_error('the input data should at least contain a label-field');            end
      if ~isfield(cfg,  'parameter'), ft_error('the configuration should at least contain a ''parameter'' field'); end
      if ~isfield(cfg,  'xparam')
        cfg.xlim = [1 1];
        xparam   = '';
        yparam   = '';
      end
  end

  % check whether rpt/subj is present and remove if necessary
  dimord = getdimord(data, cfg.parameter);
  dimtok = tokenize(dimord, '_');
  hasrpt = any(ismember(dimtok, {'rpt' 'subj'}));

  if ~hasrpt
    assert(isequal(cfg.trials, 'all') || isequal(cfg.trials, 1), 'incorrect specification of cfg.trials for data without repetitions');
  else
    assert(~isempty(cfg.trials), 'empty specification of cfg.trials for data with repetitions');
  end

  % parse cfg.channel
  if isfield(cfg, 'channel') && isfield(data, 'label')
    cfg.channel = ft_channelselection(cfg.channel, data.label);
  elseif isfield(cfg, 'channel') && isfield(data, 'labelcmb')
    cfg.channel = ft_channelselection(cfg.channel, unique(data.labelcmb(:)));
  end

  % Apply baseline correction
  if ~strcmp(cfg.baseline, 'no')
    % keep mask-parameter if it is set
    if ~isempty(cfg.maskparameter)
      tempmask = data.(cfg.maskparameter);
    end
    tmpcfg = removefields(cfg, 'inputfile');
    if strcmp(xparam, 'time') && strcmp(yparam, 'freq')
      data = ft_freqbaseline(tmpcfg, data);
    elseif strcmp(xparam, 'time') && strcmp(yparam, '')
      data = ft_timelockbaseline(tmpcfg, data);
    end
    % put mask-parameter back if it is set
    if ~isempty(cfg.maskparameter)
      data.(cfg.maskparameter) = tempmask;
    end
  end

  % time and/or frequency should NOT be selected and averaged here, since a singleplot might follow in interactive mode
  tmpcfg = keepfields(cfg, {'channel', 'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  if hasrpt
    tmpcfg.avgoverrpt = 'yes';
  else
    tmpcfg.avgoverrpt = 'no';
  end
  tmpvar = data;
  ws = ft_warning('off', 'FieldTrip:getdimord:warning_dimord_could_not_be_determined');
  [data] = ft_selectdata(tmpcfg, data);
  ft_warning(ws);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

  if isfield(tmpvar, cfg.maskparameter) && ~isfield(data, cfg.maskparameter)
    % the mask parameter is not present after ft_selectdata, because it is
    % not included in all input arguments. Make the same selection and copy
    % it over
    tmpvar = ft_selectdata(tmpcfg, tmpvar);
    data.(cfg.maskparameter) = tmpvar.(cfg.maskparameter);
  end

  clear tmpvar tmpcfg dimord dimtok hastime hasfreq hasrpt

  % ensure that the preproc specific options are located in the cfg.preproc
  % substructure, but also ensure that the field 'refchannel' remains at the
  % highest level in the structure. This is a little hack by JM because the field
  % refchannel can relate to connectivity or to an EEG reference.

  if isfield(cfg, 'refchannel'), refchannelincfg = cfg.refchannel; cfg = rmfield(cfg, 'refchannel'); end
  cfg = ft_checkconfig(cfg, 'createsubcfg',  {'preproc'});
  if exist('refchannelincfg', 'var'), cfg.refchannel  = refchannelincfg; end

  if ~isempty(cfg.preproc)
    % preprocess the data, i.e. apply filtering, baselinecorrection, etc.
    fprintf('applying preprocessing options\n');
    if ~isfield(cfg.preproc, 'feedback')
      cfg.preproc.feedback = cfg.interactive;
    end
    data = ft_preprocessing(cfg.preproc, data);
  end

  % Handle the bivariate case
  dimord = getdimord(varargin{1}, cfg.parameter);
  if startsWith(dimord, 'chan_chan_') || startsWith(dimord, 'chancmb_')
    % convert the bivariate data to univariate and call the parent plotting function again
    s = dbstack;
    cfg.originalfunction = s(2).name;
    cfg.trials = 'all'; % trial selection has been taken care off
    bivariate_common(cfg, varargin{:});
    return
  end

  % Apply channel-type specific scaling
  tmpcfg = keepfields(cfg, {'parameter', 'chanscale', 'ecgscale', 'eegscale', 'emgscale', 'eogscale', 'gradscale', 'magscale', 'megscale', 'mychan', 'mychanscale'});
  data = chanscale_common(tmpcfg, data);

  %% Section 3: select the data to be plotted and determine min/max range

  dimord = getdimord(varargin{1}, cfg.parameter);
  dimtok = tokenize(dimord, '_');

  % Create time-series of small topoplots
  if ~ischar(cfg.xlim) && length(cfg.xlim)>2 %&& any(ismember(dimtok, 'time'))
    % Switch off interactive mode:
    cfg.interactive = 'no';
    xlims = cfg.xlim;
    % Iteratively call topoplotER with different xlim values:
    nplots = numel(xlims)-1;
    nyplot = ceil(sqrt(nplots));
    nxplot = ceil(nplots./nyplot);
    tmpcfg = removefields(cfg, 'inputfile');
    for i=1:length(xlims)-1
      tmpcfg.figure = subplot(nxplot, nyplot, i);
      tmpcfg.xlim = xlims(i:i+1);
      ft_topoplotTFR(tmpcfg, data);
    end
    return
  end

  % Get physical min/max range of x
  if ~isempty(xparam)
    if strcmp(cfg.xlim, 'maxmin')
      xmin = min(data.(xparam));
      xmax = max(data.(xparam));
    else
      xmin = cfg.xlim(1);
      xmax = cfg.xlim(2);
    end
    xminindx = nearest(data.(xparam), xmin);
    xmaxindx = nearest(data.(xparam), xmax);
    xmin = data.(xparam)(xminindx);
    xmax = data.(xparam)(xmaxindx);
    selx = xminindx:xmaxindx;
  end

  % Get physical min/max range of y
  if ~isempty(yparam)
    if strcmp(cfg.ylim, 'maxmin')
      ymin = min(data.(yparam));
      ymax = max(data.(yparam));
    else
      ymin = cfg.ylim(1);
      ymax = cfg.ylim(2);
    end
    yminindx = nearest(data.(yparam), ymin);
    ymaxindx = nearest(data.(yparam), ymax);
    ymin = data.(yparam)(yminindx);
    ymax = data.(yparam)(ymaxindx);
    sely = yminindx:ymaxindx;
  end

  % Take subselection of channels, this only works if the interactive mode is switched off
  if exist('selchannel', 'var')
    sellab = match_str(data.label, selchannel);
    label  = data.label(sellab);
  else
    sellab = 1:numel(data.label);
    label  = data.label;
  end

  % Make data vector with one scalar value for each channel
  dat = data.(cfg.parameter);
  % get dimord dimensions
  ydim = find(strcmp(yparam, dimtok));
  xdim = find(strcmp(xparam, dimtok));
  zdim = setdiff(1:ndims(dat), [ydim xdim]);
  % and permute
  dat = permute(dat, [zdim(:)' ydim xdim]);

  if ~isempty(yparam)
    % time-frequency data
    dat = dat(sellab, sely, selx);
    dat = nanmean(nanmean(dat, 3), 2);
  elseif ~isempty(cfg.component)
    % component data, nothing to do
  else
    % time or frequency data
    dat = dat(sellab, selx);
    dat = nanmean(dat, 2);
  end
  dat = dat(:);

  if isfield(data, cfg.maskparameter)
    % Make mask vector with one value for each channel
    msk = data.(cfg.maskparameter);
    % get dimord dimensions
    ydim = find(strcmp(yparam, dimtok));
    xdim = find(strcmp(xparam, dimtok));
    zdim = setdiff(1:ndims(dat), [ydim xdim]);
    % and permute
    msk = permute(msk, [zdim(:)' ydim xdim]);

    if ~isempty(yparam)
      % time-frequency data
      msk = msk(sellab, sely, selx);
    elseif ~isempty(cfg.component)
      % component data, nothing to do
    else
      % time or frequency data
      msk = msk(sellab, selx);
    end

    if size(msk,2)>1 || size(msk,3)>1
      ft_warning('no masking possible for average over multiple latencies or frequencies -> cfg.maskparameter cleared')
      msk = [];
    end

  else
    msk = [];
  end

  % Select the channels in the data that match with the layout:
  [seldat, sellay] = match_str(label, cfg.layout.label);
  if isempty(seldat)
    ft_error('labels in data and labels in layout do not match');
  end

  dat = dat(seldat);
  if ~isempty(msk)
    msk = msk(seldat);
  end

  % Select x and y coordinates and labels of the channels in the data
  chanX = cfg.layout.pos(sellay,1);
  chanY = cfg.layout.pos(sellay,2);

  % Get physical min/max range of z:
  if strcmp(cfg.zlim, 'maxmin')
    zmin = min(dat);
    zmax = max(dat);
  elseif strcmp(cfg.zlim, 'maxabs')
    zmin = -max(max(abs(dat)));
    zmax = max(max(abs(dat)));
  elseif strcmp(cfg.zlim, 'zeromax')
    zmin = 0;
    zmax = max(dat);
  elseif strcmp(cfg.zlim, 'minzero')
    zmin = min(dat);
    zmax = 0;
  else
    zmin = cfg.zlim(1);
    zmax = cfg.zlim(2);
  end

  % Construct comment
  switch cfg.comment
    case 'auto'
      comment = date;
      if ~isempty(xparam)
        comment = sprintf('%0s\n%0s=[%.3g %.3g]', comment, xparam, xmin, xmax);
      end
      if ~isempty(yparam)
        comment = sprintf('%0s\n%0s=[%.3g %.3g]', comment, yparam, ymin, ymax);
      end
      if ~isempty(cfg.parameter)
        comment = sprintf('%0s\n%0s=[%.3g %.3g]', comment, cfg.parameter, zmin, zmax);
      end
    case 'xlim'
      comment = '';
      if ~isempty(xparam)
        comment = sprintf('%0s=[%.3g %.3g]', xparam, xmin, xmax);
      end
    case 'ylim'
      comment = '';
      if ~isempty(yparam)
        comment = sprintf('%0s=[%.3g %.3g]', yparam, ymin, ymax);
      end
    case 'zlim'
      comment = sprintf('%0s=[%.3g %.3g]', cfg.parameter, zmin, zmax);
    otherwise
      comment = cfg.comment; % allow custom comments (e.g., ft_clusterplot specifies custom comments)
  end % switch comment

  if ~isempty(cfg.refchannel)
    if iscell(cfg.refchannel)
      comment = sprintf('%s\nreference=%s %s', comment, cfg.refchannel{:});
    else
      comment = sprintf('%s\nreference=%s %s', comment, cfg.refchannel);
    end
  end

  % Draw topoplot
  cla
  hold on

  % check for nans
  nanInds = isnan(dat);
  if strcmp(cfg.interpolatenan, 'yes') && any(nanInds)
    ft_warning('removing NaNs from the data');
    chanX(nanInds) = [];
    chanY(nanInds) = [];
    dat(nanInds)   = [];
    if ~isempty(msk)
      msk(nanInds) = [];
    end
  end

  % Set ft_plot_topo specific options
  if strcmp(cfg.interplimits, 'head')
    interplimits = 'mask';
  else
    interplimits = cfg.interplimits;
  end
  if strcmp(cfg.style, 'both');            style = 'surfiso';     end
  if strcmp(cfg.style, 'straight');        style = 'surf';        end
  if strcmp(cfg.style, 'contour');         style = 'iso';         end
  if strcmp(cfg.style, 'fill');            style = 'isofill';     end
  if strcmp(cfg.style, 'straight_imsat');  style = 'imsat';       end
  if strcmp(cfg.style, 'both_imsat');      style = 'imsatiso';    end

  % Draw plot
  if strcmp(cfg.style, 'blank')
    ft_plot_layout(cfg.layout, 'box', 'no', 'label', 'no', 'point', 'no')
  else
    opt = {};
    opt = ft_setopt(opt, 'interpmethod',  cfg.interpolation);
    opt = ft_setopt(opt, 'interplim',     interplimits);
    opt = ft_setopt(opt, 'gridscale',     cfg.gridscale);
    opt = ft_setopt(opt, 'outline',       cfg.layout.outline);
    opt = ft_setopt(opt, 'shading',       cfg.shading);
    opt = ft_setopt(opt, 'isolines',      cfg.contournum);
    opt = ft_setopt(opt, 'mask',          cfg.layout.mask);
    opt = ft_setopt(opt, 'style',         style);
    opt = ft_setopt(opt, 'datmask',       msk);
    if strcmp(style, 'imsat') || strcmp(style, 'imsatiso')
      opt = ft_setopt(opt, 'clim',  [zmin zmax]);
    end
    ft_plot_topo(chanX, chanY, dat, opt{:});
  end

  % For Highlight (channel-selection)
  for icell = 1:length(cfg.highlight)
    if ~strcmp(cfg.highlight{icell}, 'off')
      cfg.highlightchannel{icell} = ft_channelselection(cfg.highlightchannel{icell}, data.label);
      [dum, layoutindex] = match_str(cfg.highlightchannel{icell}, cfg.layout.label);
      templay = [];
      templay.outline = cfg.layout.outline;
      templay.mask    = cfg.layout.mask;
      templay.pos     = cfg.layout.pos(layoutindex,:);
      templay.width   = cfg.layout.width(layoutindex);
      templay.height  = cfg.layout.height(layoutindex);
      templay.label   = cfg.layout.label(layoutindex);
      if strcmp(cfg.highlight{icell}, 'labels') || strcmp(cfg.highlight{icell}, 'numbers')
        labelflg = 1;
      else
        labelflg = 0;
      end
      if strcmp(cfg.highlight{icell}, 'numbers')
        for ichan = 1:length(layoutindex)
          templay.label{ichan} = num2str(match_str(data.label, templay.label{ichan}));
        end
      end

      ft_plot_layout(templay, 'box', 'no', 'label', labelflg, 'point', ~labelflg, ...
        'pointsymbol',  cfg.highlightsymbol{icell}, ...
        'pointcolor',   cfg.highlightcolor{icell}, ...
        'pointsize',    cfg.highlightsize{icell}, ...
        'fontsize',     cfg.highlightfontsize{icell}, ...
        'labeloffset',  cfg.labeloffset, ...
        'labelalignh', 'center', ...
        'labelalignv', 'middle');
    end
  end % for icell

  % For Markers (all channels)
  switch cfg.marker
    case {'off', 'no'}
      % do not show the markers
    case {'on', 'labels', 'numbers'}
      channelsToMark = 1:length(data.label);
      channelsToHighlight = [];
      for icell = 1:length(cfg.highlight)
        if ~strcmp(cfg.highlight{icell}, 'off')
          channelsToHighlight = [channelsToHighlight; match_str(data.label, cfg.highlightchannel{icell})];
        end
      end
      if strcmp(cfg.interpolatenan, 'no')
        channelsNotMark = channelsToHighlight;
      else
        channelsNotMark = union(find(isnan(dat)), channelsToHighlight);
      end
      channelsToMark(channelsNotMark) = [];
      [dum, layoutindex] = match_str(ft_channelselection(channelsToMark, data.label), cfg.layout.label);
      templay = [];
      templay.outline = cfg.layout.outline;
      templay.mask    = cfg.layout.mask;
      templay.pos     = cfg.layout.pos(layoutindex,:);
      templay.width   = cfg.layout.width(layoutindex);
      templay.height  = cfg.layout.height(layoutindex);
      templay.label   = cfg.layout.label(layoutindex);
      if strcmp(cfg.marker, 'labels') || strcmp(cfg.marker, 'numbers')
        labelflg = 1;
      else
        labelflg = 0;
      end
      if strcmp(cfg.marker, 'numbers')
        for ichan = 1:length(layoutindex)
          templay.label{ichan} = num2str(match_str(data.label,templay.label{ichan}));
        end
      end
      ft_plot_layout(templay, 'box', 'no', 'label', labelflg, 'point', ~labelflg, ...
        'pointsymbol',  cfg.markersymbol, ...
        'pointcolor',   cfg.markercolor, ...
        'pointsize',    cfg.markersize, ...
        'fontsize',     cfg.markerfontsize, ...
        'labeloffset',  cfg.labeloffset, ...
        'labelalignh', 'center', ...
        'labelalignv', 'middle');
    otherwise
      ft_error('incorrect value for cfg.marker');
  end

  if isfield(cfg, 'vector')
    % FIXME this is not documented
    vecX = nanmean(real(data.(cfg.vector)(:,selx)), 2);
    vecY = nanmean(imag(data.(cfg.vector)(:,selx)), 2);

    % scale quiver relative to largest gradiometer sample
    k = 0.15/max([max(abs(real(data.(cfg.vector)(:)))) max(abs(imag(data.(cfg.vector)(:))))]);
    quiver(chanX, chanY, k*vecX, k*vecY, 0, 'red');
  end

  % Write comment
  if strcmp(cfg.comment, 'no')
    comment_handle = [];
  elseif strcmp(cfg.commentpos, 'title')
    comment_handle = title(comment, 'FontSize', cfg.fontsize);
  elseif ~isempty(strcmp(cfg.layout.label, 'COMNT'))
    x_comment = cfg.layout.pos(strcmp(cfg.layout.label, 'COMNT'), 1);
    y_comment = cfg.layout.pos(strcmp(cfg.layout.label, 'COMNT'), 2);
    % 'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom',
    comment_handle = ft_plot_text(x_comment, y_comment, comment, 'FontSize', cfg.fontsize, 'FontWeight', cfg.fontweight);
  else
    comment_handle = [];
  end

  % Set colour axis
  if ~strcmp(cfg.style, 'blank')
    caxis([zmin zmax]);
  end

  % Plot colorbar
  if isfield(cfg, 'colorbar')
    if strcmp(cfg.colorbar, 'yes')
      c = colorbar;
      ylabel(c, cfg.colorbartext);
    elseif ~strcmp(cfg.colorbar, 'no')
      c = colorbar('location', cfg.colorbar);
      ylabel(c, cfg.colorbartext);
    end
  end

  % set the figure window title, but only if the user has not changed it
  if isempty(get(gcf, 'Name'))
    if isfield(cfg, 'funcname')
      funcname = cfg.funcname;
    else
      funcname = mfilename;
    end
    if isempty(cfg.figurename)
      dataname_str = join_str(', ', dataname);
      set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), funcname, dataname_str));
      set(gcf, 'NumberTitle', 'off');
    else
      set(gcf, 'name', cfg.figurename);
      set(gcf, 'NumberTitle', 'off');
    end
  end

  axis off
  hold off
  axis equal

  if strcmp('yes', cfg.hotkeys)
    %  Attach data and cfg to figure and attach a key listener to the figure
    set(gcf, 'KeyPressFcn', {@key_sub, zmin, zmax})
  end

  % Make the figure interactive
  if strcmp(cfg.interactive, 'yes')
    % add the cfg/data/channel information to the figure under identifier linked to this axis
    ident                    = ['axh' num2str(round(sum(clock.*1e6)))]; % unique identifier for this axis
    set(gca, 'tag',ident);
    info                     = guidata(gcf);
    info.(ident).x           = cfg.layout.pos(:, 1);
    info.(ident).y           = cfg.layout.pos(:, 2);
    info.(ident).label       = cfg.layout.label;
    info.(ident).dataname    = dataname;
    info.(ident).cfg         = cfg;
    info.(ident).commenth    = comment_handle;
    if exist('linecolor', 'var')
      info.(ident).linecolor   = linecolor;
    end
    if ~isfield(info.(ident),'datvarargin')
      info.(ident).datvarargin = varargin(1:Ndata); % add all datasets to figure
    end
    info.(ident).datvarargin{indx} = data; % update current dataset (e.g. baselined, channel selection, etc)
    guidata(gcf, info);
    if any(strcmp(dimord, {'chan_time', 'chan_freq', 'subj_chan_time', 'rpt_chan_time', 'chan_chan_freq', 'chancmb_freq', 'rpt_chancmb_freq', 'subj_chancmb_freq'}))
      set(gcf, 'WindowButtonUpFcn',     {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotER}, 'event', 'WindowButtonUpFcn'});
      set(gcf, 'WindowButtonDownFcn',   {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotER}, 'event', 'WindowButtonDownFcn'});
      set(gcf, 'WindowButtonMotionFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotER}, 'event', 'WindowButtonMotionFcn'});
    elseif any(strcmp(dimord, {'chan_freq_time', 'subj_chan_freq_time', 'rpt_chan_freq_time', 'rpttap_chan_freq_time', 'chan_chan_freq_time', 'chancmb_freq_time', 'rpt_chancmb_freq_time', 'subj_chancmb_freq_time'}))
      set(gcf, 'WindowButtonUpFcn',     {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotTFR}, 'event', 'WindowButtonUpFcn'});
      set(gcf, 'WindowButtonDownFcn',   {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotTFR}, 'event', 'WindowButtonDownFcn'});
      set(gcf, 'WindowButtonMotionFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotTFR}, 'event', 'WindowButtonMotionFcn'});
    else
      ft_warning('unsupported dimord "%s" for interactive plotting', dimord);
    end
  end

end % for numel(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which is called after selecting channels in case of cfg.interactive='yes'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function select_singleplotER(label)
ident       = get(gca, 'tag');
info        = guidata(gcf);
cfg         = info.(ident).cfg;
datvarargin = info.(ident).datvarargin;
linecolor   = ft_getopt(info.(ident), 'linecolor', lineattributes_common(cfg, datvarargin{:}));

if ~isempty(label)
  cfg = removefields(cfg, 'inputfile');       % the reading has already been done and varargin contains the data
  cfg.baseline = 'no';                        % make sure the next function does not apply a baseline correction again
  cfg.channel = label;
  cfg.dataname = info.(ident).cfg.dataname;   % put data name in here, this cannot be resolved by other means
  cfg.trials = 'all';                         % trial selection has already been taken care of
  cfg.xlim = 'maxmin';
  % if user specified a zlim, copy it over to the ylim of singleplot
  if isfield(cfg, 'zlim')
    cfg.ylim = cfg.zlim;
    cfg = rmfield(cfg, 'zlim');
  end
  fprintf('selected cfg.channel = {%s}\n', join_str(', ', cfg.channel));
  % ensure that the new figure appears at the same position
  cfg.figure = 'yes';
  cfg.position = get(gcf, 'Position');

  selchan = match_str(datvarargin{1}.label, cfg.channel);
  cfg.linecolor = linecolor(selchan, :, :); % make a subselection for the correct inheritance of the line colors
  ft_singleplotER(cfg, datvarargin{:});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which is called after selecting channels in case of cfg.interactive='yes'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function select_singleplotTFR(label)
ident       = get(gca, 'tag');
info        = guidata(gcf);
cfg         = info.(ident).cfg;
datvarargin = info.(ident).datvarargin;
if ~isempty(label)
  cfg = removefields(cfg, 'inputfile');   % the reading has already been done and varargin contains the data
  cfg.baseline = 'no';                    % make sure the next function does not apply a baseline correction again
  cfg.channel = label;
  cfg.dataname = info.(ident).dataname;   % put data name in here, this cannot be resolved by other means
  cfg.trials = 'all';                     % trial selection has already been taken care of
  cfg.xlim = 'maxmin';
  cfg.ylim = 'maxmin';
  fprintf('selected cfg.channel = {%s}\n', join_str(', ', cfg.channel));
  % ensure that the new figure appears at the same position
  cfg.figure = 'yes';
  cfg.position = get(gcf, 'Position');
  ft_singleplotTFR(cfg, datvarargin{:});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which handles hot keys in the current plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key_sub(handle, eventdata, varargin)
ident       = get(gca, 'tag');
info        = guidata(gcf);

climits = caxis;
incr_c  = abs(climits(2) - climits(1)) /10;

newz = climits;
if length(eventdata.Modifier) == 1 && strcmp(eventdata.Modifier{:}, 'control')
  % TRANSLATE by 10%
  switch eventdata.Key
    case 'pageup'
      newz = [climits(1)+incr_c climits(2)+incr_c];
    case 'pagedown'
      newz = [climits(1)-incr_c climits(2)-incr_c];
  end % switch
else
  % ZOOM by 10%
  switch eventdata.Key
    case 'pageup'
      newz = [climits(1)-incr_c climits(2)+incr_c];
    case 'pagedown'
      newz = [climits(1)+incr_c climits(2)-incr_c];
    case 'm'
      newz = [varargin{1} varargin{2}];
  end % switch
end % if

% update the color axis
caxis(newz);

if ~isempty(ident) && isfield(info.(ident), 'commenth') && ~isempty(info.(ident).commenth)
  commentstr = get(info.(ident).commenth, 'string');
  sel        = contains(commentstr, info.(ident).cfg.parameter);
  if any(sel)
    commentstr{sel} = sprintf('%0s=[%.3g %.3g]', info.(ident).cfg.parameter, newz(1), newz(2));
    set(info.(ident).commenth, 'string', commentstr);
  end
end


text:
function [P,A] = GFM(X)
% Generic front modeling

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [N,M] = size(X);
    X     = max(X,1e-12);
    P     = ones(1,M);
    A     = ones(1,M);
    lamda = 1;
	E     = sum(repmat(A,N,1).*X.^repmat(P,N,1),2) - 1;
    MSE   = mean(E.^2);
    for epoch = 1 : 1000
        % Calculate the Jacobian matrix
        J = [repmat(A,N,1).*X.^repmat(P,N,1).*log(X),X.^repmat(P,N,1)];
        % Update the value of each weight
        while true
            Delta  = -(J'*J+lamda*eye(size(J,2)))^-1*J'*E;
            newP   = P + Delta(1:M)';
            newA   = A + Delta(M+1:end)';
            newE   = sum(repmat(newA,N,1).*X.^repmat(newP,N,1),2) - 1;
            newMSE = mean(newE.^2);
            if newMSE < MSE && all(newP>1e-3) && all(newA>1e-3)
                P     = newP;
                A     = newA;
                E     = newE;
                MSE   = newMSE;
                lamda = lamda/1.1;
                break;
            elseif lamda > 1e8
                return;
            else
                lamda = lamda*1.1;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ output_args ] = net_turnOffServer( t )
%% This function is used to turn off the server on the robot
% function [ output_args ] = net_turnOffServer( t )
% t: is the TCP/IP connection object
% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='end';
fprintf(t, theCommand);
fclose(t);
%message=fgets(t);
end

repo: KST-Kuka-Sunrise-Toolbox