text:
function tapas_uniqc_new_struct(varargin)
% Creates a new script (e.g. tests/examples) including header.
% tapas_uniqc_new_struct(funname) opens the editor and pastes the content
% of a user-defined template into the file funname.m.
%
%   Example
%       tapas_uniqc_new_struct myscript
%           OR
%       tapas_uniqc_new_struct myscript author
%
%   opens the editor and pastes the following
%
% 	% script myscript
% 	%MYSCRIPT  One-line description here, please.
% 	%
% 	%   Example
% 	%   MYSCRIPT
% 	%
% 	%   See also
%
% 	% Author: Your name
% 	% Created: 2005-09-22
% 	% Copyright 2005 Your company.
%
%   See also edit, mfiletemplate

% Author: Peter (PB) Bodin
% Created: 2005-09-22
% Modified: 2019-01-05 (Saskia Bollmann and Lars Kasper)

	% See the variables repstr, repwithstr and tmpl to figure out how
	% to design your own template.
	% Edit tmpl to your liking, if you add more tokens in tmpl, make
	% sure to add them in repstr and repwithstr as well.

	% I made this function just for fun to check out some java handles to
	% the editor. It would probably be better to fprintf the template
	% to a new file and then call edit, since the java objects might change
	% names between versions.

	switch nargin
		case 0
			edit
			warning('tapas_uniqc_new_struct without argument is the same as edit')
			return;
		case 1
			fname=varargin{:};
			edit(fullfile(pwd,fname));
            authors = 'Saskia Bollmann & Lars Kasper'; %default authors, set this further down in function authors
        case 2
            fname = varargin{1};
            authors = varargin{2};
		otherwise
			error('tapas:uniqc:TooManyInputArguments', ...
                'too many input arguments')
	end

	try lasterror
		edhandle=com.mathworks.mlservices.MLEditorServices;

        % R2009a => 2009.0, R2009b = 2009.5
        vs = version('-release');
        v = str2double(vs(1:4));
        if vs(5)=='b'
            v = v + 0.5;
        end

        if v < 2009.0
            edhandle.builtinAppendDocumentText(strcat(fname,'.m'),parse(fname,authors));
        else
            edhandle.getEditorApplication.getActiveEditor.appendText(parse(fname, authors));
        end
	catch
		rethrow(lasterror)
	end

	function out = parse(func, authors)

		tmpl={ ...
  'classdef $filename'
            '% This is a predefined options/parameter structure for CLASS_OR_FUNCTION'
            '%'
            '%   outputStruct = $filename(input)'
            '%'
            '%'
            '%   See also CLASS_OR_FUNCTION'
            ' '
            '% Author:   $author'
            '% Created:  $date'
            '% Copyright (C) $year $institute'
            '%                    $company'
            '%'
            '% This file is part of the TAPAS UniQC Toolbox, which is released'
            '% under the terms of the GNU General Public License (GPL), version 3. '
            '% You can redistribute it and/or modify it under the terms of the GPL'
            '% (either version 3 or, at your option, any later version).'
            '% For further details, see the file COPYING or'
            '%  <http://www.gnu.org/licenses/>.'
            ''
            'properties'
            ''
            '   %INSERT_PROPERTY_DESCRIPTION_ABOVE_EACH_PROPERTY'
            ''
            'end'
            'methods'
            '   function this = $filename(varargin)'
            '   % Constructor of class, allows ''PropertyName'', ''PropertyValue'''
            '   % pairs and predefined shortcuts'
            '   % everything else set to default ([])'
            '       if nargin'
            '           for cnt = 1:nargin/2 % save them to object properties'
            '               this.(varargin{2*cnt-1}) = varargin{2*cnt};'
            '           end'
            '       end'
            '   end'
            'end'
            'end'
            };

		repstr={...
			'$filename'
			'$FILENAME'
			'$date'
			'$year'
			'$author'
            '$institute'
			'$company'};

		repwithstr={...
			func
			upper(func)
			datestr(now,29)
			datestr(now,10)
			authors
			'Institute for Biomedical Engineering'
            'University of Zurich and ETH Zurich'};

		for k = 1:numel(repstr)
			tmpl = strrep(tmpl,repstr{k},repwithstr{k});
		end
		out = sprintf('%s\n',tmpl{:});
	end
end

author: translationalneuromodeling

repo: tapas


text:
function [X,info] = IRrrgmres(A,b,varargin)
% IRrrgmres Range Restricted GMRES for square systems
%
% options  = IRrrgmres('defaults')
% [X,info] = IRrrgmres(A,b);
% [X,info] = IRrrgmres(A,b,K);
% [X,info] = IRrrgmres(A,b,options);
% [X,info] = IRrrgmres(A,b,K,options);
%
% This function implements the Range Restricted GMRES method for solving
% a square system.
%
% With 'defaults' as input returns the default options.  Otherwise outputs
% the iterates specified in K, using max(K) as MaxIter, and using all other
% default options.  With options as input: uses the user-specified options
% and all the other default options.
%
% Inputs:
%  A : either (a) a full or sparse matrix
%             (b) a matrix object that performs the matrix*vector operation
%             (c) user-defined function handle
%  b : right-hand side vector
%  K : (optional) integer vector that specifies which iterates are returned
%      in X; the maximum number of iterations is assumed to be max(K)
%      [ positive integer | vector of positive components ]
%  options : structure with the following fields (optional)
%      x0         - initial guess for the iterations; default = zero vector
%                   [ array | {'none'} ]
%      MaxIter    - maximum allowed number of iterations
%                   [ {100} | positive integer ]
%                   NOTE: K overrules MaxIter if both are assigned
%      x_true     - true solution; allows us to returns error norms with
%                   respect to x_true at each iteration
%                   [ array | {'none'} ]
%      NoiseLevel - norm of noise in rhs divided by norm of rhs
%                   [ {'none'} | nonnegative scalar ]
%      eta        - safety factor for the discrepancy principle
%                   [ {1.01} | scalar greater than (and close to) 1 ]
%      DecompOut  - return the decomposition to the user
%                   [ 'on' | {'off'} ]
%      IterBar    - shows the progress of the outer iterations
%                   [ {'on'} | 'off' ]
%      NoStop     - specifies weather the iterations should proceed
%                   after a stopping criterion has been satisfied
%                   [ 'on' | {'off'} ]
% Note: the options structure can be created using the function IRset.
%
% Outputs:
%   X : computed solutions, stored column-wise (at the iterations listed in K)
%   info: structure with the following fields:
%      its     - number of the last computed iteration
%      saved_iterations - iteration numbers of iterates stored in X
%      StopFlag - stringg that describes the output/stopping condition:
%                   * Performed max number of iterations
%                   * Residual tolerance satisfied (discrepancy principle)
%                   * Normal equationa residual tolerance satisfied
%      Rnrm     - relative residual norms at each iteration
%      Xnrm     - solution norms at each iteration
%      Enrm     - error norms (requires x_true) at each iteration
%      V        - matrix generated by the Arnoldi algorithm: the orthonormal
%                 columns of V are a basis for the RRGMRES solution
%      Q        - orthogonal matrix in QR factorization of Hessenberg matrix
%      T        - triangular matrix in QR factorization of Hessenberg matrix
%      StopReg  - struct containing information about the solution that
%                 satisfies the stopping criterion.  Fields:
%                   It   : iteration where stopping criterion is satisfied
%                   X    : solution satisfying the stopping criterion
%                   Enrm : the corresponding relative error (requires x_true)
%      BestReg  - struct containing information about the solution that
%                 minimizes Enrm (requires x_true). Fields:
%                   It   : iteration where the minimum is attained
%                   X    : best solution
%                   Enrm : best relative error
%
% See also: IRcgls, IRhybrid_gmres, IRget, IRset

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('x0','none', 'MaxIter',100, 'x_true','none', ...
    'NoiseLevel','none', 'eta',1.01, 'IterBar','on', ...
    'NoStop','off', 'DecompOut','off');

% If input is 'defaults,' return the default options in X.
if nargin==1 && nargout <= 1 && isequal(A,'defaults')
    X = defaultopt;
    return;
end

defaultopt.restart = 'off';
defaultopt.verbosity = 'on';

% Check for acceptable number of optional input arguments.
switch length(varargin)
    case 0
        K = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = [];
        else
            K = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = varargin{2};
        else
            K = varargin{2}; options = varargin{1};
        end
        if isfield(options, 'MaxIter') && ~isempty(options.MaxIter) && (~isempty(K) && options.MaxIter ~= max(K))
            warning('The value of MaxIter is discarded; the maximum value in K is taken as MaxIter')
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = IRset(defaultopt, options);

MaxIter    = IRget(options, 'MaxIter',    [], 'fast');
x_true     = IRget(options, 'x_true',     [], 'fast');
NoiseLevel = IRget(options, 'NoiseLevel', [], 'fast');
eta        = IRget(options, 'eta',        [], 'fast');
IterBar    = IRget(options, 'IterBar',    [], 'fast');
NoStop     = IRget(options, 'NoStop',     [], 'fast');
restart    = IRget(options, 'restart',    [], 'fast');
verbose    = IRget(options, 'verbosity',  [], 'fast');
DecompOut  = IRget(options, 'DecompOut',  [], 'fast');

restart = strcmp(restart, 'on');
verbose = strcmp(verbose, 'on');

% Setting K.
if isempty(K)
    K = MaxIter;
end
% Sorting the iterations (in case they are shuffled in input).
K = K(:); K = sort(K,'ascend'); K = unique(K);
if ~((isreal(K) && (all(K > 0)) && all(K == floor(K))))
    error('K must be a vector of positive real integers')
end
if K(end) ~= MaxIter
    MaxIter = K(end);
end

StopIt = MaxIter;

if isempty(NoiseLevel) || strcmp(NoiseLevel,'none')
    Rtol = 0;
else
    Rtol = eta*NoiseLevel;
end

n = length(b);
test_sq = ones(n,1);
try
    test_sq = A_times_vec(A, test_sq);
    if (length(test_sq)~=n)
        error('The matrix A shuold be square; check the length of b.')
    end
catch
    error('The matrix A must be square; check the length of b.')
end

% See if an initial guess is given; if not then use 0 as initial guess.
x0 = IRget(options, 'x0', [], 'fast');
Ab = A_times_vec(A, b);
nrmb = norm(b(:));
if strcmp(x0,'none')
    x0 = zeros(n,1);
    r = b;
    rr = Ab;
else
    try
        Ax = A_times_vec(A, x0);
        r = b - Ax;
        rr = Ab - Ax;
    catch
        error('Check the length of x')
    end
end
if restart
    ktotcount  = IRget(options, 'ktotcount', [],'fast');
    TotIterMax = IRget(options, 'TotIterMax',[],'fast');
    if strcmp(TotIterMax, 'none') || TotIterMax < MaxIter
        TotIterMax = MaxIter;
    end
    if strcmp(ktotcount, 'none')
        error('The total iteration counter must be assigned')
    end
    Ktot = IRget(options, 'Ktot', [], 'fast');
    % No checks on Ktot, it should be given from IRrestart.
end

% Declare matrices and prepare for the iterations.
X       = zeros(n,length(K));
Xnrm    = zeros(MaxIter,1);
Rnrm    = zeros(MaxIter,1);
V       = zeros(n,MaxIter); % Orthonormal vectors spanning Krylov subspace
h       = zeros(MaxIter,1); % New column of Hessenberg matrix H
Q       = zeros(MaxIter+1); % H = Q*T, Q orthogonal
T       = zeros(MaxIter);   % H = Q*T, T upper triangular
invT    = zeros(MaxIter);   % Inverse of the matrix T
coeff   = zeros(MaxIter,1); % Solution of the projected problem, i.e.,
                            % coefficient w.r.t. the Krylov basis.
rhs = zeros(MaxIter+1,1);   % Projected right-hand side

if strcmp(x_true,'none')
    errornorms = false;
else
    errornorms = true;
    Enrm = zeros(max(K),1);
    nrmtrue = norm(x_true(:));
    BestReg.It = [];
    BestReg.X = [];
    BestReg.Enrm = [];
    BestEnrm = 1e10;
    BestReg.Xnrm = [];
    BestReg.Rnrm = [];
end

NoStop = strcmp(NoStop,'on');

if restart
    saved_iterations = zeros(1, length(Ktot));
else
    saved_iterations = zeros(1, length(K));
end

% Prepare for the iterations.
alpharr  = norm(rr);
V(:,1) = rr/alpharr;
Q(1,1) = 1;
rhs(1) = V(:,1)'*r;
y = 0;

% Iterate.
noIterBar = strcmp(IterBar,{'off'});
if ~noIterBar
  h_wait = waitbar(0, 'Running iterations, please wait ...');
end
j = 0;
for k=1:MaxIter
    if restart, ktotcount = ktotcount + 1; end
    if ~noIterBar
        waitbar(k/MaxIter, h_wait)
    end
    w = A_times_vec(A, V(:,k));
    % Modified Gram-Schmidt on the new vector.
    for ll=1:k
        h(ll) = V(:,ll)'*w;
        w = w - V(:,ll)*h(ll);
    end
    alpha = norm(w);
    % Store new Arnoldi vector and update projected rhs.
    V(:,k+1) = w/alpha;
    rhs(k+1) = V(:,k+1)'*r;
    beta = rhs(1:k+1);
    % Apply previous rotations to h.
    T(1:k,k) = Q(1:k,1:k)'*h(1:k);
    % Compute Givens rotation parameters.
    rc = T(k,k);
    if alpha == 0
        c = 1; s = 0;
    elseif abs(alpha) > abs(rc)
        tau = -rc/alpha;
        s = 1 / sqrt(1 + abs(tau)^2);
        c = s*tau;
    else
        tau = -alpha/rc;
        c = 1 / sqrt(1 + abs(tau)^2);
        s = c*tau;
    end
    % Apply givens rotations.
    T(k,k) = c'*rc - s'*alpha;
    Q(1:k,[k,k+1]) = Q(1:k,k)*[c s];
    Q(k+1,[k,k+1]) = [-s c];
    if abs(T(k,k)) <= eps
        disp('Hessenberg matrix is (numerically) singular')
        X(:,j+1) = x;
        X = X(:,1:j+1);
        V = V(:,1:k+1);
        if restart
            saved_iterations(j+1) = ktotcount-1;
        else
            saved_iterations(j+1) = k-1;
        end
        saved_iterations = saved_iterations(1:j+1);
        % if k-1 < StopIt, StopIt = k-1; end
        if k>1
            Xnrm = Xnrm(1:k-1);
            Rnrm = Rnrm(1:k-1);
            if errornorms, Enrm = Enrm(1:k-1); end
        end
        % Stop because the Hessenberg matrix is (numerically) rank def.
        if StopIt == MaxIter
            StopFlag = 'Breakdown of the Arnoldi algorithm';
            StopReg.It = k-1;
            StopReg.X = x;
            if errornorms, StopReg.Enrm = Enrm(k-1); end
        end
        k = k-1;
        break
    end
    % Update the inverse of T.
    invT(1:k-1,k) = -(invT(1:k-1,1:k-1)*T(1:k-1,k))/T(k,k);
    invT(k,k) = 1/T(k,k);
    % Update the previous projected solution.
    coeff(1:k-1) = y;
    y = coeff(1:k) + (Q(1:k+1,k)'*beta)*invT(1:k,k);
    % Update solution.
    x = x0 + V(:,1:k)*y;
    % Compute norms.
    Xnrm(k)    = norm(x(:));
    if k==1
        Rnrm(k) = sqrt( abs( r'*r - abs(Q(1:k+1,k)'*beta)^2 ) )/nrmb;
    else
        Rnrm(k) = sqrt( abs( (Rnrm(k-1)*nrmb)^2 - abs(Q(1:k+1,k)'*beta)^2 ) )/nrmb;
    end
    if errornorms
        Enrm(k) = norm(x_true(:)-x(:))/nrmtrue;
        if Enrm(k)<BestEnrm
            BestReg.It = k;
            BestReg.X = x;
            BestEnrm = Enrm(k);
            BestReg.Enrm = BestEnrm;
            BestReg.Xnrm = Xnrm(k);
            BestReg.Rnrm = Rnrm(k);
        end
    end
    AlreadySaved = 0;
    if any(k==K)
        j = j+1;
        X(:,j) = x;
        saved_iterations(j) = k;
        if restart, saved_iterations(j) = ktotcount; end
        AlreadySaved = 1;
    end
    if restart
        if any(ktotcount == Ktot) && ~ AlreadySaved
            j = j+1;
            X(:,j) = x;
            saved_iterations(j) = ktotcount;
            AlreadySaved = 1;
        end
        if ktotcount == TotIterMax
            if ~ AlreadySaved
                j = j+1;
                saved_iterations(j) = ktotcount;
                X(:,j) = x;
                AlreadySaved = 1;
            end
            StopIt = k;
            StopReg.It = k;
            StopReg.X = x;
            if errornorms
                Enrm = Enrm(1:k);
                StopReg.Enrm = Enrm(k);
            end
            Xnrm     = Xnrm(1:k);
            Rnrm     = Rnrm(1:k);
            X = X(:,1:j);
            saved_iterations = saved_iterations(1:j);
            if verbose
                disp('Reached maximum number of iterations')
            end
            StopFlag = 'Reached maximum number of iterations';
            break
        end
    end
    if (Rnrm(k) <= Rtol) && (StopIt == MaxIter)
        if verbose
            disp('Residual tolerance satisfied')
        end
        StopFlag = 'Residual tolerance satisfied';
        if ~AlreadySaved && ~NoStop
            j = j+1;
            X(:,j) = x;
            if restart
                saved_iterations(j) = ktotcount;
            else
                saved_iterations(j) = k;
            end
            AlreadySaved = 1;
        end
        StopIt = k;
        StopReg.It = k;
        StopReg.X = x;
        if errornorms, StopReg.Enrm = Enrm(k); end
        if ~ NoStop
            Xnrm    = Xnrm(1:k);
            Rnrm    = Rnrm(1:k);
            if errornorms, Enrm = Enrm(1:k); end
            X = X(:,1:j);
            saved_iterations = saved_iterations(1:j);
            break
        end
    end
end
if k == MaxIter
  if StopIt == MaxIter
        % Stop because max number of iterations reached.
        if verbose
            disp('Reached maximum number of iterations')
        end
        StopFlag = 'Reached maximum number of iterations';
        if ~AlreadySaved
            j = j+1;
            X(:,j) = x;
            if restart
                saved_iterations(j) = ktotcount;
            else
                saved_iterations(j) = k;
            end
        end
        StopReg.It = k;
        StopReg.X = x;
        if errornorms, StopReg.Enrm = Enrm(k); end
        Xnrm    = Xnrm(1:k);
        Rnrm    = Rnrm(1:k);
        if errornorms, Enrm = Enrm(1:k); end
        X = X(:,1:j);
        saved_iterations = saved_iterations(1:j);
  end
end
if ~noIterBar, close(h_wait), end
if nargout==2
  info.its = k;
  if restart
      info.ktotcount = ktotcount;
  end
  info.saved_iterations = saved_iterations(1:j);
  info.StopFlag = StopFlag;
  info.StopReg = StopReg;
  info.Rnrm = Rnrm(1:k);
  info.Xnrm = Xnrm(1:k);
  if errornorms
    info.Enrm = Enrm(1:k);
    info.BestReg = BestReg;
  end
  if strcmp(DecompOut,'on')
    info.V = V;
    info.Q = Q;
    info.T = T;
  end
end

repo: IRtools


text:
function treepack_test005 ( )

%*****************************************************************************80
%
%% TREEPACK_TEST005 tests CATALAN and CATALAN_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TREEPACK_TEST005\n' );
  fprintf ( 1, '  CATALAN computes Catalan numbers.\n' );
  fprintf ( 1, '  CATALAN_VALUES returns some exact values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N  exact C(I)  computed C(I)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, c ] = catalan_values ( n_data );

    if ( n_data == 0 )
      break
    end

    c2 = catalan ( n );

    fprintf ( 1, '  %4d  %6d  %6d\n', n, c, c2(n+1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
close all; clear all; clc;
rng('default');

A = [3 2; 2 6];
X = [ 2 1 0 4;
     -2 2 0 4];
B = A * X;

solver = SPX_SteepestDescent(A, B);
x = solver.solve();
solver.printResults();

author: indigits

repo: sparse-plex


text:
function s = yjf_to_s_english ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_S_ENGLISH writes an English YJF date into a string.
%
%  Format:
%
%    BC OS YYYY/JJJ.FF
%    AD OS YYYY/JJJ.FF
%    AD NS YYYY/JJJ.FF
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, string S, a representation of the date.
%

%
%  Copy the input.
%
  y2 = y;
  j2 = j;
  f2 = f;
%
%  Check the input.
%
  [ y2, j2, f2, ierror ] = yjf_check_english ( y2, j2, f2 );

  if ( ierror ~= 0 )
    s = '?';
    return
  end

  if ( y2 < 0 )
    s = 'BC OS ';
    s1 = sprintf ( '%d', - y2 );
  elseif ( y < 1752 || ( y == 1752 && j < 278 ) )
    s = 'AD OS ';
    s1 = sprintf ( '%d', y2 );
  else
    s = 'AD NS ';
    s1 = sprintf ( '%d', y2 );
  end

  s = [ s, s1 ];

  s1 = sprintf ( '/%d', j2 );
  s = [ s, s1 ];

  s1 = sprintf ( '%2.2f', f );
  s = [ s, s1(2:end) ];

  return
end


text:
function Y = smoothn(X,sz,filt,std)

%SMOOTHN Smooth N-D data
%   Y = SMOOTHN(X, SIZE) smooths input data X. The smoothed data is
%       retuirned in Y. SIZE sets the size of the convolution kernel
%       such that LENGTH(SIZE) = NDIMS(X)
%
%   Y = SMOOTHN(X, SIZE, FILTER) Filter can be 'gaussian' or 'box' (default)
%       and determines the convolution kernel.
%
%   Y = SMOOTHN(X, SIZE, FILTER, STD) STD is a vector of standard deviations
%       one for each dimension, when filter is 'gaussian' (default is 0.65)

%     $Author: ganil $
%     $Date: 2001/09/17 18:54:39 $
%     $Revision: 1.1 $
%     $State: Exp $

if nargin == 2,
  filt = 'b';
elseif nargin == 3,
  std = 0.65;
elseif nargin>4 | nargin<2
  error('Wrong number of input arguments.');
end

% check the correctness of sz
if ndims(sz) > 2 | min(size(sz)) ~= 1
  error('SIZE must be a vector');
elseif length(sz) == 1
  sz = repmat(sz,ndims(X));
elseif ndims(X) ~= length(sz)
  error('SIZE must be a vector of length equal to the dimensionality of X');
end

% check the correctness of std
if filt(1) == 'g'
  if length(std) == 1
    std = std*ones(ndims(X),1);
  elseif ndims(X) ~= length(std)
    error('STD must be a vector of length equal to the dimensionality of X');
  end
  std = std(:)';
end

sz = sz(:)';

% check for appropriate size
padSize = (sz-1)/2;
if ~isequal(padSize, floor(padSize)) | any(padSize<0)
  error('All elements of SIZE must be odd integers >= 1.');
end

% generate the convolution kernel based on the choice of the filter
filt = lower(filt);
if (filt(1) == 'b')
  smooth = ones(sz)/prod(sz); % box filter in N-D
elseif (filt(1) == 'g')
  smooth = ndgaussian(padSize,std); % a gaussian filter in N-D
else
  error('Unknown filter');
end

% pad the data
X = padreplicate(X,padSize);

% perform the convolution
Y = convn(X,smooth,'valid');

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

function h = ndgaussian(siz,std)

% Calculate a non-symmetric ND gaussian. Note that STD is scaled to the
% sizes in SIZ as STD = STD.*SIZ

ndim = length(siz);
sizd = cell(ndim,1);

for i = 1:ndim
  sizd{i} = -siz(i):siz(i);
end

grid = gridnd(sizd);
std = reshape(std.*siz,[ones(1,ndim) ndim]);
std(find(siz==0)) = 1; % no smoothing along these dimensions as siz = 0
std = repmat(std,2*siz+1);

h = exp(-sum((grid.*grid)./(2*std.*std),ndim+1));
h = h/sum(h(:));

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

function argout = gridnd(argin)

% exactly the same as ndgrid but it accepts only one input argument of
% type cell and a single output array

nin = length(argin);
nout = nin;

for i=nin:-1:1,
  argin{i} = full(argin{i}); % Make sure everything is full
  siz(i) = prod(size(argin{i}));
end
if length(siz)<nout, siz = [siz ones(1,nout-length(siz))]; end

argout = [];
for i=1:nout,
  x = argin{i}(:); % Extract and reshape as a vector.
  s = siz; s(i) = []; % Remove i-th dimension
  x = reshape(x(:,ones(1,prod(s))),[length(x) s]); % Expand x
  x = permute(x,[2:i 1 i+1:nout]);% Permute to i'th dimension
  argout = cat(nin+1,argout,x);% Concatenate to the output
end

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

function b=padreplicate(a, padSize)
%Pad an array by replicating values.
numDims = length(padSize);
idx = cell(numDims,1);
for k = 1:numDims
  M = size(a,k);
  onesVector = ones(1,padSize(k));
  idx{k} = [onesVector 1:M M*onesVector];
end

b = a(idx{:});

author: Sable

repo: mcbench-benchmarks


text:
function geometry_test178 ( )

%*****************************************************************************80
%
%% TEST178 tests ROTATION_MAT2QUAT_3D, ROTATION_QUAT2MAT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  a = [ ...
    0.43301269, -0.5,        0.75; ...
    0.25,        0.86602539, 0.43301269; ...
   -0.86602539,  0.0,        0.5 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST178\n' );
  fprintf ( 1, '  ROTATION_MAT2QUAT_3D computes a rotation\n' );
  fprintf ( 1, '    quaternion from a rotation matrix.\n' );
  fprintf ( 1, '  ROTATION_QUAT2MAT_3D computes a rotation matrix\n' );
  fprintf ( 1, '    from a rotation quaternion.\n' );

  r8mat_print ( dim_num, dim_num, a, '  The rotation matrix:' );

  q = rotation_mat2quat_3d ( a );

  r8vec_print ( 4, q, '  The rotation quaternion:' );

  a = rotation_quat2mat_3d ( q );

  r8mat_print ( dim_num, dim_num, a, '  The rotation matrix:' );

  return
end


text:
function prob_test063 ( )

%*****************************************************************************80
%
%% TEST063 tests EXPONENTIAL_MEAN, EXPONENTIAL_SAMPLE, EXPONENTIAL_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST063\n' );
  fprintf ( 1, '  For the Exponential PDF:\n' );
  fprintf ( 1, '  EXPONENTIAL_MEAN computes the mean;\n' );
  fprintf ( 1, '  EXPONENTIAL_SAMPLE samples;\n' );
  fprintf ( 1, '  EXPONENTIAL_VARIANCE computes the variance.\n' );

  a = 1.0;
  b = 10.0;

  check = exponential_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST063 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = exponential_mean ( a, b );
  variance = exponential_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );
  fprintf ( 1, '  PDF mean =        %14f\n', mean );
  fprintf ( 1, '  PDF variance =    %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = exponential_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function c_nm = sphMultiplication(a_nm, b_nm, G)
%SPHMULTIPLICATION Computes coefficients of a product spherical function
%
%   The coefficients of a product function of two spherical functions can
%   be given directly as a linear combination of the coefficients of the
%   original functions, by employing Gaunt coefficients
%
%   a_nm: (N+1)^2 coefficients of spherical function a(\theta,\phi)
%   b_nm: (N+1)^2 coefficients of spherical function b(\theta,\phi)
%
%   c_nm:  (N+1)^2 coefficients of the product function
%           c(\theta,\phi) = a(\theta,\phi)*b(\theta,\phi)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Archontis Politis, 10/10/2013
%   archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% order
Na = sqrt(length(a_nm)) - 1;
Nb = sqrt(length(b_nm)) - 1;
Nc = Na+Nb;

c_nm = zeros((Nc+1)^2, 1);
% evaluate the coefficients of the product through the gaunt coefficients
if nargin<3
    G = gaunt_mtx(Na, Nb, Nc);
end
for n=0:Nc
    for m=-n:n
        q = n*(n+1)+m;
        c_nm(q+1) = a_nm.' * G(:,:,q+1) * b_nm;

    end
end

end

author: polarch

repo: Spherical-Harmonic-Transform


text:
function D = spm_eeg_convert(S)
% Convert various M/EEG formats to SPM12 format
% FORMAT D = spm_eeg_convert(S)
% S                - string (filename) or struct (see below)
%
% If S is a struct it can have the optional following fields:
% S.dataset        - file name
% S.mode           - 'header'     - only convert the header without reading data
%                    'continuous' - convert data as continuous
%                    'epoched'    - convert data as epoched (requires data that is
%                                   already epoched or a trial definition in S.trl).
% S.timewin        - for continuous mode [start end] of data segment in sec (all if empty)
%                  - for epoched mode time window in PST ms
% S.outfile        - output file name (default 'spmeeg_' + input)
% S.channels       - 'all' - convert all channels
%                    or cell array of labels
% For epoched mode:
%
% S.trl            - [N x 3] trl matrix or name of the trial definition file
%                    containing 'trl' variable with such a matrix
% S.conditionlabels- labels for the trials in the data [default: 'Undefined']
%
%   or
%
% S.trialdef       - structure array for trial definition with fields
%     S.trialdef.conditionlabel - string label for the condition
%     S.trialdef.eventtype      - string
%     S.trialdef.eventvalue     - string, numeric or empty
%
%
% S.inputformat    - data type (optional) to force the use of specific data
%                    reader
% S.eventpadding   - the additional time period around each trial for which
%                    the events are saved with the trial (to let the user
%                    keep and use for analysis events which are outside
%                    trial borders), in seconds. [default: 0]
% S.blocksize      - size of blocks used internally to split large files
%                    [default: ~100Mb]
% S.checkboundary  - 1 - check if there are breaks in the file and do not
%                        read across those breaks [default]
%                    0 - ignore breaks (not recommended).
% S.saveorigheader - 1 - save original data header with the dataset
%                    0 - do not keep the original header [default]
%
% % D              - MEEG object (also written on disk)
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_convert.m 7671 2019-10-02 12:34:50Z vladimir $

SVNrev = '$Rev: 7671 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG convert'); spm('Pointer', 'Watch');

if ischar(S)
    temp      = S;
    S         = [];
    S.dataset = temp;
end

if ~isfield(S, 'dataset')
    error('Dataset must be specified.');
end

if ~isfield(S, 'outfile'),         S.outfile = ['spmeeg_' spm_file(S.dataset,'basename')]; end
if ~isfield(S, 'channels'),        S.channels = 'all';                                     end
if ~isfield(S, 'timewin'),         S.timewin   = [];                                       end
if ~isfield(S, 'blocksize'),       S.blocksize = 3276800;                                  end %100 Mb
if ~isfield(S, 'checkboundary'),   S.checkboundary = 1;                                    end
if ~isfield(S, 'eventpadding'),    S.eventpadding = 0;                                     end
if ~isfield(S, 'saveorigheader'),  S.saveorigheader = 0;                                   end
if ~isfield(S, 'conditionlabels'), S.conditionlabels = 'Undefined' ;                       end
if ~isfield(S, 'inputformat'),     S.inputformat = [] ;                                    end

if ~iscell(S.conditionlabels)
    S.conditionlabels = {S.conditionlabels};
end

if ~isfield(S, 'mode') || ~isequal(S.mode, 'header')
    % The header is read here in a recursive call and the above if avoids reading
    % it twice which might be expensive for some formats
    S1 = [];
    S1.mode           = 'header';
    S1.dataset        = S.dataset;
    S1.outfile        = S.outfile;
    S1.inputformat    = S.inputformat;
    Dhdr              = spm_eeg_convert(S1);
    hdr               = Dhdr.hdr;
    event             = Dhdr.events;
    eventsamples      = Dhdr.events(':', 'samples');
else
    %--------- Read and check header
    hdr = ft_read_header(S.dataset, 'headerformat', S.inputformat);

    if isfield(hdr, 'label')
        [unique_label,junk,ind] = unique(hdr.label);
        if length(unique_label) ~= length(hdr.label)
            warning(['Data file contains several channels with ',...
                'the same name. These channels cannot be processed and will be disregarded']);
            % This finds the repeating labels and removes all their occurences
            sortind     = sort(ind);
            [junk,ind2] = setdiff(hdr.label, unique_label(sortind(diff(sortind)==0)));
            hdr.label   = hdr.label(ind2);
            hdr.nChans  = length(hdr.label);
        end
    end

    %--------- Read and prepare events

    try
        event = ft_read_event(S.dataset, 'detectflank', 'both', 'eventformat', S.inputformat);

        if ~isempty(strmatch('UPPT001', hdr.label))
            % This is s somewhat ugly fix to the specific problem with event
            % coding in FIL CTF. It can also be useful for other CTF systems where the
            % pulses in the event channel go downwards.
            fil_ctf_events = ft_read_event(S.dataset, 'detectflank', 'down', 'type', 'UPPT001', 'trigshift', -1, 'eventformat', S.inputformat);
            if ~isempty(fil_ctf_events)
                [fil_ctf_events(:).type] = deal('FIL_UPPT001_down');
                event = cat(1, event(:), fil_ctf_events(:));
            end
        end

        if ~isempty(strmatch('UPPT002', hdr.label))
            % This is s somewhat ugly fix to the specific problem with event
            % coding in FIL CTF. It can also be useful for other CTF systems where the
            % pulses in the event channel go downwards.
            fil_ctf_events = ft_read_event(S.dataset, 'detectflank', 'down', 'type', 'UPPT002', 'trigshift', -1, 'eventformat', S.inputformat);
            if ~isempty(fil_ctf_events)
                [fil_ctf_events(:).type] = deal('FIL_UPPT002_down');
                event = cat(1, event(:), fil_ctf_events(:));
            end
        end

        % This is another FIL-specific fix that will hopefully not affect other sites
        if isfield(hdr, 'orig') && isfield(hdr.orig, 'VERSION') && isequal(uint8(hdr.orig.VERSION),uint8([255 'BIOSEMI']))
            ind = strcmp('STATUS', {event(:).type});
            val = [event(ind).value];
            if ~isempty(val)
                bytes  = dec2bin(val, 8);
                bytes  = bytes(:, end-7:end);
                % This is a very specific criterion that assumes that
                % trigger code 1 is always used.
                if ~ismember('00000001', bytes, 'rows') && ismember('10000000', bytes, 'rows')
                    bytes = fliplr(bytes);
                end
                nval   = bin2dec(bytes);
                if (sum(val(:)>nval(:))/length(val))>0.5
                    nval    = num2cell(nval);
                    [event(ind).value] = deal(nval{:});
                end
            end
        end

    catch
        warning(['Could not read events from file ' S.dataset]);
        event = [];
    end

    % Replace samples with time
    if numel(event)>0
        for i = 1:numel(event)
            event(i).time = event(i).sample./hdr.Fs;
        end
    end
end

if ~isfield(S, 'mode')
    if (hdr.nTrials == 1)
        S.mode = 'continuous';
    else
        S.mode = 'epoched';
    end
end

%--------- Start making the header

D = [];
D.Fsample = hdr.Fs;

%--------- Select channels

if ~strcmp(S.channels, 'all') %Dhdr should be available in this case
    chansel = selectchannels(Dhdr, S.channels);
else
    if isfield(hdr, 'nChans')
        chansel = 1:hdr.nChans;
    else
        chansel = 1:length(hdr.label);
    end
end

nchan = length(chansel);

D.channels = repmat(struct('bad', 0), 1, nchan);

if isfield(hdr, 'label')
    [D.channels(:).label] = deal(hdr.label{chansel});
end
%--------- Preparations specific to reading mode (continuous/epoched)

if ismember(S.mode, {'continuous', 'header'})

    if isempty(S.timewin)
        if hdr.nTrials == 1
            segmentbounds = [1 hdr.nSamples];
        elseif ~S.checkboundary || isequal(S.mode, 'header')
            segmentbounds = [1 hdr.nSamples*hdr.nTrials];
        else
            error('The data cannot be read without ignoring trial borders');
        end
        timewindow = segmentbounds./D.Fsample;
    else
        timewindow = S.timewin;
        segmentbounds = round(timewindow.*D.Fsample);
        segmentbounds(1) = max(segmentbounds(1), 1);
    end

    %--------- Sort events and put in the trial

    if ~isempty(event)
        try
            event = rmfield(event, {'offset', 'sample'});
        end
        event = select_events(event, ...
            [timewindow(1)-S.eventpadding timewindow(2)+S.eventpadding]);
    end

    D.trials.label  = S.conditionlabels{1};
    D.trials.events = event;
    D.trials.onset  = timewindow(1);

    %--------- Break too long segments into blocks

    nblocksamples = floor(S.blocksize/nchan);
    nsampl = diff(segmentbounds)+1;

    trl = segmentbounds(1):nblocksamples:segmentbounds(2);
    if (trl(end)==segmentbounds(2))
        trl = trl(1:(end-1));
    end

    trl = [trl(:) [trl(2:end)-1 segmentbounds(2)]'];

    ntrial = size(trl, 1);

    readbytrials = 0;

    D.timeOnset = (trl(1,1)-1)./hdr.Fs;

    D.Nsamples = nsampl;
else % Read by trials
    if isfield(S, 'trl') || isfield(S, 'trialdef')
        if isfield(S, 'trl')
            if ischar(S.trl)
                trl = getfield(load(S.trl, 'trl'), 'trl');
                conditionlabels = getfield(load(S.trl, 'conditionlabels'), 'conditionlabels');
            else
                trl = S.trl;
                conditionlabels = S.conditionlabels;
            end
        else
            S1          = [];
            S1.D        = Dhdr;
            S1.timewin  = S.timewin;
            S1.trialdef = S.trialdef;
            S1.reviewtrials = 0;
            S1.save  = 0;
            [trl, conditionlabels] = spm_eeg_definetrial(S1);
        end

        trl = double(trl);

        if size(trl, 2) >= 3
            D.timeOnset = unique(trl(:, 3))./D.Fsample;
            trl = trl(:, 1:2);
        else
            D.timeOnset = 0;
        end

        if length(D.timeOnset) > 1
            error('All trials should have identical baseline');
        end

        if ~iscell(conditionlabels)
            conditionlabels = {conditionlabels};
        end

        if numel(conditionlabels) == 1
            conditionlabels = repmat(conditionlabels, 1, size(trl, 1));
        end

        readbytrials = 0;
    else
        try
            trialind = sort([strmatch('trial', {event.type}, 'exact'), ...
                strmatch('average', {event.type}, 'exact')]);
            trl = [eventsamples(trialind).sample];
            trl = double(trl(:));
            trl = [trl  trl+double([event(trialind).duration]')-1];

            try
                offset = unique([event(trialind).offset]);
            catch
                offset = [];
            end

            if length(offset) == 1 && offset~=0
                D.timeOnset = offset/D.Fsample;
            else
                D.timeOnset = -hdr.nSamplesPre/hdr.Fs;
            end
            conditionlabels = {};
            for i = 1:length(trialind)
                if isempty(event(trialind(i)).value)
                    conditionlabels{i} = S.conditionlabels{1};
                else
                    if all(ischar(event(trialind(i)).value))
                        conditionlabels{i} = event(trialind(i)).value;
                    else
                        conditionlabels{i} = num2str(event(trialind(i)).value);
                    end
                end
            end
            if  hdr.nTrials>1 && size(trl, 1)~=hdr.nTrials
                warning('Mismatch between trial definition in events and in data. Ignoring events');
                readbytrials = 1;
            else
                readbytrials = 0;
            end

            event = event(setdiff(1:numel(event), trialind));
        catch
            if hdr.nTrials == 1
                error('Could not define trials based on data. Use continuous option or trial definition file.');
            else
                readbytrials = 1;
            end
        end
    end
    if readbytrials
        nsampl = hdr.nSamples;
        ntrial = hdr.nTrials;
        trl = zeros(ntrial, 2);
        if exist('conditionlabels', 'var') ~= 1 || length(conditionlabels) ~= ntrial
            conditionlabels = repmat(S.conditionlabels, 1, ntrial);
        end
    else
        nsampl = unique(diff(trl, [], 2))+1;
        if length(nsampl) > 1
            error('All trials should have identical lengths');
        end

        inbounds = (trl(:,1)>=1 & trl(:, 2)<=hdr.nSamples*hdr.nTrials)';

        rejected = find(~inbounds);

        if ~isempty(rejected)
            trl = trl(inbounds, :);
            conditionlabels = conditionlabels(inbounds);
            warning([S.dataset ': Trials ' num2str(rejected) ' not read - out of bounds']);
        end

        ntrial = size(trl, 1);

        if ntrial == 0
            warning([S.dataset ': No trials to read. Bailing out.']);
            D = [];
            return;
        end
    end
    D.Nsamples = nsampl;
    if isfield(event, 'sample')
        event = rmfield(event, 'sample');
    end
end

%--------- Prepare for reading the data
outpath = spm_file(S.outfile,'fpath');
outfile = spm_file(S.outfile,'basename');
if isempty(outfile), outfile = 'spmeeg'; end

D.path = outpath;
D.fname = [outfile '.mat'];

if ~isequal(S.mode, 'header')

    % These are the units used for channel data in SPM
    chanunit  = units(Dhdr, chansel);
    chanunit(strcmp('T',   chanunit)) = {'fT'};
    chanunit(strcmp('T/m', chanunit)) = {'fT/mm'};
    chanunit(strcmp('V',   chanunit)) = {'uV'};

    if isequal(S.mode, 'continuous')
        D.data = file_array(fullfile(D.path, [outfile '.dat']), [nchan nsampl], 'float32-le');
    else
        D.data = file_array(fullfile(D.path, [outfile '.dat']), [nchan nsampl ntrial], 'float32-le');
    end

    % physically initialise file
    initialise(D.data);

    spm_progress_bar('Init', ntrial, 'reading and converting'); drawnow;
    if ntrial > 100, Ibar = floor(linspace(1, ntrial,100));
    else Ibar = 1:ntrial; end

    %--------- Read the data

    offset = 1;
    for i = 1:ntrial
        spm_progress_bar('Set','ylabel','reading...');
        if readbytrials
            dat = ft_read_data(S.dataset,'header',  hdr, 'begtrial', i, 'endtrial', i,...
                'chanindx', chansel, 'chanunit', chanunit, 'checkboundary', S.checkboundary, 'dataformat', S.inputformat);
        else
            dat = ft_read_data(S.dataset,'header',  hdr, 'begsample', trl(i, 1), 'endsample', trl(i, 2),...
                'chanindx', chansel, 'chanunit', chanunit, 'checkboundary', S.checkboundary, 'dataformat', S.inputformat);
        end

        % Sometimes ft_read_data returns sparse output
        dat = full(dat);

        spm_progress_bar('Set','ylabel','writing...');
        switch S.mode
            case 'continuous'
                nblocksamples = size(dat,2);

                D.data(:, offset:(offset+nblocksamples-1)) = dat;

                offset = offset+nblocksamples;
            case 'epoched'
                D.data(:, :, i) = dat;
                D.trials(i).label = conditionlabels{i};
                D.trials(i).onset = trl(i, 1)./D.Fsample;
                D.trials(i).events = select_events(event, ...
                    [ trl(i, 1)./D.Fsample-S.eventpadding  trl(i, 2)./D.Fsample+S.eventpadding]);
        end

        if ismember(i, Ibar)
            spm_progress_bar('Set', i);
        end

    end

    spm_progress_bar('Clear');

else
    D.data = [];
end

% Specify sensor positions and fiducials
if isfield(hdr, 'grad')
    D.sensors.meg = hdr.grad;
end

if isfield(hdr, 'elec')
    elec = hdr.elec;
else
    try
        elec = ft_read_sens(S.dataset, 'fileformat', S.inputformat);
        % It might be that read_sens will return the grad for MEG datasets
        if any(ismember({'ori', 'coilori', 'coilpos'}, fieldnames(elec)))
            elec = [];
        end
    catch
        elec = [];
    end
end

if ~isempty(elec)
    D.sensors.eeg = elec;
else
    sw = warning('off','backtrace');
    warning('Could not obtain electrode locations automatically.');
    warning(sw);
end

try
    D.fiducials = ft_convert_units(ft_read_headshape(S.dataset, 'fileformat', S.inputformat), 'mm');
catch
    sw = warning('off','backtrace');
    warning('Could not obtain fiducials automatically.');
    warning(sw);
end

%--------- Create meeg object
D = meeg(D);

% history
D = D.history('spm_eeg_convert', S);

if isfield(hdr, 'orig')
    if S.saveorigheader
        D.origheader = hdr.orig;
    end

    % Uses fileio function to get the information about channel types stored in
    % the original header.
    origchantypes = ft_chantype(hdr);
    [sel1, sel2] = spm_match_str(D.chanlabels, hdr.label);
    origchantypes = origchantypes(sel2);
    if length(strmatch('unknown', origchantypes, 'exact')) ~= numel(origchantypes)
        D.origchantypes = struct([]);
        D.origchantypes(1).label = hdr.label(sel2);
        D.origchantypes(1).type = origchantypes;
    end
end

S1 = [];
S1.task = 'defaulttype';
S1.D = D;
S1.updatehistory = 0;
D = spm_eeg_prep(S1);

% Assign default EEG sensor positions if possible
if ~isempty(strmatch('EEG', D.chantype, 'exact'))
    if isempty(D.sensors('EEG'))
        S1 = [];
        S1.task = 'defaulteegsens';
        S1.updatehistory = 0;
        S1.D = D;

        D = spm_eeg_prep(S1);
    else
        S1 = [];
        S1.task = 'project3D';
        S1.modality = 'EEG';
        S1.updatehistory = 0;
        S1.D = D;

        D = spm_eeg_prep(S1);
    end
end

% Create 2D positions for MEG
% by projecting the 3D positions to 2D
if ~isempty(strmatch('MEG', D.chantype)) && ~isempty(D.sensors('MEG'))
    S1 = [];
    S1.task = 'project3D';
    S1.modality = 'MEG';
    S1.updatehistory = 0;
    S1.D = D;

    D = spm_eeg_prep(S1);
end

% If channel units are available, store them.
if isequal(S.mode, 'header')
    [sel1, sel2] = spm_match_str(D.chanlabels, hdr.label);
    D = units(D, sel1, hdr.chanunit(sel2));
else
    D = units(D, ':', chanunit);
end

% The conditions will later be sorted in the original order they were defined.
if isfield(S, 'trialdef')
    D = condlist(D, {S.trialdef(:).conditionlabel});
end

% This is for the recursive call to work properly
if isequal(S.mode, 'header')
    D.hdr = hdr;
end

if ~isequal(S.mode, 'header')
    save(D);
end

%-Cleanup
%--------------------------------------------------------------------------
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','M/EEG convert: done'); spm('Pointer', 'Arrow');

%==========================================================================
% select_events
%==========================================================================
function event = select_events(event, timeseg)
% Utility function to select events according to time segment
% FORMAT event = select_events(event, timeseg)
if iscell(event)
    event = event{1};
end

if ~isempty(event)
    [time,ind] = sort([event(:).time]);

    selectind  = ind(time>=timeseg(1) & time<=timeseg(2));

    event      = event(selectind);
end

author: spm


text:
function [newModel,scores,stats] = boostingUpdate(data,labels,oldModel,binVals,bins,params)

% Learn classifier with all the data.

posEx = labels == 1;
negEx = ~posEx;
numPos = sum(posEx);
numNeg = sum(negEx);

if numPos<1 || numNeg<1,
  newModel = struct('dim',1,'error',0.5,'dir',1,'tr',0,'alpha',0);
  scores = zeros(1,numel(labels));
  return;
end

wt = ones(length(posEx),1);
posWt = numNeg/(numPos+numNeg);
negWt = numPos/(numPos+numNeg);
wt(posEx) = posWt;
wt(negEx) = negWt;
wt = wt./sum(wt);

modLabels = sign( (labels==1)-0.5);

updatedModel = updateOldWeights(data,modLabels,oldModel,wt);
[newModel,scores] = addNewRules(data,modLabels,updatedModel,binVals,bins,wt,params);

if nargout >= 3,
  % compute some statistics of how well training worked
  stats = ComputeBoostingStats(scores,labels);
end

function updatedModel = updateOldWeights(data,labels,oldModel,exWt)

wt = exWt;
updatedModel = oldModel;
scores = zeros(size(data,1),1);
for ndx = 1:numel(oldModel)
  curWkRule = oldModel(ndx);
  tr = curWkRule.tr;
  dir = curWkRule.dir;
  dim = curWkRule.dim;
  if dir>0,
    tt = ((data(:,dim)> tr)-0.5)*2;
  else
    tt = ((data(:,dim)<= tr)-0.5)*2;
  end
  curError = sum( (tt.*labels).*wt);
  if(curError<0);
    curError = 0;
  end
  updatedModel(ndx).error = 0.5-curError/2;
  updatedModel(ndx).alpha = 1-2*updatedModel(ndx).error;
  scores = scores + myBoostClassify(data,updatedModel(ndx));
  tt = scores.*labels;
  wt = exWt./(1+exp(tt));
  wt = wt./sum(wt);
end

function [newModel,scores] = addNewRules(data,labels,updatedModel,binVals,bins,exWt,params)

extraIters = params.iter_updates;
newModel = updatedModel;
scores = myBoostClassify(data,updatedModel);
tt = scores.*labels;
wt = exWt./(1+exp(tt));
wt = wt./sum(wt);
for itt = 1:extraIters
  wkRule = findWeakRuleSamples(data,labels,wt,binVals,bins,params);

  tr = wkRule.tr;
  dir = wkRule.dir;
  dim = wkRule.dim;
  if dir>0,
    tt = ((data(:,dim)> tr)-0.5)*2;
  else
    tt = ((data(:,dim)<= tr)-0.5)*2;
  end
  curError = sum( (tt.*labels).*wt);
  if(curError<0);
    curError = 0;
  end

  wkRule.error = 0.5-curError/2;
  wkRule.alpha = 1-2*wkRule.error;

  newModel(end+1) = wkRule;
  scores = scores + myBoostClassify(data,newModel(end));
  tt = scores.*labels;
  wt = exWt./(1+exp(tt));
  wt = wt./sum(wt);
end

author: kristinbranson

repo: JAABA


text:
function [data_out] = rsencodecod(data_in,rate_id,TxRx)
%%  Description: The Reed-Solomon encoder is realized according to the standard. %

if rate_id==0
    if TxRx==10
        data_out=[data_in;0;0;0;0;0;0;0;0]; %RS encoder is by passed for BPSK1/2
    elseif TxRx==01
        data_out=[data_in(1:end-8)];     %remove added zero bit
    end
else
switch (rate_id)
    case 1
         n = 32;
         k = 24;
    case 2
        n = 40;
        k = 36;
    case 3
        n = 64;
        k = 48;
    case 4
        n = 80;
        k = 72;
    case 5
        n = 108;
        k = 96;
    case 6
        n = 120;
        k = 108;
 end
m = 8;   % Number of bits per symbol
d=n-k;
if TxRx==10
    %  To realize the Reed-Solomon code, the information is needed in decimal.
    data = reshape(data_in,8,length(data_in)/8);
    data = bi2de(data.','left-msb');
    % one bytes are needed with a stuffed zero at the end of the vector:
      data=[data.' 0];
   % The Galois vector is generated, the generating polynomial of the
   % code. Then the symbols are encoded with Reed-Solomon.
     msg = gf([data],m);   % used polynomial to generate gf field array is same as describe in ieee standrds
     codeRS = rsenc(msg,n,k);  % here we used primitive polynomial as in bydefault in matlab
     out = codeRS.x ;           % convert codeRS a gf object to uni8 array
     data_out=double(out);     % convert to double array
     data_out=[data_out(end-d+1:end) data_out(1:end-d)];  % extra bits should be send befor the orignal bits
   % decimal to binary conversion for continue
   data_out=de2bi(data_out,'left-msb');
   data_out=reshape(data_out.',length(data_out)*8,1);

elseif TxRx==01
    % RS DECODER
    % binary to decimal conversion
    data = reshape(data_in,8,length(data_in)/8);
    data = bi2de(data.','left-msb');
    data=[data.'];
    data=[data(d+1:end) data(1:d)]; % putting the extra bit to last of the array
    %decoding
    msg = gf([data],m);
    decodeRS = rsdec(msg,n,k);
    out = decodeRS(1:end-1)';
    out= out.x;
    out = double (out);
    % The binary data to continue working:
    out = de2bi (out,'left-msb');
    out = reshape (out.', 1, length(out)*8);
   data_out=out';
end
end


text:
%%
X = rand(100,400);
k = 100;

tic; [U,S,V] = svdsecon(X,k); toc; norm(U*S*V' - X,'fro')
tic; [U,S,V] = svds(X,k); toc; norm(U*S*V' - X,'fro')

tic; [U,S,V] = svdecon(X); toc; norm(U*S*V' - X,'fro')
tic; [U,S,V] = svd(X,'econ'); toc; norm(U*S*V' - X,'fro')

%%
%X = rand(400,100);
X = X';
k = 100;

tic; [U,S,V] = svdsecon(X,k); toc; norm(U*S*V' - X,'fro')
tic; [U,S,V] = svds(X,k); toc; norm(U*S*V' - X,'fro')

tic; [U,S,V] = svdecon(X); toc; norm(U*S*V' - X,'fro')
tic; [U,S,V] = svd(X,'econ'); toc; norm(U*S*V' - X,'fro')

%%
%X = rand(400,100);
X = X';
k = 9;

Xm = bsxfun(@minus,X,mean(X,2));
tic; [U,T] = pcaecon(X,k); toc
norm(Xm-U*T,'fro')

tic; [U,T] = pcasecon(X,k); toc
norm(Xm-U*T,'fro')

tic; [U,S,V] = svdsecon(Xm,k); toc
norm(Xm-U*S*V','fro')

tic; [U,S,V] = svdecon(Xm); toc
norm(Xm-U(:,1:k)*S(1:k,1:k)*V(:,1:k)','fro')

tic; [U,S,V] = svds(Xm,k); toc
norm(Xm-U*S*V','fro')

tic; [U,S,V] = svd(Xm); toc
norm(Xm-U(:,1:k)*S(1:k,1:k)*V(:,1:k)','fro')

author: andrewssobral

repo: lrslibrary


text:
function ye = ubd(t, j, k, side)
%UBD	Prescribes u at inflow boundaries
% Possible values for side: 'lower', 'upper',  'left', 'right'

global geval u0 J K yu yv yseglen alpha

ye = 0;
if geval == 1		% Horizontal Poiseuille flow to the right
  if strcmp(side, 'left') == 1
    ye = u0(k*(J+1));	% Inflow profile = outflow profile
  else
    ye = 0;
  end
elseif geval == 2	% Vertical Poiseuille flow
    ye = 0;
elseif geval == 3	% Backward facing step
  if strcmp(side, 'left') == 1
    ye = 6*(yu(k) - yseglen(end)).*(yv(end) - yu(k))/(yv(end) - yseglen(end))^3;
  elseif strcmp(side, 'right') == 1
    ye = 6*(yu(k) - yv(1)).*(yv(end) - yu(k))/(yv(end) - yv(1))^3;
  else
    ye = 0;
  end
elseif geval == 4	% Driven cavity
  if strcmp(side, 'upper') == 1
    ye = 1;
  else
    ye = 0;
  end
elseif geval == 5|geval == 6	% Uniform flow under angle alpha
  ye = cos(alpha);
elseif geval == 7		% Horizontal Poiseuille flow to the left
  if strcmp(side, 'right') == 1
    ye = u0(k*(J+1));	% Inflow profile = outflow profile
  else
    ye = 0;
  end
else
  error('Wrong value in input for parameter geval')
end

repo: Aero-matlab


text:
function c = tapas_hgf_ar1_binary_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF) for AR(1) processes
% for binary inputs in the absence of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading (apart from the AR(1) process) to the
% model introduced in Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_hgf_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% The phis are estimated in 'logit space' because they are confined to the interval from 0 to 1.
% 'Logit-space' is a logistic sigmoid transformation of native space with a variable upper bound
% a>0:
%
% tapas_logit(x) = ln(x/(a-x)); x = a/(1+exp(-tapas_logit(x)))
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_binary_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.phi        row vector of phis (representing reversion slope to attractor; in ascending order of levels)
%         est.p_prc.m        row vector of ms (representing attractors; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and responses y and use
%
%   >> est = tapas_fitModel(y, u, 'tapas_hgf_binary_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'tapas_hgf_ar1_binary';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN, log(0.006), log(4)];
c.logsa_0sa = [NaN,          0,      0];

% Phis
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero (-Inf in logit space) to set to zero.
c.logitphimu = [NaN, -Inf, tapas_logit(0.1,1)];
c.logitphisa = [NaN,    0,                  2];

% ms
% Format: row vector of length n_levels.
% This should be fixed for all levels where the omega of
% the next lowest level is not fixed because that offers
% an alternative parametrization of the same model.
c.mmu = [NaN, c.mu_0mu(2), c.mu_0mu(3)];
c.msa = [NaN,           0,           1];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -2,  -6];
c.omsa = [NaN, 4^2, 4^2];

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.logitphimu,...
    c.mmu,...
    c.logkamu,...
    c.ommu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.logitphisa,...
    c.msa,...
    c.logkasa,...
    c.omsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 5*c.n_levels+(c.n_levels-1);
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_ar1_binary;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_ar1_binary_transp;

return;


text:
function D = spm_eeg_erp_correction(S)
% Applies corrections to ERPs or single trials as in DCM-ERP
% This can be used to make a sensor level analysis or source reconstruction
% consistent with DCM.
% FORMAT D = spm_eeg_erp_correction(S)
%
% S        - optional input struct
% (optional) fields of S:
% S.D        - MEEG object or filename of M/EEG mat-file with epoched data
% S.detrend  - detrending order (0 for no detrending)
% S.hanning  - apply Hanning window (true or false)
%
% Output:
% D        - MEEG object (also written on disk)
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Melanie Boly
% $Id: spm_eeg_erp_correction.m 6907 2016-10-21 09:41:59Z vladimir $

SVNrev = '$Rev: 6907 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','Correct ERPs'); spm('Pointer','Arrow');

if nargin == 0
    S = [];
end

%-Get MEEG object
%--------------------------------------------------------------------------
try
    D = S.D;
catch
    [D, sts] = spm_select(1, 'mat', 'Select M/EEG mat file');
    if ~sts, D = []; return; end
    S.D = D;
end

D = spm_eeg_load(D);

if isequal(D.type, 'continuous')
    error('Corrections can only be applied to epoched or averaged datasets.');
end

%-Get parameters
%--------------------------------------------------------------------------
if ~isfield(S, 'detrend')
    S.detrend = spm_input('Detrend order', '+1', 'n', '2', 1);
end

if ~isfield(S, 'hanning')
    S.hanning = spm_input('Apply Hanning?','+1','yes|no',[1 0], 0);
end

Ns = D.nsamples;

%-Confounds - DCT:
%--------------------------------------------------------------------------
if S.detrend == 0
    X0 = sparse(Ns, 1);
else
    X0 = spm_dctmtx(Ns, S.detrend);
end
R      = speye(Ns) - X0*X0';

%-Hanning
%--------------------------------------------------------------------------
if S.hanning
    R  = R*diag(spm_hanning(Ns))*R;
end

Dnew = clone(D, ['C' fname(D)]);

spm_progress_bar('Init', D.ntrials, 'Trials filtered'); drawnow;
if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials,100));
else Ibar = [1:D.ntrials]; end

spm('Pointer','Watch');

%-Adjust data
%--------------------------------------------------------------------------
for i = 1:D.ntrials
    Dnew(Dnew.indchantype('MEEG', 'GOOD'),:,i) = (R*spm_squeeze(D(D.indchantype('MEEG', 'GOOD'),:,i), 3)')';

    if ismember(i, Ibar), spm_progress_bar('Set', i); end
end

spm_progress_bar('Clear');

D = Dnew;

D = D.history(mfilename, S);

save(D);

%-Cleanup
%--------------------------------------------------------------------------
spm('FigName','Correct ERPs: done'); spm('Pointer','Arrow');


text:
function a = r83t_dif2 ( m, n )

%*****************************************************************************80
%
%% R83T_DIF2 returns the DIF2 matrix in R83T format.
%
%  Example:
%
%    N = 5
%
%    2 -1  .  .  .
%   -1  2 -1  .  .
%    . -1  2 -1  .
%    .  . -1  2 -1
%    .  .  . -1  2
%
%  Properties:
%
%    A is banded, with bandwidth 3.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is a special case of the TRIS or tridiagonal scalar matrix.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is positive definite.
%
%    A is an M matrix.
%
%    A is weakly diagonally dominant, but not strictly diagonally dominant.
%
%    A has an LU factorization A = L * U, without pivoting.
%
%      The matrix L is lower bidiagonal with subdiagonal elements:
%
%        L(I+1,I) = -I/(I+1)
%
%      The matrix U is upper bidiagonal, with diagonal elements
%
%        U(I,I) = (I+1)/I
%
%      and superdiagonal elements which are all -1.
%
%    A has a Cholesky factorization A = L * L', with L lower bidiagonal.
%
%      L(I,I) =    sqrt ( (I+1) / I )
%      L(I,I-1) = -sqrt ( (I-1) / I )
%
%    The eigenvalues are
%
%      LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
%                = 4 SIN^2(I*PI/(2*N+2))
%
%    The corresponding eigenvector X(I) has entries
%
%       X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
%
%    Simple linear systems:
%
%      x = (1,1,1,...,1,1),   A*x=(1,0,0,...,0,1)
%
%      x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
%
%    det ( A ) = N + 1.
%
%    The value of the determinant can be seen by induction,
%    and expanding the determinant across the first row:
%
%      det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
%                = 2 * N - (N-1)
%                = N + 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969,
%    ISBN: 0882756494,
%    LC: QA263.68
%
%    Morris Newman, John Todd,
%    Example A8,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, Number 4, pages 466-476, 1958.
%
%    John Todd,
%    Basic Numerical Mathematics,
%    Volume 2: Numerical Algebra,
%    Birkhauser, 1980,
%    ISBN: 0817608117,
%    LC: QA297.T58.
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Input, integer M, N, the order of the matrix.
%
%    Output, real A(M,3), the matrix.
%
  a = zeros(m,3);

  mn = min ( m, n );

  a(2:mn,1)   = -1.0;
  a(1:mn,2)   = +2.0;
  a(1:mn-1,3) = -1.0;

  if ( m < n )
    a(mn,3) = -1.0;
  elseif ( n < m )
    a(mn+1,1) = -1.0;
  end

  return
end


text:
function pdf = angle_pdf ( x, n )

%*****************************************************************************80
%
%% ANGLE_PDF evaluates the Angle PDF.
%
%  Discussion:
%
%    X is an angle between 0 and PI, corresponding to the angle
%    made in an N-dimensional space, between a fixed line passing
%    through the origin, and an arbitrary line that also passes
%    through the origin, which is specified by a choosing any point
%    on the N-dimensional sphere with uniform probability.
%
%  Formula:
%
%    PDF(X) = ( sin ( X ) )**(N-2) * Gamma ( N / 2 )
%             / ( sqrt ( PI ) * Gamma ( ( N - 1 ) / 2 ) )
%
%    PDF(X) = 1 / PI if N = 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Reuven Rubinstein,
%    Monte Carlo Optimization, Simulation and Sensitivity of Queueing Networks,
%    Wiley, 1986.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, integer N, the spatial dimension.
%    N must be at least 2.
%
%    Output, real PDF, the value of the PDF.
%
  if ( n < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '\ANGLE_PDF - Fatal error!\n' );
    fprintf ( 1, '\  N must be at least 2.\n' );
    fprintf ( 1, '\  The input value of N = %d\n', n );
    error ( 'ANGLE_PDF - Fatal error!' );
  end

  if ( x < 0.0 | pi < x )
    pdf = 0.0;
  elseif ( n == 2 )
    pdf = 1.0 / pi;
  else
    pdf = ( sin ( x ) )^( n - 2 ) * gamma ( n / 2.0 ) ...
      / ( sqrt ( pi ) * gamma ( ( n - 1 ) / 2.0 ) );
  end

  return
end


text:
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http ://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Slice (apply) bilateral guided upsampling.
%
% Inputs:
%
% gamma: the affine bilateral grid.
%   size: [height width depth num_output_channels num_input_channels]
% input_image is a double tensor: the (full-resolution) input image
% edge_image is a double matrix: the (full-resolution) edge image
%
function output = bguSlice(gamma, input_image, edge_image)

% Find downsampling coordinates, without rounding.
input_height = size(input_image, 1);
input_width = size(input_image, 2);
grid_height = size(gamma, 1);
grid_width = size(gamma, 2);
grid_depth = size(gamma, 3);
affine_output_size = size(gamma, 4);
affine_input_size = size(gamma, 5);

% meshgrid inputs and outputs are x, then y, with x right, y down.
[ x, y ] = meshgrid(0:(input_width - 1), 0:(input_height - 1));

% Downsample x and y to grid space (leaving them as floats).
bg_coord_x = ((x + 0.5) * (grid_width - 1) / input_width);
bg_coord_y = ((y + 0.5) * (grid_height - 1) / input_height);
bg_coord_z = edge_image * (grid_depth - 1);

% Add 1 to all coordinates for MATLAB.
bg_coord_xx = bg_coord_x + 1;
bg_coord_yy = bg_coord_y + 1;
bg_coord_zz = bg_coord_z + 1;

% interp3 takes xx, yy, zz.
affine_model = {affine_output_size, affine_input_size};
for j = 1:affine_input_size
    for i = 1:affine_output_size
            affine_model{i,j} = interp3(gamma(:,:,:,i,j), ...
                bg_coord_xx, bg_coord_yy, bg_coord_zz);
    end
end

% Iterate over each row.
% TODO: optimize this.
for i = 1:affine_output_size
    affine_model2{i,1} = cat(4, affine_model{i,:});
end
affine_model3 = cat(3, affine_model2{:,1});

input1 = cat(3, input_image, ones(size(input_image,1), size(input_image,2)));

output = 0;
for i = 1:affine_input_size
    output = output + bsxfun(@times, affine_model3(:,:,:,i), input1(:,:,i));
end

author: mahmoudnafifi


text:
function [IDX,C] = meanShiftPost( X, IDX, C, minCl, forceOutl )
% Some post processing routines for meanShift not currently being used.
%
% USAGE
%  [IDX,C] = meanShiftPost( X, IDX, C, minCl, forceOutl )
%
% INPUTS
%  see meanShift
%
% OUTPUTS
%  see meanShift
%
% EXAMPLE
%
% See also MEANSHIFT, KMEANS2
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

%%% force outliers to belong to IDX (mainly for visualization)
if( forceOutl )
  for i=find(IDX<0)'
    D = pdist2( X(i,:), C );
    [mind IDx] = min(D,[],2);
    IDX(i) = IDx;
  end
end;

%%% Delete smallest cluster, reassign points, re-sort, repeat...
k = max(IDX);
ticId = ticStatus('meanShiftPost',[],5); kinit = k;
while( 1 )
  % sort clusters [largest first]
  cnts = zeros(1,k); for i=1:k; cnts(i) = sum( IDX==i ); end
  [cnts,order] = sort( -cnts ); cnts = -cnts; C = C(order,:);
  IDX2 = IDX;  for i=1:k; IDX2(IDX==order(i))=i; end; IDX = IDX2;

  % stop if smallest cluster is big enough
  if( cnts(k)>= minCl ); break; end;

  % otherwise discard smallest [last] cluster
  C( end, : ) = [];
  for i=find(IDX==k)'
    D = pdist2( X(i,:), C );
    [mind IDx] = min(D,[],2);
    IDX(i) = IDx;
  end;
  k = k-1;

  tocStatus( ticId, (kinit-k)/kinit );
end
tocStatus( ticId, 1 );

author: zouchuhang

repo: LayoutNet


text:
function planC = save_planC(planC,optS,saveflag,optFileName)
%"save_planC"
%   Saves plan using GUI interface.  Can either behave as Save or SaveAs,
%   where Save prompts to overwrite source planC file and SaveAs raises a
%   file location menu.  saveflag is either 'save' or 'saveas'.  If save is
%   called, the filename is assumed to be contained in stateS.CERRFile.
%
%   If 'save' is called, .bz2 files are automatically compressed and
%   replaced if they existed (assuming the user agrees when prompted).
%
%   If 'saveas' is called and a .bz2 file is specified, it is also
%   automatically compressed.  If a .mat file is specified, the user is
%   asked if compression is desired.
%
%     A.I. Blanco.
%LM:  14 Oct 02, JOD.
%     26 Jun 04, JRA, Added saveas/save distinction.
%     20 Jul 04, JRA, If saveas is used, now updates current file state.
%     20 Aug 04, JRA, Automatic compression if .bz2 file specified.
%     10 Sep 04, JRA, Now check for Matlab7 and saves according to optS.
%     27 Apr 05, JRA, Added tarball handling.
%     26 Jul 06, APA, Removed tarball handling since remote variables are
%     now stored in a subdirectory ...\.._store.
%     09 Jul 07, APA, Used Matlab's built-in tar function to tar CERR plan
%     and the remotely stored files.
%     17 Dec 07, APA, Added option to compress CERR plan using zip utility.
%     11 Jan 08, APA, (incorporated KU's changes) Fixed case where file could
%                     be overwritten without a warning if filename is entered
%                     without an extension. In case of error during save attempt
%                     (eg., permission denied), added ability to re-try.
%
%Usage:
%   zipped = save_planC(planC,optS, saveflag)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS

indexS = planC{end};

%Determine if we want 'save', or 'saveas...' behavior.
if ~exist('saveflag','var')
    saveflag = 'saveas';
elseif ~any(strcmpi(saveflag, {'save','saveas','passed'}))
    error('Incorrect call to save_planC.m');
end

%Store the revision and timeSaved into header of planC
planC = updateRevision(planC);

permission = 0;
while permission == 0   %Use while statement in case permission to overwrite file is denied.   KU 12/15/07
    %Acquire save filename according to saveFlag type.
    switch upper(saveflag)
        case 'SAVE'
            button = questdlg(['Are you sure you wish to save changes made to ' stateS.CERRFile ' ?'], 'Overwrite', 'Yes', 'No', 'No');
            if ~strcmpi(button, 'yes')
                CERRStatusString(['Save cancelled. Ready.']);
                return;
            end
            saveFile = stateS.CERRFile;

        case 'SAVEAS'

            deciding = 1;
            while deciding    %Use while statement in case filename is entered without an extension.   KU 12/8/07

                whichPlan = fliplr(strtok(fliplr(fileparts(char(planC{indexS.header}.archive))),'\/'));

                try
                    directory = fileparts(stateS.CERRFile);
                catch
                    directory = cd;
                end
                if isempty(directory)
                    directory = cd;
                end
                wd = cd;
                cd(directory);
                fileTypeC = {'*.mat', 'Native Matlab/Octave (*.mat)';
                    '*.mat.bz2, *.mat.zip)','Compressed (*.mat.zip, *.mat.zip)'};
                [fname,pname] = uiputfile(fileTypeC,['Save the ' whichPlan ' data as:']);
                cd(wd);
                if isequal(fname,0) || isequal(pname,0)
                    CERRStatusString('Save cancelled. Ready.');
                    return;
                end
                if ~isempty(strfind(fname,'.'))
                    saveFile = fullfile(pname,fname);
                    deciding = 0;
                else     %in case filename is entered without an extension.   KU 12/8/07
                    fname = [fname,'.mat'];
                    list = dir(fullfile(pname, '*.*'));
                    match = 0;
                    for i = 1:size(list,1)
                        if list(i).isdir==0 && strcmpi(fname, list(i).name)
                            match = 1;
                            button = questdlg([[pname fname],' already exists! Do you want to replace it?'],'Save file as?','yes','no','no');
                            if strcmpi(button, 'yes')
                                saveFile = fullfile(pname,fname);
                                deciding = 0;
                                break
                            end
                        end
                    end

                    if match == 0
                        saveFile = fullfile(pname,fname);
                        deciding = 0;
                    end
                end
            end

            %get storage directory for remote variables in the original plan
            remotePath = '';
            if isfield(stateS,'CERRFile') && ~isempty(stateS.CERRFile)
                [pathstr, name, ext] = fileparts(stateS.CERRFile);
                remotePath = fullfile(pathstr,[name,'_store']);
            end

        case 'PASSED'
            saveFile = optFileName;
    end

    %Add .mat if required, parse other valid extensions if they exist.  If
    %invalid extensions exist, they are considered part of the root filename
    %and proper extensions are appended to them.
    [ext, roots] = getFileExtensions(saveFile);
    if ~isempty(ext)
        if strcmpi(ext(1), '.mat')
            zipFile     = 0;
            saveFile    = [roots{1} '.mat'];
        elseif length(ext) > 2 && strcmpi(ext(1), '.tar') && ...
                (strcmpi(ext(2), '.bz2') || strcmpi(ext(2), '.zip')) && strcmpi(ext(3), '.mat')
            zipFile     = 1;
            saveFile    = [roots{3} '.mat'];
        elseif length(ext) > 1 && (strcmpi(ext(1), '.bz2') || ...
                strcmpi(ext(1), '.zip')) && strcmpi(ext(2), '.mat')
            zipFile     = 1;
            saveFile    = [roots{2} '.mat'];
        else
            saveFile    = [saveFile '.mat'];
            zipFile     = 0;
        end
    else
        saveFile    = [saveFile '.mat'];
        zipFile     = 0;
    end

    %Zip if required, prompting if saveAs w/o .bz2 suffix specified.
    if zipFile
        ans = 'yes';
        %[jnk1,jnk2,ext] = fileparts(stateS.CERRFile);
        %if ~strcmpi(ext,'bz2') || ~strcmpi(ext,'zip')
        extSave = ext;
        indMat = strcmpi('.mat',extSave);
        extSave(indMat) = [];
        extSave = extSave{1};
        if ~strcmpi(extSave,'bz2') || ~strcmpi(extSave,'zip')
            if ~isfield(stateS,'optS')
                stateS.optS = opts4Exe([getCERRPath,'CERROptions.json']);
            end
            extSave = stateS.optS.CompressType;
        end
    elseif ~zipFile && strcmpi(saveflag, 'passed')
        ans = 'no';
    elseif ~zipFile && strcmpi(saveflag, 'saveas')
        ans = questdlg('Zip the .mat file using bz2/zip?');
        if ~isfield(stateS,'optS')
            stateS.optS = opts4Exe([getCERRPath,'CERROptions.json']);
        end
        extSave = stateS.optS.CompressType;
    else
        ans = 'no';
    end

    if strcmpi(ans,'yes')
        zipFile = 1;
    end
    %Save the planC to the specified file.
    try
        [saveFile, planC] = writeToDisk(planC, saveFile, zipFile);
        permission = 1;
    catch
        h = warndlg(lasterr,'Error Saving File!','modal');
        uiwait(h);
    end
end

%Remind CERR where the current file is, in case of filename change.
stateS.CERRFile = saveFile;

[pathstr, name, ext] = fileparts(stateS.CERRFile);
switch lower(ans)
    case 'yes'
        if strcmpi(extSave,'bz2')
            CERRStatusString(['Compressing ' saveFile '...']);
            outstr = gnuCERRCompression(saveFile,'compress');
            if (~ischar(outstr)) % Function returned as zero when compression complete
                errordlg('File not found by compression routine.', 'CERR compression');
                return;
            else
                stateS.CERRFile = [saveFile '.bz2'];
            end
        elseif strcmpi(extSave,'zip')
            zip(saveFile,saveFile)
            delete(saveFile)
            stateS.CERRFile = [saveFile '.zip'];
        end
        if exist(fullfile(pathstr,[name,'_store']),'dir')==7
            copyfile(fullfile(pathstr,[name,'_store']),fullfile(pathstr,[name,'.mat_store']),'f')
            rmdir(fullfile(pathstr,[name,'_store']),'s')
            [planC, stateS] = updateRemotePaths(planC, stateS, stateS.CERRFile);
        end
end

%tar the files if required
if isfield(stateS,'reqdRemoteFiles') && ~isempty(stateS.reqdRemoteFiles)
    filesToTar = stateS.reqdRemoteFiles;
    if ~isempty(filesToTar)
        dirToTar = fileparts(filesToTar{1});
    end
    %filesToTar{end+1} = stateS.CERRFile;
    tar([stateS.CERRFile,'.tar'],{dirToTar stateS.CERRFile})
    %delete the .mat or .bz2 file since it is in tarred format now
    delete(stateS.CERRFile)
else
    %delete the tar file since no remote files exist
    if exist([stateS.CERRFile,'.tar'],'file')
        %try
            delete([stateS.CERRFile,'.tar'])
        %end
    end
end

if strcmpi(saveflag,'SAVEAS')
    if exist(remotePath,'dir') == 7
        rmdir(remotePath,'s')
    end
end

CERRStatusString(['Saved ' name ext '. Ready. (' datestr(now) ')']);

try
    stateS.planMerged = 0;
end

function [saveFile, planC] = writeToDisk(planC, saveFile, zipFile)
%"writeToDisk"
%   Save the planC to disk.  If remote variables exist, handle them.

global stateS
if iscell(planC)
    indexS = planC{end};
end

%Get a list of remote variable files using the save flag to store.
% remoteFiles = listRemoteFiles(planC, 0);
remoteFiles = listRemoteScanAndDose(planC);

%Prepare a list of files that will be added to _store directory.
filesLocal = [];
remoteFullFile = {};
remotePathLocal = {};
for i = 1:length(remoteFiles)
    switch upper(remoteFiles(i).storageType)
        case {'LOCAL'}
            remotePathLocal{i} = remoteFiles(i).remotePath;
            filesLocal{i} = remoteFiles(i).filename;
            remoteFullFile{i} = fullfile(remotePathLocal{i},filesLocal{i});
    end
end

[pathstr, name, ext] = fileparts(saveFile);
CERRStatusString(['Saving ' name ext '...']);

% save Remotely stored files

if zipFile
    name = [name '.mat'];
end
if ~isempty(filesLocal) && ~exist(fullfile(pathstr,[name,'_store']),'dir')
    mkdir(fullfile(pathstr,[name,'_store']))
elseif isempty(filesLocal)
    if length(remotePathLocal) > 1
        uniqRemothPathC = unique(remotePathLocal);
        for iPath = 1:length(uniqRemothPathC)
            if exist(uniqRemothPathC{iPath},'dir') == 7
                rmdir(uniqRemothPathC{iPath},'s')
            end
        end
    end
end

for i = 1:length(filesLocal)
    if ~isequal(remoteFullFile{i},fullfile(pathstr,[name,'_store'],filesLocal{i}))
        copyfile(remoteFullFile{i},fullfile(pathstr,[name,'_store'],filesLocal{i}))
    end
end

% delete unnecessary remote files in original plan that were created
% temporarily
for i = 1:length(remoteFullFile)
    if isfield(stateS,'reqdRemoteFiles') && ~ismember(remoteFullFile{i},stateS.reqdRemoteFiles) && ...
            ~isequal(fileparts(remoteFullFile{i}),fullfile(pathstr,[name,'_store']))
        delete(remoteFullFile{i})
    end
end

%delete unnecessary remote files in original plan that changed
if ~isempty(stateS) && isfield(stateS,'reqdRemoteFiles')
    for i = 1:length(stateS.reqdRemoteFiles)
        if ~ismember(stateS.reqdRemoteFiles(i),remoteFullFile) && ...
                isequal(fileparts(stateS.reqdRemoteFiles{i}),fullfile(pathstr,[name,'_store']))
            delete(stateS.reqdRemoteFiles{i})
        end
    end
end

[planC, stateS] = updateRemotePaths(planC, stateS, saveFile, zipFile);

%Save functions... modified to work with matlab 7
varInfoS = whos('planC');
if varInfoS.bytes > 2*10^9 %Use save version 7.3 for plans > 2GB
    saveOpt = '-v7.3';
else
saveOpt = getSaveInfo;     %Default save version (defined in CERROptions.m)
end
if ~isempty(saveOpt)
    save(saveFile, 'planC', saveOpt);
else
    save(saveFile, 'planC');
end

function [ext, root] = getFileExtensions(filename)
%"getFileExtensions"
%   Returns the file extensions of filename, each in an individual cell, in
%   reverse order.  Also returns the root filename strings for each
%   iteration of the while loop.

ext  = {};
root = {};

[pathstr,fname,fext] = fileparts(filename);
while ~isempty(fext)
    ext{end+1} = fext;
    root{end+1} = fullfile(pathstr, fname);
    [pathstr,fname,fext] = fileparts(fullfile(pathstr, fname));
end

function [planC, stateS] = updateRemotePaths(planC, stateS, saveFile, zipFile)
%"updateRemotePaths"
% update remote path in planC for remote doseArrays and scanArrays. Also
% update fullFiles in stateS

if iscell(planC)
    indexS = planC{end};
else
    return;
end

if ~exist('zipFile','var')
    zipFile = 0;
end
[pathstr, name, ext] = fileparts(saveFile);
if zipFile
    name = [name '.mat'];
end

stateS.reqdRemoteFiles = {};
for i = 1:length(planC{indexS.scan})
    if ~isLocal(planC{indexS.scan}(i).scanArray)
        planC{indexS.scan}(i).scanArray.remotePath = fullfile(pathstr,[name,'_store']);
        stateS.reqdRemoteFiles{end+1} = fullfile(pathstr,[name,'_store'],planC{indexS.scan}(i).scanArray.filename);
        planC{indexS.scan}(i).scanArraySuperior.remotePath = fullfile(pathstr,[name,'_store']);
        stateS.reqdRemoteFiles{end+1} = fullfile(pathstr,[name,'_store'],planC{indexS.scan}(i).scanArraySuperior.filename);
        planC{indexS.scan}(i).scanArrayInferior.remotePath = fullfile(pathstr,[name,'_store']);
        stateS.reqdRemoteFiles{end+1} = fullfile(pathstr,[name,'_store'],planC{indexS.scan}(i).scanArrayInferior.filename);
    end
end
for i = 1:length(planC{indexS.dose})
    if ~isLocal(planC{indexS.dose}(i).doseArray)
        planC{indexS.dose}(i).doseArray.remotePath = fullfile(pathstr,[name,'_store']);
        stateS.reqdRemoteFiles{end+1} = fullfile(pathstr,[name,'_store'],planC{indexS.dose}(i).doseArray.filename);
    end
end

for i = 1:length(planC{indexS.IM})
    if ~isempty(planC{indexS.IM}(i).IMDosimetry.beams) && ~isLocal(planC{indexS.IM}(i).IMDosimetry.beams(1).beamlets)
        for iBeam = 1:length(planC{indexS.IM}(i).IMDosimetry.beams)
            planC{indexS.IM}(i).IMDosimetry.beams(iBeam).beamlets.remotePath = fullfile(pathstr,[name,'_store']);
            stateS.reqdRemoteFiles{end+1} = fullfile(pathstr,[name,'_store'],planC{indexS.IM}(i).IMDosimetry.beams(iBeam).beamlets.filename);
        end
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef Data < handle
    properties (Access = public)
        type
    end

    methods
        function obj = Data(type)
            if nargin == 1
                obj.type = type;
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
%DEMO_PASSGP Demonstration of PASS-GP routine for GP classification
%
%  Description
%    Here we demonstrate pass-gp routine for Gaussian Processes
%    classification. Data used is 2-dimensional toy data with Gaussian
%    bumbs defining classes. We demonstrate with both fixed and not fixed
%    sizes of active set for pass-gp.
%
%  See also PASSGP

% Copyright (c) 2013 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Generate toy data
prevstream=setrandstream(0);
[x1,x2]=meshgrid(-5:0.1:5,-5:0.1:5);
x=[x1(:) x2(:)]; x=x(randperm(size(x,1),3000),:);
y=2.*mnorm_pdf(x, [0 0], [0.5 0;0 0.5]) + mnorm_pdf(x, [3 3], [0.5 0;0 0.5]) + mnorm_pdf(x, [-3 -3], [0.5 0;0 0.5]);
y=y+mnorm_pdf(x, [3 -3], [0.5 0;0 0.5])+mnorm_pdf(x, [-3 3], [0.5 0;0 0.5]);
y=y+0.03.*randn(size(y,1),1);
y(y>0.15)=1; y(y<=0.15)=-1;

[xt1, xt2]=meshgrid(-5:0.23:5,-5:0.23:5);
xt=[xt1(:) xt2(:)];
yt=ones(size(xt,1),1);

[n, nin] = size(x);

% Define covariance and likelihood functions and create model
gpcf = gpcf_sexp();
lik=lik_probit();
gp=gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6);

opt=optimset('TolX',1e-3,'TolFun',1e-3,'display','on');
w0=gp_pak(gp);

% fPASS-GP with fixed size of 800 points in active set and data divided to
% 10 subsets with 4 sweeps over data
start=tic;[gp, indA]=passgp(gp, x, y, 'opt', opt, 'npass', 4, 'ninit', 800, 'nsub', 10, 'display', 'on', 'fixed', 'on', 'pexc', 0.1, 'optimn', 2);time=toc(start);
tt=time;
[Eft, Varft, lpyt, Eyt, Varyt]=gp_pred(gp, x(indA,:), y(indA,:), xt, 'yt', yt);
figure, [cc,hh]=contour(reshape(xt(:,1),size(xt1,1), size(xt1,1)), reshape(xt(:,2),size(xt1,1), size(xt1,1)), reshape(exp(lpyt),size(xt1,1), size(xt1,1)), [0.1 0.9]);
clabel(cc,hh); title('Pr(y==1) (fpass-gp)')
param=gp_pak(gp);

% PASS-GP with inclusion threshold 0.65, deletion threshold 0.99, intial
% size of 400 points in active set and 3 sweeps over data.
gp=gp_unpak(gp,w0);
start=tic;[gp, indA2]=passgp(gp, x, y, 'opt', opt, 'pinc', 0.65, 'pdel', 0.99, 'npass', 3, 'ninit', 400, 'nsub', 10, 'display', 'on', 'optimn', 2);time=toc(start);
tt2=time;
[Eft2, Varft2, lpyt2, Eyt2, Varyt2]=gp_pred(gp, x(indA2,:), y(indA2,:), xt, 'yt', yt);
figure, [cc,hh]=contour(reshape(xt(:,1),size(xt1,1), size(xt1,1)), reshape(xt(:,2),size(xt1,1), size(xt1,1)), reshape(exp(lpyt2),size(xt1,1), size(xt1,1)), [0.1 0.9]);
clabel(cc,hh); title('Pr(y==1) (pass-gp)')
param2=gp_pak(gp);

% Full Gaussian Process for comparison
gp=gp_unpak(gp,w0);
opt.Display='iter';
start=tic;gp=gp_optim(gp,x,y,'opt',opt);tt3=toc;
[Eft3, Varft3, lpyt3, Eyt3, Varyt3]=gp_pred(gp, x, y, xt, 'yt', yt);
figure, [cc,hh]=contour(reshape(xt(:,1),size(xt1,1), size(xt1,1)), reshape(xt(:,2),size(xt1,1), size(xt1,1)), reshape(exp(lpyt3),size(xt1,1), size(xt1,1)), [0.1 0.9]);
clabel(cc,hh); title('Pr(y==1) (full gp)')

% Display some statistics

mlpd_fpassgp=mean(mean(lpyt,2))
time_fpassgp=mean(tt)
mlpd_passgp=mean(mean(lpyt2,2))
time_passgp=mean(tt2)
mlpd_full=mean(lpyt3)
time_full=tt3

% Plot data and active sets for both methods
figure(4), subplot(1,2,1),  plot(x(y==1,1),x(y==1,2),'or',x(y==-1,1),x(y==-1,2),'ob');
hold all; plot(x(indA,1), x(indA,2), '.k'); title('Data and active set (fpass-gp)')
legend('y=1', 'y=-1', 'Active set for fpass-gp');
subplot(1,2,2),  plot(x(y==1,1),x(y==1,2),'or',x(y==-1,1),x(y==-1,2),'ob');
hold all; plot(x(indA2,1), x(indA2,2), '.k'); title('Data and active set (pass-gp)')
legend('y=1', 'y=-1', 'Active set for pass-gp');
setrandstream(prevstream);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Machine Learning Online Class - Exercise 1: Linear Regression

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%
% x refers to the population size in 10,000s
% y refers to the profit in $10,000s
%

%% Initialization
clear ; close all; clc

%% ==================== Part 1: Basic Function ====================
% Complete warmUpExercise.m
fprintf('Running warmUpExercise ... \n');
fprintf('5x5 Identity Matrix: \n');
warmUpExercise()

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

%% ======================= Part 2: Plotting =======================
fprintf('Plotting Data ...\n')
data = load('ex1data1.txt');
X = data(:, 1); y = data(:, 2);
m = length(y); % number of training examples

% Plot Data
% Note: You have to complete the code in plotData.m
plotData(X, y);

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

%% =================== Part 3: Gradient descent ===================
fprintf('Running Gradient Descent ...\n')

X = [ones(m, 1), data(:,1)]; % Add a column of ones to x
theta = zeros(2, 1); % initialize fitting parameters

% Some gradient descent settings
iterations = 1500;
alpha = 0.01;

% compute and display initial cost
computeCost(X, y, theta)

% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent: ');
fprintf('%f %f \n', theta(1), theta(2));

% Plot the linear fit
hold on; % keep previous plot visible
plot(X(:,2), X*theta, '-')
legend('Training data', 'Linear regression')
hold off % don't overlay any more plots on this figure

% Predict values for population sizes of 35,000 and 70,000
predict1 = [1, 3.5] *theta;
fprintf('For population = 35,000, we predict a profit of %f\n',...
    predict1*10000);
predict2 = [1, 7] * theta;
fprintf('For population = 70,000, we predict a profit of %f\n',...
    predict2*10000);

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

%% ============= Part 4: Visualizing J(theta_0, theta_1) =============
fprintf('Visualizing J(theta_0, theta_1) ...\n')

% Grid over which we will calculate J
theta0_vals = linspace(-10, 10, 100);
theta1_vals = linspace(-1, 4, 100);

% initialize J_vals to a matrix of 0's
J_vals = zeros(length(theta0_vals), length(theta1_vals));

% Fill out J_vals
for i = 1:length(theta0_vals)
    for j = 1:length(theta1_vals)
	  t = [theta0_vals(i); theta1_vals(j)];
	  J_vals(i,j) = computeCost(X, y, t);
    end
end

% Because of the way meshgrids work in the surf command, we need to
% transpose J_vals before calling surf, or else the axes will be flipped
J_vals = J_vals';
% Surface plot
figure;
surf(theta0_vals, theta1_vals, J_vals)
xlabel('\theta_0'); ylabel('\theta_1');

% Contour plot
figure;
% Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100
contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))
xlabel('\theta_0'); ylabel('\theta_1');
hold on;
plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);

author: AvaisP


text:
function done = mission1(occupiers, player)

captured = zeros(1,6);
for i = 1:6
    captured(i) = all(occupiers{i} == player);
end

done = captured(3) && captured(6) && sum(captured) > 2;