text:
function [ordered_alpha,index_alpha]=order_subsets(angles,blocksize, mode)
index_alpha=1:size(angles,2);

block_alpha=mat2cell(angles     ,size(angles,1),[repmat(blocksize,1,floor(size(angles,2)/blocksize)) mod(size(angles,2),blocksize)]);
index_alpha=mat2cell(index_alpha,1,[repmat(blocksize,1,floor(size(angles,2)/blocksize)) mod(size(angles,2),blocksize)]);

block_alpha=block_alpha(~cellfun('isempty',block_alpha));
index_alpha=index_alpha(~cellfun('isempty',index_alpha));

if strcmp(mode,'ordered')
    ordered_alpha=block_alpha;
    return;
end
if strcmp(mode,'random')
    neworder=randperm(length(block_alpha));
    ordered_alpha=block_alpha(neworder);
    index_alpha=index_alpha(neworder);
    return;
end
%% not finished

if strcmp(mode,'angularDistance')
    % we need them sorted, so we need to recompute the blocks, but sorted
    [angles,sortindex] = sortrows(angles.',1);
    angles=angles.';sortindex=sortindex.';
    index_alpha=1:size(angles,2);
    index_alpha=index_alpha(sortindex);

    block_alpha=mat2cell(angles      ,size(angles,1),[repmat(blocksize,1,floor(size(angles,2)/blocksize)) mod(size(angles,2),blocksize)]);
    index_alpha=mat2cell(index_alpha,1,[repmat(blocksize,1,floor(size(angles,2)/blocksize)) mod(size(angles,2),blocksize)]);

    block_alpha=block_alpha(~cellfun('isempty',block_alpha));
    index_alpha=index_alpha(~cellfun('isempty',index_alpha));

    avrg=cellfun(@mean,block_alpha);
    used_avrg=[];
    % start from the beggining
    ordered_alpha{1}=block_alpha{1};
    auxindex_alpha=index_alpha;
    index_alpha{1}=auxindex_alpha{1};
    used_avrg(end+1)=avrg(1);
    for ii=2:length(block_alpha)
        dist=[];
        for jj=1:length(used_avrg)
            dist(jj,:)=abs(mod((avrg- used_avrg(jj))+pi,2*pi)-pi);
        end
        dist=bsxfun(@times,dist,all(dist,1));
        [~,midx]=max(dist(:));
        [~,avrgindx]=ind2sub(size(dist),midx);
        index_alpha{ii}=auxindex_alpha{avrgindx};
        ordered_alpha{ii}=block_alpha{avrgindx};
        used_avrg(end+1)=avrg(avrgindx);
    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function [ERP erpcom] = pop_getChanInfoFromEeglab(ERP)
ERPaux = ERP;
ERP    = pop_clearerpchanloc( ERP, 'Warning', 'off');
ERP    = borrowchanloc(ERP);
erpcom = sprintf('%s = getChanInfoFromEeglab(%s)', inputname(1), inputname(1));
[ERP, issave, erpcom_save] = pop_savemyerp(ERP,'gui','erplab', 'History', 'implicit');
if issave>0
        % generate text command
        if issave==2
                erpcom = sprintf('%s\n%s', erpcom, erpcom_save);
                msgwrng = '*** Your ERPset was saved on your hard drive.***';
        else
                msgwrng = '*** Warning: Your ERPset was only saved on the workspace.***';
        end
        fprintf('\n%s\n\n', msgwrng)
        try cprintf([0 0 1], 'COMPLETE\n\n');catch,fprintf('COMPLETE\n\n');end ;
else
        ERP     = ERPaux;
        msgwrng = 'ERPLAB Warning: Your changes were not saved';
        try cprintf([1 0.52 0.2], '%s\n\n', msgwrng);catch,fprintf('%s\n\n', msgwrng);end ;
end
displayEquiComERP(erpcom);
%
% Completion statement
%
msg2end
return

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [infostring warnstring errstring] = est_checkMVARParams(EEG, params)
%
% This function performs a series of sanity checks on the parameters chosen
% for fitting a VAR model. The results are returned in cell vectors
% indicating information bulletins (notifications), warnings, and/or
% errors.
%
% Inputs:
%
%   EEG:            EEG data structure
%   params:         a parameter structure as output from est_fitMVAR
%
% Outputs:
%
%   infostring:     cell array of information bulletins (one per cell)
%   warnstring:     cell array of warning bulletins (one per cell)
%   errstring:      cell array of error bulletins (one per cell)
%
% References:
%
% [1] Korzeniewska, et al (2008). Dynamics of Event-Related Causality in
%     Brain Electrical Activity. Human brain mapping 29:1170?1192
% [2] Schlogl and Supp, (2006). Analyzing event-related EEG data with
%     multivariate AR parameters. Prog. in Brain Researc. Neuper and Klimesh, Eds.
% [3] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%     Theoretical Handbook and User Manual.
%     Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% See Also: pop_est_fitMVAR()
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD
% Email:  tim@sccn.ucsd.edu

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

srate = EEG.srate;
wlen = srate*params.winlen;  % winlen in samples
T = EEG.CAT.trials;
p = params.morder;
M = EEG.CAT.nbchan;

% addtional error checks and info
infostring = {};
warnstring = {};
errstring  = {};

% critical error checks
if params.winstep==0
    error('Step size must be greater than 0\n');
end

if floor(params.winlen*1000) > floor((EEG.xmax-EEG.xmin+1)*1000)
    error('Winlen cannot be greater than the trial length (%0.1f sec)\n',EEG.xmax-EEG.xmin);
end

if round(params.winlen*EEG.srate) <= max(params.morder)
    error('The window length must be greater than the model order. Increase your window to at least %0.2f sec\n',(params.morder(end)+1)/EEG.srate);
end

if length(p)==2
    infostring = [infostring sprintf('Two model orders specified [%d %d]\n',p(1),p(2))];
    warnstring = [warnstring sprintf('\tI assume you are providing a [min max] range for model order selection.\n\tI will use p=(%d) for the remaining checks...\n',p(2))];
    errstring  = [errstring {''}];
    p = p(2);
end

% Check if we have a reasonable ratio of datapoints to free parameters
% Optimal ratio derived from [1]
winrat = ((M^2)*p)/(wlen*T);
infostring = [infostring sprintf('Ratio of number of parameters to datapoints is %0.3f.\n',winrat)];
if winrat > 1
    reclen = ((M^2)*p/T)/srate;
    warnstring = [warnstring sprintf('\tIf using an unconstrained (unregularized) model fitting apprach: The ratio of number of parameters to datapoints must be <= 1.\n\tYour window length must be at least %0.3f sec\n',reclen)];
    errstring  = [errstring {''}];
elseif winrat > 0.1
    reclen = (10*(M^2)*p/T)/srate;
    warnstring = [warnstring sprintf('\tIf using an unconstrained (unregularized) model fitting apprach: For best results, ratio of number of parameters to datapoints should be < 0.1.\n\tI recommend using window length of at least %0.3f sec\n',reclen)];
    errstring  = [errstring {''}];
else
    warnstring = [warnstring {''}];
    errstring  = [errstring {''}];
end

% check if time-frequency principle is violated [2]
TFprod = wlen*sqrt(T);
infostring = [infostring sprintf('Time-Frequency Product is %0.3f. This should be greater than p=%d\n',TFprod,p)];
if wlen<=p/sqrt(T)
    reclen = p/sqrt(T);
    warnstring = [warnstring sprintf('\tTime window does not satisfy the Time-Frequency Uncertainty Principle.\n\tI recommend a window length of at least %0.3f sec\n', reclen)];
    errstring  = [errstring {''}];
else
    warnstring = [warnstring {''}];
    errstring  = [errstring {''}];
end

% notify user of max number of spectral peaks
infostring = [infostring sprintf('Given your model order of p=%d, a maximum of p/2=%0.1f frequency components (spectral peaks) can be estimated for each pair of variables\n',p, p/2)];
warnstring = [warnstring {''}];
errstring  = [errstring {''}];

repo: REST


text:
classdef StiffnessPlaneStressVoigtTensor < SymmetricFourthOrderPlaneStressVoigtTensor ...
                                           & StiffnessDescriptor

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function M = euclideansubspacefactory(E, proj, dim)
% Returns a manifold struct to optimize over a subspace of a linear manifold
%
% The factory produces a manifold description M for a linear subspace S of
% a linear space E.
%
% Points and tangent vectors on M are represented in memory in exactly the
% same way as they are for E. Furthermore, E is considered the embedding
% space of S. In particular, the Euclidean gradient of a function is
% understood as the gradient of that function in E, whereas the gradient of
% that function in S is here called the Riemannian gradient. Similar
% considerations hold for the Hessian.
%
% The linear space E is described by a structure obtained from a factory.
% For example, E can be the space of vectors, matrices or nD-arrays of real
% or complex numbers, or a subspace of those. E is equipped with a (real)
% inner product, making it a (real) Euclidean space.
%
% The subspace S is described by an orthogonal projector proj from E to E,
% whose image (i.e., span, range, ...) is S. Orthogonality is judged with
% respect to the inner product on E, and S inherits that inner product,
% making it a Euclidean space itself, and a Riemannian submanifold of E.
%
% The dimension of the subspace S must ideally also be specified.
%
% Inputs:
%   E is the factory for a linear space; for example:
%       E = euclideanfactory(n)
%       E = euclideanfactory(n, m)
%       E = euclideanfactory([n1, n2, n3])
%       E = euclideancomplexfactory(n) % and other dimensions too
%       ...
%
%   proj is a function handle that takes as input a point u of E and
%       returns another point of E: the orthogonal projection of u to S.
%       Orthogonality is understood with respect to the inner product on E,
%       that is, E.inner(x, u, v) = <u, v>. Like any orthogonal projector,
%       it must be linear, and it must satisfy:
%           <u, Pv> = <Pu, v>    and    PPw = Pw
%       for all points u, v, w in E.
%       To check these properties on random vectors, call M.checkproj();
%
%   dim is the dimension of the subspace S (an integer). As always in
%       Manopt, we consider linear spaces over the real numbers, hence,
%       this is the dimension of S as a real linear space. For example, the
%       dimension of R^n is n, and the dimension of C^n is 2n.
%
% Output:
%   M is a factory for optimization over the subspace S.
%
%
% Example:
%
% Among complex vector of length n, consider the subspace S of such vectors
% that can be the discrete Fourier transform of a real vector of length n
% (that is, the subspace generated by feeding all real vectors of length n
% to fft). The factory M below describes that subspace.
%
%   n = 17;
%   E = euclideancomplexfactory(n);
%   proj = @(u) (u + conj(u([1 ; (n:-1:2)'])))/2;
%   M = euclideansubspacefactory(E, proj, n);
%   M.checkproj(); % for debugging
%
%
% Note 1: this factory is designed to work with linear subspaces only: it
% does not work for affine subspaces. Explicitly: S must contain the origin
% of E. If you need support for affine subspaces, let us know on the Manopt
% forum and we can help (or share your improved code :)).
%
% Note 2: the linear space E can itself be a linear subspace of another.
% For example, we can have:
%       E = symmetricfactory(n, k)
%       E = skewfactory(n, k)
%       E = euclideansubspacefactory(...) % recursive nesting
% For these use-cases, bear in mind that the embedding space of M is E,
% hence, the Euclidean gradient and Hessian are expected to be given in E,
% not in the 'bigger' linear space E possibly lives in.
%
% See also: euclideanfactory euclideancomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 25, 2019.
% Contributors:
% Change log:

    M = E;

    M.name = @() ['Subspace of ', E.name()];

    M.proj = @(x, u) proj(u);

    M.egrad2rgrad = @(x, eg) proj(eg);

    M.ehess2rhess = @(x, eg, eh, d) proj(eh);

    M.tangent = @(x, u) proj(u);

    M.rand = @() proj(E.rand());

    M.randvec = @randvec;
    function v = randvec(x)
        v = proj(E.randvec(x));
        v = v / E.norm(x, v);
    end

    if exist('dim', 'var')
        M.dim = @() dim;
        M.typicaldist = @() sqrt(dim);
    else
        M = rmfield(M, 'dim');
        warning('manopt:subspacedim', ...
               ['Since the dimension of the subspace was not specified' ...
                ', M.dim() is not available in the returned factory.']);
    end

    M.checkproj = @() checkproj(E, proj);

end

% Tool to check that proj is indeed an orthogonal projector to a subspace
% of E, with respect to the inner product on the linear space E.
function checkproj(E, proj)
    x = E.rand();
    u = E.randvec(x);
    v = E.randvec(x);
    a = randn();
    b = randn();
    % Check that P is linear.
    Paubv = proj(E.lincomb(x, a, u, b, v));
    aPvbPv = E.lincomb(x, a, proj(u), b, proj(v));
    fprintf(['Is proj linear? This number should be zero up to ' ...
             'machine precision:\n\t%.16e\n'], ...
                E.norm(x, E.lincomb(x, 1, Paubv, -1, aPvbPv)));
    % Check that PPw = Pw.
    fprintf(['Is it a projector? This number should be zero up to ' ...
             'machine precision:\n\t%.16e\n'], ...
                E.norm(x, E.lincomb(x, 1, proj(u), -1, proj(proj(u)))));
    % Check that <u, Pv> = <Pu, v>.
    fprintf(['Is it self-adjoint? These two numbers should be equal ' ...
             'up to machine precision:\n\t%.16e\n\t%.16e\n'], ...
                E.inner(x, u, proj(v)), ...
                E.inner(x, proj(u), v));
end

author: NicolasBoumal

repo: manopt


text:
function varargout = quiver(N, varargin)
% QUIVER   Draw phase plot diagrams, based on ODEs specified with CHEBOPS
%
% Calling sequence:
%   H = QUIVER(N, AXIS, 'OPT1', VAL1, ...)
%
% Here, the inputs are:
%   N    : A chebop, whose N.op arguments specifies a first or second order
%          scalar ODE, or a coupled system of two first order ODEs.
%   AXIS : A 4-vector with elements [XMIN XMAX YMIN YMAX] that specify the
%          rectangular region shown on the phase plot. If none is passed, the
%          default values [-1 1 -1 1] are used.
%
% It is possible to pass the method various option pairs of the form
% 'OPTIONNAME', OPTIONVALUE. The options supported are:
%   'XPTS'      : An integer, specifying the resolution of the x-axis for the
%                 quiver plot. Default value: 10.
%   'YPTS'      : An integer, specifying the resolution of the y-axis for the
%                 quiver plot. Default value: 10.
%   'NORMALIZE' : A Boolean, which determines whether the arrows on the quiver
%                 plot or normalized all to have the same length. Default value:
%                 false.
%   'SCALE':      By default, quiver automatically scales the arrows to fit
%                 within the grid. By passing a SCALE argument S, the arrows are
%                 fitted within the grid and then stretched by S. Use S = 0 to
%                 plot the arrows without the automatic scaling. See the
%                 documentation for the built-in MATLAB QUIVER for more
%                 information.
%   LINESPEC      Specifies options for the quiver plot. The LINESPEC argument
%                 can either be a string supported by the MATLAB PLOT command,
%                 such as 'ro', or a parameter/value pair to specify additional
%                 properties of the quiver lines, such as
%                   quiver(N,[0 2 0 2], 'k', 'linewidth', 1.4).
%
% The optional output is
%   H   : A quivergroup handle.
%
% Note: The CHEBOP QUIVER command works by reformulating higher order problems
% as coupled first order systems, evaluating the resulting first order system at
% grid that should be interpreted as values of u and u', then calling the
% built-ing MATLAB QUIVER method on the results. In the case of first order
% scalar problems, the grid should be interpreted as values of t and u.
%
% Example 1 -- van der Pol equation (second order ODE)
%   N = chebop(0,100);
%   N.op = @(t,u) diff(u, 2) - 3*(1-u^2)*diff(u) + u;
%   quiver(N,[-2.5 2.5 -5.5 5.5], 'r', 'xpts', 40, 'ypts', 40, 'scale', .5, ...
%       'normalize', true, 'linewidth', 1.5)
%   hold on % Plot a particular solution on top of quiver plot
%   N.lbc = [2; 0];
%   u = N\0;
%   plot(u, diff(u),'linewidth',2)
%
% Example 2 -- Lotka-Volterra (first order coupled system)
%   N = chebop(@(t,u,v) [diff(u)-2.*u+u.*v; diff(v)+v-u.*v], [0 4]);
%   quiver(N, [0 2.5 0 4])
%   hold on
%   N.lbc = @(u,v) [u - 0.5; v - 1]; % Initial populations
%   [u, v] = N\0;
%   plot(u, v, 'linewidth', 2)
%   plot(0.5, 1,'m*','markersize',15) % Mark initial condition
%
% Example 3 -- Slopefield for a first order problem
%   N = chebop(@(t,u) diff(u)-sin(t)*u);
%   quiver(N,[-1.2*pi 1.2*pi -1 1])

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

% Set default values:
scale = 1;
normalize = false;
u0 = [];
xpts = 20;
ypts = 20;
axisLims = [-1 1 -1 1];

% The limits for the axes have to appear in the first entry of varargin if
% they're specified. They'll always be a 4-vector.
if ( ~isempty(varargin) && length(varargin{1}) == 4 )
    axisLims = varargin{1};
    % Throw away the axisLims:
    varargin(1) = [];
end

% Cell for linespecs
linespec = {};

% Parse VARARGIN, go through all elements:
while ( ~isempty(varargin) )
    if ( ~ischar(varargin{1}) && ~isnumeric(varargin{2}) )
        error('followpath:inputArgument','Incorrect options input arguments');
    end

    throwAwayLength = 2;
    switch lower(varargin{1})
        case 'normalize'
            normalize = varargin{2};
        case 'xpts'
            xpts = varargin{2};
        case 'ypts'
            ypts = varargin{2};
        case 'scale'
            scale = varargin{2};
        otherwise
            % Must have gotten something to do with linespec
            findNorm = strcmpi(varargin, 'normalize');
            findXpts = strcmpi(varargin, 'xpts');
            findYpts = strcmpi(varargin, 'ypts');
            findScale = strcmpi(varargin, 'scale');
            optionPos = findNorm | findXpts | findYpts | findScale;
            nextOption = find(optionPos, 1, 'first');
            % If nextOption is empty, we only got linespec options left
            if ( isempty(nextOption) )
                nextOption = length(varargin) + 1;
            end
            linespec = [linespec, varargin{1:nextOption - 1}];
            throwAwayLength = nextOption - 1;
    end

    % Throw away option name and argument and move on:
    varargin(1:throwAwayLength) = [];
end

% Extract the x and y limits:
xl = axisLims(1:2);
yl = axisLims(3:4);

% If ylim is empty, we solve the problem to obtain a range for plotting on:
if ( isempty(xl) )
    u0 = N\0;
    xl = 1.1*minandmax(u0);
end

if ( isempty(yl) )
    if ( isempty(u0) )
        u0 = N\0;
    end
    yl = 1.1*minandmax(diff(u0));
end

% Convert the operator in N to first order.
[firstOrderFun, ~, ~, ~, diffOrd] = treeVar.toFirstOrder(N.op, 0, N.domain);

% Vectors for constructing a meshgrid:
y1 = linspace(xl(1), xl(end), xpts);
y2 = linspace(yl(1), yl(end), ypts);

% Get a meshgrid for points of interests in the phase plane.
[x,y] = meshgrid(y1, y2);
u = zeros(size(x));
v = zeros(size(x));

% Phase plane portraits really only make sense for autonomous systems, which
% shouldn't depend on t, hence, we simply take t = 0 for evaluating.
t = 0;

% Check if we got passed a system of too high order:
numEquations = length(strfind(func2str(firstOrderFun),';'));

if ( numEquations > 2 )
    error('CHEBFUN:CHEBOP:quiver:tooHighOrder', ...
        ['The ODE passed to chebop/quiver must either be a scalar, second ' ...
        'order ODE or a system of two first order equations.'])
end

% Are we plotting a slope field (cf. #2238), or a standard phase plane?
if( (length(diffOrd) == 1) && (diffOrd == 1))   % Slope field
    for i = 1:numel(x)
        res = firstOrderFun(x(i), y(i));
        u(i) = 1;
        v(i) = res(1);
    end
else                                            % Phase plane
    for i = 1:numel(x)
        res = firstOrderFun(t,[x(i); y(i)]);
        u(i) = res(1);
        v(i) = res(2);
    end
end

if ( normalize )
    % Make all arrows equal length:
    nrm = sqrt(u.^2 + v.^2);
    u = u./nrm;
    v = v./nrm;
end

h = quiver(x, y, u, v, scale,linespec{:});

% Set x and y limits:
xlim(xl);
ylim(yl);
if ( nargout > 0 )
    varargout{1} = h;
else

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% (this is just an interface)
function varargout = mydatestd_inv (varargin)
    [varargout{1:nargout}] = mydatestdi (varargin{:});
end

%!test
%! % just check whether it runs:
%! mydatestd_inv(1, 2);

author: Sable

repo: mcbench-benchmarks


text:
function TF = antelope_exists()

%ANTELOPE_EXISTS check for Antelope Toolbox for Matlab
%  TF = ANTELOPE_EXISTS returns 1 if the Antelope Toolbox for Matlab is
%  present in the Matlab path. This function does not check to see whether
%  the Antelope Toolbox is functioning correctly. It merely searches the
%  path for the existence of certain functions. To test the actual validity
%  of the Antelope toolbox you might try the example suite executed by
%  dbexample_runall (an example in the Antelope Toolbox for Matlab). The
%  Antelope Toolbox for Matlab is not part of GISMO though GISMO includes
%  many utilities built on the toolbox. It is maintained by Lindquist
%  Consulting, Inc.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if exist('dbopen','file') && exist('trload_css','file');
    TF = 1;
else
    TF = 0;
end

author: geoscience-community-codes

repo: GISMO


text:
function my_clear_net_trn(net_config)

    cnn_update_net_group(net_config, @do_clear, []);

end

function one_group_info=do_clear(one_group_info, opts)

net_info=one_group_info.net_info;
if isempty(net_info)
    return;
end

for i=1:numel(net_info.ref.layers)

  if strcmp(net_info.ref.layers{i}.type,'conv')

      net_info.ref.layers{i}.filtersMomentum = zeros('like',net_info.ref.layers{i}.filters) ;
      net_info.ref.layers{i}.biasesMomentum = zeros('like',net_info.ref.layers{i}.biases) ;

      if ~isfield(net_info.ref.layers{i}, 'filtersLearningRate')
        net_info.ref.layers{i}.filtersLearningRate = 1 ;
      end
      if ~isfield(net_info.ref.layers{i}, 'biasesLearningRate')
        net_info.ref.layers{i}.biasesLearningRate = 1 ;
      end
      if ~isfield(net_info.ref.layers{i}, 'filtersWeightDecay')
        net_info.ref.layers{i}.filtersWeightDecay = 1 ;
      end
      if ~isfield(net_info.ref.layers{i}, 'biasesWeightDecay')
        net_info.ref.layers{i}.biasesWeightDecay = 1 ;
      end

  end

  if strcmp(net_info.ref.layers{i}.type,'my_custom')
      layer_update_fn=net_info.ref.layers{i}.layer_update_fn;
      if ~isempty(layer_update_fn)
          update_info=[];
          update_info.group_info=one_group_info;
          net_info.ref.layers{i}=net_info.ref.layers{i}.layer_update_fn(...
                'init_trn', net_info.ref.layers{i}, update_info);
      end
  end

end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function tissueModel = INIT(model, weights, tol, runtime, logfile, epsilon)
% Use the INIT algorithm (`Agren et al., 2012`) to extract a context
% specific model using data. INIT algorithm find the optimal trade-off
% between inluding and removing reactions based on their given weights. If
% desired, accumulation of certain metabolites can be allowed or even
% forced.
%
% USAGE:
%
%    tissueModel = INIT(model, weights, tol, runtime, logfile)
%
% INPUTS:
%    model:               input model (COBRA model structure)
%    weights:             column with positive and negative weights for each reaction
%                         positive weights are reactions with high expression, negative
%                         weigths for reaction with low expression (must be same length
%                         as model.rxns)
%
% OPTIONAL INPUTS:
%    tol:                 minimum flux threshold for "expressed" reactions
%                         (default 1e-8)
%    logfile:             name of the file to save the MILP log (string)
%    runtime:             maximum solve time for the MILP (default value - 7200s)
%    epsilon:             value added/subtracted to upper/lower bounds
%                         (default 1)
%
% OUTPUTS:
%    tissueModel:         extracted model
%
% `Agren et al. (2012). Reconstruction of genome-scale active metabolic
% networks for 69 human cell types and 16 cancer types using INIT. PLoS
% Comput. Biol. 8, e1002518.`
%
% .. Authors:  - Implementation adapted from the cobra toolbox (createTissueSpecificModel.m) by S. Opdam and A. Richelle, May 2017

if isfield(model,'C') || isfield(model,'E')
    issueConfirmationWarning('INIT does not handle the additional constraints and variables defined in the model structure (fields .C and .E.)\n It will only use the stoichiometry provided.');
end

if nargin < 6 || isempty(epsilon)
    epsilon=1;
end
if nargin < 5 || isempty(runtime)
    runtime = 7200;
end
if nargin < 4 || isempty(logfile)
    logfile = 'MILPlog';
end
if nargin < 3 || isempty(tol)
    tol = 1e-8;
end

    RHindex = find(weights > 0);
    RLindex = find(weights < 0);

    %Weights of 0 will be handled the same as in iMAT

    S = model.S;
    lb = model.lb;
    ub = model.ub;

    % Creating A matrix
    A = sparse(size(S,1)+2*length(RHindex)+2*length(RLindex),size(S,2)+2*length(RHindex)+length(RLindex));
    [m,n,s] = find(S);
    for i = 1:length(m)
        A(m(i),n(i)) = s(i);
    end

    for i = 1:length(RHindex)
        A(i+size(S,1),RHindex(i)) = 1;
        A(i+size(S,1),i+size(S,2)) = lb(RHindex(i)) - epsilon;
        A(i+size(S,1)+length(RHindex),RHindex(i)) = 1;
        A(i+size(S,1)+length(RHindex),i+size(S,2)+length(RHindex)+length(RLindex)) = ub(RHindex(i)) + epsilon;
    end

    for i = 1:length(RLindex)
        A(i+size(S,1)+2*length(RHindex),RLindex(i)) = 1;
        A(i+size(S,1)+2*length(RHindex),i+size(S,2)+length(RHindex)) = lb(RLindex(i));
        A(i+size(S,1)+2*length(RHindex)+length(RLindex),RLindex(i)) = 1;
        A(i+size(S,1)+2*length(RHindex)+length(RLindex),i+size(S,2)+length(RHindex)) = ub(RLindex(i));
    end

    % Creating csense
    csense1(1:size(S,1)) = 'E';
    csense2(1:length(RHindex)) = 'G';
    csense3(1:length(RHindex)) = 'L';
    csense4(1:length(RLindex)) = 'G';
    csense5(1:length(RLindex)) = 'L';
    csense = [csense1 csense2 csense3 csense4 csense5];

    % Creating lb and ub
    lb_y = zeros(2*length(RHindex)+length(RLindex),1);
    ub_y = ones(2*length(RHindex)+length(RLindex),1);
    lb = [lb;lb_y];
    ub = [ub;ub_y];

    % Creating c
    c_v = zeros(size(S,2),1);
    c_y = ones(2*length(RHindex)+length(RLindex),1);
    c_w = [weights(RHindex);weights(RHindex);abs(weights(RLindex))];
    c = [c_v;c_w.*c_y];

    % Creating b
    b_s = zeros(size(S,1),1);
    lb_rh = lb(RHindex);
    ub_rh = ub(RHindex);
    lb_rl = lb(RLindex);
    ub_rl = ub(RLindex);
    b = [b_s;lb_rh;ub_rh;lb_rl;ub_rl];

    % Creating vartype
    vartype1(1:size(S,2),1) = 'C';
    vartype2(1:2*length(RHindex)+length(RLindex),1) = 'B';
    vartype = [vartype1;vartype2];

    MILPproblem.A = A;
    MILPproblem.b = b;
    MILPproblem.c = c;
    MILPproblem.lb = lb;
    MILPproblem.ub = ub;
    MILPproblem.csense = csense;
    MILPproblem.vartype = vartype;
    MILPproblem.osense = -1;
    MILPproblem.x0 = [];

    solution = solveCobraMILP(MILPproblem, 'timeLimit', runtime, 'logFile', logfile, 'printLevel', 3);

    x = solution.cont;
    rxnRemList = model.rxns(abs(x) < tol);
    tissueModel = removeRxns(model,rxnRemList);
end

author: opencobra

repo: cobratoolbox


text:
function x = real( x )

%   Disciplined convex/geometric programming information for REAL:
%      REAL(X) may be freely applied to any CVX expression. However,
%      since REAL(X)=X for all real expressions (including convex,
%      concave, log-convex, and log-concave), it is only useful in
%      the complex affine case.

b = x.basis_;
if ~isreal( b ),
    x = cvx( x.size_, real( b ) );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%DEMOIMSHAPES2D Display various 2D demo shapes
%
%   output = demoShapes2d(input)
%
%   Example
%   demoImShapes
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% initialisations

% generate cubic images
lx = 1:100;
ly = 1:100;

% choose a center not aligned with the grid
center = [50+sqrt(2)-1 50+sqrt(3)-1];

angles = [...
    0; ...
    10; ...
    30; ...
    45];

%% Disc

% sphere is defined by center and radius
disc = [center 40];

% generation of 3D image
img = discreteDisc(lx, ly, disc);

% display image
f = figure;
imshow(~img);

% decorate
title('Disc, radius=40');

%% Ellipse

figure;
for i = 1:4
    % sphere is defined by center and radius
    ellipse = [center 40 20 angles(i)];

    % generation of 3D image
    img = discreteEllipse(lx, ly, ellipse);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Ellipse, Th=%02d', angles(i)));
end

%% Square

figure;
for i = 1:4
    % sphere is defined by center and radius
    square = [center 60 angles(i)];

    % generation of 3D image
    img = discreteSquare(lx, ly, square);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Square, Th=%02d', angles(i)));
end

%% Rectangle

figure;
for i = 1:4
    % sphere is defined by center and radius
    rect = [center 60 30 angles(i)];

    % generation of 3D image
    img = discreteRectangle(lx, ly, rect);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Rectangle, Th=%02d', angles(i)));
end

%% Capsule

figure;
for i = 1:4
    % sphere is defined by center and radius
    [dx dy] = pol2cart(deg2rad(angles(i)), 30);
    caps = [center-[dx dy] center+[dx dy] 15];

    % generation of 3D image
    img = discreteCapsule(lx, ly, caps);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Capsule, Th=%02d', angles(i)));
end

%% Egg

figure;
for i = 1:4
    % egg defined by a center, a size and an angle
    egg = [center 30 angles(i)];

    % generation of 3D image
    img = discreteEgg(lx, ly, egg);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Egg, Th=%02d', angles(i)));
end

%% Trefoil

figure;
for i = 1:4
    % trefoil defined by a center, two radii and an angle
    trefoil = [center 40 15 angles(i)];

    % generation of 3D image
    img = discreteTrefoil(lx, ly, trefoil);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Trefoil, Th=%02d', angles(i)));
end

%% Starfish

figure;
for i = 1:4
    % starfish defined by a center, two radii and an angle
    starfish = [center 40 20 angles(i)];

    % generation of 3D image
    img = discreteStarfish(lx, ly, starfish);

    % display image
    subplot(2, 2, i);
    imshow(~img);

    % decorate
    title(sprintf('Starfish, Th=%02d', angles(i)));
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function hh = quiver5(varargin)
% Slightly modified version of quiver3 function that plot arrows with
% true 3D arrow heads
% Bertrand Dano 05-12-08
%QUIVER5 3-D quiver plot.
%   QUIVER5(X,Y,Z,U,V,W) plots velocity vectors as arrows with components
%   (u,v,w) at the points (x,y,z).  The matrices X,Y,Z,U,V,W must all be
%   the same size and contain the corresponding position and velocity
%   components.  QUIVER3 automatically scales the arrows to fit.
%
%   QUIVER5(Z,U,V,W) plots velocity vectors at the equally spaced
%   surface points specified by the matrix Z.
%
%   QUIVER5(Z,U,V,W,S) or QUIVER3(X,Y,Z,U,V,W,S) automatically
%   scales the arrows to fit and then stretches them by S.
%   Use S=0 to plot the arrows without the automatic scaling.
%
%   QUIVER5(...,LINESPEC) uses the plot linestyle specified for
%   the velocity vectors.  Any marker in LINESPEC is drawn at the base
%   instead of an arrow on the tip.  Use a marker of '.' to specify
%   no marker at all.  See PLOT for other possibilities.
%
%   QUIVER5(...,'filled') fills any markers specified.
%
%   H = QUIVER3(...) returns a vector of line handles.
%
%   Example:
%       [x,y] = meshgrid(-2:.2:2,-1:.15:1);
%       z = x .* exp(-x.^2 - y.^2);
%       [u,v,w] = surfnorm(x,y,z);
%       quiver5(x,y,z,u,v,w);
%       axis vis3d; rotate3d on
%
%   See also QUIVER, PLOT, PLOT3, SCATTER.
%   Clay M. Thompson 3-3-94
%   Copyright 1984-2002 The MathWorks, Inc.
%   $Revision: 1.23 $  $Date: 2002/06/05 20:05:16 $
% Arrow head parameters
alpha = 0.33; % Size of arrow head relative to the length of the vector
beta = 0.33;  % Width of the base of the arrow head relative to the length
autoscale = 1; % Autoscale if ~= 0 then scale by this.
plotarrows = 1;
filled = 0;
ls = '-';
ms = '';
col = '';
nin = nargin;
% Parse the string inputs
while isstr(varargin{nin}),
  vv = varargin{nin};
  if ~isempty(vv) & strcmp(lower(vv(1)),'f')
    filled = 1;
    nin = nin-1;
  else
    [l,c,m,msg] = colstyle(vv);
    if ~isempty(msg),
      error(sprintf('Unknown option "%s".',vv));
    end
    if ~isempty(l), ls = l; end
    if ~isempty(c), col = c; end
    if ~isempty(m), ms = m; plotarrows = 0; end
    if isequal(m,'.'), ms = ''; end % Don't plot '.'
    nin = nin-1;
  end
end
error(nargchk(4,7,nin));
% Check numeric input arguments
if nin<6, % quiver3(z,u,v,w) or quiver3(z,u,v,w,s)
  [msg,x,y,z] = xyzchk(varargin{1});
  u = varargin{2};
  v = varargin{3};
  w = varargin{4};
else % quiver3(x,y,z,u,v,w) or quiver3(x,y,z,u,v,w,s)
  [msg,x,y,z] = xyzchk(varargin{1:3});
  u = varargin{4};
  v = varargin{5};
  w = varargin{6};
end
if ~isempty(msg), error(msg); end
% Scalar expand u,v,w.
if prod(size(u))==1, u = u(ones(size(x))); end
if prod(size(v))==1, v = v(ones(size(u))); end
if prod(size(w))==1, w = w(ones(size(v))); end
% Check sizes
if ~isequal(size(x),size(y),size(z),size(u),size(v),size(w))
  error('The sizes of X,Y,Z,U,V, and W must all be the same.');
end
% Get autoscale value if present
if nin==5 | nin==7, % quiver3(z,u,v,w,s) or quiver3(x,y,z,u,v,w,s)
  autoscale = varargin{nin};
end
if length(autoscale)>1,
  error('S must be a scalar.');
end
if autoscale,
  % Base autoscale value on average spacing in the x and y
  % directions.  Estimate number of points in each direction as
  % either the size of the input arrays or the effective square
  % spacing if x and y are vectors.
  if min(size(x))==1, n=sqrt(prod(size(x))); m=n; else [m,n]=size(x); end
  delx = diff([min(x(:)) max(x(:))])/n;
  dely = diff([min(y(:)) max(y(:))])/m;
  delz = diff([min(z(:)) max(y(:))])/max(m,n);
  del = sqrt(delx.^2 + dely.^2 + delz.^2);
  if del>0
    len = sqrt((u/del).^2 + (v/del).^2 + (w/del).^2);
    maxlen = max(len(:));
  else
    maxlen = 0;
  end

  if maxlen>0
    autoscale = autoscale*0.9 / maxlen;
  else
    autoscale = autoscale*0.9;
  end
  u = u*autoscale; v = v*autoscale; w = w*autoscale;
end
ax = newplot;
next = lower(get(ax,'NextPlot'));
hold_state = ishold;
% Make velocity vectors
x = x(:).'; y = y(:).'; z = z(:).';
u = u(:).'; v = v(:).'; w = w(:).';
uu = [x;x+u;repmat(NaN,size(u))];
vv = [y;y+v;repmat(NaN,size(u))];
ww = [z;z+w;repmat(NaN,size(u))];
h1 = plot3(uu(:),vv(:),ww(:),[col ls]);
if plotarrows,
  beta = beta * sqrt(u.*u + v.*v + w.*w) ./ (sqrt(u.*u + v.*v) + eps);
  uv=sqrt(u.*u + v.*v);
  % Make arrow heads and plot them
  hu = [x+u; x+u-alpha*(u+beta.*(v+eps)); ...
             x+u-alpha*(u-beta.*(v+eps)); ...
        x+u; x+u-alpha*u; x+u-alpha*u; x+u;...
        repmat(NaN,size(u))];
  hv = [y+v; y+v-alpha*(v-beta.*(u+eps)); ...
             y+v-alpha*(v+beta.*(u+eps)); ...
        y+v; y+v-alpha*v; y+v-alpha*v; y+v;...
        repmat(NaN,size(v))];
  hw = [z+w; z+w-alpha*w; z+w-alpha*w; ...
        z+w; z+w-alpha*(w+beta.*(uv+eps)); ...
             z+w-alpha*(w-beta.*(uv+eps)); z+w; ...
        repmat(NaN,size(w))];
  hold on
  h2 = plot3(hu(:),hv(:),hw(:),[col ls]);
else
  h2 = [];
end
if ~isempty(ms), % Plot marker on base
  hu = x; hv = y; hw = z;
  hold on
  h3 = plot3(hu(:),hv(:),hw(:),[col ms]);
  if filled, set(h3,'markerfacecolor',get(h1,'color')); end
else
  h3 = [];
end
if ~hold_state, hold off, view(3); grid on, set(ax,'NextPlot',next); end
if nargout>0, hh = [h1;h2;h3]; end

repo: ifem


text:
problemToSolve = 'test2d_triangle';
origin = 'P1';
variable = 'd_u';
projectDestination = 'P0';


text:
function start_LearnSimulink(varargin)
%{
This function simply opens an HTML document to guide you in this way to
learn Simulink.
Created by Roni Peer, Dec-2011.
%}

cd(fullfile(pwd, 'files'));
open(fullfile(pwd,'html','LearnSimulinkScript.html')); % Open the HTML.


text:
function mat2ades(dat, fileName, Fs, chanlabel, chanunit, dattype, datunit)

% write in the current folder ADES and DAT files from matrix in MATLAB workspace
% data = matrix of data (nbchannel * time points) - the data have to be in microVolt
% fileName = string of the output files without extension ; the ades and dat files will have the same name
% FS = sampling rate
% labels = cell-array with channel labels
% labelType : 'EEG' or 'MEG'
%
% Data are stored in a binary file which name is exactly the same than the header file except the extension: .dat
% The samples are stored as float, 4 bytes per sample, little endian. The channels are multiplexed.
%
% Sophie Chen - January 2014
% Modified by Robert Oostenveld - February 2019

%% generate the ADES file
adesFile = [fileName '.ades'];

fid = fopen_or_error(adesFile, 'wt');

fprintf(fid, '#ADES header file\r\n');
fprintf(fid, 'samplingRate = %d\r\n', Fs);
fprintf(fid, 'numberOfSamples = %d\r\n', size(dat,2));

for lab = 1:length(dattype)
  fprintf(fid, 'Unit = %s, %s\r\n', dattype{lab}, datunit{lab});
end

for lab = 1:length(chanlabel)
  fprintf(fid, '%s = %s\r\n', chanlabel{lab}, chanunit{lab});
end

fclose(fid);

%% generate the DAT file

datFile = [fileName '.dat'];

fad = fopen_or_error(datFile, 'wb');
fwrite(fad, dat, 'float32', 'l');
fclose(fad);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Vpolicy = mdp_eval_policy_matrix(P, R, discount, policy)

% mdp_eval_policy_matrix    Evaluation of the value function of a policy
% Arguments -------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA) = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount = discount rate in ]0; 1[
%   policy(S) = a policy
% Evaluation -------------------------------------------------------------
%   Vpolicy(S) = value function of the policy

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * 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.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% 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.

% check of arguments
if iscell(P); S = size(P{1},1); else S = size(P,1); end;
if discount <= 0 || discount >= 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1[')
    disp('--------------------------------------------------------')
elseif size(policy,1)~=S || any(mod(policy,1)) || any(policy<1) || any(policy>S)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: policy must be a (1xS) vector with integer from 1 to S')
    disp('--------------------------------------------------------')
else

    [Ppolicy, PRpolicy] = mdp_computePpolicyPRpolicy(P, R, policy);

    % V = PR + gPV  => (I-gP)V = PR  => V = inv(I-gP)*PR
    Vpolicy = (speye(S,S) - discount*Ppolicy) \ PRpolicy;

end


text:
function [E,C] = nonmanifold_edges(F)
  % NONMANIFOLD_EDGES List of non-manifold edges
  %
  % Inputs:
  %   F  #F by dim=3 list of facet indices
  % Outputs:
  %   E  #E by 2 list of nonmanifold edges
  %   C  #E by 1 list of unsigned Counts (>2)
  %   %pC  #E by 1 list of positive Counts
  %   %nC  #E by 1 list of negative Counts
  %

  allE = [F(:,[2 3]); F(:,[3 1]); F(:,[1 2])];
  sortallE = sort(allE,2);
  sC = sparse(sortallE(:,1),sortallE(:,2),(allE(:,1)<allE(:,2))*2-1);
  C = sparse(sortallE(:,1),sortallE(:,2),1);
%  assert(nnz(tril(C)) == 0);
% assert(nnz(tril(sC)) == 0);
  [EI,EJ,C] = find((C>2 | abs(sC)>1).* C);
  E = [EI EJ];
end

author: alecjacobson

repo: gptoolbox


text:
%this class retrieves a preprocessing algorithm from the UI
classdef FeatureExtractionConfigurator < handle

    properties (Access = private)
        algorithmConfigurator;

        uiElements FeatureExtractionConfiguratorUIElements;

        selectedFeatures;
        currentSelectedFeature;
    end

    properties (Access = public)
        defaultFeatures;
        loadedFeatureExtractionFiles;
    end

    methods (Access = public)
        function obj = FeatureExtractionConfigurator(defaultFeatures,...
                loadedFeatureExtractionFiles, uiElements)

            obj.defaultFeatures = defaultFeatures;
            obj.loadedFeatureExtractionFiles = loadedFeatureExtractionFiles;

            obj.uiElements = uiElements;

            obj.uiElements.addFeatureButton.ButtonPushedFcn = @obj.handleAddButtonClicked;
            obj.uiElements.removeFeatureButton.ButtonPushedFcn = @obj.handleRemoveButtonClicked;
            obj.uiElements.featuresSourceButtonGroup.SelectionChangedFcn = @obj.handleSelectionTypeChanged;
            obj.uiElements.selectedFeaturesList.ValueChangedFcn = @obj.handleSelectedFeatureChanged;
            obj.uiElements.featureStartRangeEditText.ValueChangedFcn = @obj.handleStartValueChanged;
            obj.uiElements.featureEndRangeEditText.ValueChangedFcn = @obj.handleEndValueChanged;
            obj.uiElements.featureFullSegmentCheckBox.ValueChangedFcn = @obj.handleFullSegmentCheckBoxChanged;
            obj.uiElements.featureAxisEditText.ValueChangedFcn = @obj.handleAxisValueChanged;

            obj.selectedFeatures = [];
            if ~isempty(obj.loadedFeatureExtractionFiles)
                obj.fillLoadedFeaturesList();
                obj.uiElements.loadedFeatureExtractorsList.Value = obj.uiElements.loadedFeatureExtractorsList.Items{1};
            end

            if ~isempty(obj.defaultFeatures)
                obj.reloadUI();
            end
        end

        function reloadUI(obj)
            obj.fillDefaultFeaturesList();
        end

        function featureExtractor = createFeatureExtractorWithUIParameters(obj)
            if obj.isManualFeatureSelectionMode()
                featureExtractor = FeatureExtractor(obj.selectedFeatures);
            else
                fileIdx = obj.getSelectedLoadedFileIdx();
                featureExtractionFile = obj.loadedFeatureExtractionFiles{fileIdx};
                featureExtractor = DataLoader.LoadAlgorithm(featureExtractionFile);
            end
        end

        function idx = getDefaultSelectedFeatureIdx(obj)
            idxStr = obj.uiElements.defaultFeaturesList.Value;
            [~,idx] = ismember(idxStr,obj.uiElements.defaultFeaturesList.Items);
        end

        function feature = getDefaultSelectedFeature(obj)
            idx = obj.getDefaultSelectedFeatureIdx();
            feature = obj.defaultFeatures{idx};
        end

        function idx = getSelectedFeatureIdx(obj)
            idxStr = obj.uiElements.selectedFeaturesList.Value;
            [~,idx] = ismember(idxStr,obj.uiElements.selectedFeaturesList.Items);
        end

        function feature = getSelectedFeature(obj)
            if isempty(obj.uiElements.selectedFeaturesList.Items)
                feature = [];
            else
                idx = obj.getSelectedFeatureIdx();
                feature = obj.selectedFeatures{idx};
            end
        end
    end

    methods(Access = private)

        function updateSelectedFeatureUI(obj)
            if ~isempty(obj.currentSelectedFeature)
                axisSelector = obj.currentSelectedFeature.root;
                rangeSelector = axisSelector.nextAlgorithms{1};

                if isempty(rangeSelector.rangeEnd)
                    obj.uiElements.featureStartRangeEditText.Visible = false;
                    obj.uiElements.featureEndRangeEditText.Visible = false;
                    obj.uiElements.featureFullSegmentCheckBox.Value = true;
                else
                    obj.uiElements.featureStartRangeEditText.Visible = true;
                    obj.uiElements.featureEndRangeEditText.Visible = true;
                    obj.uiElements.featureFullSegmentCheckBox.Value = false;
                    obj.uiElements.featureStartRangeEditText.Value = rangeSelector.rangeStart;
                    obj.uiElements.featureEndRangeEditText.Value = rangeSelector.rangeEnd;
                end
                obj.uiElements.featureAxisEditText.Value = Helper.ArrayToString(axisSelector.axes,' ');
            end
        end

        function updateCurrentSelectedFeatureName(obj)
            idx = obj.getSelectedFeatureIdx();
            featureExtractor = obj.selectedFeatures{idx};
            obj.uiElements.selectedFeaturesList.Items{idx} = FeatureExtractionConfigurator.StringForFeature(featureExtractor);
            obj.uiElements.selectedFeaturesList.Value = obj.uiElements.selectedFeaturesList.Items{idx};
        end

        function addFeatureToSelectedList(obj, algorithm)
            obj.selectedFeatures{end+1} = algorithm;
            obj.uiElements.selectedFeaturesList.Items{end+1} = FeatureExtractionConfigurator.StringForFeature(algorithm);
            obj.uiElements.selectedFeaturesList.Value = obj.uiElements.selectedFeaturesList.Items(end);
            obj.currentSelectedFeature = algorithm;
            obj.updateSelectedFeatureUI();
        end

        function idx = getSelectedLoadedFileIdx(obj)
            idxStr = obj.uiElements.loadedFeatureExtractorsList.Value;
            [~,idx] = ismember(idxStr,obj.uiElements.loadedFeatureExtractorsList.Items);
        end

        function b = isManualFeatureSelectionMode(obj)
            selectedButton = obj.uiElements.featuresSourceButtonGroup.SelectedObject;
            b = (selectedButton == obj.uiElements.featuresSourceManuallyRadio);
        end

        function updateSegmentRangeVisibility(obj)
            visible = ~obj.uiElements.featureFullSegmentCheckBox.Value;
            obj.uiElements.featureStartRangeEditText.Visible = visible;
            obj.uiElements.featureEndRangeEditText.Visible = visible;
        end

        function handleSelectionTypeChanged(obj,~,~)
            if obj.isManualFeatureSelectionMode()
                obj.uiElements.manualFeatureExtractionPanel.Visible = true;
                obj.uiElements.loadedFeatureExtractorsList.Visible = false;
            else
                obj.uiElements.manualFeatureExtractionPanel.Visible = false;
                obj.uiElements.loadedFeatureExtractorsList.Visible = true;
            end
        end

        function handleAddButtonClicked(obj,~,~)

            axisSelector = AxisSelector();
            rangeSelector = RangeSelector();
            featureExtractor = obj.getDefaultSelectedFeature();

            Algorithm.AlgorithmWithSequence({axisSelector,rangeSelector,featureExtractor});

            compositeAlgorithm = CompositeAlgorithm(axisSelector,featureExtractor);
            compositeAlgorithm.name = featureExtractor.name;

            obj.addFeatureToSelectedList(compositeAlgorithm);
        end

        function handleRemoveButtonClicked(obj,~,~)
            idx = obj.getSelectedFeatureIdx();
            obj.selectedFeatures(idx) = [];
            obj.selectedFeaturesList.Items(idx) = [];
            obj.selectFirstAddedFeature();
        end

        function selectFirstAddedFeature(obj)
            if ~isempty(obj.selectedFeatures)
                obj.uiElements.selectedFeaturesList.Value = obj.uiElements.selectedFeaturesList.Items{1};
                obj.currentSelectedFeature = obj.selectedFeatures{1};
            end
        end

        function fillLoadedFeaturesList(obj)
            obj.uiElements.loadedFeatureExtractorsList.Items = obj.loadedFeatureExtractionFiles;
        end

        function fillDefaultFeaturesList(obj)
            obj.uiElements.defaultFeaturesList.Items = Helper.GenerateAlgorithmNamesArray(obj.defaultFeatures);
        end

        function handleSelectedFeatureChanged(obj,~,~)
            obj.currentSelectedFeature = obj.getSelectedFeature();
            obj.updateSelectedFeatureUI();
        end

        function handleStartValueChanged(obj,~,~)
            obj.currentSelectedFeature = obj.getSelectedFeature();

            if ~isempty(obj.currentSelectedFeature)
                rangeSelector = obj.currentSelectedFeature.root.nextAlgorithms{1};
                rangeSelector.rangeStart = obj.uiElements.featureStartRangeEditText.Value;
                obj.updateCurrentSelectedFeatureName();
            end
        end

        function handleEndValueChanged(obj,~,~)
            obj.currentSelectedFeature = obj.getSelectedFeature();

            if ~isempty(obj.currentSelectedFeature)
                rangeSelector = obj.currentSelectedFeature.root.nextAlgorithms{1};
                rangeSelector.rangeEnd = obj.uiElements.featureEndRangeEditText.Value;
                obj.updateCurrentSelectedFeatureName();
            end
        end

        function handleAxisValueChanged(obj,~,~)
            obj.currentSelectedFeature = obj.getSelectedFeature();

            if ~isempty(obj.currentSelectedFeature)
                axisSelector = obj.currentSelectedFeature.root;
                axisSelector.axes = str2num(obj.uiElements.featureAxisEditText.Value);
                obj.updateCurrentSelectedFeatureName();
            end
        end

        function handleFullSegmentCheckBoxChanged(obj,~,~)
            obj.updateSegmentRangeVisibility();

            obj.currentSelectedFeature = obj.getSelectedFeature();
            if ~isempty(obj.currentSelectedFeature)
                rangeSelector = obj.currentSelectedFeature.root.nextAlgorithms{1};

                value = obj.uiElements.featureFullSegmentCheckBox.Value;
                if (value)
                    rangeSelector.rangeEnd = [];
                else
                    rangeSelector.rangeEnd = obj.uiElements.featureEndRangeEditText.Value;
                end
                obj.updateCurrentSelectedFeatureName();
            end
        end
    end

    methods (Access = private, Static)
        function str = StringForFeature(featureAlgorithm)
            axisSelector = featureAlgorithm.root;
            rangeSelector = axisSelector.nextAlgorithms{1};
            featureExtractor = rangeSelector.nextAlgorithms{1};

            str = sprintf('%s_%s',featureExtractor.name,Helper.ArrayToString(axisSelector.axes,' '));
            if ~isempty(rangeSelector.rangeEnd)
                str = sprintf('%s (%d-%d)',str,rangeSelector.rangeStart,rangeSelector.rangeEnd);
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function check = isExceptionStrip(grstruct, fieldname)

% initialize output
check = false;

% exception list {fieldname, val_types}
exceptions = {'color', @iscell, 'width', @(x)(length(x)>1), 'size', @(x)(length(x)>1)};

for e = 1:2:length(exceptions)

    %comparison function
    compfun = exceptions{e+1};

    % look for fieldnames of type exceptions{e} and compare the underyling data using exceptions{e+1}
    if strcmp(fieldname, exceptions{e}) && compfun(grstruct.(fieldname))
        check = true;
    end

end
end

author: plotly


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

%%
clear; close all; clc;

%% Syntax
% |cVoigt=voigtMap(c);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function joint = combine_marginals_into_joint(marginalsT, hnodes, ns)

jointT = dpot(hnodes, ns(hnodes));
for i=hnodes(:)'
  jointT = multiply_by_pot(jointT, marginalsT{i});
end
m = pot_to_marginal(jointT);
joint = m.T(:);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%DEMO_SURVIVAL_COXPH  Survival model using Cox proportional model
%
%  Description
%    Survival model using Cox proportional model with a piecewise
%    log-constant baseline hazard. The hazard rate is
%
%       h(t) = h_0(t)*exp(f),
%
%    where the baseline hazard is assumed to piecewise log-constant.
%
%    The inference is conducted via Laplace, where we find
%    Gaussian approximation for p(f| th, data), where th is the
%    maximum a posterior (MAP) estimate for the parameters.
%
%    The censoring indicator ye is
%
%      ye = 0 for uncensored event
%      ye = 1 for right censored event.
%
%    If survival times y for n observation are given as nx2 matrix with
%    entry times into follow-up in the first column and exit times from
%    follow-up in the second column, left truncated right censored
%    modelling is possible, for instance, in cases where age is wanted to
%    be set as a baseline hazard.
%
%    Example data set is leukemia survival data in Northwest England
%    presented in (Henderson, R., Shimakura, S., and Gorst, D. (2002).
%    Modeling spatial variation in leukemia survival data. Journal of the
%    American Statistical Association, 97:965–972). Data set was downloaded
%    from http://www.math.ntnu.no/%7Ehrue/r-inla.org/examples/leukemia/leuk.dat
%
%  See also  DEMO_SURVIVAL_WEIBULL

% Copyright (c) 2011 Jaakko Riihimäki

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

% First load data
S = which('demo_survival_weibull');
L = strrep(S,'demo_survival_weibull.m','demodata/leukemia.txt');
leukemiadata=load(L);

% leukemiadata consists of:
% 'time', 'cens', 'xcoord', 'ycoord', 'age', 'sex', 'wbc', 'tpi', 'district'

% survival times
y=leukemiadata(:,1);
% scale survival times
y=y/max(y);

ye=1-leukemiadata(:,2); % event indicator, ye = 0 for uncensored event
                  %                        ye = 1 for right censored event

% choose (for example) 'age', 'sex', 'wbc', and 'tpi' covariates
x0=leukemiadata(:,5:8);
x=x0;
% normalize continuous covariates
x(:,[1 3:4])=bsxfun(@rdivide,bsxfun(@minus,x0(:,[1 3:4]),mean(x0(:,[1 3:4]),1)),std(x0(:,[1 3:4]),1));

[n, nin]=size(x);

% number of time intervals
ntime=50;
% create finite partition of time axis
S=linspace(0,max(y)+0.001,ntime+1);

% Create the covariance functions
pl = prior_t('s2',1, 'nu', 4);
pm = prior_t('s2',1, 'nu', 4);

% covariance for hazard function
gpcfh = gpcf_sexp('lengthScale', 1, 'magnSigma2', 1.1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
% covariance for proportional part
gpcf = gpcf_sexp('lengthScale', ones(1,size(x,2)), 'magnSigma2', 1.2, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% Create the likelihood structure
lik = lik_coxph('S', S);

% NOTE! if multiple covariance functions per latent is used, define
% gp.comp_cf as follows:
% gp = gp_set(..., 'comp_cf' {[1 2] [5 6]};
% where [1 2] are for hazard function, and [5 6] for proportional part
gp = gp_set('lik', lik, 'cf', {gpcfh gpcf}, 'jitterSigma2', 1e-6, 'comp_cf', {1 2});

% Set the approximate inference method to Laplace
gp = gp_set(gp, 'latent_method', 'Laplace');

opt=optimset('TolFun',1e-2,'TolX',1e-2,'Display','iter');
gp=gp_optim(gp,x,y,'z',ye,'opt',opt);

% Make prediction
xt1=zeros(200,nin); xt1(:,2)=1;
xt2=zeros(200,nin); xt2(:,2)=-1;

xt1(:,1)=linspace(min(x(:,1)), max(x(:,1)), 200);
xt2(:,1)=linspace(min(x(:,1)), max(x(:,1)), 200);
xt01(:,1)=linspace(min(x0(:,1)), max(x0(:,1)), 200);
xt02(:,1)=linspace(min(x0(:,1)), max(x0(:,1)), 200);

[Ef1, Covf1] = gp_pred(gp, x, y, xt1, 'z', ye);
[Ef2, Covf2] = gp_pred(gp, x, y, xt2, 'z', ye);
Varf1 = diag(Covf1);
Varf2 = diag(Covf2);

figure, hold on, set(gcf, 'color', 'w'),
plot(gp.lik.xtime, Ef1(1:ntime),'k', 'linewidth', 3)
plot(gp.lik.xtime, Ef1(1:ntime)+1.96*sqrt(Varf1(1:ntime)),'--k', 'linewidth', 2)
plot(gp.lik.xtime, Ef1(1:ntime)-1.96*sqrt(Varf1(1:ntime)),'--k', 'linewidth', 2)
title('log-baseline hazard (follow-up time)')
xlabel('time')

% Compute posterior mean and 95% credible intervals of latent function as a
% function of age and sex
col1=ones(1,3)*0.7;
col2=ones(1,3)*0.3;
figure, hold on, set(gcf, 'color', 'w'),
plot(xt01(:,1), Ef1(ntime+1:end), 'color', col1, 'linewidth', 3)
plot(xt01(:,1), Ef1(ntime+1:end)+1.96*sqrt(Varf1(ntime+1:end)), '--', 'color', col1, 'linewidth', 2)
plot(xt01(:,1), Ef1(ntime+1:end)-1.96*sqrt(Varf1(ntime+1:end)), '--', 'color', col1, 'linewidth', 2)

plot(xt02(:,1), Ef2(ntime+1:end), 'color', col2, 'linewidth', 3)
plot(xt02(:,1), Ef2(ntime+1:end)+1.96*sqrt(Varf2(ntime+1:end)), '--', 'color', col2, 'linewidth', 2)
plot(xt02(:,1), Ef2(ntime+1:end)-1.96*sqrt(Varf2(ntime+1:end)), '--', 'color', col2, 'linewidth', 2)
xlabel('age')
title('effect of age for both sexes')

%- Age as baseline hazard

% Age in the beginning and in the end of follow-up
y=[x0(:,1) x0(:,1)+leukemiadata(:,1)/365];
y2=y;
% normalise ages
y2=y2-min(y2(:)); y2=y2./max(y(:));

x2=x;
x2(:,1)=[];

% covariance for proportional part
gpcf = gpcf_sexp('lengthScale', ones(1,size(x2,2)), 'magnSigma2', .5, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% NOTE! if multiple covariance functions per latent is used, define
% gp.comp_cf as follows:
% gp = gp_set(..., 'comp_cf', {[1 2] [5 6]});
% where [1 2] are for hazard function, and [5 6] for proportional part
gp = gp_set('lik', lik, 'cf', {gpcfh gpcf}, 'jitterSigma2', 1e-6, 'comp_cf', {[1] [2]});

% Set the approximate inference method to Laplace
gp = gp_set(gp, 'latent_method', 'Laplace');

opt=optimset('TolFun',1e-2,'TolX',1e-2,'Display','iter');
gp=gp_optim(gp,x2,y2,'z',ye,'opt',opt);

[Ef1, Covf1] = gp_pred(gp, x2, y2, x2, 'z', ye);
Varf1 = diag(Covf1);
xtmpl=linspace(min(y(:)),max(y(:)),50);

figure, hold on, set(gcf, 'color', 'w'),
plot(xtmpl, Ef1(1:ntime),'k', 'linewidth', 3)
plot(xtmpl, Ef1(1:ntime)+1.96*sqrt(Varf1(1:ntime)),'--k', 'linewidth', 2)
plot(xtmpl, Ef1(1:ntime)-1.96*sqrt(Varf1(1:ntime)),'--k', 'linewidth', 2)
title('log-baseline hazard (age)')
xlabel('age')


text:
function [n,filters] = plx_chan_filters(filename)
% plx_chan_filters(filename): read channel filter settings for each spike channel from a .plx or .pl2 file
%
% [n,filters] = plx_chan_filters(filename)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%
% OUTPUT:
%   filter - array of filter values (0 or 1)
%   n - number of channels

if nargin ~= 1
    error 'expected 1 input argument';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2 = PL2GetFileIndex(filename);
    n = numel(pl2.SpikeChannels);
    if n > 0
        filters = zeros(n,1);
    end
    return
end

[n,filters] = mexPlex(10, filename);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function dicomrt_plotgammamap2dvol(gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,ct,ct_xmesh,ct_ymesh,ct_zmesh,VOI,voi2use)
% dicomrt_plotgammamap2dvol(gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,ct,ct_xmesh,ct_ymesh,ct_zmesh,VOI,voi2use)
%
% Plot gamma maps.
%
% gamma is the ganmma matrix
% gamma_xmesh,gamma_ymesh,gamma_zmesh are x-y-z coordinates of the center of the gamma-pixel
% ct is the 3D CT dataset (OPTIONAL)
% ct_xmesh, ct_ymesh, ct_zmesh are the coordinates of the ct voxels (OPTIONAL)
% VOI is an OPTIONAL cell array which contain the patients VOIs as read by dicomrt_loadVOI
% voi2use is an OPTIONAL number pointing to the number of VOIs to be displayed
%
% Example:
%
% dicomrt_plotgammamap2dvol(gamma2dvol,xmesh,ymesh,zmesh,demo1_ct,ct_xmesh,ct_ymesh,ct_zmesh,demo1_voi,9);
%
% plot the gamma map in gamma2dvol on top of a ct slice and with contour of VOI 9.
%
% See also dicomrt_GAMMAcal2D, dicomrt_GAMMAcal2DVol, dicomrt_rendervoi
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(4,11,nargin))

if size(gamma,3)==1
    error('dicomrt_plotgammamap2dvol: Format not supported, try dicomrt_plotgammamap. Exit now !');
end

if exist('VOI')==1 & exist('voi2use')==1
    [mask_gamma]=dicomrt_mask(VOI,gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,voi2use,'n');
else
    mask_gamma=gamma;
end

if size(gamma,3)~=length(VOI{voi2use,2})
    disp('Warning: the size of gamma does not match the size of the VOI!');
    disp('There is the possibility that the gamma map was calculated for a different VOI.');
    goahead=input('dicomrt_plotgammamap2dvol: Do you want to continue (y/n) ? (n=default)','s');
    if goahead~='y' | goahead~='Y'
        error('dicomrt_plotgammamap2dvol: Exiting as requested');
    elseif isempty(goahead)==1
        error('dicomrt_plotgammamap2dvol: Exiting as requested');
    end
end

if length(size(ct))==2 & ct==0
    for k=1:size(gamma,3) % Loop over gamma slices
        % plot ct map
        figure;
        set(gcf,'Name',['dicomrt_plotgammamap2dvol: ',inputname(1),', Z= ',num2str(gamma_zmesh(k))]);
        set(gcf,'NumberTitle','off');
        if exist('VOI')==1 & exist('voi2use')==1
            % plot voi on the same slice
            p=plot3(VOI{voi2use,2}{k,1}(:,1),VOI{voi2use,2}{k,1}(:,2),...
                VOI{voi2use,2}{k,1}(:,3));
            % set VOI contour to Z=-1 to display together with gamma map
            p_handle=get(p);
            ZData_temp_2=p_handle.ZData;
            for i=1:length(ZData_temp_2)
                ZData_temp_2(i)=-1;
            end
            set(p, 'Color', 'k','LineWidth', 1.5);
            set(p,'ZData',ZData_temp_2);
        end
        hold;
        % plot gamma map
        %mesh(mask_gamma(:,:,k),'XData',gamma_xmesh,'YData',gamma_ymesh);
        % ms=surf(mask_gamma(:,:,k));
        ms=surfl(mask_gamma(:,:,k));
        shading interp;
        set(ms,'XData',gamma_xmesh,'YData',gamma_ymesh);
        % add labels
        xlabel('X axis (cm)','Fontsize',12);
        ylabel('Y axis (cm)','Fontsize',12);
        zlabel('\gamma','Fontsize',14);
        title(['\gamma map slice: ',num2str(k)],'Fontsize',18);
        axis tight;
    end
else
    % rebuilding mesh matrices
    [xmin,xmax,ymin,ymax,zmin,zmax]=dicomrt_voiboundaries(ct_xmesh,ct_ymesh, ct_zmesh,VOI,voi2use);
    [ct3d_xmesh,ct3d_ymesh,ct3d_zmesh]=dicomrt_rebuildmatrix(ct_xmesh(xmin:xmax),ct_ymesh(ymin:ymax),ct_zmesh);
    for k=1:size(gamma,3) % Loop over gamma slices
        % plot ct map
        figure;
        set(gcf,'Name',['dicomrt_plotgammamap2dvol: ',inputname(1),', Z= ',num2str(ct_zmesh(dicomrt_findsliceVECT(gamma_zmesh,k,ct_zmesh)))]);
        set(gcf,'NumberTitle','off');
        handle=slice(ct3d_xmesh,ct3d_ymesh,ct3d_zmesh,ct(ymin:ymax,xmin:xmax,:),[],[],ct_zmesh(dicomrt_findsliceVECT(gamma_zmesh,k,ct_zmesh)));
        % set ct map to Z=-1 to display together with gamma map
        set(handle,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
        gco_handle=get(handle);
        ZData_temp_1=gco_handle.ZData;
        for i=1:size(ZData_temp_1,1)
            for j=1:size(ZData_temp_1,2)
                ZData_temp_1(i,j)=-1;
            end
        end
        set(handle,'ZData',ZData_temp_1);
        hold;
        if exist('VOI')==1 & exist('voi2use')==1
            % plot voi on the same slice
            if k==23
                disp(num2str(k));
            end
            p=plot3(VOI{voi2use,2}{k,1}(:,1),VOI{voi2use,2}{k,1}(:,2),...
                VOI{voi2use,2}{k,1}(:,3));
            % set VOI contour to Z=-1 to display together with gamma map
            p_handle=get(p);
            ZData_temp_2=p_handle.ZData;
            for i=1:length(ZData_temp_2)
                ZData_temp_2(i)=-1;
            end
            set(p, 'Color', 'w','LineWidth', 1.5);
            set(p,'ZData',ZData_temp_2);
        end
        % plot gamma map
        mesh(mask_gamma(:,:,k),'XData',gamma_xmesh,'YData',gamma_ymesh);
        % add labels
        xlabel('X axis (cm)','Fontsize',12);
        ylabel('Y axis (cm)','Fontsize',12);
        zlabel('\gamma','Fontsize',14);
        title(['\gamma map slice: ',num2str(k)],'Fontsize',18);
        axis tight;
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcSurfl(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

figure;
surfl(double(handles.BMatrix));
shading interp;
colormap(handles.V.Color_map);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [ a, seed ] = d3_uniform ( n, seed )

%*****************************************************************************80
%
%% D3_UNIFORM randomizes a D3 matrix.
%
%  Discussion:
%
%    The D3 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%  Example:
%
%    Here is how a D3 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the linear system.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(3,N), the D3 matrix.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  a(1,1) = 0.0;
  [ a(1,2:n), seed ] = r8vec_uniform ( n-1, 0.0, 1.0, seed );

  [ a(2,1:n), seed ] = r8vec_uniform ( n,   0.0, 1.0, seed );

  [ a(3,1:n-1), seed ] = r8vec_uniform ( n-1, 0.0, 1.0, seed );
  a(3,n) = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [GAmp,GTime]=GySpiral(p)

global VCtl;
global VObj;
global VVar;

tStart=p.tStart;
dt = p.dt;

% 2D spiral encoding
FOV = VCtl.FOVFreq; % choose FOVFreq as FOV
Res = VCtl.ResFreq; % choose ResFreq as effective resolution

KMax  = Res/ (2*FOV);
ThetaMax = KMax / VCtl.S_Lamda;
Beta = (VObj.Gyro/(2*pi))*(VCtl.S_SlewRate / VCtl.S_Lamda);
a2 = (9*Beta/4)^(1/3);
CLamda = VCtl.S_SlewRate/VCtl.S_SlewRate0;
Ts = ((3*VObj.Gyro*VCtl.S_Gradient)/(4*pi*VCtl.S_Lamda*a2^2))^3;

ThetaTs = (0.5 * Beta * Ts^2)/(CLamda + (Beta/(2 * a2)) * Ts^(4/3));
if ThetaTs < ThetaMax
    t1 = 0:dt:Ts;
    Theta1 = (0.5 * Beta * t1.^2)./(CLamda + (Beta/(2 * a2)) * t1.^(4/3));

    Tacq = t1(end) + ((pi*VCtl.S_Lamda)/(VObj.Gyro * VCtl.S_Gradient))*(ThetaMax^2 - Theta1(end)^2);
    t2 = t1(end):dt:Tacq;
    Theta2 = sqrt(Theta1(end)^2 + (VObj.Gyro / (pi * VCtl.S_Lamda))*VCtl.S_Gradient*(t2 - t1(end)));

    Theta = [Theta1(1:end-1) Theta2];
    t     = [t1(1:end-1)     t2];

else
    Tacq = ((2*pi*FOV)/(3*VCtl.S_ShotNum))*sqrt(pi/(VObj.Gyro*VCtl.S_SlewRate*VCtl.RFreq^3));
    t1 = 0:dt:Tacq;
    Theta1 = (0.5 * Beta * t1.^2)./(CLamda + (Beta/(2 * a2)) * t1.^(4/3));

    Theta = Theta1;
    t     = t1;
end

DTheta = [0 diff(Theta)./diff(t)];

GAmp = ((2*pi)/(VObj.Gyro))*VCtl.S_Lamda*DTheta.*(sin(Theta + (VVar.PhaseCount-1)*(2*pi/VCtl.S_ShotNum)) ...
                                         + Theta.*cos(Theta + (VVar.PhaseCount-1)*(2*pi/VCtl.S_ShotNum)));
GTime = tStart + VCtl.TEAnchorTime + t;

GTime = [GTime GTime(end) + max(VCtl.MinUpdRate,abs(GAmp(end)/VCtl.S_SlewRate))];
GAmp = [GAmp 0];

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

end


text:
classdef MimPluginGroupPanel < GemPanel
    % MimPluginGroupPanel. Part of the gui for the TD MIM Toolkit.
    %
    %     This class is used internally within the TD MIM Toolkit to help
    %     build the user interface.
    %
    %     MimPluginGroupPanel is a panel containing a group of Plugin buttons
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        CachedMinY
        CachedMaxY
    end

    properties (Access = private)
        CachedPanelHeight
        CachedPanelWidth

        Category
        CurrentCategoryMap
        PluginButtonHandlesMap
        LoadManualSegmentationCallback;
    end

    methods
        function obj = MimPluginGroupPanel(parent, category, current_category_map, load_manual_segmentation_callback)
            obj = obj@GemPanel(parent);
            obj.LoadManualSegmentationCallback = load_manual_segmentation_callback;
            obj.Enabled = false;

            obj.Category = category;
            obj.CurrentCategoryMap = current_category_map;

            obj.AddPlugins(current_category_map);
        end

        function CreateGuiComponent(obj, position)
            CreateGuiComponent@GemPanel(obj, position);
            set(obj.GraphicalComponentHandle, 'Title', obj.Category, 'BorderType', 'etchedin');
        end

        function Resize(obj, new_position)
            Resize@GemPanel(obj, new_position);
            width = new_position(3);
            if isempty(obj.CachedPanelHeight) || (width ~= obj.CachedPanelWidth)
                obj.ResizePanel(width);
            end
        end

        function height = GetRequestedHeight(obj, width)
            if isempty(obj.CachedPanelHeight) || (width ~= obj.CachedPanelWidth)
                obj.ResizePanel(width);
            end
            height = obj.CachedPanelHeight;
        end

        function AddPlugins(obj, category_map)
            obj.CachedPanelHeight = [];
            obj.CachedPanelWidth = [];
            obj.PluginButtonHandlesMap = containers.Map;

            % Add the buttons to the panel
            for current_plugin_key = category_map.keys
                current_plugin = category_map(char(current_plugin_key));

                if isempty(current_plugin)
                    button_handle = MimSegmentationButton(obj, current_plugin_key{1}, obj.LoadManualSegmentationCallback);
                else
                    callback_function_handle = @current_plugin.RunPlugin;
                    button_handle = MimPluginButton(obj, callback_function_handle, current_plugin);
                end
                obj.AddChild(button_handle);

                obj.PluginButtonHandlesMap(char(current_plugin_key)) = button_handle;
            end
        end

        function UpdateForNewImage(obj, preview_fetcher, window, level)
            % Refresh the preview images for every button using the supplied dataset

            for button = obj.PluginButtonHandlesMap.values
                button{1}.AddPreviewImage(preview_fetcher, window, level);
            end
        end

        function plugin_found = AddPreviewImage(obj, plugin_name, preview_fetcher, window, level)
            % Refresh the preview image for the specified plugin

            plugin_found = obj.PluginButtonHandlesMap.isKey(plugin_name);
            if plugin_found
                plugin_button = obj.PluginButtonHandlesMap(plugin_name);
                plugin_button.AddPreviewImage(preview_fetcher, window, level);
            end
        end
    end

    methods (Access = private)
        function ResizePanel(obj, panel_width)

            category_map = obj.CurrentCategoryMap;
            button_spacing_w = 10;
            button_spacing_h = 5;
            header_height = 20;
            footer_height = 10;
            left_right_margins = 10;

            max_x = panel_width;
            position_x = left_right_margins;
            position_y = 0;

            last_y_coordinate = 0;
            row_height = 0;

            % Determine coordinates of buttons and the required panel size
            for current_plugin_key = category_map.keys
                current_plugin = category_map(char(current_plugin_key));

                button_handle = obj.PluginButtonHandlesMap(char(current_plugin_key));
                button_width = button_handle.ButtonWidth;
                button_height = button_handle.ButtonHeight;

                current_plugin.ParsedPluginInfo.X = position_x;
                current_plugin.ParsedPluginInfo.Y = position_y;
                current_plugin.ParsedPluginInfo.W = button_width;
                current_plugin.ParsedPluginInfo.H = button_height;

                last_y_coordinate = position_y;
                row_height = max(row_height, button_height);
                last_row_height = row_height;

                category_map(char(current_plugin_key)) = current_plugin;

                position_x = position_x + button_spacing_w + button_width;
                if (position_x + button_width) > (max_x - button_spacing_w)
                    position_y = position_y + button_spacing_h + row_height;
                    position_x = left_right_margins;
                    row_height = 0;
                end

            end

            obj.CachedPanelHeight = last_y_coordinate + last_row_height + header_height + footer_height;
            obj.CachedPanelWidth = panel_width;

            % Resize the buttons
            for current_plugin_key = category_map.keys
                current_plugin = category_map(char(current_plugin_key));

                position_x = current_plugin.ParsedPluginInfo.X;
                button_width = current_plugin.ParsedPluginInfo.W;
                button_height = current_plugin.ParsedPluginInfo.H;
                position_y = obj.CachedPanelHeight - button_height - header_height - current_plugin.ParsedPluginInfo.Y;

                new_position = [position_x, position_y, button_width, button_height];

                button_handle = obj.PluginButtonHandlesMap(char(current_plugin_key));
                button_handle.Resize(new_position);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function p = problemclass(F,h)

p = problemclass(lmi(F),h);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% X = canlab_glm_convolve(stimObj, dsgn, nscan, basis)
%
% This function takes an spm style conditions object (see below), a canlab
% DSGN object (see below), a frame count and a basis function type and
% returns a timeseries vector corresponding to onsets and durations
% specified in stimObj and any modulations specified in dsgn, convolved
% with the basis function specified by basis. Useful if you want to treat
% experimental factors as confounds in a multiple regressors matrix instead
% of treating them as events in SPM. Primarly desinged to allow for
% multiple different basis functions to be used for different design
% factors (e.g. spline interpolation for factors of interest and canonical
% interpolation for confounding factors; here you would implement the
% confounding factors as multiple regressor columns instead after
% convolution).
%
% Input ::
%
%   stimObj     - an SPM-stype condition object. see
%                   canlab_glm_subject_levels('dsgninfo') under
%                   DSGN.conditions for an example
%
%   dsgn        - see canlab_glm_subject_levels('dsgninfo') for how to
%                   format this object.
%
%   nscan       - number of TRs in the scan for which stimObj applies.
%
%   basis       - currently only 'hrf' is supported, although others could
%                   be implemented easily enough. See spm_get_bf.m for some
%                   example implementations of other basis functions. The
%                   implementation here follows that function closely.
% Optional Input ::
%
%   hrfParam    - a struct following the basis argument. If specifying
%                   spline or hrf you will need to specify a structure with
%                   windowlength, order, and (for splines) degree
%                   properties.
%
% Output ::
%
%   X           - A timeseries vector corresponding to the stimuli
%                   specified in stimObj.
%
% Written by Bogdan Petre, 5/19/2021

function [X, bf] = canlab_glm_convolve(stimObj, dsgn, nscan, basis, varargin)
if ~isempty(varargin)
    hrfParams = varargin{1};
end

% we're going to manually convolve stimFactor
spm('defaults','fmri')
spm_jobman('initcfg')

fMRI_T = spm_get_defaults('stats.fmri.t');
if ismember('fmri_t',fieldnames(dsgn))
	if ~isempty(dsgn.fmri_t), fMRI_T = dsgn.fmri_t; end
end

T = dsgn.tr;
dt = dsgn.tr/fMRI_T;

xBF.dt = dt;
xBF.UNITS = 'secs';
xBF.T     = fMRI_T;
xBF.T0    = dsgn.fmri_t0;

if strcmp(basis,'hrf')
    xBF.name = 'hrf';
else
    switch basis
        case 'spline'
            xBF.name   = 'Spline set';
            xBF.degree = hrfParams.degree;
        case 'fourier'
            xBF.name = 'Fourier set';
        case 'fourier_han'
            xBF.name = 'Fourier set (Hanning)';
        case 'gamma'
            xBF.name = 'Gamma functions';
        case 'fir'
            xBF.name = 'Finite Impulse Response';
    otherwise
        error('Unknown basis functions.');
    end
    xBF.length = hrfParams.windowlength;
	xBF.order  = hrfParams.order;
end

switch basis
    case 'hrf'
        [bf, p] = spm_hrf(dt,[],fMRI_T);

        % time derivative (untested)
        try
            if logical(dsgn.convolution.derivs(1))
                dp       = 1;
                p(6)     = p(6) + dp;
                D        = (bf(:,1) - spm_hrf(dt,p,fMRI_T))/dp;
                bf       = [bf D(:)];
                p(6)     = p(6) - dp;
            end
        end

        % dispersion derivative (untested)
        try
            if logical(dsgn.convolution.derivs(2))
                dp   = 0.01;
                p(3) = p(3) + dp;
                D    = (bf(:,1) - spm_hrf(dt,p,fMRI_T))/dp;
                bf   = [bf D(:)];
            end
        end

    case {'spline','fir'}
        xBF = spm_get_bf(xBF);
        bf = xBF.bf;
    otherwise
        error('Only Canonical HRF, spline or FIR convolution is currently supported');
end

if ~strcmp(xBF.name,'Spline set')
    bf = spm_orth(bf);
end

U = [struct(...
    'name','',...
    'ons',[],...
    'dur',[],...
    'orth', 1,...
    'P', struct('name','none','h',0))];

if isempty(stimObj)
    X = [];
else
    % rating and cue convolution
    U(1).name = {'manualConvFactor'};
    U(1).ons = stimObj.onset{1}(:);
    U(1).dur = stimObj.duration{1}(:);

    % this borrws from spm_fMRI_design.m which implements the
    % convolution when spm is invoked directly.
    xBF = struct('T',fMRI_T,'T0',dsgn.fmri_t0,'dt',dt,'UNITS','secs');
    this = struct('Sess',struct('U',U), 'nscan', nscan, 'xBF', xBF);

    U = spm_get_ons(this,1);

    [X,Xn,Fc] = spm_Volterra(U, bf, 1);

    if ~isempty(X)
        X = X((0:(this.nscan - 1))*fMRI_T + dsgn.fmri_t0 + 32,:);
    end
    %{
    switch xBF.name
        case {'hrf', 'fir'}
            %-Orthogonalise (within trial type)
            %----------------------------------------------------------------------
            for i = 1:length(Fc)
                if i<= numel(U) && ... % for Volterra kernels
                        (~isfield(U(i),'orth') || U(i).orth)
                    p = ones(size(Fc(i).i));
                else
                    p = Fc(i).p;
                end
                for j = 1:max(p)
                    X(:,Fc(i).i(p==j)) = spm_orth(X(:,Fc(i).i(p==j)));
                end
            end
    end
    %}

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Calculate border coverage for detected fields in a circular arena
%
% This function calculates firing map border coverage that is further used
% in calculation of a border score. This function must be used with recordings
% done in a circular environment. See USAGE for details.
%
%  USAGE
%   coverage = analyses.borderCoverageCircular(fieldsMap, <options>)
%   fieldsMap   2D binary matrix that repressents firing properties of the field.
%               fieldsMap must be a polar version of a regular firing rate map. This
%               means that x-axis should contain values in range 1-360, y-axis range
%               is 1 to radius of the circular arena.
%               fieldsMap should contain NaN for unvisited bins. Other possible values:
%               0 for zero firing rate, i.e. area that doesn't belong to a fields.
%               Any positive number marks an area of a field. For example, if there is
%               just a single field, then fieldsMap could consist of 3 values: NaNs, zeros,
%               and ones.
%   <options>   Optional list of property-value pairs (see table below)
%
%   ==============================================================================================
%    Properties    Values
%   ----------------------------------------------------------------------------------------------
%    'searchWidth'  If map is not perfect, but contains NaN values along borders, then
%                   search for border pixels can have NaNs. To mitigate this, we check
%                   searchWidth rows/columns near border and if the closest to the border pixel
%                   equals to NaN, we search for first non-NaN value in searchWidth rows/columns.
%                   This argument is optional and default value is 8 bins.
%    'walls'        Definition of walls along which the coverage is calculated. Provided by
%                   a 2D matrix Nx2. Each row defines 1 wall. Each row should have two values in degrees:
%                   1. angle at which wall starts.
%                   2. angle at which wall ends.
%                   Default value is [1 360], which represents a single wall around the whole circular
%                   environment. Note that you can not use 0 as an angle.
%   ==============================================================================================
%   coverage    Border coverage, ranges from 0 to 1.
%
function coverage = borderCoverageCircular(fieldsMap, varargin)
    coverage = 0;

    inp = inputParser;
    defaultSearchWidth = 8;
    defaultWalls = [1 360];

    checkSearchWidth = @(x) isnumeric(x) && isscalar(x) && (x > 0);
    checkWalls = @(x) size(x, 2) == 2 && all(x(:) > 0) && all(x(:) <= 360);

    addRequired(inp, 'fieldsMap', @(x) isnumeric(x) && size(x, 2) == 360);
    addParameter(inp, 'searchWidth', defaultSearchWidth, checkSearchWidth);
    addParameter(inp, 'walls', defaultWalls, checkWalls);

    parse(inp, fieldsMap, varargin{:});

    walls = inp.Results.walls;
    searchWidth = inp.Results.searchWidth;
    fieldsMap(fieldsMap > 1) = 1; % make it binary

    for i = 1:size(walls, 1)
        wall = walls(i, :);

        aux_map = fieldsMap(end-searchWidth+1:end, wall(1):wall(2));
        [covered, norm] = wall_field(aux_map);
        coverage = max([covered/norm, coverage]);
    end
end

% 'covered' pixels will have distance to border 0.
% Essentially we need to calculate number of elements,
% that equal to zero and take NaNs into account.
function [covered, norm] = wall_field(map)
    lx = size(map, 2);
    map = flipud(map);

    D = bwdist(map);
    nanIndices = find(isnan(map(1, :)));
    numNans = length(nanIndices);
    for i = 1:numNans
        testColumn = nanIndices(i);
        nonNan = find(~isnan(map(:, testColumn)), 1, 'first');
        if ~isempty(nonNan)
            numNans = numNans - 1;
            D(1, testColumn) = D(nonNan, testColumn);
        end
    end

    norm = lx - numNans;
    covered = nansum(D(1, :) == 0) - numNans;
end

author: kavli-ntnu


text:
classdef (Abstract) Attribute < matlabshared.libiio.base
    % Attribute IIO attribute function calls

    methods (Hidden)

        function setAttributeLongLong(obj,id,attr,value,isOutput,tol)
            phydev = getDev(obj, obj.phyDevName);
            chanPtr = iio_device_find_channel(obj,phydev,id,isOutput);%FIXME (INVERSION)
            status = cPtrCheck(obj,chanPtr);
            cstatus(obj,status,['Channel: ' id ' not found']);
            status = iio_channel_attr_write_longlong(obj,chanPtr,attr,value);
            cstatus(obj,status,['Attribute write failed for : ' attr ' with value ' num2str(value)]);
            % Check
            [status, rValue] = iio_channel_attr_read_longlong(obj,chanPtr,attr);
            cstatus(obj,status,['Error reading attribute: ' attr]);
            if nargin<6
                tol = sqrt(eps);
            end
            if abs(value - rValue) > tol
                status = -1;
                cstatus(obj,status,['Attribute ' attr ' return value ' num2str(rValue) ', expected ' num2str(value)]);
            end
        end

        function rValue = getAttributeLongLong(obj,id,attr,isOutput)
            phydev = getDev(obj, obj.phyDevName);
            chanPtr = iio_device_find_channel(obj,phydev,id,isOutput);%FIXME (INVERSION)
            status = cPtrCheck(obj,chanPtr);
            cstatus(obj,status,['Channel: ' id ' not found']);
            [status, rValue] = iio_channel_attr_read_longlong(obj,chanPtr,attr);
            cstatus(obj,status,['Error reading attribute: ' attr]);
        end

        function setAttributeBool(obj,id,attr,value,isOutput)
            phydev = getDev(obj, obj.phyDevName);
            chanPtr = iio_device_find_channel(obj,phydev,id,isOutput);%FIXME (INVERSION)
            status = cPtrCheck(obj,chanPtr);
            cstatus(obj,status,['Channel: ' id ' not found']);
            status = iio_channel_attr_write_bool(obj,chanPtr,attr,value);
            cstatus(obj,status,['Attribute write failed for : ' attr]);
            % Check
            [status, rValue] = iio_channel_attr_read_bool(obj,chanPtr,attr);
            cstatus(obj,status,['Error reading attribute: ' attr]);
            if value ~= rValue
                status = -1;
                cstatus(obj,status,['Attribute ' attr ' return value ' num2str(rValue) ', expected ' num2str(value)]);
            end
        end

        function rValue = getAttributeBool(obj,id,attr,isOutput)
            phydev = getDev(obj, obj.phyDevName);
            chanPtr = iio_device_find_channel(obj,phydev,id,isOutput);%FIXME (INVERSION)
            status = cPtrCheck(obj,chanPtr);
            cstatus(obj,status,['Channel: ' id ' not found']);
            [status, rValue] = iio_channel_attr_read_bool(obj,chanPtr,attr);
            cstatus(obj,status,['Error reading attribute: ' attr]);
        end

        function setAttributeRAW(obj,id,attr,value,isOutput)
            phydev = getDev(obj, obj.phyDevName);
            chanPtr = iio_device_find_channel(obj,phydev,id,isOutput);%FIXME (INVERSION)
            status = cPtrCheck(obj,chanPtr);
            cstatus(obj,status,['Channel: ' id ' not found']);
            bytes = iio_channel_attr_write(obj,chanPtr,attr,value);
            if bytes <= 0
                status = -1;
                cstatus(obj,status,['Attribute write failed for : ' attr ' with value ' value]);
            end
        end

        function rValue = getAttributeRAW(obj,id,attr,isOutput)
            phydev = getDev(obj, obj.phyDevName);
            chanPtr = iio_device_find_channel(obj,phydev,id,isOutput);%FIXME (INVERSION)
            status = cPtrCheck(obj,chanPtr);
            cstatus(obj,status,['Channel: ' id ' not found']);
            [status, rValue] = iio_channel_attr_read(obj,chanPtr,attr,1024);
            cstatus(obj,status,['Error reading attribute: ' attr]);
        end

        function setDeviceAttributeRAW(obj,attr,value)
            phydev = getDev(obj, obj.phyDevName);
            bytes = iio_device_attr_write(obj,phydev,attr,value);
            if bytes <= 0
                status = -1;
                cstatus(obj,status,['Attribute write failed for : ' attr ' with value ' value]);
            end
        end

        function rValue = getDeviceAttributeRAW(obj,attr)
            phydev = getDev(obj, obj.phyDevName);
            [status, rValue] = iio_device_attr_read(obj,phydev,attr);
            cstatus(obj,status,['Error reading attribute: ' attr]);
        end

    end
end

author: analogdevicesinc


text:
function Angles = QUADangles(vertices)

% Purpose: To find the angles of Quadrilateral when vertices are known
% Synopsis: Angles = FindAngles(vertices)
%
% Variable description:
% INPUT: vertices - coordinates of Quadrilateral in order 4X2 . The vertices
%                   are given input in anticlock wise direction.
%
% OUTPUT: Angles - Row vector of angles of Quadrilateral in degrees of order 4X1
%
% Example: vertices = [4  4;4  6;5 6; 5  4]
%          angles = FindAngles(vertices);
%          angles = [90 90 90 90]
%
% Algorithm: Calculates the angles using dot product
%            theta = arccos((a.b)/(|a||b|)
%
% Release : 1.0
% Release Date: 03/07/2012
%--------------------------------------------------------------------------
% Author : Siva Srinivas Kolukula
%          Senior Research Fellow
%          Structural Mechanics Laboratory
%          Indira Gandhi Center for Atomic Research
%          India
% E-mail : allwayzitzme@gmail.com
%          http://sites.google.com/site/kolukulasivasrinivas/
%--------------------------------------------------------------------------

Angles = zeros(1,4) ;
order = [ 1 2 4; 2 1 3; 3 4 2;4 3 1] ;
for i = 1:4
    % Form the vectors along the vertices
    data1 = [vertices(i,:); vertices(order(i,2),:)];
    data2 = [vertices(i,:); vertices(order(i,3),:)];
    % Use dor ptodict to find the angle between vector1 and vector2
    Angles(i) = acos((diff(data1)*diff(data2)')/.....
        (sqrt(sum(diff(data1,[],1).^2,2))*sqrt(sum(diff(data2,[],1).^2,2)))) ;
end
% Convert the angles in radians to degrees
Angles = Angles*180./pi ;        %

author: Sable

repo: mcbench-benchmarks


text:
function regresdiagplot3d(sdist,odist,rdist,cutoffsd,cutoffod,cutoffrd,k,class,multi,labsd,labod,labresd,labels)

%REGRESDIAGPLOT3D is a 3D-outlier map which visualizes the orthogonal distance,
% the score distance and the residual distance calculated in a PCR or PLSR analysis.
%
% I/O: regresdiagplot3D(sdist,odist,rdist,cutoffsd,cutoffrd,k,class,labsd,labod,labresd,labels)
%
% Example: out=rpcr(X,Y);
%          regresdiagplot3d(out.sd,out.od,out.resd,out.cutoff.sd,cutoff.od,...
%          out.cutoff.resd,out.k,out.class,1,out.labsd,out.labod,out.labresd,0)
%
% Uses function: putlabel
%
% Created on by S.Verboven
% Last revision: 09/04/2004
%

%INITIALIZATION%
if nargin<13
    labels=0;
end
if nargin<10
    labsd=3;
    labod=3;
    labresd=3;
    labels=0;
end
if nargin==10
    labsd=3
    labels=0;
    labresd=3;
    labod=3;
end
if nargin==11
    labels=0;
    labresd=3;
    labod=3;
end
if nargin==12
    labels=0;
    labresd=3;
end
if nargin<9
    error('A required input variable is missing!')
end

%all LTS-analysis in RPCR are intercept included!!!
% if ask==2 %multivariate analysis
%     %residual distances
%
% else  %univariate analysis
%     %standardized residuals
%     cutoffz=2.5;
% end
cutoffxx=cutoffsd;
cutoffxy=cutoffsd;
cutoffyy=cutoffod;
cutoffyx=cutoffod;
cutoffz=cutoffrd;
x=sdist;
y=odist;
z=rdist;

%%%%%%%MAIN FUNCTION%%%%%%%
set(gcf,'Name', '3D-Outlier map (regression)', 'NumberTitle', 'off')%,'Renderer','OpenGL');

%%%%%%%%Odist=0 not yet included in this standalone!!!!!!!!! included in
%%%%%%%%makeplot function!!!!
plot3(x,y,z,'ko','markerfacecolor',[0.75 0.75 0.75])

hold on
axhandle=gca;
ylen=get(axhandle, 'Ylim');
xlen=get(axhandle,'Xlim');
zlen=get(axhandle,'Zlim');
xrange=xlen(2)-xlen(1);
upLimx=max(cutoffxx,xlen(2))+xrange*0.1;
lowLimx=xlen(1)-xrange*0.1;
yrange=ylen(2)-ylen(1);
upLimy=ylen(2)+yrange*0.1;
lowLimy=ylen(1)-yrange*0.1;
zrange=zlen(2)-zlen(1);
upLimz=max(cutoffz,zlen(2))+zrange*0.1;
if cutoffz==2.5
    lowLimz=min(-2.5,zlen(1))-zrange*0.1;
else
    lowLimz=min(cutoffz,zlen(1))-zrange*0.1;
end

%axis square;
set(gca, 'Xlim',[lowLimx upLimx],'Ylim',[lowLimy upLimy],'Zlim',[lowLimz,upLimz])
hold on

xlabel('Score distance')
ylabel('Orthogonal distance')
if cutoffz~=2.5
    zlabel('Residual distance')
else
    zlabel('Standardized Residual')
end
grid on

%in XY-space
%red plane "vertical" on x axis
oppy=[upLimy:-0.01:lowLimy];
n=length(oppy);
h=(upLimz-lowLimz)/n;
oppz=[lowLimz:h:upLimz];

X1=cutoffxx*ones(n);
Y1=repmat(oppy,n,1);
Z1=repmat([oppz(1:n-1)'; upLimz],1,n);
surf(X1,Y1,Z1,'edgecolor','none','facecolor','r')
alpha(.3)

%blue "horizontal" planes orthogonal on z-axis
%in ZX and ZY-space
oppx=[lowLimx:0.05:upLimx];
n=length(oppx);
X2=repmat(oppx,n,1);
h=(upLimy-lowLimy)/n;
oppy=[lowLimy:h:upLimy];
Y2=repmat(oppy(1:n)',1,n);
Z2=cutoffz*ones(n);
surf(X2,Y2,Z2,'edgecolor','none','facecolor','b')
alpha(0.3)

% in XY-space
%green plane "vertical" on y axis
oppx=[upLimx:-0.01:lowLimx];
n=length(oppx);
h=(upLimz-lowLimz)/n;
oppz=[lowLimz:h:upLimz];
X1=repmat(oppx,n,1);
Y1=cutoffyy*ones(n);
Z1=repmat([oppz(1:n-1)'; upLimz],1,n);
surf(X1,Y1,Z1,'edgecolor','none','facecolor','g')
alpha(0.3)

if multi==0 %univariate case
    surf(X2,Y2,-Z2,'edgecolor','none','facecolor','b')
    alpha(0.3)
end

if labels~=0
    putlabel(x,y,labels,z,labels)
else
    plotnumbers(x,y,labsd,labod,5,z,labresd)
end
hold off

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit