text:
function [z,s,exitflag] = smoothn(varargin)

%SMOOTHN Robust spline smoothing for 1-D to N-D data.
%   SMOOTHN provides a fast, automatized and robust discretized smoothing
%   spline for data of any dimension.
%
%   Z = SMOOTHN(Y) automatically smoothes the uniformly-sampled array Y. Y
%   can be any N-D noisy array (time series, images, 3D data,...). Non
%   finite data (NaN or Inf) are treated as missing values.
%
%   Z = SMOOTHN(Y,S) smoothes the array Y using the smoothing parameter S.
%   S must be a real positive scalar. The larger S is, the smoother the
%   output will be. If the smoothing parameter S is omitted (see previous
%   option) or empty (i.e. S = []), it is automatically determined using
%   the generalized cross-validation (GCV) method.
%
%   Z = SMOOTHN(Y,W) or Z = SMOOTHN(Y,W,S) specifies a weighting array W of
%   real positive values, that must have the same size as Y. Note that a
%   nil weight corresponds to a missing value.
%
%   Robust smoothing
%   ----------------
%   Z = SMOOTHN(...,'robust') carries out a robust smoothing that minimizes
%   the influence of outlying data.
%
%   [Z,S] = SMOOTHN(...) also returns the calculated value for S so that
%   you can fine-tune the smoothing subsequently if needed.
%
%   An iteration process is used in the presence of weighted and/or missing
%   values. Z = SMOOTHN(...,OPTION_NAME,OPTION_VALUE) smoothes with the
%   termination parameters specified by OPTION_NAME and OPTION_VALUE. They
%   can contain the following criteria:
%       -----------------
%       TolZ:       Termination tolerance on Z (default = 1e-3)
%                   TolZ must be in ]0,1[
%       MaxIter:    Maximum number of iterations allowed (default = 100)
%       Initial:    Initial value for the iterative process (default =
%                   initial data)
%       -----------------
%   Syntax: [Z,...] = SMOOTHN(...,'MaxIter',500,'TolZ',1e-4,'Initial',Z0);
%
%   [Z,S,EXITFLAG] = SMOOTHN(...) returns a boolean value EXITFLAG that
%   describes the exit condition of SMOOTHN:
%       1       SMOOTHN converged.
%       0       Maximum number of iterations was reached.
%
%   Class Support
%   -------------
%   Input array can be numeric or logical. The returned array is of class
%   double.
%
%   Notes
%   -----
%   The N-D (inverse) discrete cosine transform functions DCTN and IDCTN are required.
%
%   To be made
%   ----------
%   Estimate the confidence bands (see Wahba 1983, Nychka 1988).
%
%   Reference
%   ---------
%   Garcia D, Robust smoothing of gridded data in one and higher dimensions
%   with missing values. Computational Statistics & Data Analysis, 2010.
%   PDF download
%
%   Examples:
%   --------
%   % 1-D example
%   x = linspace(0,100,2^8);
%   y = cos(x/10)+(x/50).^2 + randn(size(x))/10;
%   y([70 75 80]) = [5.5 5 6];
%   z = smoothn(y); % Regular smoothing
%   zr = smoothn(y,'robust'); % Robust smoothing
%   subplot(121), plot(x,y,'r.',x,z,'k','LineWidth',2)
%   axis square, title('Regular smoothing')
%   subplot(122), plot(x,y,'r.',x,zr,'k','LineWidth',2)
%   axis square, title('Robust smoothing')
%
%   % 2-D example
%   xp = 0:.02:1;
%   [x,y] = meshgrid(xp);
%   f = exp(x+y) + sin((x-2*y)*3);
%   fn = f + randn(size(f))*0.5;
%   fs = smoothn(fn);
%   subplot(121), surf(xp,xp,fn), zlim([0 8]), axis square
%   subplot(122), surf(xp,xp,fs), zlim([0 8]), axis square
%
%   % 2-D example with missing data
%   n = 256;
%   y0 = peaks(n);
%   y = y0 + rand(size(y0))*2;
%   I = randperm(n^2);
%   y(I(1:n^2*0.5)) = NaN; % lose 1/2 of data
%   y(40:90,140:190) = NaN; % create a hole
%   z = smoothn(y); % smooth data
%   subplot(2,2,1:2), imagesc(y), axis equal off
%   title('Noisy corrupt data')
%   subplot(223), imagesc(z), axis equal off
%   title('Recovered data ...')
%   subplot(224), imagesc(y0), axis equal off
%   title('... compared with original data')
%
%   % 3-D example
%   [x,y,z] = meshgrid(-2:.2:2);
%   xslice = [-0.8,1]; yslice = 2; zslice = [-2,0];
%   vn = x.*exp(-x.^2-y.^2-z.^2) + randn(size(x))*0.06;
%   subplot(121), slice(x,y,z,vn,xslice,yslice,zslice,'cubic')
%   title('Noisy data')
%   v = smoothn(vn);
%   subplot(122), slice(x,y,z,v,xslice,yslice,zslice,'cubic')
%   title('Smoothed data')
%
%   % Cardioid
%   t = linspace(0,2*pi,1000);
%   x = 2*cos(t).*(1-cos(t)) + randn(size(t))*0.1;
%   y = 2*sin(t).*(1-cos(t)) + randn(size(t))*0.1;
%   z = smoothn(complex(x,y));
%   plot(x,y,'r.',real(z),imag(z),'k','linewidth',2)
%   axis equal tight
%
%   % Cellular vortical flow
%   [x,y] = meshgrid(linspace(0,1,24));
%   Vx = cos(2*pi*x+pi/2).*cos(2*pi*y);
%   Vy = sin(2*pi*x+pi/2).*sin(2*pi*y);
%   Vx = Vx + sqrt(0.05)*randn(24,24); % adding Gaussian noise
%   Vy = Vy + sqrt(0.05)*randn(24,24); % adding Gaussian noise
%   I = randperm(numel(Vx));
%   Vx(I(1:30)) = (rand(30,1)-0.5)*5; % adding outliers
%   Vy(I(1:30)) = (rand(30,1)-0.5)*5; % adding outliers
%   Vx(I(31:60)) = NaN; % missing values
%   Vy(I(31:60)) = NaN; % missing values
%   Vs = smoothn(complex(Vx,Vy),'robust'); % automatic smoothing
%   subplot(121), quiver(x,y,Vx,Vy,2.5), axis square
%   title('Noisy velocity field')
%   subplot(122), quiver(x,y,real(Vs),imag(Vs)), axis square
%   title('Smoothed velocity field')
%
%   See also SMOOTH, GSMOOTHN, DCTN, IDCTN.
%
%   -- Damien Garcia -- 2009/03, revised 2010/06
%   Visit my website for more details about SMOOTHN

% Check input arguments
error(nargchk(1,10,nargin));

%% Test & prepare the variables
%---
k = 0;
while k<nargin && ~ischar(varargin{k+1}), k = k+1; end
%---
% y = array to be smoothed
y = double(varargin{1});
sizy = size(y);
noe = prod(sizy); % number of elements
if noe<2, z = y; return, end
%---
% Smoothness parameter and weights
W = ones(sizy);
s = [];
if k==2
    if isempty(varargin{2}) || isscalar(varargin{2}) % smoothn(y,s)
        s = varargin{2}; % smoothness parameter
    else % smoothn(y,W)
        W = varargin{2}; % weight array
    end
elseif k==3 % smoothn(y,W,s)
        W = varargin{2}; % weight array
        s = varargin{3}; % smoothness parameter
end
if ~isequal(size(W),sizy)
        error('MATLAB:smoothn:SizeMismatch',...
            'Arrays for data and weights must have same size.')
elseif ~isempty(s) && (~isscalar(s) || s<0)
    error('MATLAB:smoothn:IncorrectSmoothingParameter',...
        'The smoothing parameter must be a scalar >=0')
end
%---
% "Maximal number of iterations" criterion
I = find(strcmpi(varargin,'MaxIter'),1);
if isempty(I)
    MaxIter = 100; % default value for MaxIter
else
    try
        MaxIter = varargin{I+1};
    catch
        error('MATLAB:smoothn:IncorrectMaxIter',...
            'MaxIter must be an integer >=1')
    end
    if ~isnumeric(MaxIter) || ~isscalar(MaxIter) ||...
            MaxIter<1 || MaxIter~=round(MaxIter)
        error('MATLAB:smoothn:IncorrectMaxIter',...
            'MaxIter must be an integer >=1')
    end
end
%---
% "Tolerance on smoothed output" criterion
I = find(strcmpi(varargin,'TolZ'),1);
if isempty(I)
    TolZ = 1e-3; % default value for TolZ
else
    try
        TolZ = varargin{I+1};
    catch
        error('MATLAB:smoothn:IncorrectTolZ',...
            'TolZ must be in ]0,1[')
    end
    if ~isnumeric(TolZ) || ~isscalar(TolZ) || TolZ<=0 || TolZ>=1
        error('MATLAB:smoothn:IncorrectTolZ',...
            'TolZ must be in ]0,1[')
    end
end
%---
% "Initial Guess" criterion
I = find(strcmpi(varargin,'Initial'),1);
if isempty(I)
    isinitial = false; % default value for TolZ
else
    isinitial = true;
    try
        z0 = varargin{I+1};
    catch
        error('MATLAB:smoothn:IncorrectInitialGuess',...
            'Z0 must be a valid initial guess for Z')
    end
    if ~isnumeric(z0) || ~isequal(size(z0),sizy)
        error('MATLAB:smoothn:IncorrectTolZ',...
            'Z0 must be a valid initial guess for Z')
    end
end
%---
% Weights. Zero weights are assigned to not finite values (Inf or NaN),
% (Inf/NaN values = missing data).
IsFinite = isfinite(y);
nof = nnz(IsFinite); % number of finite elements
W = W.*IsFinite;
if any(W<0)
    error('MATLAB:smoothn:NegativeWeights',...
        'Weights must all be >=0')
else
    W = W/max(W(:));
end
%---
% Weighted or missing data?
isweighted = any(W(:)<1);
%---
% Robust smoothing?
isrobust = any(strcmpi(varargin,'robust'));
%---
% Automatic smoothing?
isauto = isempty(s);
%---
% DCTN and IDCTN are required
test4DCTNandIDCTN

%% Creation of the Lambda tensor
%---
% Lambda contains the eingenvalues of the difference matrix used in this
% penalized least squares process.
d = ndims(y);
Lambda = zeros(sizy);
for i = 1:d
    siz0 = ones(1,d);
    siz0(i) = sizy(i);
    Lambda = bsxfun(@plus,Lambda,...
        cos(pi*(reshape(1:sizy(i),siz0)-1)/sizy(i)));
end
Lambda = -2*(d-Lambda);
if ~isauto, Gamma = 1./(1+s*Lambda.^2); end

%% Upper and lower bound for the smoothness parameter
% The average leverage (h) is by definition in [0 1]. Weak smoothing occurs
% if h is close to 1, while over-smoothing appears when h is near 0. Upper
% and lower bounds for h are given to avoid under- or over-smoothing. See
% equation relating h to the smoothness parameter (Equation #12 in the
% referenced CSDA paper).
N = sum(sizy~=1); % tensor rank of the y-array
hMin = 1e-6; hMax = 0.99;
sMinBnd = (((1+sqrt(1+8*hMax.^(2/N)))/4./hMax.^(2/N)).^2-1)/16;
sMaxBnd = (((1+sqrt(1+8*hMin.^(2/N)))/4./hMin.^(2/N)).^2-1)/16;

%% Initialize before iterating
%---
Wtot = W;
%--- Initial conditions for z
if isweighted
    %--- With weighted/missing data
    % An initial guess is provided to ensure faster convergence. For that
    % purpose, a nearest neighbor interpolation followed by a coarse
    % smoothing are performed.
    %---
    if isinitial % an initial guess (z0) has been provided
        z = z0;
    else
        z = InitialGuess(y,IsFinite);
    end

else
    z = zeros(sizy);
end
%---
z0 = z;
y(~IsFinite) = 0; % arbitrary values for missing y-data
%---
tol = 1;
RobustIterativeProcess = true;
RobustStep = 1;
nit = 0;
%--- Error on p. Smoothness parameter s = 10^p
errp = 0.1;
opt = optimset('TolX',errp);
%--- Relaxation factor RF: to speedup convergence
RF = 1 + 0.75*isweighted;

%% Main iterative process
%---
while RobustIterativeProcess
    %--- "amount" of weights (see the function GCVscore)
    aow = sum(Wtot(:))/noe; % 0 < aow <= 1
    %---
    while tol>TolZ && nit<MaxIter
        nit = nit+1;
        DCTy = dctn(Wtot.*(y-z)+z);
        if isauto && ~rem(log2(nit),1)
            %---
            % The generalized cross-validation (GCV) method is used.
            % We seek the smoothing parameter s that minimizes the GCV
            % score i.e. s = Argmin(GCVscore)^.
            % Because this process is time-consuming, it is performed from
            % time to time (when nit is a power of 2)
            %---
            fminbnd(@gcv,log10(sMinBnd),log10(sMaxBnd),opt);
        end
        z = RF*idctn(Gamma.*DCTy) + (1-RF)*z;

        % if no weighted/missing data => tol=0 (no iteration)
        tol = isweighted*norm(z0(:)-z(:))/norm(z(:));

        z0 = z; % re-initialization
    end
    exitflag = nit<MaxIter;

    if isrobust %-- Robust Smoothing: iteratively re-weighted process
        %--- average leverage
        h = sqrt(1+16*s); h = sqrt(1+h)/sqrt(2)/h; h = h^N;
        %--- take robust weights into account
        Wtot = W.*RobustWeights(y-z,IsFinite,h);
        %--- re-initialize for another iterative weighted process
        isweighted = true; tol = 1; nit = 0;
        %---
        RobustStep = RobustStep+1;
        RobustIterativeProcess = RobustStep<4; % 3 robust steps are enough.
    else
        RobustIterativeProcess = false; % stop the whole process
    end
end

z = cast(z, 'like', varargin{1});

%% Warning messages
%---
if isauto
    if abs(log10(s)-log10(sMinBnd))<errp
        warning('MATLAB:smoothn:SLowerBound',...
            ['s = ' num2str(s,'%.3e') ': the lower bound for s ',...
            'has been reached. Put s as an input variable if required.'])
    elseif abs(log10(s)-log10(sMaxBnd))<errp
        warning('MATLAB:smoothn:SUpperBound',...
            ['s = ' num2str(s,'%.3e') ': the upper bound for s ',...
            'has been reached. Put s as an input variable if required.'])
    end
end
if nargout<3 && ~exitflag
    warning('MATLAB:smoothn:MaxIter',...
        ['Maximum number of iterations (' int2str(MaxIter) ') has ',...
        'been exceeded. Increase MaxIter option or decrease TolZ value.'])
end

%% GCV score
%---
function GCVscore = gcv(p)
    % Search the smoothing parameter s that minimizes the GCV score
    %---
    s = 10^p;
    Gamma = 1./(1+s*Lambda.^2);
    %--- RSS = Residual sum-of-squares
    if aow>0.9 % aow = 1 means that all of the data are equally weighted
        % very much faster: does not require any inverse DCT
        RSS = norm(DCTy(:).*(Gamma(:)-1))^2;
    else
        % take account of the weights to calculate RSS:
        yhat = idctn(Gamma.*DCTy);
        RSS = norm(sqrt(Wtot(IsFinite)).*(y(IsFinite)-yhat(IsFinite)))^2;
    end
    %---
    TrH = sum(Gamma(:));
    GCVscore = RSS/nof/(1-TrH/noe)^2;
end

end

%% Robust weights
function W = RobustWeights(r,I,h)
    % weights for robust smoothing.
    MAD = median(abs(r(I)-median(r(I)))); % median absolute deviation
    u = abs(r/(1.4826*MAD)/sqrt(1-h)); % studentized residuals
    c = 4.685; W = (1-(u/c).^2).^2.*((u/c)<1); % bisquare weights
    % c = 2.385; W = 1./(1+(u/c).^2); % Cauchy weights
    % c = 2.795; W = u<c; % Talworth weights
    W(isnan(W)) = 0;
end

%% Test for DCTN and IDCTN
function test4DCTNandIDCTN
    if ~exist('dctn','file')
        error('MATLAB:smoothn:MissingFunction',...
            ['DCTN and IDCTN are required. Download DCTN here.'])
    elseif ~exist('idctn','file')
        error('MATLAB:smoothn:MissingFunction',...
            ['DCTN and IDCTN are required. Download IDCTN here.'])
    end
end

%% Initial Guess with weighted/missing data
function z = InitialGuess(y,I)
    %-- nearest neighbor interpolation (in case of missing values)
    if any(~I(:))
        if license('test','image_toolbox')
            [z,L] = bwdist(I);
            z = y;
            z(~I) = y(L(~I));
        else
        % If BWDIST does not exist, NaN values are all replaced with the
        % same scalar. The initial guess is not optimal and a warning
        % message thus appears.
            z = y;
            z(~I) = mean(y(I));
            warning('MATLAB:smoothn:InitialGuess',...
                ['BWDIST (Image Processing Toolbox) does not exist. ',...
                'The initial guess may not be optimal; additional',...
                ' iterations can thus be required to ensure complete',...
                ' convergence. Increase ''MaxIter'' criterion if necessary.'])
        end
    else
        z = y;
    end
    %-- coarse fast smoothing using one-tenth of the DCT coefficients
    siz = size(z);
    z = dctn(z);
    for k = 1:ndims(z)
        z(ceil(siz(k)/10)+1:end,:) = 0;
        z = reshape(z,circshift(siz,[0 1-k]));
        z = shiftdim(z,1);
    end
    z = idctn(z);
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function test_BartholomewBiggs

  addpath('../lib') ;

  auxdata = {} ;

  options.lb = [ 1, 1, 1, 1 ] ; % Lower bound on the variables.
  options.ub = [ 5, 5, 5, 5 ] ; % Upper bound on the variables.

  % The constraint functions are bounded to zero
  options.cl = [0,   0]; %  constraints
  options.cu = [inf, 0];

  % Set up the auxiliary data.
  options.auxdata = auxdata ;

  % Set the IPOPT options.
  options.ipopt.jac_d_constant   = 'no';
  options.ipopt.hessian_constant = 'no';
  options.ipopt.mu_strategy      = 'adaptive';
  options.ipopt.max_iter         = 400;
  options.ipopt.tol              = 1e-10;
  %options.ipopt.linear_solver    = 'ma57';
  options.ipopt.linear_solver    = 'mumps';
  %options.ipopt.linear_solver    = 'pardiso';

  % The callback functions.
  funcs.objective         = @objective;
  funcs.constraints       = @constraints;
  funcs.gradient          = @gradient;
  funcs.jacobian          = @jacobian;
  funcs.jacobianstructure = @jacobianstructure;
  if true
    funcs.hessian           = @hessian;
    funcs.hessianstructure  = @hessianstructure;
    options.ipopt.derivative_test = 'second-order';
  else
    options.ipopt.hessian_approximation      = 'limited-memory';
    %options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
    %options.ipopt.limited_memory_update_type = 'sr1' ;
    options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
  end

  % Run IPOPT.
  x0 = [ 1, 5, 5, 1 ] ;

  tic
  [x, info] = ipopt_auxdata(x0,funcs,options);
  elapsed = toc ;

  info;

  x

end

%%
% map the indices with the corresponding index in the spase matrix
function f = objective( x, auxdata )
  x1 = x(1) ;
  x2 = x(2) ;
  x3 = x(3) ;
  x4 = x(4) ;
  f  = x1*x4*(x1+x2+x3) + x3 ;
end

%%
% map the indices with the corresponding index in the spase matrix
function g = gradient(x,auxdata)
  x1 = x(1) ;
  x2 = x(2) ;
  x3 = x(3) ;
  x4 = x(4) ;
  g  = [ x4*(2*x1+x2+x3), x1*x4, x1*x4 + 1, x1*(x1+x2+x3)] ;
end

function c = constraints(x,auxdata)
  x1 = x(1) ;
  x2 = x(2) ;
  x3 = x(3) ;
  x4 = x(4) ;
  c  = [ x1*x2*x3*x4 ; x1^2+x2^2+x3^2+x4^2 - 40 ] ;
end

function jac = jacobian(x,auxdata)
  x1 = x(1) ;
  x2 = x(2) ;
  x3 = x(3) ;
  x4 = x(4) ;
  jac = sparse([ x2*x3*x4, x1*x3*x4, x1*x2*x4, x1*x2*x3 ; ...
                 2*x1, 2*x2, 2*x3, 2*x4 ]) ;
end

function jac = jacobianstructure(auxdata)
  jac = sparse(ones(2,4)) ;
end

function H = hessian(x, sigma, lambda, auxdata)
  x1 = x(1) ;
  x2 = x(2) ;
  x3 = x(3) ;
  x4 = x(4) ;
  Hf = [ 2*x4,      x4,  x4, 2*x1+x2+x3 ; ...
         x4,         0,   0, x1         ; ...
         x4,         0,   0, x1         ; ...
         2*x1+x2+x3, x1, x1, x1         ] ;

  H1 = [ 0, x3*x4, x2*x4, x2*x3 ; ...
         x3*x4, 0, x1*x4, x1*x3 ; ...
         x2*x4, x1*x4, 0, x1*x2 ; ...
         x2*x3, x1*x3, x1*x2, 0 ] ;

     H2 = 2*eye(4,4) ;
  H  = tril(sparse(sigma * Hf + lambda(1)*H1 + lambda(2)*H2)) ;
end

function H = hessianstructure(auxdata)
  H = tril(sparse(ones(4,4))) ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [selection, value] = seldellistdlg(varargin)
%LISTDLG  List selection dialog box.
%   [SELECTION,OK] = LISTDLG('ListString',S) creates a modal dialog box
%   which allows you to select a string or multiple strings from a list.
%   SELECTION is a vector of indices of the selected strings (length 1 in
%   the single selection mode).  This will be [] when OK is 0.  OK is 1 if
%   you push the OK button, or 0 if you push the Cancel button or close the
%   figure.
%
%   Double-clicking on an item or pressing <CR> when multiple items are
%   selected has the same effect as clicking the OK button.  Pressing <CR>
%   is the same as clicking the OK button. Pressing <ESC> is the same as
%   clicking the Cancel button.
%
%   Inputs are in parameter,value pairs:
%
%   Parameter       Description
%   'ListString'    cell array of strings for the list box.
%   'SelectionMode' string; can be 'single' or 'multiple'; defaults to
%                   'multiple'.
%   'ListSize'      [width height] of listbox in pixels; defaults
%                   to [160 300].
%   'InitialValue'  vector of indices of which items of the list box
%                   are initially selected; defaults to the first item.
%   'Name'          String for the figure's title; defaults to ''.
%   'PromptString'  string matrix or cell array of strings which appears
%                   as text above the list box; defaults to {}.
%   'OKString'      string for the OK button; defaults to 'OK'.
%   'CancelString'  string for the Cancel button; defaults to 'Cancel'.
%
%   A 'Select all' button is provided in the multiple selection case.
%
%   Example:
%     d = dir;
%     str = {d.name};
%     [s,v] = listdlg('PromptString','Select a file:',...
%                     'SelectionMode','single',...
%                     'ListString',str)
 %
%  See also DIALOG, ERRORDLG, HELPDLG, INPUTDLG,
%    MSGBOX, QUESTDLG, WARNDLG.

%   Copyright 1984-2005 The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 2008/05/22 22:20:02 $

%   'uh'            uicontrol button height, in pixels; default = 22.
%   'fus'           frame/uicontrol spacing, in pixels; default = 8.
%   'ffs'           frame/figure spacing, in pixels; default = 8.

% simple test:
%
% d = dir; [s,v] = listdlg('PromptString','Select a file:','ListString',{d.name});
%
error(nargchk(1,inf,nargin))

figname = '';
smode = 2;   % (multiple)
promptstring = {};
liststring = [];
listsize = [160 300];
initialvalue = [];
okstring = 'OK';
cancelstring = 'Cancel';
fus = 8;
ffs = 8;
uh = 22;

if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error('MATLAB:listdlg:InvalidArgument', 'Arguments to LISTDLG must come param/value in pairs.')
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
     case 'name'
      figname = varargin{i+1};
     case 'promptstring'
      promptstring = varargin{i+1};
     case 'selectionmode'
      switch lower(varargin{i+1})
       case 'single'
        smode = 1;
       case 'multiple'
        smode = 2;
      end
     case 'listsize'
      listsize = varargin{i+1};
     case 'liststring'
      liststring = varargin{i+1};
     case 'initialvalue'
      initialvalue = varargin{i+1};
     case 'uh'
      uh = varargin{i+1};
     case 'fus'
      fus = varargin{i+1};
     case 'ffs'
      ffs = varargin{i+1};
     case 'okstring'
      okstring = varargin{i+1};
     case 'cancelstring'
      cancelstring = varargin{i+1};
     otherwise
      error('MATLAB:listdlg:UnknownParameter', ['Unknown parameter name passed to LISTDLG.  Name was ' varargin{i}])
    end
end

if ischar(promptstring)
    promptstring = cellstr(promptstring);
end

if isempty(initialvalue)
    initialvalue = 1;
end

if isempty(liststring)
    error('MATLAB:listdlg:NeedParameter', 'ListString parameter is required.')
end

ex = get(0,'defaultuicontrolfontsize')*1.7;  % height extent per line of uicontrol text (approx)

fp = get(0,'defaultfigureposition');
w = 2*(fus+ffs)+listsize(1);
h = 2*ffs+6*fus+ex*length(promptstring)+listsize(2)+uh+(smode==2)*(fus+uh);
fp = [fp(1) fp(2)+fp(4)-h w h];  % keep upper left corner fixed

fig_props = { ...
    'name'                   figname ...
    'color'                  get(0,'defaultUicontrolBackgroundColor') ...
    'resize'                 'off' ...
    'numbertitle'            'off' ...
    'menubar'                'none' ...
    'windowstyle'            'modal' ...
    'visible'                'off' ...
    'createfcn'              ''    ...
    'position'               fp   ...
    'closerequestfcn'        'delete(gcbf)' ...
            };

liststring=cellstr(liststring);

fig = figure(fig_props{:});

if length(promptstring)>0
    prompt_text = uicontrol('style','text','string',promptstring,...
        'horizontalalignment','left',...
        'position',[ffs+fus fp(4)-(ffs+fus+ex*length(promptstring)) ...
        listsize(1) ex*length(promptstring)]); %#ok
end

btn_wid = (fp(3)-2*(ffs+fus)-fus)/2;

listbox = uicontrol('style','listbox',...
                    'position',[ffs+fus ffs+uh+4*fus+(smode==2)*(fus+uh) listsize],...
                    'string',liststring,...
                    'backgroundcolor','w',...
                    'max',smode,...
                    'tag','listbox',...
                    'value',initialvalue, ...
                    'callback', {@doListboxClick});

ok_btn = uicontrol('style','pushbutton',...
                   'string',okstring,...
                   'position',[ffs+fus ffs+fus btn_wid uh],...
                   'callback',{@doOK,listbox});

cancel_btn = uicontrol('style','pushbutton',...
                       'string',cancelstring,...
                       'position',[ffs+2*fus+btn_wid ffs+fus btn_wid uh],...
                       'callback',{@doCancel,listbox});

if smode == 2
    selectall_btn = uicontrol('style','pushbutton',...
                              'string','Select all',...
                              'position',[ffs+fus 4*fus+ffs+uh listsize(1) uh],...
                              'tag','selectall_btn',...
                              'callback',{@doSelectAll, listbox});

    if length(initialvalue) == length(liststring)
        set(selectall_btn,'enable','off')
    end
    set(listbox,'callback',{@doListboxClick, selectall_btn})
end

set([fig, ok_btn, cancel_btn, listbox], 'keypressfcn', {@doKeypress, listbox});

set(fig,'position',getdialoglocation(fp, get(fig,'Units')));
% Make ok_btn the default button.
%setdefaultbutton(fig, ok_btn);

% make sure we are on screen
movegui(fig)
set(fig, 'visible','on'); drawnow;

try
    % Give default focus to the listbox *after* the figure is made visible
    uicontrol(listbox);
    uiwait(fig);
catch
    if ishandle(fig)
        delete(fig)
    end
end

if isappdata(0,'ListDialogAppData__')
    ad = getappdata(0,'ListDialogAppData__');
    selection = ad.selection;
    value = ad.value;
    rmappdata(0,'ListDialogAppData__')
else
    % figure was deleted
    selection = [];
    value = 0;
end

%% figure, OK and Cancel KeyPressFcn
function doKeypress(src, evd, listbox) %#ok
switch evd.Key
 case 'escape'
  doCancel([],[],listbox);
end

%% OK callback
function doOK(ok_btn, evd, listbox) %#ok
if (~isappdata(0, 'ListDialogAppData__'))
    ad.value = 1;
    ad.selection = get(listbox,'value');
    setappdata(0,'ListDialogAppData__',ad);
    delete(gcbf);
end

%% Cancel callback
function doCancel(cancel_btn, evd, listbox) %#ok
ad.value = 0;
ad.selection = get(listbox,'value');
setappdata(0,'ListDialogAppData__',ad)
delete(gcbf);

%% SelectAll callback
function doSelectAll(selectall_btn, evd, listbox) %#ok
set(selectall_btn,'enable','off')
set(listbox,'value',1:length(get(listbox,'string')));

%% Listbox callback
function doListboxClick(listbox, evd, selectall_btn) %#ok
% if this is a doubleclick, doOK
if strcmp(get(gcbf,'SelectionType'),'open')
    %doOK([],[],listbox);
elseif nargin == 3
    if length(get(listbox,'string'))==length(get(listbox,'value'))
        set(selectall_btn,'enable','off')
    else
        set(selectall_btn,'enable','on')
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [lat, lon, sphrad] = proj_on_sphere(x)
% PROJ_ON_SPHERE  Project a point configuration onto a sphere
%
% [LAT, LON, SPHRAD] = proj_on_sphere(X)
%
%   X is a 3-row matrix. Each colum has the (x,y,z)-coordinates of a point.
%
%   LAT, LON are vectors with the latitude and longitude of the points in X
%   projected onto a sphere. The centre of the sphere is the centroid of X.
%
%   SPHRAD is the median of the distance of the points in X to the
%   centroid. It can be used as the radius of the sphere that best contains
%   the projection of X in terms of distances between the the points.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 1);
nargoutchk(0, 3);

if (size(x, 2) ~= 3)
    error('X must be a matrix with 3-columns')
end

% project 3D Cartesian coordinates onto
% the surface of a sphere lon, lat given in radians, centered around 0
[lon, lat, r] = cart2sph(...
    x(:, 1) - mean(x(:, 1)), ...
    x(:, 2) - mean(x(:, 2)), ...
    x(:, 3) - mean(x(:, 3)));

% we use the median radius of the points in the configuration as the radius
% of the sphere that best can contains their projections
sphrad = median(r);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% corrmap_createclust() - creates a new cluster in the STUDY structure
%                         containing the ICs selected and stored using
%                         corrmap() function.
%
% Usage:
%          >> STUDY = corrmap_createclust(STUDY,CORRMAP)
%
% Inputs:
%  STUDY - input STUDY structure
%  ALLEEG - input ALLEEG structure
%
% Outputs:
%  STUDY - the input STUDY set structure modified with the new cluster.
%
% See also:  pop_corrmap(), corrmap()
%
% Author: F. Campos Viola, MRC-IHR, 19/10/2007

%123456789012345678901234567890123456789012345678901234567890123456789012

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

% revised by F Campos-Viola - corrmap1.01 (30/01/2009)

function STUDY = corrmap_createclust(STUDY,CORRMAP)

clusters = [];

cls = length(STUDY.cluster)+1;% number of cluster
nc  = length(STUDY.cluster); % index of last cluster

%%%%% CREATING NEW CLUSTER %%%%%

STUDY.cluster(cls).name = CORRMAP.input.clname; %cluster's name

%auxiliary variables to create matrix with indices - used to access
%correct preclustdata

% saving info int the STUDY.cluster fields
STUDY.cluster(cls).sets  = CORRMAP.corr.sets{2}(1:CORRMAP.clust.ics(2))';
STUDY.cluster(cls).comps = CORRMAP.corr.ics{2}(1:CORRMAP.clust.ics(2))';
STUDY.cluster(cls).algorithm = 'correlation (CORRMAP)';
STUDY.cluster(cls).parent = 'ParentCluster 1';
STUDY.cluster(cls).child = [];
STUDY.cluster(cls).centroid = [];
%STUDY.cluster(cls).preclust.preclustdata = STUDY.etc.preclust.preclustdata(ind,:);
STUDY.cluster(cls).preclust.preclustdata = STUDY.etc.preclust.preclustdata
STUDY.cluster(cls).preclust.preclustparams = STUDY.etc.preclust.preclustparams;
STUDY.cluster(cls).preclust.preclustcomps = STUDY.etc.preclust.preclustcomps;

%update parents clusters with cluster child indices
STUDY.cluster(1).child{cls-1} = CORRMAP.input.clname;

% clusters = [clusters 1:cls];%the new created clusters indices.

fprintf('>');
fprintf('a new cluster was created. \n'); %info for the user

author: PatternRecognition

repo: OpenBMI


text:
function vis_multicut(expidx, bRecompute, firstidx, nImgs, cidxs)

p = exp_params(expidx);
exp_dir = [p.expDir '/' p.shortName];
fprintf('exp_dir: %s\n',exp_dir);

resDir = [exp_dir '/multicut/'];
visDir = [exp_dir '/multicut/vis/'];

if (isfield(p,'testGTnopad'))
    load(p.testGTnopad,'annolist');
    bProject = true;
else
    load(p.testGT,'annolist');
    bProject = false;
end

if (nargin < 2)
    bRecompute = true;
end

if (nargin < 3)
    firstidx = 1;
end

if (nargin < 4)
    nImgs = length(annolist);
end

if (nargin < 5)
    cidxs = p.cidxs;
end

lastidx = firstidx + nImgs - 1;

colors = {'r','g','b','c','m','y'};
markers = {'+','o','s'}; %,'x','.','-','s','d'
scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4) scrsz(3)/2 scrsz(4)]);

for imgidx = firstidx:lastidx
    out_img_fname = [visDir '/imgidx_' padZeros(num2str(imgidx),4) '_cidx_1_14' '.png'];
    if (exist(out_img_fname, 'file') == 2) && ~bRecompute
        continue;
    end

    try
        templ = [resDir '/imgidx_' padZeros(num2str(imgidx),4) '_*'];
        files = dir(templ);
        fname = fullfile(resDir, files(end).name);
        %fname = [resDir '/imgidx_' padZeros(num2str(imgidx),4) '_cidx_' num2str(cidxs(1)) '_' num2str(cidxs(end))];
        load(fname,'unLab','unPos');
        fprintf('imgidx: %d\n',imgidx);
        clf;
        subplot(2,1,1);
        imgname = annolist(imgidx).image.name;
        img = imread(imgname);
        imagesc(img); axis equal; hold on;

        keypoints.imgname = imgname;
        keypoints.det = unPos;

        if (bProject)
            keypoints = projectKeypoints(keypoints,fileparts(p.testGTnopad),fileparts(p.testGT));
        end

        cb = keypoints.det;

        for j = 1:size(unLab,1)
            lp = unLab(j,1);
            cb_cur = cb(j,:);
            if (lp < 1000)
                lp = cidxs(lp+1)-1;
                if (lp <= 5)
                    plot(cb_cur(1),cb_cur(2),[colors{mod(lp,6)+1} markers{1}],'MarkerSize',5);
                elseif (lp > 5 && lp <= 11)
                    plot(cb_cur(1),cb_cur(2),[colors{mod(lp,6)+1} markers{2}],'MarkerSize',5);
                else
                    plot(cb_cur(1),cb_cur(2),[colors{mod(lp,6)+1} markers{3}],'MarkerSize',5);
                end
            else
                plot(cb_cur(1),cb_cur(2),'kx','MarkerSize',5);
            end
        end
        axis off;

        subplot(2,1,2);
        imagesc(img); axis equal; hold on;
        idxs = find(unLab(:,1) < 1000);
        lp_uniq = unique(unLab(idxs,2));
        for j = 1:length(idxs)
            lp = find(lp_uniq == unLab(idxs(j),2))-1;
            cb_cur = cb(idxs(j),:);
            if (lp <= 5)
                plot(cb_cur(1),cb_cur(2),[colors{mod(lp,6)+1} markers{1}],'MarkerSize',5);
            elseif (lp > 5 && lp <= 11)
                plot(cb_cur(1),cb_cur(2),[colors{mod(lp,6)+1} markers{2}],'MarkerSize',5);
            else
                plot(cb_cur(1),cb_cur(2),[colors{mod(lp,6)+1} markers{3}],'MarkerSize',5);
            end
        end
        axis off;
        print(gcf,'-dpng',out_img_fname);
    catch
    end
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [mcs, mcs_time] = calculateMCS(model_struct, n_mcs, max_len_mcs, varargin)
% Calculate Minimal Cut Sets (MCSs) using the warm-start strategy available
% in CPLEX, namely cplex.populate(), with or without selecting a given
% knockout, among all the reactions included in the model or a given subset
% of them. Tobalina et al., 2016 (Bioinformatics); von Kamp and Klamt, 2014
% (PLoS Computational Biology).
%
% USAGE:
%
%    [mcs, mcs_time] = populateMCS(model_struct, n_mcs, max_len_mcs, options)
%
% INPUTS:
%    model_struct:    Metabolic model structure (COBRA Toolbox format).
%    n_mcs:           Number of MCSs to calculate.
%    max_len_mcs:     Number of reactions in the largest MCS to be calculated.
%
% OPTIONAL INPUT:
% OPTIONAL INPUTS:
%    KO:                 Selected reaction knockout. (default = [])
%    rxn_set:            Cell array containing the set of reactions among which
%                        the MCSs are wanted to be calculated.
%                        taken from model.rxns
%                        (default = [], all reactions)
%    target_b:           Desired activity level of the metabolic task to be
%                        disrupted. (default = 1e-3)
%    forceLength:        1 if the constraint limiting the length of the
%                        MCSs is to be active (recommended for
%                        enumerating low order MCSs), 0 otherwise
%                        (default = 1)
%    timelimit:          Time limit for the calculation of MCSs each time
%                        the solver is called. (default = 1e75)
%    numWorkers:         Integer: is the maximun number of workers used
%                        by Cplex and GPR2models. 0 = automatic,
%                        1 = sequential, > 1 = parallel. (default = 0)
%    printLevel:         Integer. 1 if the process is wanted to be shown
%                        on the screen, 0 otherwise. (default = 1)
%
% OUTPUTS:
%    mcs:         Cell array containing the calculated MCSs.
%    mcs_time:    Calculation times of the different processes in
%                 the algorithm.
%
% EXAMPLE:
%    %With optional values
%    [mcs, mcs_time] = calculateMCS(modelR204, 100, 10, ...
%                                    'KO', 'r1651', ...
%                                    'rxn_set', {'r1652'; 'r1653'; 'r1654'}, ...
%                                    'timelimit', 300, ...
%                                    'target_b', 1e-4, ...
%                                    'forceLength', 0, ...
%                                    'printLevel', 0, ...
%
%    %Without optional values
%    [mcs, mcs_time] = calculateMCS(model, 100, 10)
%
% .. Authors:
%       - Inigo Apaolaza, 30/01/2017, University of Navarra, TECNUN School of Engineering.
%       - Luis V. Valcarcel, 19/11/2017, University of Navarra, TECNUN School of Engineering.
%       - Francisco J. Planes, 20/11/2017, University of Navarra, TECNUN School of Engineering.
% .. Revisions:
%       - Inigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering.
%       - Luis V. Valcarcel, 17/04/2018, University of Navarra, TECNUN School of Engineering.
%       - Luis V. Valcarcel, 30/06/2021, University of Navarra, TECNUN School of Engineering.
%       - Luis V. Valcarcel, 16/09/2022, University of Navarra, TECNUN School of Engineering.

% Check the installation of cplex
global SOLVERS;
global CBT_MILP_SOLVER;
if SOLVERS.ibm_cplex.installed && SOLVERS.ibm_cplex.working
    if ~strcmp(CBT_MILP_SOLVER,'ibm_cplex')
        warning('calculateMCS will use IBM CPLEX although it is not selected for MILP')
    end
else
    error('This version calculateMCS only works with IBM CPLEX. Newer versions will include more solvers included in COBRA Toolbox')
end

time_aa = tic;
% Set Parameters
p = inputParser;
% check required arguments
addRequired(p, 'model_struct');
addRequired(p, 'n_mcs', @isnumeric);
addRequired(p, 'max_len_mcs', @isnumeric);
% Add optional name-value pair argument
addParameter(p, 'KO', [], @(x)ischar(x)||isempty(x));
addParameter(p, 'rxn_set', [], @(x)iscell(x)||isempty(x));
addParameter(p, 'target_b', 1e-3, @(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'timelimit', 1e75, @(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'forceLength', true, @(x)islogical(x)||(isnumeric(x)&&isscalar(x)));
addParameter(p, 'numWorkers', 0, @(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x));
% extract variables from parser
parse(p, model_struct, n_mcs, max_len_mcs, varargin{:});
model_struct = p.Results.model_struct;
n_mcs = p.Results.n_mcs;
max_len_mcs = p.Results.max_len_mcs;
KO = p.Results.KO;
rxn_set = p.Results.rxn_set;
target_b = p.Results.target_b;
timelimit = p.Results.timelimit;
forceLength = p.Results.forceLength;
numWorkers = p.Results.numWorkers;
printLevel = p.Results.printLevel;

integrality_tolerance = 1e-5;
M = 1e3;    % Big Value
alpha = 1;  % used to relate the lower bound of v variables with z variables
c = 1e-3;   % used to activate w variable
b = 1e-3;   % used to activate KnockOut constraint
phi = 1000; % b/c;

% Build the K Matrix
[~, n_ini_rxns] = size(model_struct.S);
K = speye(n_ini_rxns);

% Splitting
S = [model_struct.S -model_struct.S(:, model_struct.lb<0)];
K = [K K(:, model_struct.lb<0)];
K_ind = model_struct.rxns;
n_K_ind = length(K_ind);
[n_mets, n_rxns] = size(S);
nbio = find(model_struct.c);
t = zeros(n_rxns, 1);
t(nbio) = 1;

% Permit only KOs in rxn_set
if ~isempty(rxn_set)
    if ~isempty(KO)
        rxn_set = [rxn_set; {KO}];
    end
    rxn_set = unique(rxn_set);
    tmp_set = cellfun(@ismember, model_struct.rxns, repmat({rxn_set}, n_ini_rxns, 1), 'UniformOutput', false);
    pos_set = find(cell2mat(tmp_set));
    K = K(pos_set, :);
    K_ind = model_struct.rxns(pos_set);
    n_K_ind = length(K_ind);
end

if isempty(KO)
% ENUMERATE MCSs
% Define variables
    var.u = 1:n_mets;
    var.vp = var.u(end)+1:var.u(end)+n_K_ind;
    var.w = var.vp(end)+1:var.vp(end)+1;
    var.zp = var.w(end)+1:var.w(end)+n_K_ind;
    var.zw = var.zp(end)+1:var.zp(end)+1;
    n_vars = var.zw(end);
    var_group.v = [var.vp var.w];
    var_group.z = [var.zp var.zw];

% Define constraints
    cons.Ndual = 1:size(S, 2);
    cons.forceBioCons = cons.Ndual(end)+1:cons.Ndual(end)+1;
    cons.forceLength = cons.forceBioCons(end)+1:cons.forceBioCons(end)+1;
    n_cons = cons.forceLength(end);

% Cplex - A matrix
    A = sparse(zeros(n_cons, n_vars));
    A(cons.Ndual, var.u) = S';
    A(cons.Ndual, var.vp) = K';
    A(cons.Ndual, var.w) = -t;
    A(cons.forceBioCons, var.w) = -target_b;
    if forceLength == 1
        A(cons.forceLength, var.zp) = 1;
    end

% Cplex - rhs and lhs vectors
    rhs = zeros(n_cons, 1);
    rhs(cons.Ndual, 1) = inf;
    rhs(cons.forceBioCons) = -c;
    if forceLength == 1
        rhs(cons.forceLength) = 1;
    end
    lhs = zeros(n_cons, 1);
    lhs(cons.Ndual, 1) = 0;
    lhs(cons.forceBioCons) = -1000;
    if forceLength == 1
        lhs(cons.forceLength) = 1;
    end

% Cplex - ub and lb vectors
    ub(var.u, 1) = inf;
    ub(var.vp) = inf;
    ub(var.w) = inf;
    ub(var.zp) = 1;
    ub(var.zw) = 1;
    lb(var.u, 1) = -inf;
    lb(var.vp) = 0;
    lb(var.w) = 0;
    lb(var.zp) = 0;
    lb(var.zw) = 0;

% Cplex - obj vector
    obj(var.u, 1) = 0;
    obj(var.vp) = 0;
    obj(var.w) = 0;
    obj(var.zp) = 1;
    obj(var.zw) = 0;

% Cplex - ctype vector
    ctype(var.u) = 'C';
    ctype(var.vp) = 'C';
    ctype(var.w) = 'C';
    ctype(var.zp) = 'B';
    ctype(var.zw) = 'B';

% Cplex - sense of the optimization
    sense = 'minimize';

% Cplex - Introduce all data in a Cplex structure
    cplex = Cplex('MCS');
    Model = struct();
    [Model.A, Model.rhs, Model.lhs, Model.ub, Model.lb, Model.obj, Model.ctype, Model.sense] = deal(A, rhs, lhs, ub, lb, obj, ctype, sense);
    cplex.Model = Model;

% Cplex Indicators
    % z = 1  -->  v >= alpha
    for ivar = 1:length(var_group.z)
        a = zeros(n_vars, 1);
        a(var_group.v(ivar)) = 1;
        cplex.addIndicators(var_group.z(ivar), 0, a, 'G', alpha);
    end

% Cplex Indicators
    % z = 0  -->  v <= 0
    for ivar = 1:length(var_group.z)
        a = zeros(n_vars, 1);
        a(var_group.v(ivar)) = 1;
        cplex.addIndicators(var_group.z(ivar), 1, a, 'L', 0);
    end

% Cplex Parameters
    sP = struct();
    [sP.mip.tolerances.integrality, sP.mip.strategy.heuristicfreq, sP.mip.strategy.rinsheur] = deal(integrality_tolerance, 1000, 50);
    [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers);
    [sP.preprocessing.aggregator, sP.preprocessing.boundstrength, ...
        sP.preprocessing.coeffreduce, sP.preprocessing.dependency, ...
        sP.preprocessing.dual, sP.preprocessing.fill,...
        sP.preprocessing.linear, sP.preprocessing.numpass, ...
        sP.preprocessing.presolve, sP.preprocessing.reduce,..., ...
        sP.preprocessing.relax, sP.preprocessing.symmetry] = deal(50, 1, 2, 1, 1, 50, 1, 50, 1, 3, 1, 1);
    cplex = setCplexParam(cplex, sP);

    if printLevel == 0
        cplex.DisplayFunc = [];
    end

% Calculation of MCSs
    mcs_time{1, 1} = '------ TIMING ------';
    mcs_time{1, 2} = '--- MCSs ---';
    i = 0;
    k = 0;
    n_time = size(mcs_time, 1);
    mcs_time{n_time+1, 1} = 'Preparation';
    mcs_time{n_time+1, 2} = toc(time_aa);
    mcs = [];
    largest_mcs = 0;
    while largest_mcs <= max_len_mcs && k < n_mcs && cplex.Model.rhs(cons.forceLength) <= max_len_mcs
        ini_mcs_time = toc(time_aa);
        cplex.Param.mip.limits.populate.Cur = 40;
        cplex.Param.mip.pool.relgap.Cur = 0.1;
        cplex.populate();
        n_pool = size(cplex.Solution.pool.solution, 1);
        if n_pool ~= 0
            solution = cplex.Solution.pool.solution;
            for j = 1:n_pool
                k = k+1;
                mcs{k, 1} = K_ind((solution(j).x(var.zp))>0.9);
                n_cons = n_cons+1;
                sol = solution(j).x(var.zp)>0.9;
                cplex.Model.A(n_cons, var.zp) = sparse(double(sol));
                cplex.Model.rhs(n_cons) = sum(sol)-1;
                cplex.Model.lhs(n_cons) = 0;
            end
            i = i+1;
            mcsi_time = toc(time_aa)-ini_mcs_time;
            n_time = size(mcs_time, 1);
            mcs_time{n_time+1, 1} = ['POPULATE_ORDER_' num2str(cplex.Model.rhs(cons.forceLength))];
            mcs_time{n_time+1, 2} = mcsi_time;
        else
            mcsi_time = toc(time_aa)-ini_mcs_time;
            n_time = size(mcs_time, 1);
            mcs_time{n_time+1, 1} = ['POPULATE_ORDER_' num2str(cplex.Model.rhs(cons.forceLength)) 'NF'];
            mcs_time{n_time+1, 2} = mcsi_time;
            if forceLength == 1
                cplex.Model.rhs(cons.forceLength) = cplex.Model.rhs(cons.forceLength)+1;
                cplex.Model.lhs(cons.forceLength) = cplex.Model.lhs(cons.forceLength)+1;
            else
                n_time = size(mcs_time, 1);
                mcs_time{n_time+1, 1} = 'TOTAL MCSs';
                mcs_time{n_time+1, 2} = toc(time_aa);
                return;
            end
        end
        try save('tmp.mat', 'mcs', 'mcs_time'); end
        try largest_mcs = max(cellfun(@length, mcs)); end
    end
else
% CALCULATE MCSs WITH A GIVEN KNOCKOUT
% Select the row(s) in K_ind related to the KO under study
    tmp = repmat({KO}, n_K_ind, 1);
    dp = cellfun(@isequal, K_ind, tmp);

% Define variables
    var.u = 1:n_mets;
    var.vp = var.u(end)+1:var.u(end)+n_K_ind;
    var.w = var.vp(end)+1:var.vp(end)+1;
    var.zp = var.w(end)+1:var.w(end)+n_K_ind;
    var.zw = var.zp(end)+1:var.zp(end)+1;
    var.epsp = var.zw(end)+1:var.zw(end)+n_K_ind;
    var.epsw = var.epsp(end)+1:var.epsp(end)+1;
    var.delp = var.epsw(end)+1:var.epsw(end)+n_K_ind;
    var.delw = var.delp(end)+1:var.delp(end)+1;
    var.x = var.delw(end)+1:var.delw(end)+n_rxns+1;
    n_vars = var.x(end);
    var_group.v = [var.vp var.w];
    var_group.z = [var.zp var.zw];
    var_group.eps = [var.epsp var.epsw];
    var_group.del = [var.delp var.delw];

% Define constraints
    cons.Ndual = 1:size(S, 2);
    cons.forceBioCons = cons.Ndual(end)+1:cons.Ndual(end)+1;
    cons.forceKO = cons.forceBioCons(end)+1:cons.forceBioCons(end)+1;
    cons.linearComb = cons.forceKO(end)+1:cons.forceKO(end)+size(S, 1)+size(K, 1)+size(t, 2);
    cons.forceLength = cons.linearComb(end)+1:cons.linearComb(end)+1;
    n_cons = cons.forceLength(end);

% Cplex - A matrix
    A = sparse(zeros(n_cons, n_vars));
    A(cons.Ndual, var.u) = S';
    A(cons.Ndual, var.vp) = K';
    A(cons.Ndual, var.w) = -t;
    A(cons.forceBioCons, var.w) = -target_b;
    A(cons.forceKO, var.vp) = dp';
    A(cons.linearComb, var.x) = [S sparse(zeros(n_mets, 1)); K sparse(zeros(n_K_ind, 1)); -t' target_b];
    A(cons.linearComb, [var.epsp var.epsw]) = [sparse(zeros(n_mets, length(var.vp)+length(var.w))); -speye(length(var.vp)+length(var.w))];
    A(cons.linearComb, [var.delp var.delw]) = -[sparse(zeros(n_mets, length(var.vp)+length(var.w))); -speye(length(var.vp)+length(var.w))];
    if forceLength == 1
        A(cons.forceLength, var.zp) = 1;
    end

% Cplex - rhs and lhs vectors
    rhs = zeros(n_cons, 1);
    rhs(cons.Ndual, 1) = inf;
    rhs(cons.forceBioCons) = -c;
    rhs(cons.forceKO) = 10000;
    rhs(cons.linearComb) = [sparse(zeros(n_mets, 1)); dp; zeros(size(t, 2), 1)];
    if forceLength == 1
        rhs(cons.forceLength) = 1;
    end
    lhs = zeros(n_cons, 1);
    lhs(cons.Ndual, 1) = 0;
    lhs(cons.forceBioCons) = -1000;
    lhs(cons.forceKO) = b*10;
    lhs(cons.linearComb) = [sparse(zeros(n_mets, 1)); dp; zeros(size(t, 2), 1)];
    if forceLength == 1
        lhs(cons.forceLength) = 1;
    end

% Cplex - ub and lb vectors
    ub(var.u, 1) = inf;
    ub(var.vp) = inf;
    ub(var.w) = inf;
    ub(var.zp) = 1;
    ub(var.zw) = 1;
    ub(var.epsp) = inf;
    ub(var.epsw) = 0;
    ub(var.delp) = inf;
    ub(var.delw) = 0;
    ub(var.x) = inf;
    lb(var.u, 1) = -inf;
    lb(var.vp) = 0;
    lb(var.w) = 0;
    lb(var.zp) = 0;
    lb(var.zw) = 0;
    lb(var.epsp) = 0;
    lb(var.epsw) = 0;
    lb(var.delp) = 0;
    lb(var.delw) = 0;
    lb(var.x) = 0;
    lb(var.x(end)) = phi;

% Cplex - obj vector
    obj(var.u, 1) = 0;
    obj(var.vp) = 0;
    obj(var.w) = 0;
    obj(var.zp) = 1;
    obj(var.zw) = 0;
    obj(var.epsp) = 0;
    obj(var.epsw) = 0;
    obj(var.delp) = 0;
    obj(var.delw) = 0;
    obj(var.x) = 0;

% Cplex - ctype vector
    ctype(var.u) = 'C';
    ctype(var.vp) = 'C';
    ctype(var.w) = 'C';
    ctype(var.zp) = 'B';
    ctype(var.zw) = 'B';
    ctype(var.epsp) = 'C';
    ctype(var.epsw) = 'C';
    ctype(var.delp) = 'C';
    ctype(var.delw) = 'C';
    ctype(var.x) = 'C';

% Cplex - sense of the optimization
    sense = 'minimize';

% Cplex - Introduce all data in a Cplex structure
    cplex = Cplex('MCS');
    Model = struct();
    [Model.A, Model.rhs, Model.lhs, Model.ub, Model.lb, Model.obj, Model.ctype, Model.sense] = deal(A, rhs, lhs, ub, lb, obj, ctype, sense);
    cplex.Model = Model;

% Cplex Indicators
    % z = 1  -->  v >= alpha
    for ivar = 1:length(var_group.z)
        a = zeros(var.x(end), 1);
        a(var_group.v(ivar)) = 1;
        cplex.addIndicators(var_group.z(ivar), 0, a, 'G', alpha);
    end

% Cplex Indicators
    % z = 0  -->  v <= 0
    for ivar = 1:length(var_group.z)
        a = zeros(var.x(end), 1);
        a(var_group.v(ivar)) = 1;
        cplex.addIndicators(var_group.z(ivar), 1, a, 'L', 0);
    end

% Cplex Indicators
    % z = 1  -->  epsilon <= 0
    for ivar = 1:length(var_group.z)
        a = zeros(var.x(end), 1);
        a(var_group.eps(ivar)) = 1;
        cplex.addIndicators(var_group.z(ivar), 0, a, 'L', 0);
    end

% Cplex Indicators
    % z = 0  -->  epsilon <= M
    for ivar = 1:length(var_group.z)
        a = zeros(var.x(end), 1);
        a(var_group.eps(ivar)) = 1;
        cplex.addIndicators(var_group.z(ivar), 1, a, 'L', M);
    end

% Cplex Parameters
    sP = struct();
    [sP.mip.tolerances.integrality, sP.mip.strategy.heuristicfreq, sP.mip.strategy.rinsheur] = deal(integrality_tolerance, 1000, 50);
    [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers);
    [sP.preprocessing.aggregator, sP.preprocessing.boundstrength, ...
        sP.preprocessing.coeffreduce, sP.preprocessing.dependency, ...
        sP.preprocessing.dual, sP.preprocessing.fill,...
        sP.preprocessing.linear, sP.preprocessing.numpass, ...
        sP.preprocessing.presolve, sP.preprocessing.reduce,..., ...
        sP.preprocessing.relax, sP.preprocessing.symmetry] = deal(50, 1, 2, 1, 1, 50, 1, 50, 1, 3, 1, 1);
    cplex = setCplexParam(cplex, sP);
    if printLevel == 0
        cplex.DisplayFunc = [];
    end

% Calculation of MCSs
    mcs_time{1, 1} = '------ TIMING ------';
    mcs_time{1, 2} = '--- MCSs ---';
    i = 0;
    k = 0;
    n_time = size(mcs_time, 1);
    mcs_time{n_time+1, 1} = 'Preparation';
    mcs_time{n_time+1, 2} = toc(time_aa);
    mcs = [];
    largest_mcs = 0;
    while largest_mcs <= max_len_mcs && k < n_mcs && cplex.Model.rhs(cons.forceLength) <= max_len_mcs
        ini_mcs_time = toc(time_aa);
        cplex.Param.mip.limits.populate.Cur = 40;
        cplex.Param.mip.pool.relgap.Cur = 0.1;
        cplex.populate();
        n_pool = size(cplex.Solution.pool.solution, 1);
        if n_pool ~= 0
            solution = cplex.Solution.pool.solution;
            for j = 1:n_pool
                k = k+1;
                mcs{k, 1} = K_ind((solution(j).x(var.zp))>0.9);
                n_cons = n_cons+1;
                sol = solution(j).x(var.zp)>0.9;
                cplex.Model.A(n_cons, var.zp) = sparse(double(sol));
                cplex.Model.rhs(n_cons) = sum(sol)-1;
                cplex.Model.lhs(n_cons) = 0;
            end
            i = i+1;
            mcsi_time = toc(time_aa)-ini_mcs_time;
            n_time = size(mcs_time, 1);
            mcs_time{n_time+1, 1} = ['POPULATE_ORDER_' num2str(cplex.Model.rhs(cons.forceLength))];
            mcs_time{n_time+1, 2} = mcsi_time;
        else
            mcsi_time = toc(time_aa)-ini_mcs_time;
            n_time = size(mcs_time, 1);
            mcs_time{n_time+1, 1} = ['POPULATE_ORDER_' num2str(cplex.Model.rhs(cons.forceLength)) 'NF'];
            mcs_time{n_time+1, 2} = mcsi_time;
            if forceLength == 1
                cplex.Model.rhs(cons.forceLength) = cplex.Model.rhs(cons.forceLength)+1;
                cplex.Model.lhs(cons.forceLength) = cplex.Model.lhs(cons.forceLength)+1;
            else
                n_time = size(mcs_time, 1);
                mcs_time{n_time+1, 1} = 'TOTAL MCSs';
                mcs_time{n_time+1, 2} = toc(time_aa);
                return;
            end
        end
        try save('tmp.mat', 'mcs', 'mcs_time'); end
        try largest_mcs = max(cellfun(@length, mcs)); end
    end
end
n_time = size(mcs_time, 1);
mcs_time{n_time+1, 1} = 'TOTAL MCSs';
mcs_time{n_time+1, 2} = toc(time_aa);
end

author: opencobra

repo: cobratoolbox


text:
function str = vector2str(num)
% Converts a vector to string with colon operators
str = num2str(num(1));
incrementStart = true;
for a = 2:length(num)
    if num(a) - num(a-1) == 1
        incrementStart = handleConsecutive(incrementStart);
    else
        incrementStart = handleNonConsecutive(num(a-1), num(a), ...
            incrementStart);
    end
end
handleLastIndex(num(length(num)), incrementStart);

    function incrementStart = handleNonConsecutive(previous, current, ...
            incrementStart)
        % Handles a number that is not a consecutive number
        if ~incrementStart
            str = [str num2str(previous)];
        end
        str = [str ' ' num2str(current)];
        incrementStart = true;
    end % handleNonConsecutive

    function handleLastIndex(last, incrementStart)
        % Handles the last number in the vecor
        if ~incrementStart
            str = [str num2str(last)];
        end
    end % handleLastIndex

    function incrementStart = handleConsecutive(incrementStart)
        % Handles a number that is a consecutive number
        if incrementStart
            str = [str ':'];
            incrementStart = false;
        end
    end % handleConsecutive

end % vector2str

author: VisLab

repo: EEG-Clean-Tools


text:
function [r,s,t,w] = threed_gauss(rule)
%-----------------------------------------------------------------------
%  threed_gauss.m - calculate Gauss integration points for tetrahedral
%                   elements
%
%  Copyright (c) 2002, Jeff Borggaard, Virginia Tech
%  Version: 1.0
%
%  Usage:    [r,s,t,w] = threed_gauss(rule)
%
%  Variables:     rule
%                        Number of Gauss points:
%                 r
%                        xi coordinate of Gauss points
%                 s
%                        eta coordinate of Gauss points
%                 t
%                        zeta coordinate of Gauss points
%                 w
%                        Gauss weights corresponding to (r,s,t)
%-----------------------------------------------------------------------

  if (rule == 1)
    % The following points correspond to a 1 point rule
    % Polynomials of degree 1 are integrated exactly
    r(1) = 0.25000000;  s(1) = 0.25000000;  t(1) = 0.25000000;

    w(1) = 0.16666666;  % (the volume of a regular tet is 1/6)

  elseif (rule == 4)
    % The following points correspond to a 4 point rule
    % Polynomials of degree 2 are integrated exactly
    r    = zeros(4,1);  s    = zeros(4,1);  t    = zeros(4,1);
    r(1) = 0.13819660;  s(1) = 0.13819660;  t(1) = 0.13819660;
    r(2) = 0.58541020;  s(2) = 0.13819660;  t(2) = 0.13819660;
    r(3) = 0.13819660;  s(3) = 0.58541020;  t(3) = 0.13819660;
    r(4) = 0.13819660;  s(4) = 0.13819660;  t(4) = 0.58541020;

    w    = 0.041666666667*ones(4,1);

  elseif (rule == 5)
    % The following points correspond to a 5 point rule
    % Polynomials of degree 3 are integrated exactly
    r    = zeros(5,1);  s    = zeros(5,1);  t    = zeros(5,1);
    r(1) = 0.25000000000;  s(1) = 0.25000000000;  t(1) = 0.25000000000;
    r(2) = 0.16666666666;  s(2) = 0.16666666666;  t(2) = 0.16666666666;
    r(3) = 0.50000000000;  s(3) = 0.16666666666;  t(3) = 0.16666666666;
    r(4) = 0.16666666666;  s(4) = 0.50000000000;  t(4) = 0.16666666666;
    r(5) = 0.16666666666;  s(5) = 0.16666666666;  t(5) = 0.50000000000;

    w    = zeros(5,1);
    w(1) =-0.13333333333;
    w(2) = 0.07500000000;
    w(3) = 0.07500000000;
    w(4) = 0.07500000000;
    w(5) = 0.07500000000;

  else
    error('quadrature rules other than 1, 4 or 5 are not supported\n');
    keyboard
  end

author: johannesgerer

repo: jburkardt-m


text:
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

function [gt_data, gt_freqs, gt_levels] = prescription_minimalistic(freqs, thresholds_left, thresholds_right, offset, marginfactor, rolloff, center, focus)
  gt_freqs = [177 297 500 841 1414 2378 4000 6727 11314];
  gt_levels = -10:1:110;
  reference_freqs  = [125 250 500 1000 2000 4000 8000 16000];
  reference_levels = abs(polyval([focus 0],log2(reference_freqs./1000)))+offset;
  thresholds_left_ext = [0 0 thresholds_left 0 0];
  thresholds_right_ext = [0 0 thresholds_right 0 0];
  freqs_ext = [0 50 freqs 16000 48000];
  gt_data_left = zeros(length(gt_levels),length(gt_freqs));
  gt_data_right = zeros(length(gt_levels),length(gt_freqs));
  maxgain = 40;
  maxlevel = 110;
  for i=1:length(gt_freqs)
    reference_level = interp1(reference_freqs,reference_levels,gt_freqs(i),'extrap');
    threshold_level_left = interp1(freqs_ext,thresholds_left_ext,gt_freqs(i),'extrap');
    threshold_level_right = interp1(freqs_ext,thresholds_right_ext,gt_freqs(i),'extrap');
    low_level_gain_left = max(0,threshold_level_left-reference_level);
    low_level_gain_right = max(0,threshold_level_right-reference_level);
    margin_left = marginfactor.*(center-(reference_level+low_level_gain_left));
    margin_right = marginfactor.*(center-(reference_level+low_level_gain_right));
    gt_data_left(:,i) = interp1([gt_levels(1);reference_level+margin_left;reference_level+margin_left+low_level_gain_left.*rolloff;gt_levels(end)],[low_level_gain_left;low_level_gain_left;0;0],gt_levels);
    gt_data_right(:,i) = interp1([gt_levels(1);reference_level+margin_right;reference_level+margin_right+low_level_gain_right.*rolloff;gt_levels(end)],[low_level_gain_right;low_level_gain_right;0;0],gt_levels);
  end
  gt_data = [gt_data_left.';gt_data_right.'];
  gt_data = min(maxgain,gt_data);
  gt_data = gt_data + min(0,maxlevel - (gt_data+gt_levels));
end

author: m-r-s

repo: hearingaid-prototype


text:
classdef nme_shunt_acp_nln < mp.nme_shunt_acp & mp.nme_wrapper_ac_nln

%   MATPOWER
%   Copyright (c) 2019-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        nme_class = @mp.nme_shunt_acp;
    end

    methods
        function obj = nme_shunt_acp_nln()
            obj@mp.nme_shunt_acp();
            obj.nme_wrapper_ac_nln_init();
        end

        function build_params(obj, nm, dm)
            build_params@mp.nme_shunt_acp(obj, nm, dm);
            obj.build_nln_params(nm, dm);
        end

        function nk = count(obj, dm)
            obj.nme.count(dm);
            nk = count@mp.nme_shunt_acp(obj, dm);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Simple Low Rank Tensor Completion (SiLRTC)
% Time: 03/11/2012
% Reference: "Tensor Completion for Estimating Missing Values
% in Visual Data", PAMI, 2012.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X, errList] = SiLRTC(T, Omega, alpha, gamma, maxIter, epsilon, X)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% min(X, M1, M2, M3,... Mn): (\gamma1||X_(1)-M1||^2 + \gamma2||X_(2)-T2||^2 + \gamma3||X_(3)-T3||^2 + ...)/2 +
%               \alpha1||M1||_* + \alpha2||M2||_* + \alpha3||M3||_* + ....
%         s.t.  X_\Omega = T_\Omega
%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 7
    X = T;
    X(logical(1-Omega)) = mean(T(Omega));
end

errList = zeros(maxIter, 1);
normT = norm(T(:));
%L = errList;
dim = size(T);
M = cell(ndims(T), 1);
gammasum = sum(gamma);
tau = alpha./ gamma;

%normT = norm(T(:));
for k = 1:maxIter
    if mod(k, 20) == 0
        fprintf('SiLRTC: iterations = %d   difference=%f\n', k, errList(k-1));
    end
    Xsum = 0;
    for i = 1:ndims(T)
        M{i} = Fold(Pro2TraceNorm(Unfold(X, dim, i), tau(i)), dim, i);
        Xsum = Xsum + gamma(i) * M{i};
    end
    Xlast = X;
    X = Xsum / gammasum;
    X(Omega) = T(Omega);
    errList(k) = norm(X(:)-Xlast(:)) / normT;
    if (errList(k) < epsilon)
        errList = errList(1:k);
        break;
    end
    %L(k) = norm(X(:)-T(:)) / normT;
end
fprintf('SiLRTC ends: total iterations = %d   difference=%f\n\n', k, errList(k));

author: andrewssobral


text:
report_this_filefun(mfilename('fullpath'));

bz2 = []; z = [];
p1 = []; p1b=[];p99b=[];p50b=[];pmib=[];pmab=[];
p99 = []; zr =[];
p50 = [];pma=[]; pmi=[];
niv = 1:0.2:4;

tdiff = round((teb - t0b)*365/par1);

def = {'100','300'};
tit ='Random Zmax calculation';
prompt={'Number of events in each window?', 'Number of random samples drawn? (Grid-size)'};

ni2 = inputdlg(prompt,tit,1,def);
l = ni2{1};
nil = str2double(l);
l = ni2{2};
no = str2double(l);

%ni = str2double(prmptdlg('Number of events in each window?','100'));
%no = str2double(prmptdlg('Number of random samples drawn ?','30'));
iwl = iwl*365/par1;
zr = (-15:0.1:15)*0;
wai = waitbar(0,' Please Wait ...  ');
set(wai,'NumberTitle','off','Name','Makegrid  -Percent done');;
n0 = no;
na = no;

for iwl = 1.0:0.2:4
    iwl
    iwl = iwl*365/par1;
    zr = 0;
    na = n0; no = n0;
    while na+1000 > 1000;
        if na > 1000; na = 1000; end
        l = ceil(rand([ni na])*a.Count);
        [cumu, xt] = hist(reshape(a(l,3),ni,na),(t0b:par1/365:teb));
        for j = 2:tdiff-iwl
            cu = [cumu(1:j-1,:) ; cumu(j+iwl+1:length(cumu(:,1)),:)];
            mean1 = mean(cu);
            mean2 = mean(cumu(j:j+iwl,:));
            var1 = (std(cu)).^2;;
            var2 = (std(cumu(j:j+iwl,:)).^2);
            as = (mean1 - mean2)./(sqrt(var1/(length(cumu(:,1))-iwl)+var2/iwl));
            p1 = [p1 prctile2(as,1)];
            p99 = [p99 prctile2(as,99)];
            p50 = [p50 prctile2(as,50)];
            si =  [ si std(as)];
            pma = [pma max(as)];
            pmi = [pmi min(as)];
            [tmp, tmp2] = hist(as,-15:0.1:15);
            zr = [zr + tmp];
        end     % for j
        p1 = [p1 prctile2(as,1)];
        p99 = [p99 prctile2(as,99)];
        p50 = [p50 prctile2(as,50)];
        pma = [pma max(as)];
        pmi = [pmi min(as)];
        no = no - 1000;
        na = no;
    end % while na

    p1b = [p1b mean(p1) ];
    p99b = [p99b mean(p99) ];
    p50b = [p50b mean(p50) ];
    pmab = [pmab max(pma) ];
    pmib = [pmib min(pmi) ];
    z = [z , zr'];
    p1 = []; p99 = []; p50 = [];pma=[]; pmi=[];
    waitbar((iwl*par1/365)/max(niv));

end  % for iwl

close(wai)
figure
pl =plot(niv,p1b,'b')
hold on
set(pl,'LineWidth',2.0)
pl =plot(niv,p99b,'b')
set(pl,'LineWidth',2.0)
pl =plot(niv,p50b,'r')
pl =plot(niv,pmab,'g--')
set(pl,'LineWidth',2.0)
pl =plot(niv,pmib,'g--')
set(pl,'LineWidth',2.0)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
grid
xlabel('Windowlength in [years]')
ylabel('Range of z')
title(['ni  =  ' num2str(ni) ' events, ' num2str(n0) ' random samples'])

matdraw

figure
pcolor(z)
shading flat
colormap(jet)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
xlabel('Windowlength in [years]')
ylabel('Range of z')
title(['ni  =  ' num2str(ni) 'events, ' num2str(n0) ' random samples'])

matdraw

author: CelsoReyes


text:
function [J] = Nm2J(Nm)
% Convert energy or work from newton-meters to joules.
% Chad A. Greene 2012
J = Nm;
% You say tomato, I say tomatermort.

author: Sable

repo: mcbench-benchmarks


text:
function [Int, p, mask, w3] = estPolIntGrad(inp, N, logFlag);
% ESTPOLINTGRAD - Estimates and corrects the intensity gradient, using a polynomial model
%
% [Int pol] = estPolIntGrad(inp, N, logFlag);
%
% Inputs:
%  inp - input inplanes affected by the intensity gradient
%  N - order of the polynomial [Nx Ny Nz]
%  logFlag - operates in the logarithm of the intensity
%

if ~exist('N')
   N = [3,3,3];
end
if ~exist('logFlag')
   logFlag = 1;
end

[Ny, Nx, Nz] = size(inp);
inp(find(inp==0)) = NaN;

% chosing appropriate thresholds
II = find(~isnan(inp));
% fitting a GMM
inpMu = mean(inp(II));
inpStd = std(inp(II));
outlim = [min(inp(II)) max(inp(II))];
initCov(1,1,1) = inpStd^2/4;
size(initCov)
[Mu Cov W] = fitGMM(inp(II)', 1, inpMu+2.5*inpStd, initCov, [0.9 0.1], outlim);
%initmu = [inpMu-2*inpStd inpMu+2.5*inpStd];
%initCov(:,:,1:2) = 2*inpStd^2;
%initW = [0.45 0.45 0.1];
%outlim = [min(inp(II)) max(inp(II))];
%[Mu Cov W] = fitGMM(inp(II)', 2, initmu, initCov, initW, outlim);

[hh xx] = hist(inp(II),256);
hh = hh/sum(hh)/(xx(2)-xx(1));
p1 = W(1)*pgauss1d(xx,Mu(1),Cov(:,:,1));
%p2 = W(2)*pgauss1d(xx,Mu(2),Cov(:,:,2));
p3 = W(2)/(max(inp(II))-min(inp(II))) * ones(size(p1));
figure(20); clf
plot(xx,hh,xx,p1+p3)

% chosing as thresholds mu2+-2std2
LoT = Mu(1) - 3*sqrt(Cov(:,:,1))
UpT = Mu(1) + 3*sqrt(Cov(:,:,1))
inp(find( (inp<LoT) | (inp>UpT) )) = NaN;

figure(20)
hold on
stem([LoT UpT], [max(hh),max(hh)]);
hold off

% Change polynomial order to be 2 less than number of non-NaN slices
NzEff = 0;
for k=1:Nz
    tmp = inp(:,:,k);
    if ~all(isnan(tmp(:)))
        NzEff = NzEff + 1;
    end
end
N(3) = min(max(1,NzEff-2), N(3));

% taking logarithm
if logFlag
   inp = log(inp);
end

% selecting valid set of indices
II = find(~isnan(inp));

% building A matrix for the appropriate polynomial order
[x y z] = meshgrid(1:Nx,1:Ny,1:Nz);
A = [];
for m=0:N(1)
  for n=0:max(0,(N(2)-m))
    for o=0:max(0,(N(3)-m-n))
         A = [A x(II).^m .* y(II).^n .* z(II).^o];
      end
   end
end
clear x y z

% robust estimation of the polynomial coefficients
% NOTE - The values of the last two parameters in this function can affect the
% performance of the solution. In general, the recommended value for CB = 4.685
% is very conservative with outliers (almost nothing is considered outlier).
% For this reason I choose a value of 2.5. I also set a bit lower the value of SC
% (1.2 instead of 1.4), so that more outliers are rejected.
[p w]= robustMest(A, inp(II), 3, 1.2);
w3 = zeros(size(inp));
w3(II) = w;

% computation of the intensity correction
il = zeros(Ny,Nx,Nz); k=1;
[x y z] = meshgrid(1:Nx, 1:Ny, 1:Nz);
for m=0:N(1)
  for n=0:max(0,(N(2)-m))
    for o=0:max(0,(N(3)-m-n))
         il = il + p(k)*x.^m.*y.^n.*z.^o;
         k = k+1;
      end
   end
end

if logFlag
   Int = exp(il);
else
   Int = il;
end

mask = zeros(size(inp));
mask(II) = 1;

return

clear
load testdata
[Int, p, mask, w3] = estPolIntGrad(inp, [3 3 1], 1);
inpc = inp./Int;
sinpc = reshape(inpc, [size(inp,1) size(inp,2)*size(inp,3)]);
figure(1)
subplot(4,1,1)
vis(inp)
subplot(4,1,2)
imshow(sinpc, [0.5 1.5], 'notruesize')
subplot(4,1,3)
vis(Int.*mask)
subplot(4,1,4)
vis(w3)
figure(2)
II = find((inpc~=0)&(~isnan(inpc)));
hist(inpc(II),256)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function element_node = grid_element ( code, element_order, nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_ELEMENT returns the grid associated with any available element.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CODE(*), identifies the element desired.
%    Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3',
%    'T6' and 'T10'.
%
%    Input, integer ELEMENT_ORDER, the order of the element.
%
%    Input, integer NELEMX, NELEMY, the number of quadrilaterals along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY for quadrilaterals, or 2 * NELEMX * NELEMY for
%    triangles.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the nodes
%    that form each element.
%
  if ( s_eqi ( code, 'Q4' ) )
    element_node = grid_q4_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q8' ) )
    element_node = grid_q8_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q9' ) )
    element_node = grid_q9_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q12' ) )
    element_node = grid_q12_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q16' ) )
    element_node = grid_q16_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'QL' ) )
    element_node = grid_ql_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T3' ) )
    element_node = grid_t3_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T4' ) )
    element_node = grid_t4_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T6' ) )
    element_node = grid_t6_element ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T10' ) )
    element_node = grid_t10_element ( nelemx, nelemy );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'GRID_ELEMENT - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of CODE = %s\n', code );
    error ( 'GRID_ELEMENT - Fatal error!' );

  end

  return
end


text:
classdef PTKReportingDefault < CoreReportingDefault
    % PTKReportingDefault. Legacy support class. Replaced by CoreReportingDefault.
    %
    %
    %     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.
    %

    methods
        function obj = PTKReportingDefault(varargin)
            obj = obj@CoreReportingDefault(varargin{:});
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function y = circle_characteristic ( m, n, x )

%*****************************************************************************80
%
%% CIRCLE_CHARACTERISTIC evaluates the characteristic function of a circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points to check.
%
%    Input, real X(M,N), the coordinates of points to be checked.
%
%    Output, real Y(N,1), is 1 if the point is inside, 0 otherwise.
%
  y = sqrt ( sum ( x.^2, 1 ) );

  y = y';

  y = ( y < 1.0 );

  return
end


text:
%computes the spectral spread from the magnitude spectrum
%> called by ::ComputeFeature
%>
%> @param X: spectrogram (dimension FFTLength X Observations)
%> @param f_s: sample rate of audio data
%>
%> @retval vss spectral spread (in Hz)
% ======================================================================
function [vss] = FeatureSpectralSpread (X, f_s)

    % get spectral centroid as index
    vsc = FeatureSpectralCentroid(X, f_s) * 2 / f_s * (size(X, 1)-1);

    % allocate memory
    vss = zeros(size(vsc));

    % compute spread
    for n = 1:size(X,2)
        vss(n) = (((0:size(X, 1)-1)-vsc(n)).^2 * X(:, n)) ./ sum(X(:, n));
    end
    vss = sqrt(vss);

    % convert from index to Hz
    vss = vss / (size(X, 1)-1) * f_s / 2;

    % avoid NaN for silence frames
    vss (sum(X, 1) == 0) = 0;
end

author: alexanderlerch

repo: ACA-Code


text:
function chanlist = ie_readchans(chanfile)
%IE_READCHANS Read and parse station/channel list from file
% CHANLIST = IE_READCHANS(CHANFILE) returns a string cell array of stations
% and their associated channels, read from comma-delimited input filename CHANFILE.
% Each line of CHANFILE should contain up to four fields: a station name in
% the first column, then either one or three channel names in the
% subsequent columns.  Lines beginning with the percent sign % are ignored.
%
% IE_READCHANS performs some basic error checking on the input file, but
% does not compare against any database or other listing of station data.
% This is left up to the user.
%
% EXAMPLE
%
% Suppose CHANFILE points to the following file:
%
% % Stations and channels of interest in Uturuncu data set
% UTCA,SHZ,SHN,SHE
% UTCM,SHZ,,
% %UTKH,BHZ,BHN,BHE
% %UTLA2,SHZ,SHN,SHE
% UTSW,SHZ,SHN,SHE
%
% Then IE_READCHANS will return:
%
% ans =
%     'UTCA'    'SHZ'    'SHN'    'SHE'
%     'UTCM'    'SHZ'       ''       ''
%     'UTSW'    'SHZ'    'SHN'    'SHE'
%
% Author: Christopher Bruton, Geophysical Institute, University of Alaska Fairbanks
% $Date$
% $Revision$

	[fp, message] = fopen(chanfile,'r'); % Open the file

	if fp < 0 % Could not open the file
		error(['Could not open file ' chanfile ': ' message]);
	end

	% The following seems to be the best way to read a CSV file of text.
	% csvread() can only handle numeric data.
	filecontents = textscan(fp,'%s %s %s %s','Delimiter',',','ReturnOnError',false,'CommentStyle','%');

	fclose(fp); % Close file handle

	% Concatenate into one cell array
	chanlist = cat(2,filecontents{:});

	% Run some mechanical checks

	for n = 1:size(chanlist,1)
		% Empty station?
		if isempty(chanlist{n,1})
			error(['Empty station on line ' num2str(n) ' of input file ' chanfile]);
		end

		% Number of empty channels
		n_empty = 0;
		for c = 2:4
			if isempty(chanlist{n,c})
				n_empty = n_empty + 1;
			end
		end

		% Error if we have wrong number of empty channel names
		if (n_empty ~= 2) && (n_empty ~= 0)
			error(['Incorrect number of channels specified on line ' num2str(n) ' of input file ' chanfile '.  Please specify either 1 or 3 channel names.']);
		end

		% All channels must have unique names
		if length(unique(chanlist(n,2:4))) ~= (3-(n_empty>0))
			error(['Nonunique channel names specified on line ' num2str(n) ' of input file ' chanfile]);
		end
	end

	% All stations must have unique names
	if length(unique(chanlist(:,1))) ~= length(chanlist(:,1))
		error(['Nonunique station names specified in input file ' chanfile]);
	end

end

author: geoscience-community-codes

repo: GISMO


text:
function [n, npw, ts, wave] = plx_waves_v(filename, channel, unit)
% plx_waves_v(filename, channel, unit): read waveform data from a .plx or .pl2 file
%
% [n, npw, ts, wave] = plx_waves_v(filename, channel, unit)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   channel - 1-based channel number or channel name
%   unit  - unit number (0- unsorted, 1-4 units a-d)
%
% OUTPUT:
%   n - number of waveforms
%   npw - number of points in each waveform
%   ts - array of timestamps (in seconds)
%   wave - array of waveforms [npw, n] converted to mV

n = 0;
npw = 32;
ts = -1;
wave = -1;

if nargin ~= 3
    error 'expected 3 input arguments';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    w = PL2Waves( filename, channel, unit );
    if numel(w.Ts) > 0
        n = numel(w.Ts);
        ts = w.Ts;
        npw = w.NumPointsWave;
        wave = w.Waves;
    else
		% we need to return number of points per wave the way it is done by mexPlex
		if channel > -1
			npw = 0;
			pl2 = PL2GetFileIndex(filename);
			% find the first channel with data and get its number of points per wave
			for i=1:numel(pl2.SpikeChannels)
				if pl2.SpikeChannels{i}.SamplesPerSpike > 0
					npw = pl2.SpikeChannels{i}.SamplesPerSpike;
					break;
				end
			end
		end
    end
    return;
end

channelNumber = plx_resolve_channel(filename, channel);
if channelNumber == -1
    %fprintf('\n plx_waves_v: no header for the specified spike channel.\n');
    return
end

[n, npw, ts, wave] = mexPlex(19, filename, channelNumber, unit);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

Z = X * U(:, 1:K);

% =============================================================

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function input = handle_edit_input(text)

% HANDLE_EDIT_INPUT deals with user-entered input in the GUI. This is used to select
% channels and/or trials in FT_REJECTVISUAL and to select channels in FT_DATABROWSER
%
% The input text can consist of a string such as
%   1 2 3 4
%   1:4
%   [1 2 3 4]
%   [1:4]
% This is converted in a list of numbers.
%
% The input text can also consist of a single non-numeric string or a string that
% represents a cell-array of strings such as
%   all
%   {'MEG', '-MR*'}

% Copyright (C) 2022, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isempty(text)
  input = [];
elseif all(ismember(text, '0123456789 \t'))
  input = regexp(text, '\s+', 'split');
  input = cellfun(@str2num, input);
elseif ismember(text(1), '0123456789')
  input = eval(text); % it could be a single number or a string like 1:10
elseif startsWith(text, '[') && endsWith(text, ']')
  input = eval(text);
elseif startsWith(text, '{') && endsWith(text, '}')
  input = eval(text); % this allows for {'MEG', '-MR*'}
else
  input = text;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function run_writedataset(varargin)
% Output a raw online stream into an EEGLAB dataset.
% run__writedataset(SourceStream,FileName,UpdateFrequency,StartDelay)
%
% This function does not do any processing, but just saves a stream to a file (possibly in parallel
% to some other operation processing it).
%
% In:
%   SourceStreamName : Optional name of the stream data structure in the MATLAB base workspace to
%                      take as the data source (previously created with onl_newstream).
%                      (default: 'laststream')
%
%   FileName : File name to write to (default: 'lastdata.set')
%
%   UpdateFrequency : The are at which new chunks of data are appended to the file, in Hz (default: 1)
%
%   StartDelay : Start-up delay before real-time operation begins; grace period until file is being
%                written to, in s. (default: 3)
%
% Examples:
%   % write an input stream (named 'mystream') to a file named 'recording.set' (EEGLAB dataset)
%   run_writedataset('mystream','recording.set')
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-19

declare_properties('name','File');

% define arguments
arg_define(varargin, ...
    arg({'in_stream','SourceStreamNames','SourceStream'}, 'laststream',[],'Input Matlab stream name(s). Optional names of stream data structures in the MATLAB base workspace to consider as possible data sources (previously created with onl_newstream); if a stream contains all channels that are needed by the predictor, or alternatively has the right number and type of channels it will be considered as a potential source stream unless ambiguous.'), ...
    arg({'out_filename','FileName'},'lastdata.set',[],'The file name to write to.'), ...
    arg({'update_freq','UpdateFrequency'},1,[0 Inf],'Update frequency. This is the rate at which new chunks of data are appended to the file.'), ...
    arg({'start_delay','StartDelay'}, 3, [0 Inf],'Start-up delay. Delay before real-time operation begins; grace period until file is written.'));

if ~any(out_filename=='/' | out_filename=='\') %#ok<NODEF>
    out_filename = ['bcilab:/userdata/' out_filename]; end
out_filename = env_translatepath(out_filename);

% open the stream and write the initial set file header...
stream = evalin('base',in_stream);
% create missing fields
stream.data = randn(stream.nbchan,1024);
stream.pnts = size(stream.data,2);
stream.xmax = stream.xmin + (stream.pnts-1)/stream.srate;
stream.timestamp_at_beginning = toc(uint64(0));

% prepare .set file for saving
[fp,fn,fe] = fileparts(out_filename); %#ok<ASGLU>
EEG = rmfield(stream,{'buffer','smax','buffer_len','timestamps','timestamps_len','tmax','streamid','timestamp_at_beginning'});
[EEG.data,EEG.datfile] = deal([fn,fe]);
io_save(out_filename,'-mat','-makedirs','-attributes','''+w'',''a''','EEG');

% create the .fdt file...
fid = fopen(fullfile(eeg.filepath, eeg.datfile),'wb','ieee-le');
if fid == -1
    error('Cannot write output file, check permission and space.'); end;

% create a temporary text file to store stream event markers
event_filename = env_translatepath(['bcilab:/userdata/', fn, '_events.txt']);
marker_fid = fopen(event_filename, 'w+');
if marker_fid == -1
    error('Cannot write marker output file, check permission and space.');
end
% write a basic header to the event marker file
fprintf(marker_fid, '%s\t%s\n', 'type', 'latency');

% create timer (which periodically writes to the stream)
t = timer('ExecutionMode','fixedRate', 'Name',[in_stream '_write_timer'], 'Period',1/update_freq, ...
    'StartDelay',start_delay, 'TimerFcn',@(obj,varargin) append_data(in_stream,fid,marker_fid,stream.streamid,obj,eeg));

% start timer
start(t);

% timer callback: visualization
function append_data(stream,fid,marker_fid,streamid,timerhandle,eeg)
try
    % check if the stream and the predictor are still there
    s = evalin('base',stream);
    if s.streamid ~= streamid
        error('Stream changed.'); end

    % get an updated chunk of data
    samples_to_get = min(s.buffer_len, s.smax-ftell(fid)/(4*s.nbchan));
    chunk = s.buffer(:, 1+mod(s.smax-samples_to_get:s.smax-1,s.buffer_len));

    % and write it into the file
    fwrite(fid,chunk,'float');

    % check if this data chunk contains any event markers
    marker_chunk = s.marker_pos(:, 1+mod(s.smax-samples_to_get:s.smax-1,s.buffer_len));

    % if so, write them into the marker file
    if nnz(marker_chunk) > 0
        % find the sample(s) in this chunk with events
        [marker_pos_in_sample, marker_pos_in_chunk] = find(marker_chunk);

        % EEG samples which correspond to positions in the chunk
        chunk_samples = s.smax - (samples_to_get - 1) : s.smax;

        for m = 1:length(marker_pos_in_chunk)
            marker_idx = marker_chunk(marker_pos_in_sample(m), marker_pos_in_chunk(m));
            marker_type = s.marker_buffer(marker_idx).type;
            marker_sample_whole = chunk_samples(marker_pos_in_chunk(m));
            marker_sample_fractional = s.marker_buffer(marker_idx).latency;
            marker_sample = marker_sample_whole + marker_sample_fractional;
            fprintf(marker_fid, '%s\t%i\n', marker_type, marker_sample);
        end
    end

catch e
     if ~strcmp(e.identifier,'MATLAB:UndefinedFunction')
        env_handleerror(e); end
    finalize_dataset(fid,marker_fid,eeg);
    % interrupted: make sure that the file gets closed
    stop(timerhandle);
    delete(timerhandle);
end

function finalize_dataset(fid,marker_fid,EEG)
samples = ftell(fid)/(4*EEG.nbchan);
fclose(fid);
EEG.pnts = samples;
EEG.data = EEG.datfile;
EEG.xmax = EEG.xmin + (EEG.pnts-1)/EEG.srate;
EEG.timestamp_at_end = toc(uint64(0));

% load the event marker file from disk and include the markers in the EEG .set
marker_file = [EEG.filepath, filesep, regexprep(EEG.filename, '.set', '_events.txt')];
frewind(marker_fid);
marker_data = textscan(marker_fid, '%s%d', 'Delimiter', '\t', 'HeaderLines', 1);
if ~isempty(marker_data{1})
   EEG.event = struct('type', marker_data{1}', 'latency', num2cell(marker_data{2}'), ...
        'urevent', num2cell((1:length(marker_data{1}))));
end
fclose(marker_fid);
delete(marker_file);

save(fullfile(EEG.filepath, EEG.filename), '-v6', '-mat', 'EEG');

repo: REST


text:
classdef (TestTags = {'SPGR', 'Demo', 'Integration'}) SimDemo_SPGR_Test < matlab.unittest.TestCase

    properties
       qmrlabPath = cell2mat(regexp(cd, '.*qMRLab/', 'match'));
    end
    methods (TestClassSetup)
        function addqMRLabToPath(testCase)
            addpath(genpath(testCase.qmrlabPath));
        end
    end

    methods (TestClassTeardown)
        function removeqMRILabFromPath(testCase)
            clear all, close all
        end
    end

    methods (Test)
        function testFittedParamsNearInputValues(testCase)
            run([testCase.qmrlabPath, '/Models_Functions/SPGRfun/SimDemo_SPGR.m'])

            inputParams  = Sim.Param;
            outputParams = SimCurveResults;

            inputArr  = [inputParams.F  inputParams.kf  inputParams.R1f  inputParams.T2f  inputParams.T2r];
            outputArr = [outputParams.F outputParams.kf outputParams.R1f outputParams.T2f outputParams.T2r];

            %                                                 , # percent
            %                                                 . [F  kf R1f T2f T2r]
            testCase.verifyLessThan(pDiff(inputArr, outputArr), [30 30 30  30  30]);
        end
    end

end

function value = pDiff(inputVal, outputVal)
    value = abs((outputVal-inputVal)./inputVal).*100;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [unfoldMesh, nFaces] = mfmBuildSubMesh(mesh, params, perimeterEdges, insideNodes, ...
    orderedUniquePerimeterPoints)
%
% [unfoldMesh, nFaces] = mfmBuildSubMesh(mesh, params, perimeterEdges, insideNodes, ...
%    orderedUniquePerimeterPoints)
%
% Author: Winawer
% Purpose:
%   This routine begins with the original large mesh and extracts a
%   topologically correct sub-mesh based on the perimeter edges and inside nodes.
%
%   Sub-routine derived from Alex's unfoldMeshFromGUI code.
%
% See Also:  unfoldMeshFromGUI
%

% Get all the variables we may need
meshFileName     = params.meshFileName;
grayFileName     = params.grayFileName;
flatFileName     = params.flatFileName;
startCoords      = params.startCoords;
scaleFactor      = params.scaleFactor;
perimDist        = params.perimDist;
statusHandle     = params.statusHandle;
busyHandle       = params.busyHandle;
spacingMethod    = params.spacingMethod;
adjustSpacing    = params.adjustSpacing;
gridSpacing      = params.gridSpacing;
showFigures      = params.showFigures;
saveExtra        = params.saveExtra;
truePerimDist    = params.truePerimDist;
hemi             = params.hemi;
nperims          = params.NPERIMS;
saveIntermeidate = params.SAVE_INTERMEDIATE;
numberOfSteps    = params.NUMBEROFSTEPS;

% internal points are the ones we want
unfoldMesh.connectionMatrix=mesh.connectionMatrix(insideNodes,insideNodes);
unfoldMesh.normal=mesh.normal(insideNodes,:);

unfoldMesh.uniqueVertices=mesh.uniqueVertices(insideNodes,:);
unfoldMesh.uniqueCols=mesh.uniqueCols(insideNodes,:);
unfoldMesh.dist=mesh.dist(insideNodes);

% We need to get uniqueFaceIndexList for the unfoldMesh
indicesOfFacesInSubGroup=findFacesInGroup(mesh,insideNodes);
subGroupFaces=mesh.uniqueFaceIndexList(indicesOfFacesInSubGroup,:);
nFaces=size(subGroupFaces,1);

statusStringAdd(statusHandle,'Computing sub-mesh lookup table');

% Get a lookup table for converting indices into the full node array into indices to the unfold mesh nodes.
lookupTable=zeros(length(mesh.uniqueVertices),1);
lookupTable(insideNodes)=1:length(insideNodes);

% Use the lookup table to convert our list of face indices so that they index into unfoldMesh.uniqueVertices.
sgf=lookupTable(subGroupFaces(:));
unfoldMesh.uniqueFaceIndexList=reshape(sgf,nFaces,3);

% Convert the edges to feed into orderMeshPerimeterPoints
fullEdgePointList=perimeterEdges(:);

% How many edges do we have?
[numEdges,x]=size(perimeterEdges);

newEdges=zeros((numEdges*2),1);
statusStringAdd(statusHandle,'Finding sub-mesh edges.');

for t=1:(numEdges*2)
    if (~mod(t,100) & ~isempty(busyHandle))
        updateBusybar(busyHandle,t);
    end
    newEdges(t)=find(insideNodes==fullEdgePointList(t));
end

newEdges=reshape(newEdges,numEdges,2);
statusStringAdd(statusHandle,'Finding sub-mesh perim.');

% Find the perimeter points in the sub mesh.
unfoldMesh.orderedUniquePerimeterPoints=zeros(length(orderedUniquePerimeterPoints),1);

for t=1:length(orderedUniquePerimeterPoints)
    f1=find(insideNodes==orderedUniquePerimeterPoints(t));

    %orderMeshPerimeterPoints(newEdges);
    unfoldMesh.orderedUniquePerimeterPoints(t)=f1;
end

return;


text:
% add the path of RBM code
addpath('..');
addpath('~/work/Algorithms/liblinear-1.7/matlab');

% load MNIST
load 'mnist_14x14.mat';

% shuffle the training data
perm_idx = randperm (size(X,1));

n_all = size(X, 1);
n_train = ceil(n_all * 3 / 4);
n_valid = floor(n_all /4);

X_valid = X(perm_idx(n_train+1:end), :);
X_valid_labels = X_labels(perm_idx(n_train+1:end));
X = X(perm_idx(1:n_train), :);
X_labels = X_labels(perm_idx(1:n_train));

layers = [size(X,2), 200, 100, 50, 2];
n_layers = length(layers);
blayers = [1, 1, 1, 1, 0];

use_tanh = 0;
do_pretrain = 1;

if do_pretrain
    Ds = cell(n_layers - 1, 1);
    H = X;
    H_valid = X_valid;

    for l = 1:n_layers-1
        % construct DAE and use default configurations
        D = default_dae (layers(l), layers(l+1));

        D.data.binary = blayers(l);
        D.hidden.binary = blayers(l+1);

        if use_tanh
            if l > 1
                D.visible.use_tanh = 1;
            end
            D.hidden.use_tanh = 1;
        else
            if D.data.binary
                mH = mean(H, 1)';
                D.vbias = min(max(log(mH./(1 - mH)), -4), 4);
            else
                D.vbias = mean(H, 1)';
            end
        end

        D.learning.lrate = 1e-1;
        D.learning.lrate0 = 5000;
        D.learning.weight_decay = 0.0001;
        D.learning.minibatch_sz = 128;

        D.valid_min_epochs = 10;

        D.noise.drop = 0.2;
        D.noise.level = 0;

        %D.adagrad.use = 1;
        %D.adagrad.epsilon = 1e-8;
        D.adagrad.use = 0;
        D.adadelta.use = 1;
        D.adadelta.epsilon = 1e-8;
        D.adadelta.momentum = 0.99;

        D.iteration.n_epochs = 500;

        % save the intermediate data after every epoch
        D.hook.per_epoch = {@save_intermediate, {sprintf('dae_mnist_%d.mat', l)}};

        % print learining process
        D.verbose = 0;
        % display the progress
        D.debug.do_display = 0;

        % train RBM
        fprintf(1, 'Training DAE (%d)\n', l);
        tic;
        D = dae (D, H, H_valid, 0.1);
        fprintf(1, 'Training is done after %f seconds\n', toc);

        H = dae_get_hidden(H, D);
        H_valid = dae_get_hidden(H_valid, D);

        Ds{l} = D;
    end
end

S = default_sdae (layers);

S.data.binary = blayers(1);
S.bottleneck.binary = blayers(end);
S.hidden.use_tanh = use_tanh;

S.hook.per_epoch = {@save_intermediate, {'sdae_mnist.mat'}};

S.learning.lrate = 1e-1;
S.learning.lrate0 = 5000;
%S.learning.momentum = 0.9;
S.learning.weight_decay = 0.0001;
S.learning.minibatch_sz = 128;

%S.noise.drop = 0.2;
%S.noise.level = 0;
S.adadelta.use = 1;
S.adadelta.epsilon = 1e-8;
S.adadelta.momentum = 0.99;

%S.adagrad.use = 1;
%S.adagrad.epsilon = 1e-8;
S.valid_min_epochs = 10;

S.iteration.n_epochs = 100;

if do_pretrain
    for l = 1:n_layers-1
        S.biases{l+1} = Ds{l}.hbias;
        S.W{l} = Ds{l}.W;
    end
else
    if S.data.binary
        mH = mean(X, 1)';
        S.biases{1} = min(max(log(mH./(1 - mH)), -4), 4);
    else
        S.biases{1} = mean(X, 1)';
    end
end

fprintf(1, 'Training sDAE\n');
tic;
S = sdae (S, X, X_valid, 0.1);
fprintf(1, 'Training is done after %f seconds\n', toc);

H = sdae_get_hidden (X, S);
save 'sdae_mnist_vis.mat' H X_labels;

vis_mnist;

author: kyunghyuncho

repo: deepmat


text:
function rgb = skin_light

% This returns a predefined color as [red green blue] values
%   red               = [255   0   0]/255;
%   green             = [  0 192   0]/255;
%   blue              = [  0   0 255]/255;
%   magenta           = [255 255   0]/255;
%   cyan              = [  0 255 255]/255;
%   yellow            = [255 255   0]/255;
%   white             = [255 255 255]/255;
%   black             = [  0   0   0]/255;
%
%   skull             = [140  85  85]/255
%   cortex            = [255 213 119]/255;
%   cortex_light      = [199 194 169]/255;
%   cortex_dark       = [100  97  85]/255;
%   skin              = [249 223 192]/255;
%   skin_light        = [249 223 192]/255;
%   skin_medium_light = [225 194 158]/255;
%   skin_medium       = [188 142 106]/255;
%   skin_medium_dark  = [155 102	65]/255;
%   skin_dark         = [ 91  71  61]/255;
%
% The different skin-based colors follow the Fitzpatrick scale with type I and II
% combined, and return RGB values that approximate those used by Apple in the emoji
% skin tones. See also https://emojipedia.org/emoji-modifier-sequence/
%
% If no specific skin tone is specified, this function returns a light skin color.
% This corresponds with that of one of the developers who approximated his own skin
% color more than 15 years ago upon the first implementation of this function.

% returns a predefined color as [red green blue] values
rgb = [249	223	192]/255;


text:
function fgWrite(fg,name,type)
% Write out a fiber group structure for use with Quench. Pdb format.
%
%  fgWrite(fg,[name],[type])
%
% INPUTS:
%     fg - Fiber group structure.
%   name - Name used to save the fiber group. Can bee a full path
%          terminating in a name. Defaults to the current directory with
%          the name in the fg.name field.
%   type - File type to save as.
%            Options are:
%             'pdb'    - quench file format [DEFAULT]
%             'quench' - same as 'pdb'
%             'mat'    - matlab file format
%             'tck'    - MRTRIX file format (also 'mrtrix'/'.tck')
%
%          * If "name" ends in either .pdb or .mat the correct file type
%            will be used.
%
% WEB RESOURCES:
%   mrvBrowseSVN('fgWrite');
%
% EXAMPLE USAGE:
%   fgWrite(fg,fg.name,'pdb');
%
% See Also:
%   fgRead.m
%
%
% (C) Stanford VISTA, 2011
%
%     Brent McPherson and Franco Pestilli
%     Indiana University 2017

%% Check inputs

% Check that fg is a structure
if ~isstruct(fg) || ~isfield(fg, 'fibers')
    error('fg must be a fiber group structure with a fibers field.');
end

% Should we do this?
% This will set fg.params = [] and fg.pathwayInfo = [].
% fg = dtiClearQuenchStats(fg);

% Check for name input variable and use fg.name if empty
if ~exist('name','var') || isempty(name)
    name = fgGet(fg,'name');
end

% Check type and set to pdb by default
if ~exist('type','var') || isempty(type)
    type = 'pdb';
end

% If the name ends in '.mat' set type to 'mat'
[p n e] = fileparts(name);
if strcmp(e,'.mat')
    type = 'mat';
end

%% Write out the fibers
switch type
    case {'pdb', 'quench'}
        % Arguments are:
        % fg, fileName, xformToAcPc, xform, fiberToStrFlag, version
        mtrExportFibers(fg, name, [], [], [], 3);
    case 'mat'
        dtiWriteFiberGroup(fg, name);
    case {'tck', 'mrtrix', '.tck'}
        dtiExportFibersMrtrix(fg, name);
end

return