text:
function t = bartlett_sample ( m, df, sigma )

%*****************************************************************************80
%
%% BARTLETT_SAMPLE samples the Bartlett distribution.
%
%  Discussion:
%
%    If the matrix T is sampled from the Bartlett distribution, then
%    the matrix W = T' * T is a sample from the Wishart distribution.
%
%    This function requires functions from the PDFLIB and RNGLIB libraries.
%
%    The "initialize()" function from RNGLIB must be called before using
%    this function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Odell, Alan Feiveson,
%    A numerical procedure to generate a sample covariance matrix,
%    Journal of the American Statistical Association,
%    Volume 61, Number 313, March 1966, pages 199-203.
%
%    Stanley Sawyer,
%    Wishart Distributions and Inverse-Wishart Sampling,
%    Washington University,
%    30 April 2007, 12 pages.
%
%  Parameters:
%
%    Input, integer M, the order of the matrix.
%
%    Input, integer DF, the number of degrees of freedom.
%    M <= DF.
%
%    Input, real SIGMA(M,M), the covariance matrix, which should be
%    a symmetric positive definite matrix.
%
%    Output, real T(M,M), the sample matrix from the Bartlett distribution.
%
  if ( df < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BARTLETT_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  DF = %d < M = %d.\n', df, m );
    error ( 'BARTLETT_SAMPLE - Fatal error!\n' );
  end
%
%  Get the upper triangular Cholesky factor of SIGMA.
%
  r = chol ( sigma );
%
%  Sample the unit Bartlett distribution.
%
  tu = bartlett_unit_sample ( m, df );
%
%  Construct the matrix.
%
  t = tu * r;

  return
end

author: johannesgerer


text:
function [f, df] = myFunNew(X0, x, y, sn)
if nargin == 0
    clc; close all;
    x = [-20:0.1:60]';
    y = rand(length(x), 1);
    n = length(x);
    sn = 0.1;
    X0 = [0.1 0.1 log([2.0 2.0])]';
end
if isrow(x)
    x = x';
end
if isrow(y)
    y = y';
end
a = X0(1);
b = X0(2);
l = exp(X0(3));
l2 = l*l;
sf = exp(X0(4));
sf2 = sf*sf;
sn2 = sn^2;

n = length(x);
[K0, dK] = CalKFun(x, x, X0(3:4));
K = K0 + sn2*eye(n);
scale = 1000;
tt = eig(scale*K)/scale;
%%%%%%%%% make sure K is positive definite.
if min(tt) < 0.0
    K = K - min(tt)*eye(n);
end
if rank(K) < n
    K = K + 1e-6*eye(n);
end
% save('./K.mat');
L = chol(K, 'lower');
iL = inv(L);
iK = iL'*iL;
logDetK = 2*sum(log(diag(L)));

mx = a*x + b;
f = -0.5*(y-mx)'*iK*(y-mx) - 0.5*logDetK-0.5*n*log(2*pi);
%%%%%%%%%% calculate gradient.

Alpha = iK*(y-mx);
df2m = [x'; ones(1, n)] * Alpha;
A = Alpha * Alpha' - iK;
df2c = zeros(2, 1);
for i = 1 : 1 : 2
     nSum = 0;
     for id = 1 : 1 : n
         nSum = nSum + A(id, :) * dK{i}(:, id);
     end
     df2c(i) = 0.5*nSum;
end
df = [df2m; df2c];
f = -f;
df = -df;
bTest = 1;
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function d = vl_numder(func, x, varargin)
% VL_NUMDER  Numerical derivative
%   D = VL_NUMDER(FUNC, X) computes the numerical derivative of the
%   function FUNC at point X. X is a real array and is passed as first
%   argument of FUNC.
%
%   D = VL_NUMDER(FUNC, X, ARG1, ARG2, ...) passes ARG1, ARG2, ...  as
%   additional arguments to the function FUNC.
%
%   See also: VL_NUMDER2(), VL_HELP().

% TODO: uniform sacaling of axis is not a good idea

dx = 1e-7 ;
N  = numel(x) ;
fx = feval(func, x, varargin{:}) ;
d  = zeros(length(fx(:)),N) ;

for n=1:N
  e   = zeros(size(x)) ; e(n) = 1 ;
  fxn = feval(func, x+dx*e, varargin{:}) ;
  di  = (fxn - fx) / dx ;
  d(:,n) = di(:) ;
end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function prior = priorExpandParam(prior, params)

% PRIOREXPANDPARAM Expand the prior model's parameters from params vector.
% FORMAT
% DESC returns a prior structure filled with the parameters in the
% given vector. This is used as a helper function to enable
% parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG prior : the prior structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% prior structure.
% RETURN prior : prior structure with the given parameters in the
% relevant locations.
%
% As well as extracting the parameters, some transformation of
% parameters is also undertaken in this file. If the field
% transforms is not empty, it dictactes how the prior parameters
% are to be transformed (for example by a exponential to keep them
% positive).
%
% SEEALSO : priorExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004

% PRIOR

if isfield(prior, 'transforms')
  for i = 1:length(prior.transforms)
    index = prior.transforms(i).index;
    fhandle = str2func([prior.transforms(i).type 'Transform']);
    params(index) = fhandle(params(index), 'atox');
  end
end

fhandle = str2func([prior.type 'PriorExpandParam']);
prior = fhandle(prior, params);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef (Abstract) MetricGeneratorX < BaseX
% MetricGeneratorX Abstract class
%
% Summary of MetricGeneratorX:
% This is the base class for all TrackingX metric generators.
% Any custom defined metric generator should be derived from this MetricGeneratorX
% base class.
%
% MetricGeneratorX Properties:
%   None
%
% MetricGeneratorX Methods:
%   + MetricGeneratorX - Constructor method
%   ~ evaluate - Evaluates and returns a given metric
%
% (+) denotes puplic properties/methods
%
% February 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
    end

    methods (Abstract)
        % evaluate Evaluate the correctness of a set of tracks compared to
        % the groundtrouth.
        %   METRIC = EVALUATE(THIS, TRACKLIST, GROUNDTRUTH) evaluates the
        %   correctness of the TRACKLIST compared to the GROUNDTRUTH.
        %   TRACKLIST and GROUNDTRUTH can be structured based on the needs
        %   of the respective metric generator.

        metric = evaluate(~,TrackList,GroundTruth);
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function c=dctiii(f,L,dim)
%DCTIII  Discrete Consine Transform type III
%   Usage:  c=dctiii(f);
%           c=dctiii(f,L);
%           c=dctiii(f,[],dim);
%           c=dctiii(f,L,dim);
%
%   `dctiii(f)` computes the discrete cosine transform of type III of the
%   input signal *f*. If *f* is multi-dimensional, the transformation is
%   applied along the first non-singleton dimension.
%
%   `dctiii(f,L)` zero-pads or truncates *f* to length *L* before doing the
%   transformation.
%
%   `dctiii(f,[],dim)` or `dctiii(f,L,dim)` applies the transformation along
%   dimension *dim*.
%
%   The transform is real (output is real if input is real) and orthonormal.
%
%   This is the inverse of |dctii|.
%
%   Let f be a signal of length *L*, let `c=dctiii(f)` and define the vector
%   *w* of length *L* by
%
%   ..  w = [1/sqrt(2) 1 1 1 1 ...]
%
%   .. math:: w\left(n\right)=\begin{cases}\frac{1}{\sqrt{2}} & \text{if }n=0\\1 & \text{otherwise}\end{cases}
%
%   Then
%
%   ..                     L-1
%     c(n+1) = sqrt(2/L) * sum w(m+1)*f(m+1)*cos(pi*(n+.5)*m/L)
%                          m=0
%
%   .. math:: c\left(n+1\right)=\sqrt{\frac{2}{L}}\sum_{m=0}^{L-1}w\left(m\right)f\left(m+1\right)\cos\left(\frac{\pi}{L}\left(n+\frac{1}{2}\right)m\right)
%
%   Examples:
%   ---------
%
%   The following figures show the first 4 basis functions of the DCTIII of
%   length 20:::
%
%     % The dctii is the adjoint of dctiii.
%     F=dctii(eye(20));
%
%     for ii=1:4
%       subplot(4,1,ii);
%       stem(F(:,ii));
%     end;
%
%   See also:  dctii, dctiv, dstii
%
%   References: rayi90 wi94

%   AUTHOR: Peter L. Søndergaard
%   TESTING: TEST_PUREFREQ
%   REFERENCE: REF_DCTIII

complainif_argnonotinrange(nargin,1,3,mfilename);

if nargin<3
  dim=[];
end;

if nargin<2
  L=[];
end;

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,'DCTIII');

if ~isempty(L)
  f=postpad(f,L);
end;
c=comp_dct(f,3);
% c=zeros(2*L,W,assert_classname(f));
%
% m1=1/sqrt(2)*exp(-(0:L-1)*pi*i/(2*L)).';
% m1(1)=1;
%
% m2=1/sqrt(2)*exp((L-1:-1:1)*pi*i/(2*L)).';
%
% for w=1:W
%   c(:,w)=[m1.*f(:,w);0;m2.*f(L:-1:2,w)];
% end;
%
% c=fft(c)/sqrt(L);
%
% c=c(1:L,:);
%
% if isreal(f)
%   c=real(c);
% end;

c=assert_sigreshape_post(c,dim,permutedsize,order);

% This is a slow, but convenient way of expressing the above algorithm.
%R=1/sqrt(2)*[diag(exp(-(0:L-1)*pi*i/(2*L)));...
%	     zeros(1,L); ...
%	     [zeros(L-1,1),flipud(diag(exp((1:L-1)*pi*i/(2*L))))]];

%R(1,1)=1;

%c=fft(R*f)/sqrt(L);

%c=c(1:L,:);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [returnVal, txtReturnVal] = lvmTwoDPlot(X, lbl, symbol, fhandle)

% LVMTWODPLOT Helper function for plotting the labels in 2-D.
% FORMAT
% DESC helper function for plotting an embedding in 2-D with symbols.
% ARG X : the data to plot.
% ARG lbl : the labels of the data point.
% ARG symbol : the symbols to use for the different labels.
%
% SEEALSO : lvmScatterPlot, lvmVisualise
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2008

% MLTOOLS

if nargin < 2
    lbl = [];
end
if(strcmp(lbl, 'connect'))
    connect = true;
    lbl = [];
else
    connect = false;
end

if iscell(lbl)
    lblStr = lbl{2};
    lbl = lbl{1};
    labelsString = true;
else
    labelsString = false;
end

if nargin < 3 || isempty(symbol)
    if isempty(lbl)
        symbol = getSymbols(1);
    else
        symbol = getSymbols(size(lbl,2));
    end
end
if nargin > 3 && ~isempty(fhandle)
    axisHand = fhandle;
else
    axisHand = gca;
end
returnVal = [];
textReturnVal = [];
nextPlot = get(axisHand, 'nextplot');
for i = 1:size(X, 1)
    if i == 2
        set(axisHand, 'nextplot', 'add');
    end
    if ~isempty(lbl)
        labelNo = find(lbl(i, :));
    else
        labelNo = 1;
    end
    %try
    returnVal = [returnVal; plot(X(i, 1), X(i, 2), symbol{labelNo})];
    if labelsString
        textReturnVal = [textReturnVal; text(X(i, 1), X(i, 2), ['   ' lblStr{i}])];
    end
    if connect
        if i>1
            line([X(i-1, 1) X(i, 1)], [X(i-1, 2) X(i, 2)]);
        end
    end
    %catch
    % if strcmp(lasterr, 'Index exceeds matrix dimensions.')
    %  error(['Only ' num2str(length(symbol)) ' labels supported (it''s easy to add more!)'])
    % end
    %end
end
set(axisHand, 'nextplot', nextPlot);
set(returnVal, 'markersize', 10);
set(returnVal, 'linewidth', 2);
end


text:
function msh = meshStoreSettings(msh, name);
% Store Settings in a mesh, updating any GUIs that reflect this.
%
% msh = meshStoreSettings(msh, <name='Setting#'>);
%
% This code will use the meshSettings command to grab mesh view settings,
% and append the result to the msh.settings field, creating it if it
% doesn't exist. In addition, if a uicontrol is found with the tag
% 'MeshSettingsList', will update the string in this code to reflect the
% new setting.
%
% Also saves it in a meshDir/'MeshSettings.mat' file, since the mesh
% saving/loading has gotten complicated.
%
% ras, 03/14/06.
if notDefined('msh'), msh = viewGet(getSelectedGray, 'selectedMesh');  end

% get the saved settings
meshDir = fileparts(msh.path);
if isempty(meshDir)
    meshDir = viewGet(getSelectedGray, 'meshDir');
end

settingsFile = fullfile(meshDir, 'MeshSettings.mat');
if exist(settingsFile, 'file')
    load(settingsFile, 'settings');
    N = length(settings) + 1; %#ok<NODEF>
else
    N = 1;
end

if notDefined('name')
    name = sprintf('Setting%i',N);
end

S = meshSettings(msh);
S.name = name;
settings(N) = S;

% look for a uicontrol from a mesh setting GUI
h = findobj('Tag', 'MeshSettingsList');

if length(h) > 1
    % more than one menu/list open: loop through lists
    for ii = 1:length(h)
       if isequal( get(h(ii), 'Type'), 'uimenu' ) && ...
           get( get(h(ii), 'Parent'), 'Parent' )==gcf
            h = h(ii);
            break
       end
    end
end

if ~isempty(h)
	% ras 11/07: this can either be a uicontrol listbox (as in the 3D
	% Window and mrViewer) or a uimenu (in the new gray menu)
	if isequal( get(h, 'Type'), 'uicontrol' )
	    set(h, 'String', {settings.name});

	elseif isequal( get(h, 'Type'), 'uimenu' )
		%% make a submenu

		% (prompt for a name now, since renaming later is a little less
		% convenient)
		prompt = {'Enter the new name for the mesh settings'};
	    name = inputdlg(prompt, mfilename, 1, {settings(N).name});
		if isempty(name), return; end   % exit gracefully if cancel

		% the callback should retrieve the selected settings
		cb = sprintf('meshRetrieveSettings([], %i); ', N);

		% if this is the first stored setting, we want a separator between
		% this and the other options (store, rename, delete); otherwise, no
		% separator:
		if N==1, sep = 'on'; else sep = 'off'; end

		uimenu(h, 'Label', name{1}, 'Callback', cb, 'Separator', sep);

		% set this as the name in the settings
		settings(N).name = name{1}; %#ok<NASGU>
	end
end

% save the settings
save(settingsFile, 'settings');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function samples=cosmo_sample_unique(k,n,count,varargin)
% sample without replacement from subsets of integers in balanced manner
%
% function samples=cosmo_sample_unique(k,n[,count,varargin])
%
% Inputs:
%   k           number of elements to return in each subset
%   n           size of integer range from which to sample
%   count       number of subsets to select (default: 1)
%   'seed', s   Use seed s for pseudo-random sampling (optional). If this
%               option is omitted, then different calls to this function
%               may (usually: will) return different results
%
%
% Output:
%   samples     k x count indices, all in the range 1:n, with the following
%               properties:
%               - each value is randomly sampled from the range 1:n
%               - each column forms a subset of 1:n (without repeats)
%               - across the entire matrix, each value in the range 1:n
%                 occurs approximately equally often
%
% Example:
%     % get 4 random subsets of 3 elements in range 1:7
%     % (in this example a seed is used to get the same result upon every
%     % function call)
%     cosmo_sample_unique(3,6,4,'seed',3)
%     %||      1     2     1     2
%     %||      3     5     3     4
%     %||      4     6     5     6
%
% Notes:
%   - this is a utility function; it does not work on dataset structures.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<3 || isempty(count)
        count=1;
    end

    ensure_scalar_nat({k,n,count});
    if k>n
        error(['First argument (%d) cannot be greater than '...
                    'second argument (%d)'],k,n);
    end

    opt=cosmo_structjoin(varargin{:});

    % random elements, each column is a permutation of 1:count
    rs_mat=random_permutations(n,count+1,opt);

    % vectorize
    rs=rs_mat(:);

    % allocate space for output
    samples=zeros(k,count);

    % keep track of which elements in rs were visited
    visited=false((k+1)*count,1);

    first_non_visited_pos=1;
    in_bin=false(n,k); % re-use this for each column
    for col=1:count
        % no elements added so far
        in_bin(:)=false;

        pos=first_non_visited_pos;
        for row=1:k
            % avoid duplicates in in_bin
            while visited(pos) || in_bin(rs(pos))
                pos=pos+1;
            end

            r=rs(pos);
            in_bin(r)=true;
            samples(row,col)=r;
            visited(pos)=true;
        end

        % update first_non_visited pos
        while visited(first_non_visited_pos)
            first_non_visited_pos=first_non_visited_pos+1;
        end
    end

    samples=sort(samples,1);

function r=random_permutations(n,count,opt)
    % output r has in each column the values 1:count in randomly permuted
    % order
    if isfield(opt,'seed') && ~isempty(opt.seed)
        args={'seed',opt.seed};
    else
        args={};
    end

    v=cosmo_rand(n,count,args{:});
    [unused,r]=sort(v,1);

function ensure_scalar_nat(vs)
    for k=1:numel(vs)
        v=vs{k};

        if ~(isnumeric(v) && isscalar(v) && v>0 && round(v)==v)
            error('Argument %d must be positive scalar integer',k);
        end
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function img = discreteCube(varargin)
%DISCRETECUBE discretize a 3D cube
%
%   IMG = discreteCube(LX, LY, LZ, CUBE)
%   Creates a 3D image of a cube, given coordinates of voxel centers in LX,
%   LY and LZ, and parameters of the shape in CUBE.
%
%   CUBE has the format:
%   [XC YC ZC SIDE], with (XC YC ZC) being the coordinates of the center,
%   and SIDE being the side length of the cube.
%
%   Oriented cubes are represented as follow:
%   [XC YC ZC SIDE THETA PHI PSI], with (THETA, PHI) being the angle of the
%   normal of the cube (THETA is colatitude, PHI is longitude), and PSI is
%   rotation of the cube around the normal. All angles are in degrees.
%
%   IMG = discreteCube(DIM, ...)
%   DIM is the size of image, with the format [x0 dx x1;y0 dy y1;z0 dz z1]
%
%   IMG = discreteCube(DIM, CENTER, SIDE, THETA, PHI)
%   Also specify spherical angle of the normal of a face of the cube. THETA
%   is the angle from the vertical, and PHI is the angle from the 0x axis,
%   both between 0 and 360.
%
%   IMG = discreteCube(DIM, CENTER, SIDE, THETA, PHI, PSI)
%   Also consider 'roll' of the cube around its normal axis. PSI is an
%   angle given in degrees between 0 and 360.
%
%   IMG = discreteCube(DIM, CUBE)
%   send parameters in a row vector, where CUBE contains at least the
%   center coordinate, and possibly the other parameters.
%
%   Example
%   img = discreteCube(1:100, 1:100, 1:100, [50 50 50 30 10 60 45]);
%   img = discreteCube([1 1 100;1 1 100;1 1 100], [50 50 50], 30, 10);
%   img = discreteCube([1 1 100;1 1 100;1 1 100], [50 50 50 30 10]);
%
%   See Also
%   imShapes, discreteBall, discreteEllipsoid, discreteCylinder, drawCube
%

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

%   HISTORY
%   03/11/2006: concatenate transforms, add doc for PSI
%   04/03/2009: use meshgrid
%   30/04/2009: update transforms
%   29/05/2009: use more possibilities for specifying grid
%   22/01/2010: fix auto center with odd image size
%   2011-06-21 use degrees for angles

% compute coordinate of image voxels
[lx, ly, lz, varargin] = parseGridArgs3d(varargin{:});
[x, y, z] = meshgrid(lx, ly, lz);

% default parameters
center  = [lx(ceil(end/2)) ly(ceil(end/2)) lz(ceil(end/2))];
side    = center;
theta   = 0; phi = 0; psi=0;

% process input parameters
if length(varargin)==1
    var = varargin{1};
    center = var(:,1:3);
    if size(var, 2)>3
        side = var(:,4);
    end
    if size(var, 2)>4
        theta = var(:,5);
    end
    if size(var, 2)>5
        phi = var(:,6);
    end
    if size(var, 2)>6
        psi = var(:,7);
    end

elseif ~isempty(varargin)
    center = varargin{1};
    if length(varargin)>1
        side = varargin{2};
    end
    if length(varargin)>2
        theta = varargin{3};
    end
    if length(varargin)>3
        phi = varargin{4};
    end
    if length(varargin)>4
        psi = varargin{5};
    end
end

if length(side) == 1
    side = [side side side];
end

% compute coordinate of image voxels in cube reference system
trans = composeTransforms3d(...
    createTranslation3d(-center),...
    createRotationOz(-deg2rad(phi)),...
    createRotationOy(-deg2rad(theta)), ...
    createRotationOz(-deg2rad(psi)));
[x, y, z] = transformPoint3d(x, y, z, trans);

% create image: simple threshold over 3 dimensions
img = abs(x) <= side(1)/2 & abs(y) <= side(2)/2 & abs(z) <= side(3)/2;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ x_scaled ] = ml_scale(x,a,b )
%ML_SCALE Summary of this function goes here
%   Detailed explanation goes here

x_scaled = (x-min(x))*(b-a)/(max(x)-min(x)) + a;
end

author: epfl-lasa


text:
function pass = test_iszero(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% Test scalars:
f = chebfun(0, pref);
pass(1) = iszero(f);

f = chebfun([], pref);
pass(2) = iszero(f);

f = chebfun(2, pref);
pass(3) = ~iszero(f);

% Test piecewise domains:
f = chebfun(0, [-1, 0, 1], pref);
pass(4) = iszero(f);

f = chebfun([], [-1, 0, 1], pref);
pass(5) = iszero(f);

f = chebfun({0, 1}, [-1, 0, 1], pref);
pass(6) = ~iszero(f);

f = chebfun({1, 0}, [-1, 0, 1], pref);
pass(7) = ~iszero(f);

% Test arrays:
f = chebfun([0 0], pref);
pass(8) = all(iszero(f));

f = chebfun([0 1], pref);
pass(9) = all(iszero(f) == [1 0]);

%% Test on singular function:
dom = [-2 7];
pow = -1.64;
f = chebfun(@(x) sin(100*x).*(x-dom(1)).^pow, dom, 'exps', [pow 0], ...
    'splitting', 'on');
pass(10) = ~iszero(f);

%% Test for functions defined on unbounded domain:

% Set the domain:
dom = [-Inf Inf];

% Blow-up function:
op = @(x) x.^2.*(1-exp(-x.^2));
f = chebfun(op, dom, 'exps', [2 2]);
pass(11) = ~iszero(f);

% Function defined on [0 Inf]:

% Specify the domain:
dom = [0 Inf];

op = @(x) 0.75+sin(10*x)./exp(x);
f = chebfun(op, dom, 'splitting', 'on');
pass(12) = ~iszero(f);

% Function defined on [0 Inf]:

% Set the domain:
dom = [-Inf -3*pi];

% Blow-up function:
op = @(x) x.*(5+exp(x.^3))./(dom(2)-x);
f = chebfun(op, dom, 'exps', [0 -1]);
pass(13) = ~iszero(f);

% Zero function:
f = chebfun(@(x) 0*x, dom);
pass(14) = iszero(f);

% [TODO]: Add these once SUBSREF is implemented.
% f = chebfun(0, pref);
% f(0) = 1;
% pass(10) = iszero(f);
%
% f = chebfun([0 0], pref);
% f(0) = [0, 1];
% pass(11) = all(iszero(f) == [1, 0]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef printingVigdergauzMicroStructures < handle

    properties (Access = private)
        mesh
        outputName
        pdim
        postprocess
    end

    methods (Access = public)

        function init(obj)
            obj.ouputName = 'VigergauzMicroStructure';
        end

        function obj = printingVigdergauzMicroStructures()
            settings = Settings('VigergauzLevelSetInput');

            settings.superEllipseRatio = tan(pi/6);%tan(pi/5);
            settings.volumeMicro = 0.8;

            translator = SettingsTranslator();
            translator.translate(settings);
            fileName = translator.fileName;
            settingsTopOpt = SettingsTopOptProblem(fileName);
            topOpt = TopOpt_Problem(settingsTopOpt);

            dV = topOpt.designVariable;

        end

    end

    methods (Access = private)

        function createPostProcess(obj)
            dB = obj.createPostProcessDataBase();
            dB.printers = obj.printers;
            postCase = 'NumericalHomogenizer';
            obj.postProcess = Postprocess(postCase,dB);
        end

        function dB = createPostProcessDataBase(obj)
            dI.mesh            = obj.mesh;
            dI.outName         = obj.outputName;
            dI.pdim            = obj.pdim;
            dI.ptype           = 'MICRO';
            ps = PostProcessDataBaseCreator(dI);
            dB = ps.getValue();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
clear
clc
close all
addpath('./tensor_fit/');
addpath('./RiceOptVST/');
addpath('./GL-HOSVD/');
addpath('./HOSVD/');
addpath('./data/in_vivo_data/');

load '6dir_b1000_NSA10.mat';
load 'bvalue.mat'
load 'gradient.mat'
load 'mask_slice_4th.mat'
[N1,N2,numDWI] = size(im_r);

 % ---- noise estimation ---- %
idata  = im_r.^2; % assumption of the rician distribution
tdata  = idata(1:30,1:30,2:end);
for i  = 1:numDWI-1
    tbg = tdata(:,:,i);
    sigmas(i) = sqrt(mean(tbg(:))/2);
end
sigmah = mean(sigmas(2:end)); % estimate of sigma

 % ---- Setting up parameters ---- %
kglobal=0.4;
klocal=0.5; %the denoising effect can be improved by adjusting the parameter klocal

% ---- start image denoising---- %
rimavst = riceVST(im_r,sigmah,'A');
ims_denoised =glhosvd(rimavst,1,kglobal,klocal);
ims_denoised = riceVST_EUI(ims_denoised ,sigmah,'A');

% ---- FA estimation---- %
display   = 0; bVal=bvalue';   bacq      = bvalue(2);
[FA_denoised, RGB_denoised, tensors_denoised, MD_denoised] = tensor_est(ims_denoised,gradientDirections,bVal,bacq,display,mask);
[FA_noisy, RGB_noisy, tensors_noisy,MD_noisy] = tensor_est(im_r,gradientDirections,bVal,bacq,display,mask);
[FA_reference, RGB_reference, tensors_reference,MD_reference] = tensor_est(im_ref,gradientDirections,bVal,bacq,display,mask);

% ---- Display results---- %
figure,imshow(FA_noisy,[0 1]);figure,imshow(FA_denoised,[0 1]);figure,imshow(FA_reference,[0 1])
figure,imshow(RGB_noisy,[]);figure,imshow(RGB_denoised,[]);figure,imshow(RGB_reference,[])
figure,imshow(im_r(:,:,5),[]);figure,imshow(ims_denoised(:,:,5),[]);figure,imshow(im_ref(:,:,5),[])

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [label, energy] = knKmeansPred(model, Xt)
% Prediction for kernel kmeans clusterng
% Input:
%   model: trained model structure
%   Xt: d x n testing data
% Ouput:
%   label: 1 x n predict label
%   engery: optimization target value
% Written by Mo Chen (sth4nth@gmail.com).
X = model.X;
t = model.label;
kn = model.kn;

n = size(X,2);
k = max(t);
E = sparse(t,1:n,1,k,n,n);
E = E./sum(E,2);
Z = E*kn(X,Xt)-dot(E*kn(X,X),E,2)/2;
[val, label] = max(Z,[],1);
energy = sum(kn(Xt))-2*sum(val);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function cS = magnetite

%
cs_Mag = crystalSymmetry('m3m');
N = Miller({1,0,0},{1,1,1},cs_Mag);
dist = [0.95, 1];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
function treepack_test16 ( )

%*****************************************************************************80
%
%% TREEPACK_TEST16 tests TREE_ROOTED_DEPTH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
  nnode = 12;

  parent = [ 0, 1, 1, 2, 2, 2, 3, 3, 5, 5, 6, 10 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TREEPACK_TEST16\n' );
  fprintf ( 1, '  TREE_ROOTED_DEPTH: depth of a rooted tree.\n' );

  i4vec_print ( nnode, parent, '  Parent vector for tree:' );

  [ depth, depth_node ] = tree_rooted_depth ( nnode, parent );

  i4vec_print ( nnode, depth_node, '  Individual node depths:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Overall rooted tree depth: %d\n', depth );

  return
end

repo: jburkardt-m


text:
function [data] = ft_megrealign(cfg, data)

% FT_MEGREALIGN interpolates MEG data towards standard gradiometer locations by
% projecting the individual timelocked data towards a coarse source reconstructed
% representation and computing the magnetic field on the standard gradiometer
% locations.
%
% Use as
%   [interp] = ft_megrealign(cfg, data)
% where the input data corresponds to the output from FT_PREPROCESSING.
%
% Required configuration options are
%   cfg.template
%   cfg.inwardshift
%
% The new gradiometer definition is obtained from a template dataset,
% or can be constructed by averaging the gradiometer positions over
% multiple datasets.
%   cfg.template       = single dataset that serves as template
%   cfg.template(1..N) = datasets that are averaged into the standard
%
% The realignment is done by computing a minumum norm estimate using a
% large number of dipoles that are placed in the upper layer of the brain
% surface, followed by a forward computation towards the template
% gradiometer array. This requires the specification of a volume conduction
% model of the head and of a source model.
%
% A volume conduction model of the head should be specified with
%   cfg.headmodel   = structure, see FT_PREPARE_HEADMODEL
%
% A source model (i.e. a superficial layer with distributed sources) can be
% constructed from a headshape file, or from inner surface of the volume conduction
% model using FT_PREPARE_SOURCEMODEL using the following options
%   cfg.spheremesh  = number of dipoles in the source layer (default = 642)
%   cfg.inwardshift = depth of the source layer relative to the headshape
%                     surface or volume conduction model (no default
%                     supplied, see below)
%   cfg.headshape   = a filename containing headshape, a structure containing a
%                     single triangulated boundary, or a Nx3 matrix with surface
%                     points
%
% If you specify a headshape and it describes the skin surface, you should specify an
% inward shift of 2.5 cm.
%
% For a single-sphere or a local-spheres volume conduction model based on the skin
% surface, an inward shift of 2.5 cm is reasonable.
%
% For a single-sphere or a local-spheres volume conduction model based on the brain
% surface, you should probably use an inward shift of about 1 cm.
%
% For a realistic single-shell volume conduction model based on the brain surface, you
% should probably use an inward shift of about 1 cm.
%
% Other configuration options are
%   cfg.tolerance  = tolerance ratio for leadfield matrix inverse based on a truncated svd,
%                    reflects the relative magnitude of the largest singular value
%                    to retain (default =s 1e-3)
%   cfg.verify     = 'yes' or 'no', show the percentage difference (default = 'yes')
%   cfg.feedback   = 'yes' or 'no' (default = 'no')
%   cfg.channel    =  Nx1 cell-array with selection of channels (default = 'MEG'),
%                     see FT_CHANNELSELECTION for details
%   cfg.trials     = 'all' or a selection given as a 1xN vector (default = 'all')
%
% This implements the method described by T.R. Knosche, Transformation
% of whole-head MEG recordings between different sensor positions.
% Biomed Tech (Berl). 2002 Mar;47(3):59-62. For more information and
% related methods, see Stolk et al., Online and offline tools for head
% movement compensation in MEG. NeuroImage, 2012.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_PREPARE_LOCALSPHERES, FT_PREPARE_SINGLESHELL

% Copyright (C) 2004-2014, Robert Oostenveld
%
% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% store the original datatype
dtype = ft_datatype(data);

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes', 'hassampleinfo', 'yes', 'ismeg', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'plot3d',      'feedback'});
cfg = ft_checkconfig(cfg, 'renamedval', {'headshape',   'headmodel', []});
cfg = ft_checkconfig(cfg, 'required',   {'inwardshift', 'template'});
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',         'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'grid',        'sourcemodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'pruneratio',  'tolerance'});

% set the default configuration
cfg.headshape  = ft_getopt(cfg, 'headshape',  []);
cfg.pruneratio = ft_getopt(cfg, 'tolerance',  1e-3);
cfg.spheremesh = ft_getopt(cfg, 'spheremesh', 642);
cfg.verify     = ft_getopt(cfg, 'verify',     'yes');
cfg.feedback   = ft_getopt(cfg, 'feedback',   'yes');
cfg.trials     = ft_getopt(cfg, 'trials',     'all', 1);
cfg.channel    = ft_getopt(cfg, 'channel',    'MEG');
cfg.topoparam  = ft_getopt(cfg, 'topoparam',  'rms');

% do realignment per trial
pertrial = all(ismember({'nasX';'nasY';'nasZ';'lpaX';'lpaY';'lpaZ';'rpaX';'rpaY';'rpaZ'}, data.label));

% put the low-level options pertaining to the dipole grid in their own field
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'}); % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

if isstruct(cfg.template)
  % this should be a cell-array
  cfg.template = {cfg.template};
end

% retain only the MEG channels and temporarily store the rest of the channels
% elsewhere, these will be added back to the transformed data later.

% select trials and channels of interest, first of the non-MEG channels, then of the MEG channels
tmpcfg = keepfields(cfg, {'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'}); % don't keep tolerance, it is used differently here
tmpcfg.channel = setdiff(data.label, ft_channelselection(cfg.channel, data.label), 'stable');
rest = ft_selectdata(tmpcfg, data);
tmpcfg.channel = ft_channelselection(cfg.channel, data.label);
data = ft_selectdata(tmpcfg, data);

% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

Ntrials = length(data.trial);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% construct the average template gradiometer array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template = struct([]); % initialize as 0x0 empty struct array with no fields
for i=1:length(cfg.template)
  if ischar(cfg.template{i})
    ft_info('reading template sensor position from %s\n', cfg.template{i});
    tmp = ft_read_sens(cfg.template{i}, 'senstype', 'meg');
  elseif isstruct(cfg.template{i}) && isfield(cfg.template{i}, 'coilpos') && isfield(cfg.template{i}, 'coilori') && isfield(cfg.template{i}, 'tra')
    tmp = cfg.template{i};
  elseif isstruct(cfg.template{i}) && isfield(cfg.template{i}, 'pnt') && isfield(cfg.template{i}, 'ori') && isfield(cfg.template{i}, 'tra')
    % it seems to be a pre-2011v1 type gradiometer structure, update it
    tmp = ft_datatype_sens(cfg.template{i});
  else
    ft_error('unrecognized template input');
  end
  % prevent "Subscripted assignment between dissimilar structures" error
  template = appendstruct(template, tmp); clear tmp
end

grad = ft_average_sens(template);

% construct the final template gradiometer definition
template = [];
template.grad = grad;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PREPARE_HEADMODEL will match the data labels, the gradiometer labels and the
% volume model labels (in case of a localspheres model) and result in a gradiometer
% definition that only contains the gradiometers that are present in the data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

volcfg = [];
volcfg.headmodel = cfg.headmodel;
volcfg.grad      = data.grad;
volcfg.channel   = data.label; % this might be a subset of the MEG channels

% FIXME As of yet the next steps might not entirely correct, because it does not keep
% track of the balancing of the gradiometer array. This may require some thought
% because the leadfields are computed with low level functions and do not easily
% accommodate for matching the correct channels with each other (in order to compute
% the projection matrix).

% PREPARE_HEADMODEL will match the data labels, the gradiometer labels and the
% volume model labels (in case of a localspheres model) and result in a gradiometer
% definition that only contains the gradiometers that are present in the data.
[volold, data.grad] = prepare_headmodel(volcfg, []);

% Note that it is necessary to keep the two volume conduction models separate, since
% the single-shell Nolte model contains gradiometer specific precomputed parameters.
% Also note that this is not guaranteed to result in a good projection for local
% sphere models.
volcfg.grad    = template.grad;
volcfg.channel = 'MEG'; % include all MEG channels
[volnew, template.grad] = prepare_headmodel(volcfg, []);

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

if strcmp(ft_senstype(data.grad), ft_senstype(template.grad))
  [id, it] = match_str(data.grad.label, template.grad.label);
  ft_info('mean distance towards template gradiometers is %.2f %s\n', mean(sum((data.grad.chanpos(id,:)-template.grad.chanpos(it,:)).^2, 2).^0.5), template.grad.unit);
else
  % the projection is from one MEG system to another MEG system, which makes a comparison of the data difficult
  cfg.feedback = 'no';
  cfg.verify = 'no';
end

% copy all options that are potentially used in ft_prepare_sourcemodel
tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
tmpcfg.headmodel = volold;
tmpcfg.grad      = data.grad;
% create the source positions on which the data will be projected
sourcemodel = ft_prepare_sourcemodel(tmpcfg);

% compute the forward model for the new gradiometer positions
ft_info('computing forward model for %d dipoles\n', size(sourcemodel.pos,1));
lfnew = ft_compute_leadfield(sourcemodel.pos, template.grad, volnew, leadfieldopt{:});
if ~pertrial
  % this needs to be done only once
  lfold = ft_compute_leadfield(sourcemodel.pos, data.grad, volold, leadfieldopt{:});
  [realign, noalign, bkalign] = computeprojection(lfold, lfnew, cfg.pruneratio, cfg.verify);
end

% interpolate the data towards the template gradiometers
for i=1:Ntrials
  ft_info('realigning trial %d\n', i);
  if pertrial
    %warp the gradiometer array according to the motiontracking data
    sel   = match_str(rest.label, {'nasX';'nasY';'nasZ';'lpaX';'lpaY';'lpaZ';'rpaX';'rpaY';'rpaZ'});
    hmdat = rest.trial{i}(sel,:);
    if ~all(hmdat==repmat(hmdat(:,1),[1 size(hmdat,2)]))
      ft_error('only one position per trial is at present allowed');
    else
      M    = ft_headcoordinates(hmdat(1:3,1),hmdat(4:6,1),hmdat(7:9,1));
      grad = ft_transform_geometry(M, data.grad);
    end

    volcfg.grad = grad;
    % compute volume conductor
    [volold, grad] = prepare_headmodel(volcfg, []);
    % compute forward model
    lfold = ft_compute_leadfield(sourcemodel.pos, grad, volold, leadfieldopt{:});
    % compute projection matrix
    [realign, noalign, bkalign] = computeprojection(lfold, lfnew, cfg.pruneratio, cfg.verify);
  end
  data.realign{i} = realign * data.trial{i};
  if strcmp(cfg.verify, 'yes')
    % also compute the residual variance when interpolating
    [id,it]   = match_str(data.grad.label, template.grad.label);
    rvrealign = rv(data.trial{i}(id,:), data.realign{i}(it,:));
    ft_info('original -> template             RV %.2f %%\n', 100 * mean(rvrealign));
    datnoalign = noalign * data.trial{i};
    datbkalign = bkalign * data.trial{i};
    rvnoalign = rv(data.trial{i}, datnoalign);
    rvbkalign = rv(data.trial{i}, datbkalign);
    ft_info('original             -> original RV %.2f %%\n', 100 * mean(rvnoalign));
    ft_info('original -> template -> original RV %.2f %%\n', 100 * mean(rvbkalign));
  end
end

% plot the topography before and after the realignment
if strcmp(cfg.feedback, 'yes')

  ft_warning('showing MEG topography (RMS value over time) in the first trial only');
  Nchan = length(data.grad.label);
  [id,it]   = match_str(data.grad.label, template.grad.label);
  pos1 = data.grad.chanpos(id,:);
  pos2 = template.grad.chanpos(it,:);
  prj1 = elproj(pos1); tri1 = delaunay(prj1(:,1), prj1(:,2));
  prj2 = elproj(pos2); tri2 = delaunay(prj2(:,1), prj2(:,2));

  switch cfg.topoparam
    case 'rms'
      p1 = sqrt(mean(data.trial{1}(id,:).^2, 2));
      p2 = sqrt(mean(data.realign{1}(it,:).^2, 2));
    case 'svd'
      [u, s, v] = svd(data.trial{1}(id,:)); p1 = u(:,1);
      [u, s, v] = svd(data.realign{1}(it,:)); p2 = u(:,1);
    otherwise
      ft_error('unsupported cfg.topoparam');
  end

  X = [pos1(:,1) pos2(:,1)]';
  Y = [pos1(:,2) pos2(:,2)]';
  Z = [pos1(:,3) pos2(:,3)]';

  % show figure with old an new helmets, volume model and source positions
  figure
  hold on
  ft_plot_headmodel(volold);
  plot3(sourcemodel.pos(:,1),sourcemodel.pos(:,2),sourcemodel.pos(:,3),'b.');
  plot3(pos1(:,1), pos1(:,2), pos1(:,3), 'r.') % original positions
  plot3(pos2(:,1), pos2(:,2), pos2(:,3), 'g.') % template positions
  line(X,Y,Z, 'color', 'black');
  view(-90, 90);

  % show figure with data on old helmet location
  figure
  hold on
  plot3(pos1(:,1), pos1(:,2), pos1(:,3), 'r.') % original positions
  plot3(pos2(:,1), pos2(:,2), pos2(:,3), 'g.') % template positions
  line(X,Y,Z, 'color', 'black');
  axis equal; axis vis3d
  bnd1 = [];
  bnd1.pos = pos1;
  bnd1.tri = tri1;
  ft_plot_mesh(bnd1,'vertexcolor',p1,'edgecolor','none')
  title('RMS, before realignment')
  view(-90, 90)

  % show figure with data on new helmet location
  figure
  hold on
  plot3(pos1(:,1), pos1(:,2), pos1(:,3), 'r.') % original positions
  plot3(pos2(:,1), pos2(:,2), pos2(:,3), 'g.') % template positions
  line(X,Y,Z, 'color', 'black');
  axis equal; axis vis3d
  bnd2 = [];
  bnd2.pos = pos2;
  bnd2.tri = tri2;
  ft_plot_mesh(bnd2,'vertexcolor',p2,'edgecolor','none')
  title('RMS, after realignment')
  view(-90, 90)
end

% store the realigned data in a new structure
interp.label   = template.grad.label;
interp.grad    = template.grad;   % replace with the template gradiometer array
interp.trial   = data.realign;    % remember the processed data
interp.fsample = data.fsample;
interp.time    = data.time;

% add the rest channels back to the data, these were not interpolated
if ~isempty(rest.label)
  ft_info('adding %d non-MEG channels back to the data (', length(rest.label));
  ft_info('%s, ', rest.label{1:end-1});
  ft_info('%s)\n', rest.label{end});
  for trial=1:length(rest.trial)
    interp.trial{trial} = [interp.trial{trial}; rest.trial{trial}];
  end
  interp.label = [interp.label; rest.label];
end

% copy the trial specific information into the output
if isfield(data, 'trialinfo')
  interp.trialinfo = data.trialinfo;
end

% copy the sampleinfo field as well
if isfield(data, 'sampleinfo')
  interp.sampleinfo = data.sampleinfo;
end

% convert back to input type if necessary
switch dtype
  case 'timelock'
    interp = ft_checkdata(interp, 'datatype', 'timelock');
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data

% rename the output variable to accomodate the savevar postamble
data = interp;

ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction that computes the projection matrix(ces)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [realign, noalign, bkalign] = computeprojection(lfold, lfnew, tolerance, verify)

% compute this inverse only once, although it is used twice
tmp = ft_inv(lfold, 'method', 'tsvd', 'tolerance', tolerance);
% compute the three interpolation matrices
ft_info('computing interpolation matrix #1\n');
realign = lfnew * tmp;
if strcmp(verify, 'yes')
  ft_info('computing interpolation matrix #2\n');
  noalign = lfold * tmp;
  ft_info('computing interpolation matrix #3\n');
  bkalign = (lfold * ft_inv(lfnew, 'method', 'tsvd', 'tolerance', tolerance)) * realign;
else
  noalign = [];
  bkalign = [];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
    %% GRADIENTDESCENT Performs gradient descent to learn theta
    %   theta = GRADIENTDESENT(X, y, theta, alpha, num_iters) updates theta by
    %   taking num_iters gradient steps with learning rate alpha

    % Initialize some useful values
    n = length(y); % number of training examples
    J_history = zeros(num_iters, 1);

    for iter = 1 : num_iters

        % ====================== YOUR CODE HERE ======================
        % Instructions: Perform a single gradient step on the parameter vector
        %               theta.
        %
        % Hint: While debugging, it can be useful to print out the values
        %       of the cost function (computeCost) and gradient here.
        %

        % compute the gradient
        grad = zeros(size(theta));
        error = X * theta - y;
        for j = 1 : length(grad)
            grad(j) = alpha * sum(error .* X(:, j)) / n;
        end

        % simultaneously update all theta entries
        theta = theta - grad;

        % Save the cost J in every iteration
        J_history(iter) = computeCost(X, y, theta);
    end
end

author: worldveil

repo: coursera-ml


text:
function clear(logger)
% Removes all logfiles from the list.  Files that were added using
% their name (i.e. not their fid) are closed.

logger.fidlist = [];
logger.fidlevels = [];
for ii=1:length(logger.filenamefids)
    fclose(logger.filenamefids(ii));
end
logger.filenamelist = {};
logger.filenamefids = [];
logger.filenamelevels = [];

assert(logger.validate())

author: nesl


text:
function val = get(P, prop_name)
%
%  GET Get psf properties from specified object and return
%      the value.
%
%     val = get(P, prop_name);
%
%  Valid choices for prop_name are:
%    'image' and 'center'
%

%  J. Nagy 5/1/01

switch prop_name
case 'image'
  val = P.image;
case 'center'
  val = P.center;
otherwise
  error([prop_name, 'Is not valid psf property'])
end

repo: IRtools


text:
% findcircle - returns the coordinates of a circle in an image using the Hough transform
% and Canny edge detection to create the edge map.
%
% Usage:
% [row, col, r] = findcircle(image,lradius,uradius,scaling, sigma, hithres, lowthres, vert, horz)
%
% Arguments:
%	image		    - the image in which to find circles
%	lradius		    - lower radius to search for
%	uradius		    - upper radius to search for
%	scaling		    - scaling factor for speeding up the
%			          Hough transform
%	sigma		    - amount of Gaussian smoothing to
%			          apply for creating edge map.
%	hithres		    - threshold for creating edge map
%	lowthres	    - threshold for connected edges
%	vert		    - vertical edge contribution (0-1)
%	horz		    - horizontal edge contribution (0-1)
%
% Output:
%	circleiris	    - centre coordinates and radius
%			          of the detected iris boundary
%	circlepupil	    - centre coordinates and radius
%			          of the detected pupil boundary
%	imagewithnoise	- original eye image, but with
%			          location of noise marked with
%			          NaN values
%
% Author:
% Libor Masek
% masekl01@csse.uwa.edu.au
% School of Computer Science & Software Engineering
% The University of Western Australia
% November 2003

function [row, col, r] = findcircle(image,lradius,uradius,scaling, sigma, hithres, lowthres, vert, horz)

lradsc = round(lradius*scaling);
uradsc = round(uradius*scaling);
rd = round(uradius*scaling - lradius*scaling);

% generate the edge image
[I2 or] = canny(image, sigma, scaling, vert, horz);
I3 = adjgamma(I2, 1.9);
I4 = nonmaxsup(I3, or, 1.5);
edgeimage = hysthresh(I4, hithres, lowthres);

% perform the circular Hough transform
h = houghcircle(edgeimage, lradsc, uradsc);

maxtotal = 0;

% find the maximum in the Hough space, and hence
% the parameters of the circle
for i=1:rd

    layer = h(:,:,i);
    [maxlayer] = max(max(layer));

    if maxlayer > maxtotal

        maxtotal = maxlayer;

        r = int32((lradsc+i) / scaling);

        [row,col] = ( find(layer == maxlayer) );

        row = int32(row(1) / scaling); % returns only first max value
        col = int32(col(1) / scaling);

    end

end

author: Qingbao

repo: iris

save_path: github-repos/MATLAB/Qingbao-iris


text:
function opt = TSEMO_options
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

% Description: Create TSEMO options structure.
opt.maxeval = 35;                   % Maximum number of function evaluations

for i = 1:3                         % this is only valid for up to 3 objectives. Number needs to be increased for more than 3.
opt.GP(i).nSpectralpoints = 4000;   % Number of spectral sampling points for objective i
opt.GP(i).matern = 1;               % Matern type 1 / 3 / 5 / inf for objective i
opt.GP(i).fun_eval = 200;           % Function evaluations by direct algorithm per input dimension for objective i
end

opt.pop = 100;                      % Genetic algorithm population size
opt.Generation = 100;               % Genetic algorithm number of generations
opt.NoOfBachSequential = 1;         % Batch size
end

author: Eric-Bradford

repo: TS-EMO


text:
%% Setting of the problem
global s
pde = checkboarddata; % f = (x1^2-1)*(x2^2-1); Kellogg problem
pde.L = 1;
option.theta = 0.3;
option.estType = 'star';
option.maxIt = 25;
option.maxN = 2e4;
option.solver = 'mg';
option.tol = 1e-6;
[node,elem] = squaremesh([-1,1,-1,1],0.5);
bdFlag = setboundary(node,elem,'Dirichlet');

%% s = 0.2
s = 0.2; %#ok<*NASGU>
err1 = afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.4
s = 0.4;
err2 = afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
err3 = afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
err4 = afemfracLap(node,elem,pde,bdFlag,option);

%% save data and plot the table
save Kelloggcf err1 err2 err3 err4
load Kelloggcf
figure;
plot_error_table(err1.N,err1.energyError,err2.N,err2.energyError,...
                 err3.N,err3.energyError,err4.N,err4.energyError);
%Save figure
saveas(gcf, 'error_Kellogg_cf', 'pdf')
figure;
plot_error_table(err1.N,err1.eta,err2.N,err2.eta,err3.N,err3.eta,err4.N,err4.eta);
ylabel('Estimator','interpreter','latex', 'FontSize', 22)
saveas(gcf, 'eta_Kellogg_cf', 'pdf')

repo: ifem


text:
classdef MimContextHierarchy < CoreBaseClass
    % MimContextHierarchy. Part of the internal framework of the TD MIM Toolkit.
    %
    %     This class is used to switch beteen different contexts.
    %
    %     There are a number of contexts, which each represent particular
    %     regions of the lung. For example, the OriginalImage context is the
    %     entire image, whereas the LungROI is the parrallipiped region
    %     containing the lungs and airways. The LeftLung and RightLung comprise
    %     just the voumes of the left and right lung respectively.
    %
    %     Context sets describe a collection of related contexts. For example,
    %     the SingleLung context set contains LeftLung and RightLung.
    %
    %     Each plugin specifies its context set, which is the domain of the
    %     results produced by the plugin. Some plugins operate over the whole
    %     LungROI region, while some operate on individual lungs.
    %
    %     This class manages the situations where a result is requested for a
    %     particualar context (e.g. LungROI) but the plugin defines a different
    %     context set (e.g. SingleLung). The LungROI can be built from the two
    %     contexts in the SingleLung set. Therefore in this case, the plugin is run twice,
    %     once for the left and once for the right lung. Then the resutls are
    %     combined to produce the result for the LungROI context.
    %
    %     Thie class manages this heirarchy of contexts and context sets. In
    %     this way plugins can operate on whichever context is appropriate,
    %     while results can be requested for any context, and the conversions
    %     are handled automatically.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ContextSets
        Contexts
        ContextDef

        DatasetDiskCache
        ImageTemplates
        Pipelines
        FrameworkAppDef % Framework configuration
        PluginCache % Used to check plugin version numbers
    end

    methods
        function obj = MimContextHierarchy(context_def, dataset_disk_cache, image_templates, pipelines, framework_app_def, plugin_cache)
            obj.DatasetDiskCache = dataset_disk_cache;
            obj.ImageTemplates = image_templates;
            obj.Pipelines = pipelines;
            obj.FrameworkAppDef = framework_app_def;
            obj.PluginCache = plugin_cache;

            obj.ContextDef = context_def;
            obj.Contexts = context_def.GetContexts;
            obj.ContextSets = context_def.GetContextSets;
        end

        function context_list = GetContextList(obj, output_context, plugin_info, reporting)
            % Parses a list of Contexts, expands out any ContextSets into
            % Contexts, and returns an error if a Context is not known

            % Determines the context and context type requested by the calling function
            if isempty(output_context)
                output_context = obj.ContextDef.ChooseOutputContext(plugin_info.Context);
            end

            context_list = [];

            for next_output_context_set = CoreContainerUtilities.ConvertToSet(output_context)
                next_output_context = next_output_context_set{1};
                if obj.Contexts.isKey(char(next_output_context))
                    % Note that if the enum name of the context set matches
                    % the enum name of a context, we need to be sure the
                    % contet gets added, not the context set
                    context_list{end + 1} = obj.Contexts(char(next_output_context)).Context;
                elseif obj.ContextSets.isKey(char(next_output_context))
                    context_set_mapping = obj.ContextSets(char(next_output_context));
                    context_mapping_list = context_set_mapping.ContextList;
                    context_list = [context_list, CoreContainerUtilities.GetFieldValuesFromSet(context_mapping_list, 'Context')];
                else
                    if obj.DatasetDiskCache.ManualSegmentationExists(char(next_output_context), reporting)
                        context_list{end + 1} = char(next_output_context);
                    else
                        % Allow contexts with specific label indices
                        [context_prefix, context_suffix] = CoreTextUtilities.SplitAtLastDelimiter(char(next_output_context), '.');
                        if obj.DatasetDiskCache.ManualSegmentationExists(context_prefix, reporting)
                            context_list{end + 1} = char(next_output_context);
                        else
                            reporting.Error('MimContextHierarchy:UnknownOutputContext', 'I do not understand the requested output context.');
                        end
                    end
                end
            end
        end

        function value = GetParameter(obj, name, dataset_stack, reporting)
            value = dataset_stack.GetParameterAndAddDependenciesToPluginsInStack(name, reporting);
        end

        function combined_result = GetResultRecursive(obj, plugin_name, output_context, parameters, linked_dataset_chooser, plugin_info, plugin_class, dataset_uid, dataset_stack, force_generate_image, memory_cache_policy, disk_cache_policy, reporting)

            % Records that there has been an attempt to call this plugin,
            % so that we can detect if a plugin call failed
            obj.ImageTemplates.NoteAttemptToRunPlugin(plugin_name, output_context, reporting);

            [plugin_context_set, plugin_context_set_mapping] = obj.GetContextSetMappings(plugin_info.Context);
            [output_context, output_context_set, output_context_mapping, output_context_set_mapping] = obj.GetContextMaps(output_context);

            % If the input and output contexts are of the same type, or if the
            % plugin context is of type 'Any', then proceed to call the plugin
            if isempty(output_context_set) || obj.ContextDef.ContextSetMatches(plugin_context_set, output_context_set)

                % Code moved from PluginDependencyTracker

                % Fetch plugin result from the disk cache
                result = [];

                edited_result_exists = obj.DatasetDiskCache.EditedResultExists(plugin_name, output_context, reporting);

                % We can skip fetching the result if an edited result exists
                % and does not depend on the automated result
                if ~edited_result_exists || plugin_info.EditRequiresPluginResult

                    if ~plugin_info.AlwaysRunPlugin

                        [result, cache_info] = obj.DatasetDiskCache.LoadPluginResult(plugin_name, output_context, memory_cache_policy, reporting);

                        % Check dependencies of the result. If they are invalid, set the
                        % result to null to force a re-run of the plugin
                        if ~isempty(cache_info)
                            dependencies = cache_info.DependencyList;
                            if ~obj.CheckPluginVersion(cache_info.InstanceIdentifier, reporting)
                                reporting.ShowWarning('MimContextHierarchy:InvalidDependency', ['The plugin ' plugin_name ' has changed since the cache was generated. I am forcing this plugin to re-run to generate new results.'], []);
                                result = [];
                            end

                            if ~obj.CheckDependenciesValid(linked_dataset_chooser, dependencies, dataset_stack, parameters, reporting)
                                reporting.ShowWarning('MimContextHierarchy:InvalidDependency', ['The cached value for plugin ' plugin_name ' is no longer valid since some of its dependencies have changed. I am forcing this plugin to re-run to generate new results.'], []);
                                result = [];
                            end
                        end

                        % Add the dependencies of the cached result to any other
                        % plugins in the callstack
                        if ~isempty(result) && ~isempty(cache_info)
                            dependencies = cache_info.DependencyList;
                            dataset_stack.AddDependenciesToAllPluginsInStack(dependencies, reporting);

                            dependency = cache_info.InstanceIdentifier;
                            dependency_list_for_this_plugin = PTKDependencyList();
                            dependency_list_for_this_plugin.AddDependency(dependency, reporting);
                            dataset_stack.AddDependenciesToAllPluginsInStack(dependency_list_for_this_plugin, reporting);
                        end

                    end

                    % Run the plugin
                    if isempty(result)
                        plugin_has_been_run = true;

                        % At present we ignore dependency checks if the results are
                        % not permanently cached to disk. This is to ensure we
                        % don't get dependency check errors when a plugin is
                        % fetched twice (by two different plugins) and has to be
                        % generated twice becase the result is not cached.
                        ignore_dependency_checks = plugin_info.DiskCachePolicy ~= MimCachePolicy.Permanent;

                        % Pause the self-timing of the current plugin
                        dataset_stack.PauseTiming();

                        % Create a new dependency. The dependency relates to the plugin
                        % being called (plugin_name) and the UID of the dataset the
                        % result is being requested from; however, the stack belongs to
                        % the primary dataset
                        plugin_version = plugin_info.PluginVersion;
                        dataset_stack.CreateAndPush(plugin_name, output_context, parameters, dataset_uid, ignore_dependency_checks, false, obj.FrameworkAppDef.TimeFunctions, plugin_version, reporting);

                        dataset_callback = MimDatasetCallback(linked_dataset_chooser, dataset_stack, output_context, reporting);

                        try
                            % This is the actual call which runs the plugin
                            if strcmp(plugin_info.PluginInterfaceVersion, '1')
                                result = plugin_class.RunPlugin(dataset_callback, reporting);
                            else
                                result = plugin_class.RunPlugin(dataset_callback, output_context, reporting);
                            end

                            new_cache_info = dataset_stack.Pop();

                            if obj.FrameworkAppDef.TimeFunctions()
                                dataset_stack.ResumeTiming();
                            end

                            if ~strcmp(plugin_name, new_cache_info.InstanceIdentifier.PluginName)
                                reporting.Error('MimContextHierarchy:GetResult', 'Inconsistency in plugin call stack. To resolve this error, try deleting the cache for this dataset.');
                            end

                            % Get the newly calculated list of dependencies for this
                            % plugin
                            dependencies = new_cache_info.DependencyList;

                            % Cache the plugin result according to the specified cache policies
                            obj.DatasetDiskCache.SavePluginResult(plugin_name, result, new_cache_info, output_context, disk_cache_policy, memory_cache_policy, reporting);

                            dataset_stack.AddDependenciesToAllPluginsInStack(dependencies, reporting);

                            dependency = new_cache_info.InstanceIdentifier;
                            dependency_list_for_this_plugin = PTKDependencyList();
                            dependency_list_for_this_plugin.AddDependency(dependency, reporting);
                            dataset_stack.AddDependenciesToAllPluginsInStack(dependency_list_for_this_plugin, reporting);

                            cache_info = new_cache_info;
                        catch ex
                                % For certain plugins we throw a special exception
                                % which indicates that the user should be offered
                                % the ability to create a manual edit
                                if plugin_info.SuggestManualEditOnFailure
                                    throw(MimSuggestEditException(plugin_name, output_context, ex, CoreTextUtilities.RemoveHtml(plugin_class.ButtonText)));
                                else
                                    rethrow(ex);
                                end
                        end

                    else
                        plugin_has_been_run = false;
                    end

                else
                    cache_info = [];
                    plugin_has_been_run = false;
                end

                % Fetch the edited result, if it exists
                if edited_result_exists

                    [edited_result, edited_cache_info] = obj.DatasetDiskCache.LoadEditedPluginResult(plugin_name, output_context, reporting);

                    % If the edited result does not depend on the automated
                    % result, we won't have (and don't need) cache info for the
                    % plugin result so just use the edited result
                    if isempty(cache_info)
                        cache_info = edited_cache_info;
                    end

                    % In case the cache is out of sync with the existance of
                    % the edited result, this will update the cache
                    obj.DatasetDiskCache.UpdateEditedResults(plugin_name, edited_cache_info, output_context, reporting);

                    % Call the plugin to create an edited output
                    result = plugin_class.GetEditedResult(result, edited_result, reporting);

                    % Get the dependency for this edit and add to calling functions
                    edited_dependency = edited_cache_info.InstanceIdentifier;
                    dependency_list_for_edit = PTKDependencyList();
                    dependency_list_for_edit.AddDependency(edited_dependency, reporting);
                    dataset_stack.AddDependenciesToAllPluginsInStack(dependency_list_for_edit, reporting);
                    cache_info.MarkEdited();
                else
                    % In case the cache is out of sync with the existance of
                    % the edited result, this will delete the edited result entry from the cache
                    obj.DatasetDiskCache.UpdateEditedResults(plugin_name, [], output_context, reporting);

                end

                % End of code moved from PluginDependencyTracker

                % If the plugin has been re-run, then we will generate an output
                % image, in order to create a new preview image
                if (plugin_info.GeneratePreview && plugin_has_been_run)
                    force_generate_image = true;
                end

                % We generate an output image if requested, or if the plugin has been re-run (indictaing that we will need to generate a new preview image)
                if force_generate_image
                    output_image = obj.GenerateImageFromResults(result, plugin_class, linked_dataset_chooser, dataset_stack, reporting);
                else
                    output_image = [];
                end

                combined_result = MimCombinedPluginResult(result, output_image, plugin_has_been_run, cache_info);

            % Otherwise, if the plugin's context set is higher in the hierarchy
            % than the requested output, then get the result for the higher
            % context and then extract out the desired context
            elseif output_context_set_mapping.IsOtherContextSetHigher(plugin_context_set_mapping)
                combined_result = MimContextResultCombiner(obj.ImageTemplates);
                higher_context_mapping = output_context_mapping.Parent;
                if isempty(higher_context_mapping)
                    reporting.Error('MimContextHierarchy:UnexpectedSituation', 'The requested plugin call cannot be made as I am unable to determine the relationship between the plugin context and the requested result context.');
                end
                parent_context = higher_context_mapping.Context;
                this_context_results = obj.GetResultRecursive(plugin_name, parent_context, parameters, linked_dataset_chooser, plugin_info, plugin_class, dataset_uid, dataset_stack, force_generate_image, memory_cache_policy, disk_cache_policy, reporting);
                combined_result.AddParentResult(parent_context, this_context_results, output_context, output_context_mapping, dataset_stack, reporting);

            % If the plugin's context set is lower in the hierarchy, then get
            % the plugin result for all lower contexts and concatenate the results
            elseif plugin_context_set_mapping.IsOtherContextSetHigher(output_context_set_mapping)
                combined_result = MimContextResultCombiner(obj.ImageTemplates);
                child_context_mappings = output_context_mapping.Children;

                for child_mapping = child_context_mappings
                    child_context = child_mapping{1}.Context;
                    this_context_results = obj.GetResultRecursive(plugin_name, child_context, parameters, linked_dataset_chooser, plugin_info, plugin_class, dataset_uid, dataset_stack, force_generate_image, memory_cache_policy, disk_cache_policy, reporting);

                    combined_result.AddChildResult(child_context, this_context_results, output_context, output_context_mapping, dataset_stack, reporting);
                end

            else
                reporting.Error('MimContextHierarchy:UnexpectedSituation', 'The requested plugin call cannot be made as I am unable to determine the relationship between the plugin context and the requested result context.');
            end

            % Note that the plugin completed successfully
            obj.ImageTemplates.NoteSuccessRunPlugin(plugin_name, output_context, reporting);

            % Allow pipelines to be run if required
            obj.Pipelines.RunPipelines(plugin_name, output_context, parameters, combined_result.GetPluginHasBeenRun, dataset_stack, dataset_uid, reporting);
        end

        function SaveEditedResultRecursive(obj, plugin_name, input_context, edited_result_image, plugin_info, dataset_stack, dataset_uid, reporting)
            % This function saves a manually edited result which plugins
            % may use to modify their results

            % Determines the context and context type requested by the calling function
            if isempty(input_context)
                reporting.Error('MimContextHierarchy:NoContextSpecified', 'When calling SaveEditedResult(), the contex of the input image must be specified.');
                input_context = obj.ContextDef.GetDefaultContext;
            end

            [plugin_context_set, plugin_context_set_mapping] = obj.GetContextSetMappings(plugin_info.Context);

            [input_context, input_context_set, input_context_mapping, input_context_set_mapping] = obj.GetContextMaps(input_context);

            % If the input and output contexts are of the same type, or if the
            % plugin context is of type 'Any', then proceed to save the results
            % for this context.
            if obj.ContextDef.ContextSetMatches(plugin_context_set, input_context_set)
                plugin_version = plugin_info.PluginVersion;
                obj.SaveEditedResult(plugin_name, input_context, edited_result_image, dataset_uid, plugin_version, reporting);

            % Otherwise, if the input's context set is lower in the hierarchy
            % than that of the plugin, then resize the input to match the plugin
            elseif input_context_set_mapping.IsOtherContextSetHigher(plugin_context_set_mapping)

                higher_context_mapping = input_context_mapping.Parent;
                if isempty(higher_context_mapping)
                    reporting.Error('MimContextHierarchy:UnexpectedSituation', 'The requested plugin call cannot be made as I am unable to determine the relationship between the plugin context and the requested result context.');
                end
                parent_context = higher_context_mapping.Context;
                parent_image_template = obj.ImageTemplates.GetTemplateImage(parent_context, dataset_stack, reporting);
                new_edited_image_for_context_image = edited_result_image.Copy;
                new_edited_image_for_context_image.ResizeToMatch(parent_image_template);
                obj.SaveEditedResultRecursive(plugin_name, parent_context, new_edited_image_for_context_image, plugin_info, dataset_stack, dataset_uid, reporting);

            % If the plugin's context set is lower in the hierarchy, then get
            % reduce the edited image input to each lower context in turn and save to each
            elseif plugin_context_set_mapping.IsOtherContextSetHigher(input_context_set_mapping)

                child_context_mappings = input_context_mapping.Children;

                for child_mapping = child_context_mappings
                    child_context = child_mapping{1}.Context;

                    % Create a new image from the edited result, reduced to this
                    % context
                    new_edited_image_for_context_image = MimReduceResultToContext(edited_result_image, child_mapping{1}.Context, obj.ImageTemplates, dataset_stack, reporting);

                    % Save this new image
                    obj.SaveEditedResultRecursive(plugin_name, child_context, new_edited_image_for_context_image, plugin_info, dataset_stack, dataset_uid, reporting);
                end

            else
                reporting.Error('MimContextHierarchy:UnexpectedSituation', 'The requested plugin call cannot be made as I am unable to determine the relationship between the plugin context and the requested result context.');
            end
        end

        function [valid, edited_result_exists] = CheckDependencyValid(obj, next_dependency, reporting)
            % Checks a given dependency against the cached values to
            % determine if it is valid (ie it depends on the most recent
            % computed results)

            plugin_results_info = obj.DatasetDiskCache.GetPluginResultsInfo();

            if isfield(next_dependency.Attributes, 'IsParameter')
                is_parameter = next_dependency.Attributes.IsParameter;
            else
                is_parameter = false;
            end
            % ToDo: process dependency values

            if isfield(next_dependency.Attributes, 'IsEditedResult')
                is_edited_result = next_dependency.Attributes.IsEditedResult;
            else
                is_edited_result = false;
            end

            if isfield(next_dependency.Attributes, 'IsManualSegmentation')
                is_manual = next_dependency.Attributes.IsManualSegmentation;
            else
                is_manual = false;
            end

            if isfield(next_dependency.Attributes, 'IsMarkerSet')
                is_marker = next_dependency.Attributes.IsMarkerSet;
            else
                is_marker = false;
            end

            if is_manual
                type_of_dependency = 'manual segmentation';
            else
                type_of_dependency = 'plugin';
            end

            edited_result_exists = plugin_results_info.CachedInfoExists(next_dependency.PluginName, next_dependency.Context, MimCacheType.Edited);

            if is_edited_result
                cache_type = MimCacheType.Edited;
            elseif is_manual
                cache_type = MimCacheType.Manual;
            elseif is_marker
                cache_type = MimCacheType.Markers;
            else
                cache_type = MimCacheType.Results;
            end

            % The full list should always contain the most recent dependency
            % uid, unless the dependencies file was deleted
            if ~plugin_results_info.CachedInfoExists(next_dependency.PluginName, next_dependency.Context, cache_type)
                reporting.Log(['No dependency record for this ' type_of_dependency ' - forcing re-run.']);
                valid = false;
                return;
            end

            current_info = plugin_results_info.GetCachedInfo(next_dependency.PluginName, next_dependency.Context, cache_type);
            current_dependency = current_info.InstanceIdentifier;

            if current_info.IgnoreDependencyChecks
                reporting.LogVerbose(['Ignoring dependency checks for ' type_of_dependency ' ' next_dependency.PluginName '(' char(next_dependency.Context) ')']);
            else
                % Sanity check - this case should never occur
                if ~strcmp(next_dependency.DatasetUid, current_dependency.DatasetUid)
                    reporting.Error('MimPluginResultsInfo:DatsetUidError', 'Code error - not matching dataset UID during dependency check');
                end

                if ~strcmp(next_dependency.Uid, current_dependency.Uid)
                    reporting.Log('Mismatch in dependency version uids - forcing re-run');
                    valid = false;
                    return;
                else
                    reporting.LogVerbose(['Dependencies Ok for ' type_of_dependency ' ' next_dependency.PluginName]);
                end
            end

            valid = true;

        end

    end

    methods (Access = private)

        function valid = CheckDependenciesValid(obj, linked_dataset_chooser, dependencies, dataset_stack, parameters_for_next_plugin_call, reporting)
            % Checks the dependencies in this result with the current dependency
            % list, and determine if the dependencies are still valid

            dependency_list = dependencies.DependencyList;

            % Separate dependencies into parameters and plugin
            % dependencies, and also create list of dependencies which are
            % edited values
            plugin_dependencies = {};
            known_edited_values = containers.Map();
            parameter_list = {};
            for index = 1 : length(dependency_list)
                next_dependency = dependency_list(index);
                dataset_uid = next_dependency.DatasetUid;
                if ~known_edited_values.isKey(dataset_uid)
                    known_edited_values(dataset_uid) = {};
                end
                if isfield(next_dependency.Attributes, 'IsParameter') && (next_dependency.Attributes.IsParameter)
                    parameter_list{end + 1} = next_dependency;
                else
                    plugin_dependencies{end + 1} = next_dependency;
                    if ~obj.CheckPluginVersion(next_dependency, reporting)
                        valid = false;
                        reporting.Log(['A newer version of plugin ' next_dependency.PluginName ' has been found. This result must be regenerated.']);
                        return;
                    end
                    if isfield(next_dependency.Attributes, 'IsEditedResult') && (next_dependency.Attributes.IsEditedResult)
                        edited_list = known_edited_values(dataset_uid);
                        edited_list{end + 1} = next_dependency.PluginName;
                        known_edited_values(dataset_uid) = edited_list;
                    end
                end
            end

            if ~dataset_stack.CheckParameterDependencies(parameter_list, parameters_for_next_plugin_call, reporting)
                valid = false;
                return;
            end

            % Iterate through plugin dependencies
            for index = 1 : length(plugin_dependencies)
                next_dependency = plugin_dependencies{index};

                dataset_uid = next_dependency.DatasetUid;
                [valid, edited_result_exists] = linked_dataset_chooser.GetDataset(reporting, dataset_uid).CheckDependencyValid(next_dependency, reporting);
                if ~valid
                    valid = false;
                    return;
                end

                % If the dependency is based on a non-edited result, but we have an edited
                % result in the cache, then this plugin result is invalid
                if edited_result_exists && ~ismember(next_dependency.PluginName, known_edited_values(dataset_uid))
                    reporting.Log(['The result for dependency ' next_dependency.PluginName '(' char(next_dependency.Context) ') has been edited - forcing re-run for plugin.']);
                    valid = false;
                    return;
                end
            end

            valid = true;
        end

        function valid = CheckPluginVersion(obj, next_dependency, reporting)
            if isfield(next_dependency.Attributes, 'PluginVersion')
                dependency_version = next_dependency.Attributes.PluginVersion;
            else
                dependency_version = 1;
            end

            plugin_info = obj.PluginCache.GetPluginInfo(next_dependency.PluginName, [], reporting);
            valid = plugin_info.PluginVersion == dependency_version;
        end

        function [context_set, context_set_mapping] = GetContextSetMappings(obj, context)
            % Determines the type of context supported by the plugin

            context_set = context;
            if isempty(context_set)
                    context_set = obj.ContextDef.GetDefaultContextSet();
            end
            context_set_mapping = obj.ContextSets(char(context_set));
        end

        function [context, context_set, context_mapping, context_set_mapping] = GetContextMaps(obj, context)
            % Determines the context and context type requested by the calling function

            if isempty(context)
                context = obj.ContextDef.GetDefaultContext();
            end
            if obj.Contexts.isKey(char(context))
                context_mapping = obj.Contexts(char(context));
                context_set_mapping = context_mapping.ContextSet;
                context_set = context_set_mapping.ContextSet;
            else
                context_mapping = [];
                context_set_mapping = [];
                context_set = [];
            end
        end

        function output_image = GenerateImageFromResults(obj, result, plugin_class, linked_dataset_chooser, dataset_stack, reporting)
            template_callback = MimTemplateCallback(linked_dataset_chooser, dataset_stack, reporting);

            if isa(result, 'PTKImage')
                output_image = plugin_class.GenerateImageFromResults(result.Copy, template_callback, reporting);
            else
                output_image = plugin_class.GenerateImageFromResults(result, template_callback, reporting);
            end
        end

        function edited_cached_info = SaveEditedResult(obj, plugin_name, context, edited_result, dataset_uid, plugin_version, reporting)
            % Saves the result of a plugin after semi-automatic editing

            attributes = [];
            attributes.IgnoreDependencyChecks = false;
            attributes.IsEditedResult = true;
            attributes.PluginVersion = plugin_version;
            instance_identifier = PTKDependency(plugin_name, context, CoreSystemUtilities.GenerateUid, dataset_uid, attributes);
            edited_cached_info = obj.FrameworkAppDef.GetClassFactory.CreateDatasetStackItem(instance_identifier, PTKDependencyList(), false, false, reporting);
            edited_cached_info.MarkEdited();

            obj.DatasetDiskCache.SaveEditedResult(plugin_name, context, edited_result, edited_cached_info, reporting);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% REJKURT  - calculation of kutosis of a 1D, 2D or 3D array and
%              rejection of outliers values of the input data array
%              using the discrete kutosis of the values in that dimension.
%
% Usage:
%   >>  [kurtosis rej] = rejkurt( signal, threshold, kurtosis, normalize);
%
% Inputs:
%   signal     - one dimensional column vector of data values, two
%                dimensional column vector of values of size
%                sweeps x frames or three dimensional array of size
%                component x sweeps x frames. If three dimensional,
%                all components are treated independently.
%   threshold  - Absolute threshold. If normalization is used then the
%                threshold is expressed in standard deviation of the
%                mean. 0 means no threshold.
%   kurtosis   - pre-computed kurtosis (only perform thresholding). Default
%                is the empty array [].
%   normalize  - 0 = do not not normalize kurtosis. 1 = normalize kurtosis.
%                2 is 20% trimming (10% low and 10% high) kurtosis before
%                normalizing. Default is 0.
%
% Outputs:
%   kurtosis    - normalized joint probability  of the single trials
%                (same size as signal without the last dimension)
%   rej         - rejected matrix (0 and 1, size: 1 x sweeps)
%
% Remarks:
%   The exact values of kurtosis depend on the size of a time
%   step and thus cannot be considered as absolute.
%   This function uses the kurtosis function from the statistival
%   matlab toolbox. If the statistical toolbox is not installed,
%   it uses the 'kurt' function of the ICA/EEG toolbox.
%
% See also: KURT, KURTOSIS

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

function [kurto, rej] = rejkurt( signal, threshold, oldkurtosis, normalize);

if nargin < 1
	help rejkurt;
	return;
end;
if nargin < 2
	threshold = 0;
end;
if nargin < 4
	normalize = 0;
end;
if nargin < 3
	oldkurtosis = [];
end;

if size(signal,2) == 1 % transpose if necessary
	signal = signal';
end

nbchan = size(signal,1);
pnts = size(signal,2);
sweeps = size(signal,3);
kurto = zeros(nbchan,sweeps);

if ~isempty( oldkurtosis ) % speed up the computation
	kurto = oldkurtosis;
else
	for rc = 1:nbchan
		% compute all kurtosis
		% --------------------
		for index=1:sweeps
			try
			    kurto(rc, index) = kurtosis(signal(rc,:,index));
			catch
				kurto(rc, index) = kurt(signal(rc,:,index));
			end;
		end
	end

	% normalize the last dimension
	% ----------------------------
	if normalize
        tmpkurt = kurto;
        if normalize == 2,
            tmpkurt = sort(tmpkurt);
            minind  = max(round(length(tmpkurt)*0.1),1);
            maxind  = round(length(tmpkurt)-round(length(tmpkurt)*0.1));
            if size(tmpkurt,2) == 1
                 tmpkurt = tmpkurt(minind:maxind);
            else tmpkurt = tmpkurt(:,minind:maxind);
            end
        end
	    switch ndims( signal )
	    	case 2,	kurto = (kurto-mean(tmpkurt)) / std(tmpkurt);
	    	case 3,	kurto = (kurto-mean(tmpkurt,2)*ones(1,size(kurto,2)))./ ...
				        (std(tmpkurt,0,2)*ones(1,size(kurto,2)));
		end
	end
end

% reject
% ------
if threshold(1) ~= 0
    if length(threshold) > 1
    	rej = (threshold(1) > kurto) | (kurto > threshold(2));
    else
    	rej = abs(kurto) > threshold;
    end
else
	rej = zeros(size(kurto));
end;

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [prob,sol,fmin] = milp_prob(varargin)
%MILP_PROB  Return an OPTI MILP
%
%   prob = milp_prob(no) return a pre-built optiprob of a saved MILP.
%
%   [prob,sol,fmin] = milp_prob(no) returns the optimum solution and
%   function eval at the optimum
%
%   no = milp_prob() returns the number of problems available for testing.

%   (C) 2011 Jonathan Currie (I2C2)

%Check if just returning no problems
if(nargin < 1)
    prob = 10; sol = []; fmin = [];
    return;
else
    no = varargin{1};
end

%Big switch yard
switch(no)
    case 1
        f = -[-1, 2]';
        A = [2, 1;-4, 4];
        b = [5, 5]';
        e = -[1, 1];
        xint = 'II';
        prob = optiprob('f',f,'mix',A,b,e,'int',xint);
        sol = [1;2];
        fmin = -3;

    case 2
        f = -[50, 100];
        A = [10, 5;4, 10; 1, 1.5];
        b = [2500, 2000, 450]';
        e = [-1, -1, -1];
        xint = 'II';
        prob = optiprob('f',f,'mix',A,b,e,'int',xint);
        sol = [187;125];
        fmin = -21850;

    case 3
        f = [40, 36];
        A = [5, 3];
        b = 45;
        e = 1;
        ub = [8, 10];
        xint = 'II';
        prob = optiprob('f',f,'mix',A,b,e,'ub',ub,'int',xint);
        sol = [8;2];
        fmin = 392;

    case 4
        f = [3, -7, -12];
        A = [-3, 6, 8;6, -3, 7;-6, 3, 3];
        b = [12, 8, 5]';
        e = [-1, -1, -1];
        xint = 'III';
        prob = optiprob('f',f,'mix',A,b,e,'int',xint);
        sol = [2;3;0];
        fmin = -15;

    case 5
        f = [2, 3, 7, 7];
        A = [1, 1, -2, -5;-1, 2, 1, 4];
        b = [2, 3]';
        e = [1, 1];
        lb = zeros(4,1);
        ub = [30 100 20 1]';
        xint = 'CICI';
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'int',xint);
        sol = [0;2;0;0];
        fmin = 6;

    case 6
        f = [1, 2, 3, 7, 8, 8];
        A = [5, -3, 2, -3, -1, 2; -1, 0, 2, 1, 3, -3;1, 2, -1, 0, 5, -1];
        b = [-5, -1, 3]';
        e = [1, 1, 1];
        lb = zeros(6,1);
        ub = 10*ones(6,1);
        xint = 'IIIIII';
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'int',xint);
        sol = [1;1;0;0;0;0];
        fmin = 3;

    case 7
        n = 40;
        t = (0:n-1)';
        y = 3.5 -.2*t;
        b = y + 0.5*ones(size(y));
        m = [ones(n,1),t(:)];
        A = [m,-m,eye(n)];
        f = [sum(m),sum(-m),2*ones(1,n)];
        e = ones(n,1);
        lb = zeros(n+4,1);
        ub = [10, 10, 10, 10, 5*ones(1,n)];
        xint = [1 3];
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'int',xint);
        sol = [4;0;0;0.2;zeros(40,1)];
        fmin = 4.0000000000000044;

    case 8
        f = -[8, 15];
        A = [10, 21;2, 1];
        b = [156, 22]';
        e = [-1, -1];
        xint = 'IC';
        prob = optiprob('f',f,'mix',A,b,e,'int',xint);
        sol = [9;3.14285714285714];
        fmin = -119.14285714285714;

    case 9
        f = -[3, 13];
        A = [2, 9;11, -8];
        b = [40, 82]';
        e = [-1, -1];
        xint = 'IC';
        prob = optiprob('f',f,'mix',A,b,e,'int',xint);
        sol = [9;2.44444444444444];
        fmin = -58.77777777777778;

    case 10
        f = -[592, 381, 273, 55, 48, 37, 23];
        A = [3534, 2356, 1767, 589, 528, 451, 304];
        b = 119567;
        e = -1;
        lb = zeros(7,1);
        ub = [100 50 33 20 77 44 20]';
        xint = 'IIIIIII';
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'int',xint);
        sol = [32;2;1;0;0;0;0];
        fmin = -19979;

    otherwise
        error('Problem not available or not implemented yet');
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ indx_extract, n, indx ] = r8vec_indexed_heap_d_extract ( n, a, indx )

%*****************************************************************************80
%
%% R8VEC_INDEXED_HEAP_D_EXTRACT: extract from heap descending indexed R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    An indexed R8VEC is an R8VEC of data values, and an R8VEC of N indices,
%    each referencing an entry of the data vector.
%
%    The routine finds the maximum value in the heap, returns that value to the
%    user, deletes that value from the heap, and restores the heap to its
%    proper form.
%
%    Note that the argument N must be a variable, which will be decremented
%    before return, and that INDX will hold one less value on output than it
%    held on input.
%
%    This is one of three functions needed to model a priority queue.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Thomas Cormen, Charles Leiserson, Ronald Rivest,
%    Introduction to Algorithms,
%    MIT Press, 2001,
%    ISBN: 0262032937,
%    LC: QA76.C662.
%
%  Parameters:
%
%    Input, integer N, the number of items in the index vector.
%
%    Input, real A(*), the data vector.
%
%    Input, integer INDX(N), the index vector.
%
%    Output, integer INDX_EXTRACT, the index in A of the item of
%    maximum value, which has now been removed from the heap.
%
%    Output, integer N, the number of items in the revised index vector.
%
%    Output, integer INDX(N), the revised index vector.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_INDEXED_HEAP_D_EXTRACT - Fatal error!\n' );
    fprintf ( 1, '  The heap is empty.\n' );
    error ( 'R8VEC_INDEXED_HEAP_D_EXTRACT - Fatal error!' )
  end
%
%  Get the index of the maximum value.
%
  indx_extract = indx(1);

  if ( n == 1 )
    n = 0;
    indx = [];
    return
  end
%
%  Shift the last index down.
%
  indx(1) = indx(n);
%
%  Restore the heap structure.
%
  n = n - 1;
  indx = r8vec_indexed_heap_d ( n, a, indx );

  return
end

author: johannesgerer