text:
function [EDF]=closedf(EDF)

% EDF=closedf(EDF)
% Opens an EDF File (European Data Format for Biosignals) into MATLAB
% About EDF
%
% EDF   struct of EDF-Header of a EDF-File

%   Version 2.0
%   15.12.1997
%   Copyright (c) 1997 by  Alois Schloegl
%   a.schloegl@ieee.org
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
%

EDF.FILE.OPEN=0;
fclose(EDF.FILE.FID);
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x = l_polynomial_zeros ( n )

%*****************************************************************************80
%
%% L_POLYNOMIAL_ZEROS: zeros of the Laguerre polynomial L(n,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2012
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the order of the polynomial.
%
%    Output, real X(N), the zeros.
%

%
%  Define the zero-th moment.
%
  zemu = 1.0;
%
%  Define the Jacobi matrix.
%
  b = zeros ( n, 1 );
  for i = 1 : n
    bj(i) = i;
  end

  x = zeros ( n, 1 );
  for i = 1 : n
    x(i) = 2 * i - 1;
  end

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ model ] = initFBAFields( model, printLevel )
% This function initializes all fields in a model that are required for
% downstream FBA analysis. It does so if and only if a Stoichiometric
% matrix S is provided.
%
% USAGE:
%    model = convertOldStyleModel(model, printLevel)
%
% INPUT:
%    model:         a COBRA Model structure with at least the model.S field.
%                   All fields already present are retained and absent
%                   fields are initialized with their defaults.
%
% OPTIONAL INPUT:
%    printLevel:    indicates whether warnings and messages are given (default, 1).
%
% OUTPUT:
%    model:         a COBRA model struct with the following fields:
%                   .S (same as input)
%                   .rxns (default: a vector of strings R1 .. R size(S,2)
%                   .mets (default: a vector of strings M1 .. M size(S,1)
%                   .lb   (default: -1000 * ones(size(S,2),1) );
%                   .ub   (default: 1000 * ones(size(S,2),1) );
%                   .genes   (default: cell(0,1));
%                   .rules   (default: repmat({''},size(S,2),1))
%                   .osense (default: -1)
%                   .csense (default: a char vector of 'E' of the size size(S,1) x 1)
%
% .. Author: - Thomas Pfau Nov 2017

if ~isfield(model,'S')
    error('Cannot create FBA fields without stoichiometric matrix');
else
    %Init basic sizes
    rxns = size(model.S,2);
    mets = size(model.S,1);
    genes = 0;
end

structfields = fieldnames(model);
optionalFields = getDefinedFieldProperties();
fluxConsistencyFields = optionalFields(cellfun(@(x) x, optionalFields(:,6)),:);

for field = 1:size(fluxConsistencyFields,1)
    cfieldName = fluxConsistencyFields{field,1};
    if ~any(ismember(structfields,cfieldName))

        csize1 = fluxConsistencyFields{field,2};
        if ischar(csize1)
            eval(['csize1 = ' fluxConsistencyFields{field,2} ';']);
        end
        csize2 = fluxConsistencyFields{field,3};
        if ischar(csize2)
            eval(['csize2 = ' fluxConsistencyFields{field,3} ';']);
        end

        i = 0;
        defaultVal = fluxConsistencyFields{field,5};
        isCellEntry = false;
        if ischar(defaultVal)
            if strfind('iscell',fluxConsistencyFields(field,4))
                eval(['defaultVal = ' fluxConsistencyFields{field,5} ';']);
            else
                defaultVal = fluxConsistencyFields{field,5};
            end
            isCellEntry = true;
        end
        if ~isCellEntry
            cfield = zeros(csize1,csize2);
            for i = 1:csize1
                cfield(i) = fluxConsistencyFields{field,5};
            end
            model.(cfieldName) = cfield;
        else
            %Lets also check for this being a char.
            if isequal(fluxConsistencyFields{field,4},'ischar(x)')
                cfield = repmat(defaultVal,csize1,csize2);
                model.(cfieldName) = cfield;
            else

                cfield = cell(csize1,csize2);
                for i = 1:csize1
                    eval(['cfield{i} = ' fluxConsistencyFields{field,5} ';']);
                end
                model.(cfieldName) = cfield;
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function r8mat_house_form_test ( )

%*****************************************************************************80
%
%% R8MAT_HOUSE_FORM_TEST tests R8MAT_HOUSE_FORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  v = [ 0.0, 0.0, 1.0, 2.0, 3.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_HOUSE_FORM_TEST\n' );
  fprintf ( 1, '  R8MAT_HOUSE_FORM forms a Householder\n' );
  fprintf ( 1, '  matrix from its compact form.\n' );

  r8vec_print ( n, v, '  Compact vector form V:' );

  h = r8mat_house_form ( n, v );

  r8mat_print ( n, n, h, '  Householder matrix H:' );

  return
end


text:
function kernel = SC_ASGkernel(sccalib, geo, dus, dvs)
%% Anti-sctter grid response function
% Reference: Improved scatter correction using adaptive scatter kernel superposition
% Input:
%           geo: geometry structure
%           dus: downsampled u vector
%           dvs: downsampled v vector
% Output:
%           kernel: anti-scatter grid
%
% Date: 2021-05-04
% Author: Yi Du (yi.du@hotmail.com)

%% Geometry
% unit mm
DSD = geo.DSD;

%% Anti-Scatter Grid along X direction
% vs dimension
gamma = abs(rad2deg(atan(dvs'./DSD)));

%% Transmission modelling
k = -0.15;
b = 1;

t_ratio = k.*abs(dvs'/10) + b;

%% Kernel: [nv, nu]
kernel = repmat(t_ratio, [1,length(dus)]);
efficiency = str2double(sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{1}.GridEfficiency.LamellaTransmission.Text);
kernel(kernel < efficiency) = efficiency;

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function A = mycol2im(v,par1, par2)
% FUNCTION K = mycol2im(v,par1,par2)
%   AUTHOR:    Makoto Yamada
%              (myamada@ism.ac.jp)
%   DATE:       02/16/08
%
%  DESCRIPTION:
%
%   This function is for converting 1D signal to 2D signal.
%

count = 1;
A =zeros(par2(1), par2(2));

for ii = 1:par2(1)
    for jj = 1:par2(2)
        A(jj,ii) = v(count);
        count = count + 1;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
% POP_STUDYERP - create a simple design for ERP analysis
%
% Usage:
%       >> [STUDY ALLEEG] = pop_studyerp; % pop up interface
%
% Outputs:
%   STUDY        - an EEGLAB STUDY set of loaded EEG structures
%   ALLEEG       - ALLEEG vector of one or more loaded EEG dataset structures
%
% Author: Arnaud Delorme, SCCN, UCSD, 2011-
%
% See also: EEG_CHECKSET

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

function [STUDY ALLEEG com ] = pop_studyerp;

% first GUI, get the number of conditions and subjects
% ----------------------------------------------------
textinfo = [ 'This interface creates a simple STUDY and ' 10 ...
             'computes its condition grand average ERPs.' 10 ...
             'For each subject, trials for each condition' 10 ...
             'must first be stored in a separate dataset.' 10 ...
             'Create other STUDY using the standard editor.' ];
guispec = { ...
    {'style' 'text' 'string' 'Create simple ERP STUDY' ...
     'FontWeight' 'Bold' 'fontsize', 12} ...
    { 'style' 'text' 'string' textinfo } ...
    {'style' 'text' 'string' 'Number of conditions:' } ...
    {'style' 'edit' 'string' '1' 'tag' 'cond' } { } ...
    {'style' 'text' 'string' 'Number of subjects:' } ...
    {'style' 'edit' 'string' '15' 'tag' 'subjects' } { } };
guigeom = { [1] [1] [1 0.3 0.4] [1 0.3 0.4] };

optiongui = { 'geometry', guigeom, 'uilist'  , guispec, ...
              'geomvert', [ 1 4 1 1], ...
              'helpcom' , 'pophelp(''pop_studyerp'')', ...
              'title'   , 'Create a new STUDY set -- pop_studyerp()' };
[result, userdat2, strhalt, outstruct] = inputgui(optiongui{:});
STUDY  = [];
ALLEEG = [];
com = '';
if isempty(result), return; end

nSubjects = str2num(outstruct.subjects);
nConds    = str2num(outstruct.cond);

% second GUI, enter the datasets
% ------------------------------
guispec = { ...
    {'style' 'text' 'string' 'Create simple ERP STUDY' 'FontWeight' 'Bold' 'fontsize', 12} ...
    {} ...
    {} {'style' 'text' 'string' 'STUDY set name:' } { 'style' 'edit' 'string' '' 'tag' 'study_name' } ...
    {} };
guigeom = { [1]  [1] [0.2 1 3.5] [1] };

% define conditions
% -----------------
guigeom{end+1} = [];
for icond = 1:nConds
    if icond == 1, guigeom{end} = [ guigeom{end} 1 0.2];
    else           guigeom{end} = [ guigeom{end} 0.1 1 0.2];
    end
    if icond > 1, guispec{end+1} = {}; end
    guispec = { guispec{:}, {'style' 'text' 'string' [ 'Condition ' num2str(icond) ' name'] } {} };
end

% edit boxes for conditions
% -------------------------
guigeom{end+1} = [];
for icond = 1:nConds
    if icond == 1, guigeom{end} = [ guigeom{end} 1 0.2];
    else           guigeom{end} = [ guigeom{end} 0.1 1 0.2];
    end
    if icond > 1, guispec{end+1} = {}; end
    guispec = { guispec{:}, {'style' 'edit' 'string' '' 'tag' [ 'cond' num2str(icond) ] } {} };
end
guispec{end+1} = {};
guigeom{end+1} = [1];

% define dataset headers
% ----------------------
guigeom{end+1} = [];
for icond = 1:nConds
    if icond == 1, guigeom{end} = [ guigeom{end} 1 0.2];
    else           guigeom{end} = [ guigeom{end} 0.1 1 0.2];
    end
    if icond > 1, guispec{end+1} = {}; end
    guispec = { guispec{:}, {'style' 'text' 'string' ['Condition ' num2str(icond) ' datasets' ] } {} };
end

% create edit boxes
% -----------------
for index = 1:nSubjects
    guigeom{end+1} = [];
    for icond = 1:nConds
        if icond == 1, guigeom{end} = [ guigeom{end} 1 0.2];
        else           guigeom{end} = [ guigeom{end} 0.1 1 0.2];
        end
        select_com = ['[inputname, inputpath] = uigetfile2(''*.set;*.SET'', ''Choose dataset to add to STUDY -- pop_study()'');'...
                      'if inputname ~= 0,' ...
                      '   guiind = findobj(''parent'', gcbf, ''tag'', ''set' int2str(icond) '_' int2str(index) ''');' ...
                      '   set( guiind,''string'', fullfile(inputpath, inputname));' ...
                      'end; clear inputname inputpath;'];
        if icond > 1, guispec{end+1} = {}; end
        guispec = { guispec{:}, ...
                {'style' 'edit'       'string' ''    'tag' [ 'set' int2str(icond) '_' int2str(index) ] }, ...
                {'style' 'pushbutton' 'string' '...' 'Callback' select_com } };
    end
end

% last text
% ---------
textinfo = [  'When using more than 1 condition, datasets on each line must correspond to the same subject.' ];
guispec = { guispec{:}, {}, {'style' 'text' 'string' textinfo } };
guigeom = { guigeom{:} [1] [1] };

optiongui = { 'geometry', guigeom, ...
              'uilist'  , guispec, ...
              'helpcom' , 'pophelp(''pop_studyerp'')', ...
              'title'   , 'Create a new STUDY set -- pop_studyerp()' };
[result, userdat2, strhalt, outstruct] = inputgui(optiongui{:});
if isempty(result), return; end

% decode outstruct and build call to std_editset
% ----------------------------------------------
options  = { 'name' outstruct.study_name 'updatedat' 'off' };
commands = {};
for icond = 1:nConds

    % check that condition name is defined
    tagCond = ['cond' int2str(icond) ];
    if isempty(outstruct.(tagCond))
        outstruct.(tagCond) = [ 'condition ' int2str(icond) ];
    end

    for index = 1:nSubjects
        tagSet  = [ 'set' int2str(icond) '_' int2str(index) ];
        subject  = sprintf('S%2.2d', index);

        if ~isempty(outstruct.(tagSet))
            commands = { commands{:}, {'index' nConds*index+icond-1 'load' outstruct.(tagSet)  'subject' subject 'condition' outstruct.(tagCond) } };
        end
    end
end
options = { options{:}, 'commands', commands };

% call std_editset to create the STUDY
% ------------------------------------
com1 = sprintf( '[STUDY ALLEEG] = std_editset( STUDY, ALLEEG, %s );', vararg2str(options) );
[STUDY ALLEEG] = std_editset(STUDY, ALLEEG, options{:});
if exist(fullfile(ALLEEG(1).filepath, [ALLEEG(1).filename(1:end-4) '.daterp' ]))
    textmsg = [ 'WARNING: SOME ERP DATAFILES ALREADY EXIST, OVERWRITE THEM?' 10 ...
                '(the recomputed files will be identical to the existing ones, unless you modified the EEG' 10 ...
                'datasets or set a different ERP baseline using menu item STUDY > Precompute channel measures).' ];
    res = questdlg2(textmsg, 'Precomputed datafiles already present on disk', 'No', 'Yes', 'Yes');
    if strcmpi(res, 'No')
        error('User aborded precomputing ERPs');
    end;
end;

% call std_precomp for ERP (channels)
% -----------------------------------
com2 = '[STUDY ALLEEG] = std_precomp(STUDY, ALLEEG, ''channels'', ''interp'', ''on'', ''recompute'',''on'',''erp'',''on'');';
[STUDY ALLEEG] = std_precomp(STUDY, ALLEEG, 'channels','interp', 'on', 'recompute','on','erp','on');

% call std_erpplot to plot ERPs (channels)
% ----------------------------------------
com3 = 'tmpchanlocs = ALLEEG(1).chanlocs; STUDY = std_erpplot(STUDY, ALLEEG, ''channels'', { tmpchanlocs.labels }, ''plotconditions'', ''together'');';
tmpchanlocs = ALLEEG(1).chanlocs;
STUDY = std_erpplot(STUDY, ALLEEG, 'channels', { tmpchanlocs.labels }, 'plotconditions', 'together');
pos = get(gcf, 'position');
set(gcf, 'position', [10 pos(2) pos(3)*2 pos(4)*2]);

% call the STUDY plotting interface
% ---------------------------------
disp('Press OK to close plotting interface and save the STUDY');
disp('If you press CANCEL, the whole STUDY will be lost.');
[STUDY com4] = pop_chanplot(STUDY, ALLEEG);

com = sprintf('%s\n%s\n%s\n%s', com1, com2, com3, com4);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [zs] = besselJRobinzeros(c,k,N,freqCenter)

if ~exist('freqCenter','var')
    freqCenter = 0;
end

xmin = max(freqCenter - pi*N,0);
xmax = freqCenter + 2*pi*N; % Using zn ~ pi*n
if c == Inf
    % Returns the N first zeros of the function
    % Jk(x)
    zs = AllZeros(@(x)(besselj(k,x)),xmin,xmax,5*(N+freqCenter));

else
    % Returns the N first zeros of the function
    % c Jk(x) + xJk'(x)

    switch N
        case 0
            derBess = @(x)(-besselj(1,x));
        otherwise
            derBess = @(x)(0.5*(besselj(k-1,x)-besselj(k+1,x)));
    end

    zs = AllZeros(@(x)(c*besselj(k,x)+x*derBess(x)),xmin,xmax,5*(N+freqCenter));

end

[~,I] = sort(abs(zs-freqCenter));
zs = zs(I);

if zs(1)==0
    zs = zs(2:(N+1));
else
    zs = zs(1:N);
end

zs = zs(:);

end

function z=AllZeros(f,xmin,xmax,N)
% Inputs :
% f : function of one variable
% [xmin - xmax] : range where f is continuous containing zeros
% N : control of the minimum distance (xmax-xmin)/N between two zeros
if (nargin<4)
    N=100;
end
options=optimset('Display','off');
dx=(xmax-xmin)/N;
x2=xmin;
y2=f(x2);
z=[];
for i=1:N
    x1=x2;
    y1=y2;
    x2=xmin+i*dx;
    y2=f(x2);
    if (y1*y2<=0)                              % Rolle's theorem : one zeros (or more) present
        z=[z,fsolve(f,(x2*y1-x1*y2)/(y1-y2),options)]; % Linear approximation to guess the initial value in the [x1,x2] range.
    end
end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function mappedX = hlle(X, no_dims, k, eig_impl)
%HLLE Runs the Hessian LLE algorithm
%
%   mappedX = hlle(X, no_dims, k, eig_impl)
%
% Runs the Hessian LLE algorithm on dataset X to reduce its dimensionality
% to no_dims. The variable k specifies the number of nearest negihbors that
% is used.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('k', 'var')
        k = 12;
    end
    if ~exist('eig_impl', 'var')
        eig_impl = 'Matlab';
    end

    % Compute nearest neighbors
    if ischar(k)
        warning('Adaptive neighborhood selection often leads to problems in Hessian LLE.');
    end
    disp('Finding nearest neighbors...');
    [D, nind] = find_nn(X, k);
    max_k = size(nind, 2);

    % Size of original data
    n = size(X, 1);

    % Extra term count for quadratic form
    dp = no_dims * (no_dims + 1) / 2;
    W = sparse([], [], [], dp * n, n, dp * n * max_k);

    % For all datapoints
    disp('Building Hessian estimator for neighboring points...');
    for i=1:n
        % Center datapoints by substracting their mean
        tmp_ind = nind(i,:);
        tmp_ind = tmp_ind(tmp_ind ~= 0);
        kt = length(tmp_ind);
        thisx = X(tmp_ind,:);
        thisx = (thisx - repmat(mean(thisx, 1), kt, 1))';

        % Compute local coordinates (using SVD)
        [U, D, Vpr] = svd(thisx);
        if size(Vpr, 2) < no_dims
            no_dims = size(Vpr, 2);
            dp = no_dims * (no_dims + 1) / 2;
            warning(['Target dimensionality reduced to ' num2str(no_dims) '...']);
        end
        V = Vpr(:,1:no_dims);
		% Basically, the above is applying PCA to the neighborhood of Xi. The PCA mapping that is found
		% (and that is contained in V) is an approximation for the tangent space at Xi.

        % Build Hessian estimator
        clear Yi; clear Pii;
        ct = 0;
        for mm=1:no_dims
            startp = V(:,mm);
            for nn=1:length(mm:no_dims)
                indles = mm:no_dims;
                Yi(:,ct+nn) = startp .* (V(:,indles(nn)));
            end
            ct = ct + length(mm:no_dims);
        end
        Yi = [repmat(1, kt, 1) V Yi];

        % Gram-Schmidt orthogonalization (works different from Matlab QR function)
        [Yt, Orig] = mgs(Yi);
        Pii = Yt(:,no_dims + 2:end)';

        % Double check weights sum to 1
        for j=1:dp
            if sum(Pii(j,:)) > 0.0001
                tpp = Pii(j,:) ./ sum(Pii(j,:));
            else
                tpp = Pii(j,:);
            end

            % Fill weight matrix
            W((i - 1) * dp + j, tmp_ind) = tpp;
        end
    end

    % The weight matrix W is now entirely filled, perform eigenanalysis of W
    disp('Computing HLLE embedding (eigenanalysis)...');

    % Make sparse matrix that is inproduct of W
    G = W' * W;
	G(isnan(G)) = 0;
    G = sparse(G);

    % Clear some memory
    clear X thisx W D nind U D Vpr;

    % Perform eigendecomposition
    tol = 0;
    if strcmp(eig_impl, 'JDQR')
        options.Disp = 0;
        options.LSolver = 'bicgstab';
        [mappedX, eigenvals] = jdqr(G, no_dims + 1, tol, options);
    else
        options.disp = 0;
        options.issym = 1;
        options.isreal = 1;
        [mappedX, eigenvals] = eigs(G, no_dims + 1, tol, options);
    end

    % Sort eigenvalues and eigenvectors
    [eigenvals, ind] = sort(diag(eigenvals), 'ascend');
    if size(mappedX, 2) < no_dims + 1, no_dims = size(mappedX, 2) - 1; end
    mappedX = mappedX(:,ind(2:no_dims + 1));

    % Extract nonzero coordinates
    mappedX = mappedX(:,1:no_dims)' * sqrt(n);
    mappedX = mappedX';

repo: cluster


text:
function S = sum(d)
%SUM       Integration functional.
%   This function is deprecated. Use FUNCTIONALBLOCK.SUM.

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

S = linop( functionalBlock.sum(double(d)) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function mask = create_pad_mask(imsz, pad)
    if numel(pad) == 1
        pad = pad * ones(imsz);
    end

    if ndims(imsz) == 2
        mask = ones(imsz);
        mask(1:pad(1), :) = 0;
        mask(end-pad(1)+1:end,:) = 0;
        mask(:, 1:pad(2)) = 0;
        mask(:, end-pad(2)+1:end) = 0;
    end
end

repo: pTVreg


%% Quiver calculations
% These are calculations for the quiver dimensions as implemented in MATLAB
% (HG1) as in the |quiver.m| function.
%
% For HG2 and Octave, the situation might be different.
%
% A single quiver is defined as:
%
%                   C
%                    \
%                     \
%  A ----------------- B
%                     /
%                    /
%                   D
%
% To know the dimensions of the arrow head, MATLAB defines the quantities
%  alpha = beta = 0.33 that determine the coordinates of C and D as given below.

clc;
clear variables;
close all;

%% Parameters
try
    syms x y z u v w alpha beta epsilon real
catch
    warning('Symbolic toolbox not found. Interpret the values with care!');
    x = randn(); y = randn(); z = randn();
    u = randn(); v = randn(); w = randn();
end
alpha = 0.33;
beta = alpha;
epsilon = 0;
is2D = true;

%% Coordinates as defined in MATLAB
% Note that in 3D, the arrow head is oriented in a weird way. Let' just ignore
% that and only focus on 2D and use the same in 3D. Due to the lack
% of [u,v,w]-symmetry in those equations, the angle is bound to depend on the
% length of |delta|, i.e. something we don't know beforehand.
A = [x y z].';
delta = [u v w].';
B = A + delta;
C = B - alpha*[u+beta*(v+epsilon);
               v-beta*(u+epsilon)
               w];
D = B - alpha*[u-beta*(v+epsilon);
               v+beta*(u+epsilon)
               w];

if is2D
    A = A(1:2);
    B = B(1:2);
    C = C(1:2);
    D = D(1:2);
    delta = delta(1:2);
end

%% Calculating the angle of the arrowhead
% Calculate the cos(angle) using the inner product
unitVector = @(v) v/norm(v);
cosAngleBetween = @(a,b,c) unitVector(a-b).' * unitVector(c-b);

cosTwiceTheta = cosAngleBetween(C,B,D);
if isa(cosTwiceTheta, 'sym')
    cosTwiceTheta = simplify(cosTwiceTheta);
end

theta = acos(cosTwiceTheta) / 2

radToDeg = @(rads) (rads * 180 / pi);

thetaVal = radToDeg(theta)
try
    thetaVal = double(thetaVal)
end

% For the MATLAB parameters alpha=beta=0.33, we get theta = 18.263 degrees.


text:
function test_example_SAE
load mnist_uint8;

train_x = double(train_x)/255;
test_x  = double(test_x)/255;
train_y = double(train_y);
test_y  = double(test_y);

%%  ex1 train a 100 hidden unit SDAE and use it to initialize a FFNN
%  Setup and train a stacked denoising autoencoder (SDAE)
rand('state',0)
sae = saesetup([784 100]);
sae.ae{1}.activation_function       = 'sigm';
sae.ae{1}.learningRate              = 1;
sae.ae{1}.inputZeroMaskedFraction   = 0.5;
opts.numepochs =   1;
opts.batchsize = 100;
sae = saetrain(sae, train_x, opts);
visualize(sae.ae{1}.W{1}(:,2:end)')

% Use the SDAE to initialize a FFNN
nn = nnsetup([784 100 10]);
nn.activation_function              = 'sigm';
nn.learningRate                     = 1;
nn.W{1} = sae.ae{1}.W{1};

% Train the FFNN
opts.numepochs =   1;
opts.batchsize = 100;
nn = nntrain(nn, train_x, train_y, opts);
[er, bad] = nntest(nn, test_x, test_y);
assert(er < 0.16, 'Too big error');

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function ignore = ignorefields(purpose)

% IGNOREFIELDS returns a list of fields that can be present in the cfg structure that
% should be ignored at various places in the code, e.g. for provenance, history,
% size-checking, etc.

switch purpose

  case 'appendtimelock'
    ignore = {
      'cfg'
      'dimord'
      'elec'
      'fsample'
      'grad'
      'label'
      'negclusters'
      'negdistribution'
      'opto'
      'posclusters'
      'posdistribution'
      'sampleinfo' % this is dealt with explicitly
      'time'
      'topo'
      'topodimord'
      'topolabel'
      'trialinfo'  % this is dealt with explicitly
      'unmixing'
      'unmixingdimord'
      'brainordinate'
      };

  case 'appendfreq'
    ignore = {
      'cfg'
      'cumsumcnt'  % this is dealt with explicitly
      'cumtapcnt'  % this is dealt with explicitly
      'dimord'
      'elec'
      'freq'
      'grad'
      'label'
      'negclusters'
      'negdistribution'
      'opto'
      'posclusters'
      'posdistribution'
      'sampleinfo' % this is dealt with explicitly
      'time'
      'trialinfo'  % this is dealt with explicitly
      };

  case 'deface'
    ignore = {
      % some fields should be dealt with explicitly
      'dim'
      'hex'
      'pos'
      'tet'
      'transform'
      'tri'
      % some fields are irrelevant
      'cfg'
      'coordsys'
      'fid'
      'unit'
      };

  case 'pipeline'
    ignore = {
      % some fields that are always allowed to be present in the configuration
      'cfg'
      'inside'
      'leadfield'
      'previous'
      };

  case 'allowed'
    ignore = {
      % some fields that are always allowed to be present in the configuration
      'callinfo'
      'checkconfig'
      'checkpath'
      'checksize'
      'debug'
      'notification'
      'outputfilepresent'
      'previous'
      'progress'
      'showcallinfo'
      'spmversion'
      'toolbox'
      'trackcallinfo'
      'trackdatainfo'
      'trackusage'
      'version'
      'warning'
      };

  case {'rollback'}
    ignore = {
      % these should not be updated in rollback_provenance
      'callinfo'
      'checkconfig'
      'checksize'
      'debug'
      'notification'
      'previous'
      'showcallinfo'
      'trackcallinfo'
      'trackdatainfo'
      'trackusage'
      'version'
      'warning'
      };

  case {'provenance', 'history'}
    ignore = {
      % these should not be included in the provenance or history
      'checkconfig'
      'checksize'
      'debug'
      'notification'
      'reproducescript'
      'showcallinfo'
      'trackcallinfo'
      'trackdatainfo'
      'trackusage'
      'warning'
      };

  case {'reproducescript'}
    ignore = {
      % these should not be included in the output script
      'callinfo'
      'checkconfig'
      'checkpath'
      'checksize'
      'debug'
      'notification'
      'outputfilepresent'
      'progress'
      'reproducescript'
      'showcallinfo'
      'toolbox'
      'trackcallinfo'
      'trackdatainfo'
      'trackusage'
      'version'
      'warning'
      'event'
      };

  case 'checksize'
    ignore = {
      % the size of these fields should not be checked
      'artfctdef'
      'artifact'
      'checksize'
      'event'
      'hastoolbox'
      'previous'
      'trl'
      'trlold'
      };

  case 'makessense'
    ignore = {
      % these fields should not be used to check whether the trialinfo and sampleinfo make sense
      'cfg'
      'dimord'
      'elec'
      'freq'
      'fsample'
      'grad'
      'hdr'
      'label'
      'negclusters'
      'negdistribution'
      'opto'
      'posclusters'
      'posdistribution'
      'sampleinfo'
      'time'
      'trialinfo'
      'brainordinate'
      };

  case 'html'
    ignore = {
      % when generating a html-formatted pipeline, ignore data-like fields and fields that probably were not added by the user himself
      'checkconfig'
      'checksize'
      'debug'
      'event'
      'headmodel'
      'outputfilepresent'
      'previous'
      'progress'
      'showcallinfo'
      'sourcemodel'
      'trackcallinfo'
      'trackdatainfo'
      'trackusage'
      'warning'
      };

  case 'selectdata'
    ignore = {
      % these fields do not contain data and should be excluded
      'cfg'
      'coordsys'
      'dim'
      'elec'
      'fsample'
      'fsampleorig'
      'grad'
      'hdr'
      'negclusters'
      'negdistribution'
      'opto'
      'posclusters'
      'posdistribution'
      'topolabel'
      'transform'
      'unit'
      'brainordinate'
      'chantype'
      'chanunit'
      };

  case 'recursesize'
    ignore = {
      % these fields should not recursively be checked on their size
      'elec'
      'event'
      'grad'
      'grid'
      'headmodel'
      'headshape'
      'layout'
      'matrix'
      'mri'
      'neighbours'
      'opto'
      'sourcemodel'
      'vol'
      'trl'
      };

  otherwise
    ft_error('invalid purpose');
end % switch purpose


text:
function anonymize_dicm(src, rst, subID)
% Replace PatientName in dicom header with an ID without altering ANYTHING else.
%
%  anonymize_dicm('sourceFolder', 'anonymizedFolder', 3) % 'sub003' as ID
%
% All three input argument are optional:
%  1. A dicom file name or a folder containing dicom files to anonymize.
%  2. A file name (single source file) or a folder to save the result file(s).
%  3. Starting subject ID to replace PatientName, default 'sub001'. This can be
%     numeric, like 3 which will be converted into 'sub003', or full string like
%     'subj_04' specifying both the prefix and starting ID number.
%
% To avoid subject name confusion, this function replaces PatientName with an ID
% (third input argument). In case of multiple patients, the ID will be increased
% by 1 for each patient, and a file 'subjIDs_PatientNames.mat' will be saved
% into the source folder to indicate the correspondence between the original
% PatientName and assigned subject IDs.
%
% See also DICM_HDR, SORT_DICM, RENAME_DICM, DICM2NII

% 161226 Wrote it (xiangrui.li@gmail.com)
% 170103 Take care of confusion of multi-subjects; safer by using dicm tag.

% 3rd input can be either a number, or string like 'subj_001'
if nargin<3 || isempty(subID), subID = 'sub001'; end
if isnumeric(subID), subID = sprintf('sub%03g', subID); end
if ~isstrprop(subID(end), 'digit')
    error('subject ID must be a number or string in format of ''sub001''');
end
a = subID;
while ~all(isstrprop(a, 'digit'))
    ind = find(isstrprop(a, 'digit'), 1);
    a = a(ind:end);
end
iPN = str2double(a); % start ID number
n = numel(a);
idStr = [subID(1:numel(subID)-n) '%0' num2str(n) 'g']; % for sprintf

if nargin<1 || isempty(src)
    src = uigetdir(pwd, 'Select the folder containing dicom files');
    if isnumeric(src), return; end
end

if nargin<2 || isempty(rst)
    if isdir(src), def = src; else, def = fileparts(src); end
    if ~isdir(src) && (exist(src, 'file') || (iscell(src) && numel(src)==1))
        [rst, pth] = uiputfile([def filesep '*.dcm'], ...
            'Input file name to save the anonymized file');
        if isnumeric(rst), return; end
        rst = fullfile(pth, rst);
    else
        rst = uigetdir(def, 'Select or create a folder to save anonymized files');
        if isnumeric(rst), return; end
    end
end

if isdir(src)
    nams = dir(src);
    nams([nams.isdir]) = [];
    nams = {nams.name};
    nams = strcat(src, filesep, nams);
elseif ischar(src)
    nams = cellstr(src);
end
nFile = numel(nams);

dict = dicm_dict('', 'PatientName');
if nFile>1 && ~exist(rst, 'dir'), mkdir(rst); end
subjIDs = {}; PatientNames = {}; % for record keeping
for i = 1:nFile
    s = dicm_hdr(nams{i}, dict);
    if isempty(s), continue; end % not dicom

    if nFile == 1 % 2nd arg is file name
        nam = rst;
    else % 2nd arg is dir
        [~, nam, ext] = fileparts(s.Filename);
        nam = fullfile(rst, [nam ext]);
    end

    try
        pn = s.PatientName;
    catch % only a copy if no PatientName
        copyfile(s.Filename, nam, 'f');
        continue;
    end

    ip = strcmp(PatientNames, pn);
    if any(ip)
        subID = subjIDs{ip};
    else
        ip = numel(PatientNames) + 1;
        PatientNames{ip} = pn; %#ok<*AGROW>
        subID = sprintf(idStr, iPN+ip-1);
        if mod(numel(subID),2), subID(end+1) = char(0); end
        subjIDs{ip} = subID;
    end

    fid = fopen(s.Filename, 'r', 'l');
    b8 = fread(fid, inf, 'uint8=>uint8')';
    fclose(fid);

    try
        tsUID = s.TransferSyntaxUID;
        be = strcmp(tsUID, '1.2.840.10008.1.2.2');
        expl = ~strcmp(tsUID, '1.2.840.10008.1.2');
    catch
        be = false;
        expl = true;
    end
    if be, ed = 'b'; tag = char([0 16 0 16]); % '0010' '0010' 'PN' 'PatientName'
    else,  ed = 'l'; tag = char([16 0 16 0]);
    end
    if expl, tag = [tag 'PN']; end
    try n = s.PixelData.Start; catch, n = s.FileSize; end
    i0 = strfind(char(b8(1:n)), tag);
    if isempty(i0) % almost impossible
        fprintf(2, 'Failed to locate PatitentName in %s\n', s.Filename);
        copyfile(s.Filename, nam, 'f');
        continue;
    end
    i0 = i0(1) + numel(tag);
    n = double(b8(i0+(0:1)));
    if be, n = n(1)*256 + n(2); else, n = n(2)*256 + n(1); end % uint16

    fid = fopen(nam, 'w', ed);
    fwrite(fid, b8(1:i0-1), 'uint8'); % till len
    fwrite(fid, numel(subID), 'uint16'); % len: endian related
    fwrite(fid, subID, 'char'); % ID
    fwrite(fid, b8(i0+2+n:end), 'uint8'); % skip len & PatientName
    fclose(fid);
end

if numel(subjIDs)<2, return; end
subjIDs = deblank(subjIDs); %#ok remove padded null
matNam = fullfile(fileparts(nams{1}), 'subjIDs_PatientNames.mat');
try
    save(matNam, 'subjIDs', 'PatientNames');
catch me
    fprintf(2, '%s\n subjIDs_PatientNames saved into result folder.\n', me.message);
    matNam = fullfile(fileparts(nam), 'subjIDs_PatientNames.mat');
    save(matNam, 'subjIDs', 'PatientNames');
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [X,mP,spmP] = tor_get_physio(varargin)
% :Usage:
% ::
%
%     [X,mP,spmP] = tor_get_physio([mP],[spmP],[nvoxels],[doortho])
% arguments are optional, but you must enter them in this order.
%
% Tor Wager 10/21/02
%
% Get nuisance covariates likely to be related to physiological noise and head motion
% The algorithm:
%
% The program extracts raw/preprocessed image data from the ventricles (CSF space), as
% defined by a mask denoting which voxels are CSF for that subject.
% Either all voxels or a randomly selected subset [nvoxels] is subjected to
% principal components analysis, to determine regular patters of drift over time
% and across voxels.  Those patterns are expected to be related to global signal drift,
% head movement, and physiological noise, and are assumed to be UNrelated to the task
% of interest, by virtue of the fact that they occur in the ventricles.
%
% PCA is done twice on the timeseries' of CSF voxels.  The first time, PCA is done
% on the sums of squared values (not the correlations) of voxel timeserieses across
% the entire experiment, mean-centered based on the whole experiment.  Most of the
% coherent variation in this case is expected to be due to head movement and changes
% in shims/gradients/etc. from run to run.  The SS values are used because we want to
% weight the voxels with the highest variation most heavily, as they are presumably
% picking up most of this signal.  The first 3 eigenvariates (canonical timeseries)
% are saved.
%
% Following, a separate, second PCA is done on the correlation matrix of data
% within each session.  Session data for each voxel are mean-centered and scaled
% relative to the session (variance of each voxel = 1).  We do this because
% physiological noise-related signals may produce periodic signals of different
% magnitudes in different voxels, and we want to extract the most coherent signals
% we can within each session.  So these eigenvariates are expected to reflect
% primarily noise related to physiology (heart rate, respiration).  Up to 5 eigenvariates
% for each session are saved (nothing with eigenvalue < 1 is saved).
%
% Next, the CSF-related nuisance covariates (eigenvariates from PCA) are combined
% with existing nuisance covariates and intercept columns from the existing
% design matrix (SPMcfg xX).  The proportion of variance in each predictor of interest
% explained by this nuisance basis set is calculated using regression, and the
% nuisance covariates are orthogonalized with respect to each predictor of interest.
% There are good and bad results of this step.  The bad is that any signal that
% tracks the predictors is attributed to the task, not to noise, even if it's actually
% caused by physiological artifact.  So the orthogonalized basis set does not
% protect you from physiology or movement-related false positives.  However,
% the nuisance covariates are also unlikely to reduce power in estimating you effects
% of interest.  More importantly, it avoids false positives created when one
% predictor (A) is more highly correlated with the nuisance covariates than another
% (B).  In practice, betas for A will tend to be smaller than B, given the same
% actual response to both, and a random effects analysis on A-B will produce
% false positive activations.  Orthogonalization of the nuisance set precludes this.
%
% :Inputs:
%
%   **mP:**
%        CSF mask image file.  *_seg3.img output from SPM is appropriate
%        should be in same space and have same dims as functionals
%        but automatic reslicing is done if necessary.
%
%   **spmP:**
%        name (full path name preferred) of SPMcfg.mat file to use
%        This contains the design matrix and raw/preproc image file names to use.
%
%   **nvoxels:**
%        Number of CSF voxels to use in PCA analysis
%        More than 100 can be very slow and memory intensive.
%        Fewer than 100 voxels loads a different way, and may be slower.
%        Best is probably between 100 - 1000.  800 runs pretty fast.
%
%   **doortho:**
%        Orthogonalize nuisance covariates with respect to regs of interest
%        This assumes that any signal that covaries with the task is, in fact,
%        due to the task, so it gives you some bias towards finding positive results.
%        However, the alternative is that nuisance covariates may soak up variance
%        related to the task, and you'll miss activations.
%        In addition, if some regressors are more colinear with the nuisance set,
%        you can create false "activations" when comparing these regressors to other
%        ones.  This problem exists whether or not we choose to model nuisance
%        covariates.  One solution is to use the ortho when doing random effects analyses,
%        as the sign and magnitude of nuisance-related activations would not be expected to be
%        the same across subjects unless the variance was really task-related.
%        Default is 1, or "yes, do orthogonalization."
%
% for functions called, see this .m file.
%
% :Examples:
% ::
%
%    % get filenames for SPMcfg files and CSF mask for each subject
%    cd C:\Tor_Documents\CurrentExperiments\intext2\RESULTS\model1
%    spmP = get_filename('sub*','SPMcfg.mat');
%    cd C:\Tor_Documents\CurrentExperiments\intext2\
%    mP = get_filename('sub*','anatomy/nscalped_f*seg3.img');
%    % Now run:
%    for i = 1:size(mP,1)
%        tor_get_physio(mP(i,:),spmP(i,:),300);  % 300 voxels
%        pause(10); close all
%    end

% :Functions called:
%   - spm functions: spm_get, etc.
%   - timeseries2.m	(for < 100 voxels)
%   - read_hdr.m	(big-little endian dependent; validate for your data)
%   - timeseries3.m	(for > 100 voxels; uses SPM's image reading)
%   - reslice_imgs.m
%   - mask2voxel.m 	(only if ind2sub.m from Matlab is not found)

CSFprob = .95;      % this is the value a voxel in the mask img must have to be considered
                    % an 'on' value.
                    % if using SPM segmentation output (e.g., *_seg3) for the mask,
                    % this is something like the prob. of being in CSF
                    % low thresholds will result in HUGE eigenvalue problems
                    % and memory difficulties.

mypwd = pwd;
t1 = clock;

% ----------------------------------------------------------------------------------
% * set up input arguments
% ----------------------------------------------------------------------------------

if length(varargin) > 0
    mP = varargin{1};
else
    % mask filename
    mP = spm_get(1,'*img','Select CSF mask for this subject.');
end
if isempty(mP), mP = spm_get(1,'*img','Select CSF mask for this subject.');, end

if length(varargin) > 1
    spmP = varargin{2};
else
    spmP = spm_get(Inf,'SPMcfg.mat','Choose SPMcfg.mat file for this subject or Done to skip.');
end

if length(varargin) > 2, srand = varargin{3};, else, srand = 0;, end

if length(varargin) > 3, doortho = varargin{4};, else, doortho = 1;, end

% ----------------------------------------------------------------------------------
% * cd to SPM results directory so we can write output file there
% ----------------------------------------------------------------------------------
d = fileparts(spmP);
eval(['cd ' d])

diary physio_nuisance_covariates.out

% ----------------------------------------------------------------------------------
% * load SPMcfg.mat file, which contains all relevant info except mask image
% ----------------------------------------------------------------------------------

        load(spmP)
        nsess = length(Sess);
        P = str2mat(VY.fname);
	if ~(exist(deblank(P(1,:))) == 2)
		disp(['Looking for: ' P(1,:)])
		disp(['Can''t find original img files!! Please specify.'])
		P = spm_get(Inf,'*.img','Select raw image files.');
		VY = spm_vol(P);
		disp(['VY has been modified!!! Using: ' P(1,:)])
	end
        for i = 1:length(Sess), nimgs(i) = size(Sess{i}.row,2);,end

%  [nsess] = spm_input_ui('Enter number of runs, or 0 to choose SPM.mat file [recommended]',.1,'i',[],1);

% ----------------------------------------------------------------------------------
% * reslice the CSF mask if necessary, to be in space of functionals!
% ----------------------------------------------------------------------------------

mV = spm_vol(mP);
pV = spm_vol(P(1,:));
if any(mV.dim(1:3) - pV.dim(1:3)) | any(any(pV.mat(1:3,1:3) - mV.mat(1:3,1:3)))
    [d f e] = fileparts(mP); reslice_imgs(P(1,:),mP,0);
    mP = fullfile(d,['r' f e]);
    disp(['Resliced mask to space of functionals: ' mP])
    mV = spm_vol(mP);
end

% ----------------------------------------------------------------------------------
% * load and check mask
% ----------------------------------------------------------------------------------

mv = spm_read_vols(mV); mv = mv > CSFprob;
fprintf(1,'\nCSF mask has %3.0f voxels out of %3.0f total.\t',sum(mv(:)),prod(size(mv)))
try
    % if we have the right toolbox...
    [x y z] = ind2sub(size(mv),find(mv)); [XYZ] = [x y z];
catch
    disp('elmat toolbox function ind2sub not found; using SLOWER version mask2voxel.m')
    XYZ = mask2voxel(mv);
end

tmp = sum(sum(mv));
if ~findobj('Tag','Graphics'), spm fmri; figure(findobj('Tag','Graphics'));,end
if ~findobj('Tag','Interactive'), spm fmri;,end
%imagesc(mv(:,:,find(tmp == max(tmp)))); colormap gray; warning off; title([mP]), warning on,;drawnow
spm_check_registration(str2mat(mP,P(1,:)));

if sum(tmp) > 100,
    disp(['More than 100 CSF voxels in mask - this may be computationally intensive!'])
    % if srand is not entered as an input argument, prompt.
    if ~srand,
        figure(findobj('Tag','Interactive'))
        srand = spm_input_ui('Enter n vox to use, or 0 for all',.1,'i',[],1);,
    end
    if srand > sum(tmp),
	disp(['More voxels requested than available at threshold ' num2str(CSFprob) ': using ' num2str(max(tmp))])
	srand = sum(tmp);
    end

else
    srand = 0;
end

if srand
    wh = rand(size(XYZ,1),1) * size(XYZ,1);
    wh2 = sort(wh); wh2 = wh2(1:srand);
    for i = 1:srand, XYZ2(i,:) = XYZ(find(wh == wh2(i)),:);, end
    XYZ = XYZ2;
end

% ----------------------------------------------------------------------------------
% * Load images and mask with CSF
% ----------------------------------------------------------------------------------
% M1 is unscaled timeseries over all sessions; eigs computed over whole experiment
% M2 is timeseries scaled within each session; (cell array)
%   eigs are computed within session

fprintf(1,'\nEigenvariates based on %3.0f voxels \t',size(XYZ,1))
fprintf(1,'\nLoading volumes, extracting voxels, and scaling \t')

wh = [0 cumsum(nimgs)];
M1 = [];
M2 = [];

if size(XYZ,1) < 101
    % this is super slow for large n!
    ts = timeseries2('multi',P,struct('coords',XYZ));
    M1 = ts.indiv;

    % adjust M2 to mean 0 var 1
    for j = 1:size(M1,2)
        for i = 1:nsess
            ind = wh(i)+1:wh(i+1);
            M2{i}(:,j) = (M1(ind,j) - mean(M1(ind,j))) ./ std(M1(ind,j));
        end
    end

else

for i = 1:nsess
    subP = P(wh(i)+1:wh(i+1),:);
    fprintf(1,'.')
    ts = timeseries3(XYZ,subP);

    % save 2 matrices, one for session-mean centered and scaled and one uncentered (until overall mean is known)
    M1 = [M1;ts.all_data];

    for j = 1:size(ts.all_data,2)
        % center and scale to variance = 1
        ts.all_data(:,j) = (ts.all_data(:,j) - mean(ts.all_data(:,j))) ./ std(ts.all_data(:,j));
    end
    M2{i} = [ts.all_data];
end

end

fprintf(1,'Done.\n')

% ----------------------------------------------------------------------------------
% * Mean-center first (unscaled) timeseries and check for NaN or Inf values
% ----------------------------------------------------------------------------------
for i = 1:size(M1,2)
    M1(:,i) = M1(:,i) - mean(M1(:,i));
end

ex = any(isnan(M1) | isinf(M1)) | all(M1 == 0);
if any(ex)
	disp(['WARNING! NaN or Inf values for ' num2str(sum(ex)) ' voxels in timeseries!!  Mis-registration of funct and anatomy?'])
	M1(:,find(ex)) = [];
end

figure;imagesc(M1); title('Overall timeseries (y) for all voxels (x)')

% ----------------------------------------------------------------------------------
% * Find Principal Components overall
% ----------------------------------------------------------------------------------
fprintf(1,'\nPrincipal components overall ')
[eigvec,eigval] = eig(M1'*M1);
eigvec = eigvec(:,end-2:end);
X = M1 * eigvec;
for i = 1:size(X,2), X(:,i) = (X(:,i) - mean(X(:,i))) ./ std(X(:,i));, end
figure;subplot 221; plot(eigvec); title('Weights on original variables (=voxels) == eigenvectors')
subplot 222; plot(X),title('Components (X * v)')
legend({'Comp 3' 'Comp 2' 'Comp 1'},0)
    subplot 223; bar(diag(eigval));, title('Scree plot for eigenvalues')

    xx = abs(fft(X)); xxx = (1:size(xx,1)) ./ (xX.RT * size(xx,1));
    subplot 224; plot(xxx(1:round(length(xx)./2)),xx(1:round(length(xx)./2),:));, title('FFT of components')
    xlabel('Frequency (Hz)')
drawnow

% ----------------------------------------------------------------------------------
% * Find Principal Components for each session
% ----------------------------------------------------------------------------------
wh = [0 cumsum(nimgs)];

for i = 1:nsess
    fprintf(1,'\nPrincipal components for session %1.0f',i)

    ex = any(isnan(M2{i}) | isinf(M2{i}));
    if any(ex)
	disp('')
	disp(['WARNING! NaN or Inf values for ' num2str(sum(ex)) ' voxels in sess ' num2str(i) '!!  Mis-registration of funct and anatomy?'])
	M2{i}(:,find(ex)) = [];
    end

    [eigvec,eigval] = eig(corrcoef(M2{i}));
    fprintf(1,'\t%3.0f eigenvalues > 1\t',sum(diag(eigval) > 1))

    num2save = min(sum(diag(eigval)>1),5);    % save at most 5 eigenvalues from this session
    eigvec = eigvec(:,end-num2save+1:end);
    X2 = M2{i} * eigvec;
    for j = 1:size(X2,2), X2(:,j) = (X2(:,j) - mean(X2(:,j))) ./ std(X2(:,j));, end

    figure;subplot 221; plot(eigvec); title(['Weights (eigenvectors) for session ' num2str(i)])
    subplot 222; plot(X2),title('Components (X * v)')
    subplot 223; bar(diag(eigval));, title('Scree plot for eigenvalues')

    xx = abs(fft(X2)); xxx = (1:size(xx,1)) ./ (xX.RT * size(xx,1));
    subplot 224; plot(xxx(1:round(length(xx)./2)),xx(1:round(length(xx)./2),:));, title('FFT of components')
    xlabel('Frequency (Hz)')

    drawnow

    % pad with zeros to get in the right session
    zbef = zeros(wh(i),size(X2,2));
    X2 = [zbef; X2];
    zaft = zeros(size(X,1) - size(X2,1),size(X2,2));
    X2 = [X2; zaft];

    X = [X X2];
end

X(:,end+1:end+length(xX.iB)) = xX.X(:,xX.iB); % add the nuisance covariates already in xX

% ----------------------------------------------------------------------------------
% * Correlate / Regress design vectors on components
% ----------------------------------------------------------------------------------
px = X * pinv(X);
fprintf(1,'\n Variation in predictors explained by nuisance covariates before orthogonalization')
fprintf(1,'\n Differences among predictors could create false activations without orthogonalization.')

for i = 1:length(xX.iC)
    r = xX.X(:,xX.iC(i)) - px * xX.X(:,xX.iC(i));       % residuals
    pve = 1 - ((r' * r) ./ (xX.X(:,xX.iC(i))' * xX.X(:,xX.iC(i)))); % percentage of variance explained
    fprintf(1,'\n%s\t%3.2f%%',xX.Xnames{xX.iC(i)},100*pve)
end

% ----------------------------------------------------------------------------------
% * Orthogonalize components (?)
% This isn't a good idea from the standpoint of misattributing noise variance to signal
% but it will prevent artifactual activations based on differential correlations with
% nuisance subspace among conditions.  Lesser of two evils?
% ----------------------------------------------------------------------------------

if doortho
	fprintf(1,'\n Orthogonalizing nuisance covariates wrt model and scaling')
	mX = xX.X(xX.iC); mX(:,end+1) = 1;
	px = mX * pinv(mX);

	for i = 1:size(X,2)
    	X(:,i) = X(:,i) - px * X(:,i);                      % residuals
    	if ~(std(X(:,i)) == 0)
    		X(:,i) = (X(:,i) - mean(X(:,i))) ./ std(X(:,i));    % re-scale
    	end
	end

	figure; imagesc(X); colormap gray; title('Found and orthogonalized nuisance covariates')
	Xn = X; save Nuisance_covariates Xn
else
	figure; imagesc(X); colormap gray; title('Found nuisance covariates')
	Xn = X; save Nuisance_covariates_ortho Xn
end

% modify SPMcfg.mat
% just don't do this, it gets messy.
% Better to just add them by loading Nuisance_covariates
% and entering as user-spec regressors.
% add_nuisance_to_SPMcfg(Xn);

fprintf(1,'\n Total running time is %3.2f s.\n',etime(clock,t1))
eval(['cd ' mypwd])

diary off
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function img_patch = init_point_vis(img_patch, first_x, third_x, fourth_x, sixth_x, first_y, third_y, fourth_y, sixth_y, xmin, ymin)
    %%normalize 'em bitch
    first_x = max(first_x, xmin) - xmin;
    third_x = max(third_x, xmin) - xmin;
    fourth_x = max(fourth_x, xmin) - xmin;
    sixth_x = max(sixth_x, xmin) - xmin;
    first_y = max(first_y, ymin) - ymin;
    third_y = max(third_y, ymin) - ymin;
    fourth_y = max(fourth_y, ymin) - ymin;
    sixth_y = max(sixth_y, ymin) - ymin;
    img_patch = insertShape(img_patch, 'Circle', [first_x first_y 1], 'Color', 'red', 'LineWidth', 3);
    img_patch = insertShape(img_patch, 'Circle', [third_x third_y 1], 'Color', 'green', 'LineWidth', 3);
    img_patch = insertShape(img_patch, 'Circle', [fourth_x fourth_y 1], 'Color', 'blue', 'LineWidth', 3);
    img_patch = insertShape(img_patch, 'Circle', [sixth_x sixth_y 1], 'Color', 'yellow', 'LineWidth', 3);
end

author: cs-chan

repo: Total-Text-Dataset


text:
function [MDP] = spm_MDP_VB_sleep(MDP,BMR)
% Bayesian model reduction (sleep) for MDP models
% FORMAT [MDP] = spm_MDP_VB_sleep(MDP,BMR)
%
% MDP  - (inverted) MDP structure
%
% BMR.g - modality [default: 1]
% BMR.o - outcomes - that induce REM [default: {}]
% BMR.x - increase in concentration parameters for BMR [default: 8]
% BMR.f - hearing factors to sum over [default: 0]
% BMR.T - log Bayes factor threshold [default: 1/4]
% BMR.m - indicator function to enable BMR [@(i,i1,i2,i3,i4)1]
%
%
% MDP  - (reduced) model structure: with reduced MDP.a
%
% This routine optimises the hyperparameters of a NDP model (i.e.,
% concentration parameters encoding probabilities. It uses Bayesian model
% reduction to evaluate the evidence for models with and without a
% particular parameter in the columns of MDP.a (c.f., SWS)
%
% If specified, the scheme will then recompute posterior beliefs about the
% model parameters based upon (fictive) outcomes generated under its
% (reduced) generative model.(c.f., REM sleep)
%
% See also: spm_MDP_log_evidence.m, spm_MDP_VB and spm_MDP_VB_X.m
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_VB_sleep.m 7679 2019-10-24 15:54:07Z spm $

% deal with a sequence of trials
%==========================================================================

% BMR options
%--------------------------------------------------------------------------
try, g   = BMR.g; catch, g = 1;   end
try, o   = BMR.o; catch, o = {};  end
try, x   = BMR.x; catch, x = 8;   end
try, f   = BMR.f; catch, f = 0;   end
try, T   = BMR.T; catch, T = 1/4; end

% model selection function
%--------------------------------------------------------------------------
if isfield(BMR,'m')
    m = BMR.m;
else
    m = @(i,i1,i2,i3,i4)1;
end

% Baysian model reduction - parameters
%--------------------------------------------------------------------------
if isfield(MDP,'a')
    [sa,ra] = spm_MDP_VB_prune(MDP(end).a{g},MDP(1).a0{g},f,x,T,m);
end

% reiterate expectation maximisation (or rapid eye movement sleep)
%--------------------------------------------------------------------------
N  = numel(o);
if N

    % remove previous experience
    %----------------------------------------------------------------------
    REM  = MDP;
    try, REM = rmfield(REM,'s'); end
    try, REM = rmfield(REM,'o'); end
    try, REM = rmfield(REM,'u'); end

    % and install a generative process and reset priors
    %----------------------------------------------------------------------
    REM.a{g}  = ra;
    REM.a0{g} = ra;
    REM.o = o{1};
    for i = 1:N
        REM(i)   = REM(1);
        REM(i).o = o{i};
    end

    % Bayesian updating and updated parameters
    %----------------------------------------------------------------------
    REM    = spm_MDP_VB_X(REM);
    MDP.a  = REM(N).a;
    MDP.a0 = REM(N).a0;

else

    % otherwise, use reduced posteriors and priors
    %----------------------------------------------------------------------
    MDP.a{g}  = sa;
    MDP.a0{g} = ra;

end

function [sA,nA] = spm_MDP_VB_prune(qA,pA,f,x,T,m)
% FORMAT [sA,nA] = spm_MDP_VB_prune(qA,pA,f,x,T,m)
% qA - posterior expectations
% pA - prior expectations
% f  - hidden factor to integrate over [defult: 0]
% x  - prior counts [default: 8]
% T  - threshold for Bayesian model reduction [default: three]
%
% sA - reduced posterior expectations
% rA - reduced prior expectations
%__________________________________________________________________________

% defaults
%--------------------------------------------------------------------------
if nargin < 4, T = 3; end
if nargin < 3, x = 8; end

% identify noninformative subspaces
%--------------------------------------------------------------------------
ndim  = ndims(pA);
for i = 1:ndim
    d    = 1:ndim;
    d(i) = [];
    s    = pA < 4 & pA > 0;
    for j = 1:(ndim - 1)
        s = sum(s,d(j));
    end
    ind{i} = find(s(:));
end

% motifs of salient (low free free energy) abilities
%--------------------------------------------------------------------------
mot{1} = [1 0 0;0 1 0;0 0 1];
mot{2} = [1 0 0;0 0 1;0 1 0];
mot{3} = [0 1 0;1 0 0;0 0 1];
mot{4} = [0 1 0;0 0 1;1 0 0];
mot{5} = [0 0 1;1 0 0;0 1 0];
mot{6} = [0 0 1;0 1 0;1 0 0];

% model space: additional concentration parameters (i.e., precision)
%--------------------------------------------------------------------------
nmot  = numel(mot);
for i = 1:nmot*nmot;
    rA{i} = spm_zeros(pA);
end

for m1 = 1:numel(mot)
    [i1,i2] = find(mot{m1});
    for m2 = 1:numel(mot)
        i  = nmot*(m1 - 1) + m2;
        for j = 1:numel(i1)
            dA    = spm_zeros(pA);
            dA(1:3,i1(j),:,i2(j),4) = mot{m2};
            rA{i} = rA{i} + dA*x;
        end
    end
end

% score models using Bayesian model reduction
%--------------------------------------------------------------------------
for i = 1:numel(rA);
    G    = spm_MDP_log_evidence(qA,pA,pA + rA{i});
    F(i) = sum(G(isfinite(G)));
end

% find any model that has greater evidence than the parent model
%--------------------------------------------------------------------------
[F,j] = min(F);
if (F + T) < 0
    rA = rA{j};
else
    rA = spm_zeros(pA);
end

% reduced posterior and prior
%--------------------------------------------------------------------------
sA     = qA;
nA     = pA;
for i1 = 1:size(qA,2)
    for i2 = 1:size(qA,3)
        for i3 = 1:size(qA,4)
            for i4 = 1:size(qA,5)

                % get posteriors, priors and reduced priors
                %----------------------------------------------------------
                p  = pA(:,i1,i2,i3,i4);
                q  = qA(:,i1,i2,i3,i4);
                r  = rA(:,i1,i2,i3,i4);
                j  = find(p);
                p  = p(j);
                q  = q(j);
                r  = j(find(r(j)));

                % redistribute concentration parameters if indicated
                %----------------------------------------------------------
                if numel(r);
                    sA(:,i1,i2,i3,i4) = 0;
                    nA(:,i1,i2,i3,i4) = 0;
                    sA(r,i1,i2,i3,i4) = sum(q);
                    nA(r,i1,i2,i3,i4) = sum(p);
                else
                    sA(j,i1,i2,i3,i4) = q;
                    nA(j,i1,i2,i3,i4) = p;
                end
            end
        end
    end
end

return

% alternative formulation with specified indicator functions
%==========================================================================

% defaults
%--------------------------------------------------------------------------
if nargin < 5, m = @(i,i1,i2,i3,i4)1; end
if nargin < 4, T = 3; end
if nargin < 3, x = 8; end
if nargin < 2, f = 0; end

% column-wise model comparison
%--------------------------------------------------------------------------
for i1 = 1:size(qA,2)
    for i2 = 1:size(qA,3)
        for i3 = 1:size(qA,4)
            for i4 = 1:size(qA,5)

                % get posteriors, priors and cycle over reduced priors
                %----------------------------------------------------------
                p  = pA(:,i1,i2,i3,i4);
                q  = qA(:,i1,i2,i3,i4);
                j  = find(p);
                p  = p(j);
                q  = q(j);

                % informative state?
                %----------------------------------------------------------
                F  = 0;
                if length(j) > 1
                    for i = 1:length(j);
                        if m(i,i1,i2,i3,i4)
                            r    = p;
                            r(i) = r(i) + x;
                            F(i) = spm_MDP_log_evidence(q,p,r);
                        else
                            F(i) = 16;
                        end
                    end
                end

                % eliminate parameter
                %----------------------------------------------------------
                [F,i] = min(F);
                mF(i1,i2,i3,i4) = F;
                iF(i1,i2,i3,i4) = j(i);
            end
        end
    end
end

% pool over s{f}
%---------------------------------------------------------------------
if f, sF = sum(mF,f); end

% column-wise reduction
%--------------------------------------------------------------------------
sA     = qA;
rA     = pA;
for i1 = 1:size(qA,2)
    for i2 = 1:size(qA,3)
        for i3 = 1:size(qA,4)
            for i4 = 1:size(qA,5)

                % get posteriors, priors and cycle over reduced priors
                %----------------------------------------------------------
                p  = pA(:,i1,i2,i3,i4);
                q  = qA(:,i1,i2,i3,i4);
                i  = iF(i1,i2,i3,i4);
                j  = find(p);
                p  = p(j);
                q  = q(j);

                % BMC
                %----------------------------------------------------------
                if f == 0
                    F = mF(i1,i2,i3,i4);
                elseif f == 1
                    F = sF( 1,i2,i3,i4);
                elseif f == 2
                    F = sF(i1, 1,i3,i4);
                elseif f == 3
                    F = sF(i1,i2, 1,i4);
                elseif f == 4
                    F = sF(i1,i2,i3, 1);
                end

                % eliminate parameter
                %----------------------------------------------------------
                if F < - T;
                    sA(:,i1,i2,i3,i4) = 0;
                    rA(:,i1,i2,i3,i4) = 0;
                    sA(i,i1,i2,i3,i4) = sum(q);
                    rA(i,i1,i2,i3,i4) = sum(p);
                else
                    sA(j,i1,i2,i3,i4) = q;
                    rA(j,i1,i2,i3,i4) = p;
                end

            end
        end
    end
end

author: spm


text:
clear all
clc

% s - number of sparse coefficients
% N - length of each input vector
% n - length of output vector
% r - number of measurements

s = 15; N = 100; n = 60; r = 10;

t = randperm(N);
X0 = zeros(N,r);

for i = 1:r
    X0(t(1:s),i) = randn(s,1);
end

M = randn(n,N); Mop = opMatrix(M);
Sop = opDCT(N);

S = max(svd(M));

%%
F = opFoG(Mop, Sop);

for i = 1:r
    Y(:,i) = F(X0(:,i),1);
end

X = UnconSynthMMV(Y, F, 1e-6, 10, S);
norm(X0-X,'fro')/norm(X0,'fro')

X = SynthMMV(Y, F, 1e-6, S);
norm(X0-X,'fro')/norm(X0,'fro')

%%
for i = 1:r
    Z0(:,i) = Sop(X0(:,i),2);
    Y(:,i) = Mop(Z0(:,i),1);
end

Z = UnconAnaMMV(Y, M, Sop, 1e-6, 10, S, 1);
norm(Z0-Z,'fro')/norm(Z0,'fro')

Z = AnaMMV(Y, M, Sop, 1e-6, S, 1);
norm(Z0-Z,'fro')/norm(Z0,'fro')


text:
function perf=mae(e)
%
% calculate the absolute error of the given errors
%
%  'perf = mae(E);'
%
% see also:
%    mse, linf
%

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

perf = sum(sum(abs(e))) / numel(e);

author: peterhcharlton

repo: RRest


text:
function log_evidence = ldae_is_variants(words, topics, topic_prior, num_samples, variant, variant_iters)
%LDAE_IS_VARIANTS approximate evidence of LDA model, importance sampling with choice of q-distributions
%
% log_evidence = ldae_is_variants(words, topics, topic_prior[, num_samples=1000[, variant=3[, variant_iters=1]]]);
%
% Inputs:
%             words 1xNd
%            topics TxV each row is a distribution over a vocabulary of size V
%       topic_prior 1xT parameters of Dirichlet from which document topic vector is drawn
%       num_samples 1x1 default 1000
%           variant 1x1 1: prior, 2: q(z) = \prod p(z_n|w_n), 3: q(z) = \prod(z_n|w_n,hacky_pseudo_counts)
%                       Default is 3, as that seems to work best.
%     variant_iters 1x1 If variant needs iterative updates, use this number
%
% Outputs:
%     log_evidence  1x1

% Iain Murray, January 2009

[T, V] = size(topics);
Nd = length(words);

if ~exist('num_samples', 'var')
    num_samples = 1000;
end
if ~exist('variant', 'var')
    variant = 3;
end
if ~exist('variant_iters', 'var')
    variant_iters = 1;
end

topic_prior = topic_prior(:)';
topic_alpha = sum(topic_prior);
%topic_mean = topic_prior / topic_alpha;

if variant == 1
    % Importance sample from prior
    qstar = repmat(topic_prior', 1, Nd); % T x Nd
    qq = bsxfun(@rdivide, qstar, sum(qstar, 1));
else
    % Take w_n into account when picking z_n
    qstar = bsxfun(@times, topic_prior', topics(:, words)); % T x Nd
    qq = bsxfun(@rdivide, qstar, sum(qstar, 1));

    if variant == 3
        for i = 1:variant_iters
            % Now create pseudo-counts from qq and recompute qq using them
            pseudo_counts = bsxfun(@minus, topic_prior' + sum(qq, 2), qq);
            qstar = bsxfun(@times, pseudo_counts, topics(:, words)); % T x Nd
            qq = bsxfun(@rdivide, qstar, sum(qstar, 1));
        end
    end
end

% Draw samples from the q-distribution
samples = zeros(Nd, num_samples);
for n = 1:Nd
    samples(n, :) = discreternd(num_samples, qq(:, n))'; % Nd x num_samples
end

% Evaluate P(z, v) at samples and compare to q-distribution
Nk = histc(samples, 1:T, 1); % T x num_samples
log_pz = sum(gammaln(bsxfun(@plus, Nk, topic_prior')), 1) ...
        + gammaln(topic_alpha) - sum(gammaln(topic_prior)) ...
        - gammaln(Nd + topic_alpha); % 1 x num_samples
log_w_given_z = zeros(1, num_samples);
for n = 1:Nd
    log_w_given_z = log_w_given_z + log(topics(samples(n,:), words(n)))';
end
log_joint = log_pz + log_w_given_z;
%
log_qq = zeros(1, num_samples);
for n = 1:Nd
    log_qq = log_qq + log(qq(samples(n,:), n))';
end
log_weights = log_joint - log_qq;
log_evidence = logsumexp(log_weights(:)) - log(length(log_weights));

author: jacobeisenstein

repo: SAGE


text:
% This script train an LSTM or DNN based clean speech log spectrogram
% predictor, using simulated data of Reverb challenge.
%
% Xiong Xiao, Nanyang Technological University, Singapore
% Last modified: 08 Feb 2016.
%
%clear
function TrainEnhanceNet_Regression(modelType, hiddenLayerSize, hiddenLayerSizeFF, DeltaGenerationType, learning_rate, useGPU)
nUtt4Iteration = 10000;
nSequencePerMinibatch = 40;
useFileName = 1;

addpath('local', '../beamforming/lib', '../dereverb/local');     % remember to run addMyPath.m in SignalGraph root directory first.
% define SignalGraph settings. Type ParseOptions2() in command line to see the options.
para.IO.nStream = 2;        % two data streams, one is input signal, one is clean target log spectrogram
para.IO.mode = 'dynamicDistortion';     % distorted speech is generated dynamically during training. So each iteration uses different data.
para.NET.sequential = 1;    % do not randomize at frame level, randomize at sentence level
para.NET.variableLengthMinibatch = 1;   % a minibatch may contain multiple sentences of different lengths
para.NET.nSequencePerMinibatch = nSequencePerMinibatch;    % number of sentences per minibatch
para.NET.nSequencePerMinibatchCV = nSequencePerMinibatch; % number of sentences per minibatch for CV
para.NET.maxNumSentInBlock = 100;       % maximum number of sentences in a block
para.NET.L2weight = 3e-4;
para.NET.learning_rate = learning_rate;
para.NET.momentum = [0.9];
para.NET.learning_rate_decay_rate = 0.999;
para.NET.reduceLearnRateSpeed = 0.9;
para.NET.gradientClipThreshold = 1;
para.NET.weight_clip = 3;
para.useGPU     = useGPU;                  % don't use GPU for single sentence minibatch for LSTM. It's even slower than CPU.
para.minItr     = 100;
para.displayGPUstatus = 1;
para.displayInterval = ceil(para.NET.maxNumSentInBlock / para.NET.nSequencePerMinibatch / 10);
para.skipInitialEval = 1;

para.local.clean_wav_root = ChoosePath4OS({'E:\Data\LibriSpeech\train-clean-100', '/media/xiaoxiong/OS/data1/G/Libri/LibriSpeech/train-clean-100'});
para.local.clean_wav_ext = 'flac';
para.local.rir_wav_root = ChoosePath4OS({'E:\Data\ReverbEstimation\RIR_T60_1ch', '/media/xiaoxiong/DATA1/data1/ReverbEstimation/RIR_T60_1ch'});
para.local.rir_wav_ext = 'wav';
para.local.noise_wav_root = ChoosePath4OS({'E:\Data\musan\noise', '/media/xiaoxiong/DATA1/data1/musan/noise'});
para.local.noise_wav_ext = 'wav';
para.local.useFileName = useFileName;      % if set to 0, load all training data to memory. otherwise, only load file names.
para.local.seglen = 100;
para.local.segshift = 100;

para.topology.useChannel = 1;
para.topology.useWav = 1;
para.topology.RegressionNetType = modelType;
para.topology.useCMN = 0;
para.topology.DeltaGenerationType = DeltaGenerationType;   % Choose 'DeltaByEqn' or 'DeltaByAffine'
para.topology.MSECostWeightSDA = [1 4.5 10];
para.topology.hiddenLayerSize = hiddenLayerSize;
para.topology.hiddenLayerSizeFF = hiddenLayerSizeFF;
para = ConfigDereverbNet_Regression(para);

para.IO.DynamicDistortion.fs = para.topology.fs;
para.IO.DynamicDistortion.nUtt4Iteration = nUtt4Iteration;   % dynamically generate 10000 distorted sentences for each training iteration
para.IO.DynamicDistortion.nHours4Iteration = 10;    % dynamically generate 10 hours of distorted speech for each training iteration
para.IO.DynamicDistortion.SNR_PDF = 'uniform';      % distribution of SNR, choose [uniform|normal]
para.IO.DynamicDistortion.SNR_para = [-20 30];      % parameters of SNR PDF. If use uniform, it is the lowest and highest SNR allowed.
                                                    % if use normal, it is the mean and variance of SNR.
para.IO.DynamicDistortion.seglen = 100;             % cut distorted sentences into equal length segments.
para.IO.DynamicDistortion.segshift = 100;
para.IO.DynamicDistortion.frame_len = para.topology.frame_len;
para.IO.DynamicDistortion.frame_shift = para.topology.frame_shift;

[Data_tr, para] = LoadWavRIRNoise_Libri(para, 1);

[layer, para] = Build_EnhanceNet_Regression(Data_tr, para);

% load the cv data
para.local.cv_wav_root = ChoosePath4OS({'D:\Data\NoiseData\Libri\LibriSpeech\dev-parallel', '/media/xiaoxiong/DATA1/data1/Libri/LibriSpeech/dev-parallel'});
para.local.useFileName = 0;
[Data_cv, para] = LoadParallelWav_Libri(para, 2);

% generate directory and file names to store the networks.
if para.topology.useCMN; para.output = 'nnet/EnhanceRegression.CMN';
else;     para.output = 'nnet/EnhanceRegression.noCMN';    end
para.output = sprintf('%s.%s.MbSize%d.U%d.%d-%s', para.output, para.topology.DeltaGenerationType, ...
    para.NET.nSequencePerMinibatch, length(Data_tr(1).data),  3*(para.topology.fft_len/2+1), para.topology.RegressionNetType);
for i=1:length(para.topology.hiddenLayerSize)
    para.output = sprintf('%s-%d', para.output, para.topology.hiddenLayerSize(i));
end
if ~isempty(para.topology.hiddenLayerSizeFF)
    para.output = [para.output '-DNN'];
    for i=1:length(para.topology.hiddenLayerSizeFF)
        para.output = sprintf('%s-%d', para.output, para.topology.hiddenLayerSizeFF(i));
    end
end
para.output = sprintf('%s-%d.L2_%s.LR_%s/nnet', para.output, 3*(para.topology.fft_len/2+1), FormatFloat4Name(para.NET.L2weight),FormatFloat4Name(para.NET.learning_rate));
LOG = [];
para.displayTag = para.output(1:min(90, length(para.output)));

% show the configurations
para
para.NET
para.IO
para.topology

dnn = load('nnet/EnhanceRegression.noCMN.DeltaByEqn.MbSize20.U28539.771-LSTM-2048-771.L2_3E-4.LR_1E-4/nnet.itr62.LR4.97E-8.CV2439.245.mat');
layer = dnn.layer;

% train the network using SGD
trainGraph_SGD(layer, Data_tr, Data_cv, para, LOG);

author: singaxiong

repo: SignalGraph


text:
function dtiQuenchSaveFibersState(fiberGroupInds, stateFile, fg, pdbFile)
%
% dtiQuenchSaveFibersState(fiberGroupInds, stateFile, [fg], [pdbFile])
%
% Saves an array of fiber groups in QUENCH pdb/state format.
%
% HISTORY:
% 2009.09.22 ER: wrote it.
%

if(nargin>=4) %Dont save pdb if fg/pdb file name not supplied.
    mtrExportFibers(fg, pdbFile);
    if(numel(fiberGroupInds)~=numel(fg.fibers))
        error('fiberGroupInds does not match pdb file!');
    end
end

fid = fopen(stateFile,'wt');
% 5 header lines:
fprintf(fid,'Camera Position: -290.0,-21.0,76.0\n');
fprintf(fid,'Camera View Up: 0.2,0.02,1.0\n');
fprintf(fid,'Camera Focal Point: 0,-15,15\n');
fprintf(fid,'Volume Section''s Position: 40,40,20\n');
fprintf(fid,'Volume Section''s Visibility: 1,1,1\n');
fprintf(fid,'%d\n',length(fiberGroupInds));
fprintf(fid,'%d\n',fiberGroupInds);
fclose(fid);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug2040

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY nansum

[ftver, ftpath] = ft_version;

p  = tempname;
p1 = fullfile(p, 'mexfile');
p2 = fullfile(p, 'mfile');
p3 = fullfile(p, 'stats');

mkdir(p);
mkdir(fullfile(p1));
mkdir(fullfile(p2));
mkdir(fullfile(p3));

% separate the three implementations
cd(p1);
copyfile(fullfile(ftpath, 'src', ['nansum.' mexext]), '.');
cd(p2);
copyfile(fullfile(ftpath, 'src', 'nansum.m'), '.');
cd(p3);
% don't copy anything here, assume that the Mathworks signal processing toolbox is on the path

% make some data
n = 10;
x{1} = randn(n,1);
x{2} = x{1}';
x{3} = randn(n,n);
x{4} = randn(n,n,n);
x{5} = randn(n,n,n,n);

originalpath = path;
restoredefaultpath

% compute the solution for the three implementations
cd(p1);
s1 = cellfun(@nansum, x, 'UniformOutput', false);
cd(p2);
s2 = cellfun(@nansum, x, 'UniformOutput', false);
cd(p3);
s3 = cellfun(@nansum, x, 'UniformOutput', false);

% the solutions should all be equal
assert(isequal(s1, s2));
assert(isequal(s2, s3));

path(originalpath);

author: fieldtrip

repo: fieldtrip


text:
function x = p01_start ( option, nvar )

%*****************************************************************************80
%
%% P01_START returns a starting point for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  x = zeros ( nvar, 1 );

  if ( option == 1 | option == 2 | option == 3 )
    x(1:3) = [ 15.0, -2.0, 0.0 ]';
  elseif ( option == 4 | option == 5 | option == 6 )
    x(1:3) = [ 4.0, 3.0, 0.0 ]';
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P01_START - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized option  = %d\n', option );
    error ( 'P01_START - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef AlignMTB < handle
    %ALIGNMTB  Aligns images of the same scene with different exposures
    %
    % This algorithm converts images to median threshold bitmaps (1 for pixels
    % brighter than median luminance and 0 otherwise) and than aligns the
    % resulting bitmaps using bit operations.
    %
    % It is invariant to exposure, so exposure values and camera response are
    % not necessary.
    %
    % In this implementation new image regions are filled with zeros.
    %
    % For more information see [GW03].
    %
    % ## References
    % [GW03]:
    % > Greg Ward. "Fast, robust image registration for compositing high
    % > dynamic range photographs from hand-held exposures".
    % > Journal of graphics tools, 8(2):17-30, 2003.
    %
    % See also: cv.AlignMTB.AlignMTB
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % logarithm to the base 2 of maximal shift in each dimension.
        %
        % Values of 5 and 6 are usually good enough (31 and 63 pixels shift
        % respectively).
        MaxBits
        % range for exclusion bitmap that is constructed to suppress noise
        % around the median value.
        ExcludeRange
        % if true cuts images, otherwise fills the new regions with zeros.
        Cut
    end

    %% AlignMTB
    methods
        function this = AlignMTB(varargin)
            %ALIGNMTB  Creates AlignMTB object
            %
            %     obj = cv.AlignMTB()
            %     obj = cv.AlignMTB('OptionName',optionValue, ...)
            %
            % ## Options
            % * __MaxBits__ default 6
            % * __ExcludeRange__ default 4
            % * __Cut__ default true
            %
            % See also: cv.AlignMTB.process
            %
            this.id = AlignMTB_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.AlignMTB
            %
            if isempty(this.id), return; end
            AlignMTB_(this.id, 'delete');
        end

        function shift = calculateShift(this, img0, img1)
            %CALCULATESHIFT  Calculates shift between two images
            %
            %     shift = obj.calculateShift(img0, img1)
            %
            % ## Input
            % * __img0__ first image (`uint8` grayscale).
            % * __img1__ second image of same size and type as `img0`.
            %
            % ## Output
            % * __shift__ calculated shift `[x,y]`.
            %
            % Calculates shift between two images, i. e. how to shift the
            % second image to correspond it with the first.
            %
            % See also: cv.AlignMTB.shiftMat
            %
            shift = AlignMTB_(this.id, 'calculateShift', img0, img1);
        end

        function dst = shiftMat(this, src, shift)
            %SHIFTMAT  Helper function, that shift Mat filling new regions with zeros
            %
            %     dst = obj.shiftMat(src, shift)
            %
            % ## Input
            % * __src__ input image.
            % * __shift__ shift value `[x,y]`.
            %
            % ## Output
            % * __dst__ result image, same size and type as `src`.
            %
            % See also: cv.AlignMTB.calculateShift
            %
            dst = AlignMTB_(this.id, 'shiftMat', src, shift);
        end

        function [tb, eb] = computeBitmaps(this, img)
            %COMPUTEBITMAPS  Computes median threshold and exclude bitmaps of given image
            %
            %     [tb, eb] = obj.computeBitmaps(img)
            %
            % ## Input
            % * __img__ input image (`uint8` grayscale).
            %
            % ## Output
            % * __tb__ median threshold bitmap, of same size as `img` and
            %   `uint8` type.
            % * __eb__ exclude bitmap, of same size as `img` and `uint8` type.
            %
            % See also: cv.AlignMTB.process
            %
            [tb, eb] = AlignMTB_(this.id, 'computeBitmaps', img);
        end
    end

    %% AlignExposures
    methods
        function dst = process(this, src)
            %PROCESS  Aligns images
            %
            %     dst = obj.process(src)
            %
            % ## Input
            % * __src__ cell array of input images (RGB), all of the same size
            %   and `uint8` type.
            %
            % ## Output
            % * __dst__ cell array of aligned images, of same length as `src`.
            %
            % See also: cv.AlignMTB.AlignMTB
            %
            dst = AlignMTB_(this.id, 'process', src);
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.AlignMTB.empty, cv.AlignMTB.load
            %
            AlignMTB_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.AlignMTB.clear, cv.AlignMTB.load
            %
            b = AlignMTB_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.AlignMTB.save, cv.AlignMTB.load
            %
            name = AlignMTB_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.AlignMTB.load
            %
            AlignMTB_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.AlignMTB.save
            %
            AlignMTB_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

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

        function value = get.ExcludeRange(this)
            value = AlignMTB_(this.id, 'get', 'ExcludeRange');
        end
        function set.ExcludeRange(this, value)
            AlignMTB_(this.id, 'set', 'ExcludeRange', value);
        end

        function value = get.Cut(this)
            value = AlignMTB_(this.id, 'get', 'Cut');
        end
        function set.Cut(this, value)
            AlignMTB_(this.id, 'set', 'Cut', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function printmark(i, resolution, num_iterations)
    if ~mod(i, resolution)
        fprintf('.');
    end

    newline_thresh = 50*resolution;

    if ~mod(i, newline_thresh)
        fprintf(' %d\n', i);
    end

    if (i == num_iterations)
        fprintf(' %d\n', i);
    end

end

author: layumi

repo: Image-Text-Embedding


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function startup(in)
  % ocl.utils.startup(workingDirLocation)
  %
  % Startup script for OpenOCL
  % Adds required directories to the path. Sets up a folder for the results
  % of tests and a folder for autogenerated code.
  %
  % inputs:
  %   workingDirLocation - path to location where the working directory
  %                        should be created.

  ocl_dir  = fileparts(which('ocl'));

  if isempty(ocl_dir)
    ocl.utils.error('Can not find OpenOCL. Add root directory of OpenOCL to the path.')
  end

  workspaceLocation = fullfile(ocl_dir, 'Workspace');

  if nargin == 1 && ischar(in)
    workspaceLocation = in;
  elseif nargin == 1
    ocl.utils.error('Invalid argument.')
  end

  % create folders for tests and autogenerated code
  testDir     = fullfile(workspaceLocation,'test');
  exportDir   = fullfile(workspaceLocation,'export');
  [~,~] = mkdir(testDir);
  [~,~] = mkdir(exportDir);

  % set environment variables for directories
  setenv('OPENOCL_PATH', ocl_dir)
  setenv('OPENOCL_TEST', testDir)
  setenv('OPENOCL_EXPORT', exportDir)
  setenv('OPENOCL_WORK', workspaceLocation)

  % setup directories
  addpath(ocl_dir)
  addpath(exportDir)

  addpath(fullfile(ocl_dir,'doc'))

  if ~exist(fullfile(ocl_dir,'Lib','casadi'), 'dir')
    r = mkdir(fullfile(ocl_dir,'Lib','casadi'));
    ocl.utils.assert(r, 'Could not create directory in Lib/casadi');
    casadiFound = false;
  else
    % check if casadi is already installed
    addpath(fullfile(ocl_dir,'Lib'))
    addpath(fullfile(ocl_dir,'Lib','casadi'))
    casadiFound = checkCasadi(fullfile(ocl_dir,'Lib','casadi'));
  end

  % install casadi into Lib folder
  if ~casadiFound

    if ispc && verAtLeast('matlab','9.0')
      % Windows, >=Matlab 2016a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-windows-matlabR2016a-v3.4.5.zip';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    elseif ispc && verAtLeast('matlab','8.4')
      % Windows, >=Matlab 2014b
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-windows-matlabR2014b-v3.4.5.zip';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    elseif ispc && verAtLeast('matlab','8.3')
      % Windows, >=Matlab 2014a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-windows-matlabR2014a-v3.4.5.zip';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    elseif ispc && verAtLeast('matlab','8.1')
      % Windows, >=Matlab 2013a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-windows-matlabR2013a-v3.4.5.zip';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));

    elseif isunix && ~ismac && verAtLeast('matlab','8.4')
      % Linux, >=Matlab 2014b
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-linux-matlabR2014b-v3.4.5.tar.gz';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    elseif isunix && ~ismac && verAtLeast('matlab','8.3')
      % Linux, >=Matlab 2014a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-linux-matlabR2014a-v3.4.5.tar.gz';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));

    elseif ismac && verAtLeast('matlab','8.5')
      % Mac, >=Matlab 2015a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-osx-matlabR2015a-v3.4.5.tar.gz';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    elseif ismac && verAtLeast('matlab','8.4')
      % Mac, >=Matlab 2015a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-osx-matlabR2014b-v3.4.5.tar.gz';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    elseif ismac && verAtLeast('matlab','8.3')
      % Mac, >=Matlab 2015a
      path = 'https://github.com/casadi/casadi/releases/download/3.4.5/';
      filename = 'casadi-osx-matlabR2014a-v3.4.5.tar.gz';
      downloadCasadi(ocl_dir, path, filename, fullfile(ocl_dir,'Lib','casadi'));
    else
      ocl.utils.info(['Could not set up CasADi for you system.', ...
               'You need to install CasADi yourself and add it to your path.'])
    end
    % add Lib and Lib/casadi to path
    addpath(fullfile(ocl_dir,'Lib'))
    addpath(fullfile(ocl_dir,'Lib','casadi'))
  end

  casadiFound = checkCasadiWorking();
  if casadiFound
    ocl.utils.info('CasADi is up and running!')
  else
    ocl.utils.error('Go to https://web.casadi.org/get/ and setup CasADi.');
  end

  ocl.utils.info('OpenOCL startup procedure finished successfully.')

end

function downloadCasadi(ocl_dir, path, filename, dest)

  if exist(fullfile(dest, 'CUSTOM_CASADI'), 'file') > 0
    fprintf(['You chose to your use your custom CasADi installation. If you changed \n', ...
            'your mind delete the CUSTOM_CASADI file in %s\n'], dest);
    return;
  end

  fprintf(2,'\nYour input is required! Please read below:\n')

  confirmation = [ '\n', 'Dear User, if you continue, CasADi will be downloaded from \n', ...
                   path, filename, ' \n', ...
                   'and saved to the Workspace folder. The archive will be extracted \n', ...
                   'to the Lib folder. This will take a few minutes. \n\n', ...
                   'Hit [enter] to continue! \n\n', ...
                   'Advanced users: \n', ...
                   'If you have set-up your own CasADi version and you want to use that, you \n', ...
                   'can type `n` and hit [enter]. We will then perform some basic checks if \n', ...
                   'your version is campatible. We strongly recommend you to let OpenOCL install \n', ...
                   'the required version (if you do not save the path, the CasADi version of OpenOCL \n', ...
                   'will not be on your path at startup, and not conflict with you current CasADi \n', ...
                   'installation): '];
  m = input(confirmation,'s');

  if strcmp(m, 'n') || strcmp(m, 'no')
    try
      ocl.tests.variable;
      fid = fopen(fullfile(dest, 'CUSTOM_CASADI'),'w');
      fclose(fid);
      return;
    catch e
      warning(e.message)
      ocl.utils.error('You did not agree to download CasADi and your version is not compatible. Either run again or set-up a compatible CasADi version manually.');
    end
  end

  archive_destination = fullfile(ocl_dir, 'Workspace', filename);

  if ~exist(archive_destination, 'file')
    ocl.utils.info('Downloading...')
    websave(archive_destination, [path,filename]);
  end
  ocl.utils.info('Extracting...')
  [~,~,ending] = fileparts(archive_destination);
  if strcmp(ending, '.zip')
    unzip(archive_destination, dest)
  else
    untar(archive_destination, dest)
  end
end

function r = checkCasadi(path)

  cur_path = pwd;
  cd(path)

  if ~exist(fullfile(path,'+casadi','SX.m'),'file') > 0
    r = false;
  else
    try
      casadi.SX.sym('x');
      r = true;
    catch e
      cd(cur_path)
      ocl.utils.warning(e.message);
      casadiNotWorkingError();
    end
  end
  cd(cur_path)
end

function r = checkCasadiWorking()
  try
    casadi.SX.sym('x');
    r = true;
  catch e
    ocl.utils.warning(e.message);
    casadiNotWorkingError();
  end
end

function r = verAtLeast(software, version_number)
  r = ~verLessThan(software,version_number);
end

function casadiNotWorkingError
  ocl.utils.error(['Casadi installation not found or it does not ', ...
            'work properly. Try restarting Matlab. Remove all ', ...
            'casadi installations from your path. Run ocl.utils.clean. OpenOCL will ', ...
            'then install the correct casadi version for you.']);
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL