text:
function out = PP_ModelFit(y,model,order,randomSeed)
% PP_ModelFit   Investigates whether AR model fit improves with different preprocessings.
%
% After performing the range of transformations to the time series, returns the
% in-sample root-mean-square (RMS) prediction errors for an AR model on each
% transformed time series as a ratio of the RMS prediction error of the original
% time series.
%
% PP_PreProcess.m is used to perform the preprocessings
%
% The AR model is fitted using the function ar and pe from Matlab's System
% Identification Toolbox
%
% Transformations performed include:
% (i) incremental differencing,
% (ii) filtering of the power spectral density function,
% (iii) removal of piece-wise polynomial trends, and
% (iv) rank mapping the values of the time series to a Gaussian distribution.
%
%---INPUTS:
%
% y, the input time series
% model, the time-series model to fit to the transformed time series (currently
%           'ar' is the only option)
%
% order, the order of the AR model to fit to the data
%
% randomSeed, whether (and how) to reset the random seed, using BF_ResetSeed

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Preliminaries
% ------------------------------------------------------------------------------
N = length(y); % length of the time series

% ------------------------------------------------------------------------------
%% Inputs
% ------------------------------------------------------------------------------
% Model: the model to fit preprocessed time series to
if nargin < 2 || isempty(model)
    model = 'ar';
end

% order: the order of model to fit
if nargin < 3 || isempty(order)
    order = 2;
end

% randomSeed: how to treat the randomization
if nargin < 4
    randomSeed = []; % default
end

% ------------------------------------------------------------------------------
%% Apply a range of preprocessings
% ------------------------------------------------------------------------------
yp = PP_PreProcess(y,'',[],[],[],randomSeed);
% Returns a structure, yp, with a range of time series in it, each a different
% transformation of the original, y.
%% ____________________FIT MODEL TO ALL:_______________________ %%

fields = fieldnames(yp);
numFields = length(fields);
% statstore = struct('fpes',{});

for i = 1:numFields
    % for each preprocessing, fit the model
    data = yp.(fields{i});
    % data is the current preprocessed data

    switch model % SO MANY OPTIONS! ;-)
        case 'ar'
            %% Check that a System Identification Toolbox license is available
            BF_CheckToolbox('identification_toolbox')

            data = zscore(data); % zscore the data from this preprocessing
            m = ar(data,order); % fit the model

            % Get statistics on fit
            %     () FPE
            statstore.fpe(i) = m.EstimationInfo.FPE;
            %     () in-sample prediction error
            e = pe(m,data);
            statstore.rmserr(i) = sqrt(mean(e.^2));
            statstore.mabserr(i) = mean(abs(e));
            statstore.ac1(i) = CO_AutoCorr(e,1,'Fourier');

        otherwise
            error('Unknown model ''%s''',model);
    end
end

% ------------------------------------------------------------------------------
%% Return statistics on statistics
% ------------------------------------------------------------------------------
% actually often as you make more stationary and remove trends it becomes
% harder to predict because these trends are very easy to predict, and
% making the series whiter will obviously decrease its predictability.

% (1) ratio of fpe of preprocessed to unprocessed time series
% I think just this is ok.
% for i=2:numFields
%     eval(['out.fperat_' fields{i} ' = ' num2str(statstore.fpe(i)/statstore.fpe(1)) ';']);
% end

% No, I'll just do in-sample rms error, for a single model no point fpeing
for i = 2:numFields
    out.(sprintf('rmserrrat_%s',fields{i})) = statstore.rmserr(i)/statstore.rmserr(1);
end
% In fact, greater error in this case means a better detrending in some
% sense -- it's remobed more of the 'obvious' linear structure (assuming
% that's the aim).

% could also return statistics on other things like prediction error, but
% not alot of point, I think.

%
%     function ydt =  SUB_remps(y,n,method)
%         % Removes the first n (proportion) of power spectrum
%         % Based on my deseasonalize1.m code
%
%
%         %% Take the Fourier Transform
%
%         Ny = length(y); % number of samples in y
% %         t = linspace(0,1,Ny); % time vector
%         NFFT = 2^nextpow2(Ny); % next power of 2
%         Fy = fft(y,NFFT); % fast fourier transform of y
%         Fy1 = Fy(1:NFFT/2+1);
% %         f = 1/2*linspace(0,1,NFFT/2+1); % frequency vector
%
%         %% Remove this range
%         % set it to (mean of the rest) across this range
%         switch method
%             case 'lf'
%                 cullr = 1:floor(length(Fy1)*n);
%             case 'biggest'
%                 cullr = find(abs(Fy1)>quantile(abs(Fy1),n));
%         end
%
%         meanrest = mean(abs(Fy1(setxor(1:end,cullr))));
% %         meanrest = 0;
%         FyF = Fy;
%         FyF(cullr)=meanrest;
%         FyF(end-cullr+2)=meanrest;
%
%
%         % PLOT
% %         plot(abs(Fy)),hold on; plot(abs(FyF),'--r'); hold off
% %         input('Here''s the filtered one...')
% %         plot(abs(FyF),'k');
% %         input('Again on it''s own...')
%
%
%         %% Inverse Fourier Transform
%         ydt = ifft(FyF,NFFT);
%         ydt = zscore(ydt(1:Ny)); % crop to desired length
%
%
%         % PLOT
% %         plot(zscore(ydt),'b'); hold on; plot(y,'r'); hold off;
% %         input(['Mean difference is ' num2str(mean(y-ydt))])
%
%     end
%
%     function ydt = SUB_rempt(y,order,nbits)
%         N = length(y);
%         ydt = zeros(N,1);
%         bits = round(linspace(0,N,nbits+1));
%         for k=1:nbits
%             r = bits(k)+1 : bits(k+1); % range defined by adjacent 'bits'
%             x = (1:length(r))'; % faux x-range
%             ybit = y(r); % y-range
%             p = polyfit(x,ybit,order);
%             ydt(r) = ybit-polyval(p,x);
%         end
%         ydt = zscore(ydt);
% %         plot(y,'b'); hold on; plot(ydt,'r');
% %         input('here we are')
%     end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function ds = spm_uw_estimate(P,par)
% Estimation of partial derivatives of EPI deformation fields
%
% FORMAT [ds] = spm_uw_estimate((P),(par))
%
% P            - List of file names or headers.
% par          - Structure containing parameters governing the specifics
%                of how to estimate the fields.
% .M           - When performing multi-session realignment and Unwarp we
%                want to realign everything to the space of the first
%                image in the first time-series. M defines the space of
%                that.
% .order       - Number of basis functions to use for each dimension.
%                If the third dimension is left out, the order for
%                that dimension is calculated to yield a roughly
%                equal spatial cut-off in all directions.
%                Default: [12 12 *]
% .sfP         - Static field supplied by the user. It should be a
%                filename or handle to a voxel-displacement map in
%                the same space as the first EPI image of the time-
%                series. If using the FieldMap toolbox, realignment
%                should (if necessary) have been performed as part of
%                the process of creating the VDM. Note also that the
%                VDM must be in undistorted space, i.e. if it is
%                calculated from an EPI based field-map sequence
%                it should have been inverted before passing it to
%                spm_uw_estimate. Again, the FieldMap toolbox will
%                do this for you.
% .regorder    - Regularisation of derivative fields is based on the
%                regorder'th (spatial) derivative of the field.
%                Default: 1
% .lambda      - Fudge factor used to decide relative weights of
%                data and regularisation.
%                Default: 1e5
% .jm          - Jacobian Modulation. If set, intensity (Jacobian)
%                deformations are included in the model. If zero,
%                intensity deformations are not considered.
% .fot         - List of indexes for first order terms to model
%                derivatives for. Order of parameters as defined
%                by spm_imatrix.
%                Default: [4 5]
% .sot         - List of second order terms to model second
%                derivatives of. Should be an nx2 matrix where
%                e.g. [4 4; 4 5; 5 5] means that second partial
%                derivatives of rotation around x- and y-axis
%                should be modelled.
%                Default: []
% .fwhm        - FWHM (mm) of smoothing filter applied to images prior
%                to estimation of deformation fields.
%                Default: 6
% .rem         - Re-Estimation of Movement parameters. Set to unity means
%                that movement-parameters should be re-estimated at each
%                iteration.
%                Default: 0
% .noi         - Maximum number of Iterations.
%                Default: 5
% .exp_round   - Point in position space to do Taylor expansion around.
%                'First', 'Last' or 'Average'.
%                Default: 'Average'.
% ds           - The returned structure contains the following fields
% .P           - Copy of P on input.
% .sfP         - Copy of sfP on input (if non-empty).
% .order       - Copy of order on input, or default.
% .regorder    - Copy of regorder on input, or default.
% .lambda      - Copy of lambda on input, or default.
% .fot         - Copy of fot on input, or default.
% .sot         - Copy of sot on input, or default.
% .fwhm        - Copy of fwhm on input, or default.
% .rem         - Copy of rem on input, or default.
% .p0          - Average position vector (three translations in mm
%                and three rotations in degrees) of scans in P.
% .q           - Deviations from mean position vector of modelled
%                effects. Corresponds to deviations (and deviations
%                squared) of a Taylor expansion of deformation fields.
% .beta        - Coeffeicents of DCT basis functions for partial
%                derivatives of deformation fields w.r.t. modelled
%                effects. Scaled such that resulting deformation
%                fields have units mm^-1 or deg^-1 (and squares
%                thereof).
% .SS          - Sum of squared errors for each iteration.
%
%__________________________________________________________________________
% Copyright (C) 2003-2011 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: spm_uw_estimate.m 6824 2016-06-27 16:19:51Z guillaume $

% This is a major rewrite which uses some new ideas to speed up
% the estimation of the field. The time consuming part is the
% evaluation of A'*A where A is a matrix with the partial
% derivatives for each scan with respect to the parameters
% describing the warp-fields. If we denote the derivative
% matrix for a single scan by Ai, then the estimation of A'*A
% is A'*A = A1'*A1 + A2'*A2 + ... +An'*An where n is the number
% of scans in the time-series. If we model the partial-derivative
% fields w.r.t. two movement parameters (e.g. pitch and roll), each
% by [8 8 8] basis-functions and the image dimensions are
% 64x64x64 then each Ai is a 262144x1024 matrix and we need to
% evaluate and add n of these 1024x1024 matrices. It takes a while.
%
% The new idea is based on the realisation that each of these
% matrices is the kroneceker-product of the relevant movement
% parameters, our basis set and a linear combination (given by
% one column of the inverse of the rotation matrix) of the image
% gradients in the x-, y- and z-directions. This means that
% they really aren't all that unique, and that the amount of
% information in these matrices doesn't really warrant all
% those calculations. After a lot of head-scratching I arrived
% at the following
%
% First some definitions
%
% n: no. of voxels
% m: no. of scans
% l: no. of effects to model
% order: [xorder yorder zorder] no. of basis functions
% q: mxl matrix of scaled realignment parameters for effects to model.
% T{i}: inv(inv(P(i).mat)*P(1).mat);
% t(i,:) = T{i}(1:3,2)';
% B: kron(Bz,By,Bx);
% Ax = repmat(dx,1,prod(order)).*B;
% Ay = repmat(dy,1,prod(order)).*B;
% Az = repmat(dz,1,prod(order)).*B;
%
% Now, my hypothesis is that A'*A is given by
%
% AtA = kron((q.*kron(ones(1,l),t(:,1)))'*(q.*kron(ones(1,l),t(:,1))),Ax'*Ax) +...
% kron((q.*kron(ones(1,l),t(:,2)))'*(q.*kron(ones(1,l),t(:,2))),Ay'*Ay) +...
% kron((q.*kron(ones(1,l),t(:,3)))'*(q.*kron(ones(1,l),t(:,3))),Az'*Az) +...
% 2*kron((q.*kron(ones(1,l),t(:,1)))'*(q.*kron(ones(1,l),t(:,2))),Ax'*Ay) +...
% 2*kron((q.*kron(ones(1,l),t(:,1)))'*(q.*kron(ones(1,l),t(:,3))),Ax'*Az) +...
% 2*kron((q.*kron(ones(1,l),t(:,2)))'*(q.*kron(ones(1,l),t(:,3))),Ay'*Az);
%
% Which turns out to be true. This means that regardless of how many
% scans we have we will always be able to create AtA as a sum of
% six individual AtAs. It has been tested against the previous
% implementation and yields exactly identical results.
%
% I know this isn't much of a derivation, but there will be a paper soon. Sorry.
%
% Other things that are new for the rewrite is
% 1. We have removed the possibility to try and estimate the "static" field.
%    There simply isn't enough information about that field, and for a
%    given time series the estimation is just too likely to fail.
% 2. New option to pass a static field (e.g. estimated from a dual
%    echo-time measurement) as an input parameter.
% 3. Re-estimation of the movement parameters at each iteration.
%    Let us say we have a nodding motion in the time series, and
%    at each nod the brain appear to shrink in the phase-encode
%    direction. The realignment will find the nods, but might in
%    addition mistake the "shrinks" for translations, leading to
%    biased estimates. We attempt to correct that by, at iteration,
%    updating both parameters pertaining to distortions and to
%    movement. In order to do so in an unbiased fashion we have
%    also switched from sampling of images (and deformation fields)
%    on a grid centered on the voxel-centers, to a grid whith a
%    different sampling frequency.
% 4. Inclusion of a regularisation term. The speed-up has facilitated
%    using more basis-functions, which makes it neccessary to impose
%    regularisation (i.e. punisihing some order derivative of the
%    deformation field).
% 5. Change of interpolation model from tri-linear/Sinc to
%    tri-linear/B-spline.
% 6. Option to include the Jacobian compression/stretching effects
%    in the model for the estimation of the displacement fields.
%    Our tests have indicated that this is NOT a good idea though.

SVNid = '$Rev: 6824 $';

%-Say hello
%--------------------------------------------------------------------------
SPMid = spm('FnBanner',mfilename,SVNid);

%-Parameters
%--------------------------------------------------------------------------
if nargin < 1 || isempty(P), P = spm_select(Inf,'image'); end
if ~isstruct(P), P = spm_vol(P); end

%
% Hardcoded default input parameters.
%
defpar = struct('sfP',             [],...
                'M',               P(1).mat,...
                'fot',             [4 5],...
                'sot',             [],...
                'hold',            [1 1 1 0 1 0]);

%
% Merge hardcoded defaults and spm_defaults. Translate spm_defaults
% settings to internal defaults.
%
ud = spm_get_defaults('unwarp.estimate');
if isfield(ud,'basfcn'),    defpar.order = ud.basfcn; end
if isfield(ud,'regorder'),  defpar.regorder = ud.regorder; end
if isfield(ud,'regwgt'),    defpar.lambda = ud.regwgt; end
if isfield(ud,'jm'),        defpar.jm = ud.jm; end
if isfield(ud,'fwhm'),      defpar.fwhm = ud.fwhm; end
if isfield(ud,'rem'),       defpar.rem = ud.rem; end
if isfield(ud,'noi'),       defpar.noi = ud.noi; end
if isfield(ud,'expround'),  defpar.exp_round = ud.expround; end

defnames = fieldnames(defpar);

%
% Go through input parameters, chosing the default
% for any parameters that are missing, warning the
% user if there are "unknown" parameters (probably
% reflecting a misspelling).
%

if nargin < 2 || isempty(par)
   par = defpar;
end
ds = [];
for i=1:length(defnames)
   if isfield(par,defnames{i}) && ~isempty(par.(defnames{i}))
      ds.(defnames{i}) = par.(defnames{i});
   else
      ds.(defnames{i}) = defpar.(defnames{i});
   end
end
parnames = fieldnames(par);
for i=1:length(parnames)
   if ~isfield(defpar,parnames{i})
      warning('Unknown par field %s',parnames{i});
   end
end

%
% Resolve ambiguities.
%
if length(ds.order) == 2
   mm          = sqrt(sum(P(1).mat(1:3,1:3).^2)).*P(1).dim(1:3);
   ds.order(3) = round(ds.order(1)*mm(3)/mm(1));
end
if isfield(ds,'sfP') && ~isempty(ds.sfP)
   if ~isstruct(ds.sfP)
      ds.sfP = spm_vol(ds.sfP);
   end
end
nscan = length(P);
nof   = prod(size(ds.fot)) + size(ds.sot,1);
ds.P  = P;

%
% Get matrix of 'expansion point'-corrected movement parameters
% for which we seek partial derivative fields.
%
[ds.q,ds.p0] = make_q(P,ds.fot,ds.sot,ds.exp_round);

%
% Create matrix for regularisation of warps.
%
H = make_H(P,ds.order,ds.regorder);

%
% Create a temporary smooth time series to use for
% the estimation.
%
old_P = P;
if ds.fwhm ~= 0
   spm_uw_show('SmoothStart',length(P));
   for i=1:length(old_P)
      spm_uw_show('SmoothUpdate',i);
      sfname(i,:) = [tempname spm_file_ext ',1,1'];
      to_smooth = sprintf('%s,%d,%d',old_P(i).fname,old_P(i).n);
      spm_smooth(to_smooth,sfname(i,:),ds.fwhm);
   end
   P = spm_vol(sfname);
   spm_uw_show('SmoothEnd');
end

% Now that we have littered the disk with smooth
% temporary files we  should use a try-catch
% block for the rest of the function, to ensure files
% get deleted in the event of an error.

try        % Try block starts here

%
% Initialize some stuff.
%
beta0    = zeros(nof*prod(ds.order),10);
beta     = zeros(nof*prod(ds.order),1);
old_beta = zeros(nof*prod(ds.order),1);

%
% Sample images on irregular grid to avoid biasing
% re-estimated movement parameters with smoothing
% effect from voxel-interpolation (See Andersson ,
% EJNM (1998), 25:575-586.).
%

ss         = [1.1 1.1 0.9];
xs         = 1:ss(1):P(1).dim(1); xs = xs + (P(1).dim(1)-xs(end)) / 2;
ys         = 1:ss(2):P(1).dim(2); ys = ys + (P(1).dim(2)-ys(end)) / 2;
zs         = 1:ss(3):P(1).dim(3); zs = zs + (P(1).dim(3)-zs(end)) / 2;
M          = spm_matrix([-xs(1)/ss(1)+1 -ys(1)/ss(2)+1 -zs(1)/ss(3)+1])*inv(diag([ss 1]))*eye(4);
nx         = length(xs);
ny         = length(ys);
nz         = length(zs);
[x,y,z]    = ndgrid(xs,ys,zs);
Bx         = spm_dctmtx(P(1).dim(1),ds.order(1),x(:,1,1));
By         = spm_dctmtx(P(1).dim(2),ds.order(2),y(1,:,1));
Bz         = spm_dctmtx(P(1).dim(3),ds.order(3),z(1,1,:));
dBy        = spm_dctmtx(P(1).dim(2),ds.order(2),y(1,:,1),'diff');
xyz        = [x(:) y(:) z(:) ones(length(x(:)),1)]; clear x y z;
def        = zeros(size(xyz,1),nof);
ddefa      = zeros(size(xyz,1),nof);

%
% Create file struct for use with spm_orthviews to draw
% representations of the field.
%
dispP       = P(1);
dispP       = rmfield(dispP,{'fname','descrip','n','private'});
dispP.dim   = [nx ny nz];
dispP.dt    = [64 spm_platform('bigend')];
dispP.pinfo = [1 0]';
p          = spm_imatrix(dispP.mat); p = p.*[zeros(1,6) ss 0 0 0];
p(1)       = -mean(1:nx)*p(7);
p(2)       = -mean(1:ny)*p(8);
p(3)       = -mean(1:nz)*p(9);
dispP.mat  = spm_matrix(p); clear p;

%
% We will need to resample the static field (if one was supplied)
% on the same grid (given by xs, ys and zs) as we are going to use
% for the time series. We will assume that the fieldmap has been
% realigned to the space of the first EPI image in the time-series.
%
if isfield(ds,'sfP') && ~isempty(ds.sfP)
   T = ds.sfP.mat\ds.M;
   txyz = xyz*T(1:3,:)';
   c = spm_bsplinc(ds.sfP,ds.hold);
   ds.sfield = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds.hold);
   ds.sfield = ds.sfield(:);
   clear c txyz;
else
   ds.sfield = [];
end

msk = get_mask(P,xyz,ds,[nx ny nz]);
ssq = [];
%
% Here starts iterative search for deformation fields.
%
for iter=1:ds.noi
   spm_uw_show('NewIter',iter);

   [ref,dx,dy,dz,P,ds,sf,dispP] = make_ref(ds,def,ddefa,P,xyz,M,msk,beta,dispP);
   AtA       = build_AtA(ref,dx,dy,dz,Bx,By,Bz,dBy,ds,P);
   [Aty,yty] = build_Aty(ref,dx,dy,dz,Bx,By,Bz,dBy,ds,P,xyz,beta,sf,def,ddefa,msk);

   % Clean up a bit, cause inverting AtA may use a lot of memory
   clear ref dx dy dz

   % Check that residual error still decreases.
   if iter > 1 && yty > ssq(iter-1)
      %
      % This means previous iteration was no good,
      % and we should go back to old_beta.
      %
      beta = old_beta;
      break;
   else
      ssq(iter) = yty;

      spm_uw_show('StartInv',1);

      % Solve for beta
      Aty = Aty + AtA*beta;
      AtA = AtA + ds.lambda * kron(eye(nof),diag(H)) * ssq(iter)/(nscan*sum(msk));

      try % Fastest if it works
         beta0(:,iter) = AtA\Aty;
      catch % Sometimes necessary
         beta0(:,iter) = pinv(AtA)*Aty;
      end

      old_beta = beta;
      beta     = beta0(:,iter);

      for i=1:nof
         def(:,i)   = spm_get_def(Bx, By,Bz,beta((i-1)*prod(ds.order)+1:i*prod(ds.order)));
         ddefa(:,i) = spm_get_def(Bx,dBy,Bz,beta((i-1)*prod(ds.order)+1:i*prod(ds.order)));
      end

      % If we are re-estimating the movement parameters, remove any DC
      % components from the deformation fields, so that they end up in
      % the movement-parameters instead. It shouldn't make any difference
      % to the variance reduction, but might potentially lead to a better
      % explanation of the variance components.
      % Note that since we sub-sample the images (and the DCT basis set)
      % it is NOT sufficient to reset beta(1) (and beta(prod(order)+1) etc,
      % instead we explicitly subtract the DC component. Note that a DC
      % component does NOT affect the Jacobian.
      %
      if ds.rem ~= 0
         def = def - repmat(mean(def),length(def),1);
      end
      spm_uw_show('EndInv');
      tmp = dispP.mat;
      dispP.mat = P(1).mat;
      spm_uw_show('FinIter',ssq,def,ds.fot,ds.sot,dispP,ds.q);
      dispP.mat = tmp;
   end
   clear AtA
end

ds.P = old_P;
for i=1:length(ds.P)
   ds.P(i).mat = P(i).mat;    % Save P with new movement parameters.
end
ds.beta = reshape(refit(P,dispP,ds,def),prod(ds.order),nof);
ds.SS   = ssq;
if isfield(ds,'sfield');
   ds = rmfield(ds,'sfield');
end

cleanup(P,ds)
spm_uw_show('FinTot');

% Document outcome
F = spm_figure('FindWin','Graphics');
if ~isempty(F), spm_figure('Focus', F), spm_print; end

catch
   cleanup(P,ds)
   spm_uw_show('FinTot');
   fprintf('Unwarp terminated abnormally.\n');
   rethrow(lasterror);
end

fprintf('%-40s: %30s\n','Completed',spm('time'))                        %-#

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Utility functions.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [q,ep] = make_q(P,fot,sot,exp_round)
%
% P         : Array of file-handles
% fot       : nx1 array of indicies for first order effect.
% sot       : nx2 matrix of indicies of second order effects.
% exp_round : Point in position space to perform Taylor-expansion
%             around ('First','Last' or 'Average'). 'Average' should
%             (in principle) give the best variance reduction. If a
%             field-map acquired before the time-series is supplied
%             then expansion around the 'First' MIGHT give a slightly
%             better average geometric fidelity.

if strcmpi(exp_round,'average');
   %
   % Get geometric mean of all transformation matrices.
   % This will be used as the zero-point in the space
   % of object position vectors (i.e. the point at which
   % the partial derivatives are estimated). This is presumably
   % a little more correct than taking the average of the
   % parameter estimates.
   %
   mT = zeros(4);
   for i=1:length(P)
      mT = mT+logm(inv(P(i).mat) * P(1).mat);
   end
   mT = real(expm(mT/length(P)));

elseif strcmpi(exp_round,'first');
   mT = eye(4);
elseif strcmpi(exp_round,'last');
   mT = inv(P(end).mat) * P(1).mat;
else
   warning(sprintf('Unknown expansion point %s',exp_round));
end

%
% Rescaling to degrees makes translations and rotations
% roughly equally scaled. Since the scaling influences
% strongly the effects of regularisation, it would surely
% be nice with a more principled approach. Suggestions
% anyone?
%
ep = [1 1 1 180/pi 180/pi 180/pi zeros(1,6)] .* spm_imatrix(mT);

%
% Now, get a nscan-by-nof matrix containing
% mean (expansion point) corrected values for each effect we
% model.
%
q = zeros(length(P),prod(size(fot)) + size(sot,1));
for i=1:length(P)
   p = [1 1 1 180/pi 180/pi 180/pi zeros(1,6)] .*...
      spm_imatrix(inv(P(i).mat) * P(1).mat);
   for j=1:prod(size(fot))
      q(i,j) = p(fot(j)) - ep(fot(j));
   end
   for j=1:size(sot,1)
      q(i,prod(size(fot))+j) = (p(sot(j,1)) - ep(sot(j,1))) * (p(sot(j,2)) - ep(sot(j,2)));
   end
end
return
%_______________________________________________________________________

%_______________________________________________________________________
function H = make_H(P,order,regorder)
% Utility Function to create Regularisation term of AtA.

mm = sqrt(sum(P(1).mat(1:3,1:3).^2)).*P(1).dim(1:3);
kx=(pi*((1:order(1))'-1)/mm(1)).^2;
ky=(pi*((1:order(2))'-1)/mm(2)).^2;
kz=(pi*((1:order(3))'-1)/mm(3)).^2;

if regorder == 0
   %
   % Cost function based on sum of squares
   %
   H = (1*kron(kz.^0,kron(ky.^0,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^0,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^0,kx.^0)) );
elseif regorder == 1
   %
   % Cost function based on sum of squared 1st derivatives
   %
   H  = (1*kron(kz.^1,kron(ky.^0,kx.^0)) +...
         1*kron(kz.^0,kron(ky.^1,kx.^0)) +...
         1*kron(kz.^0,kron(ky.^0,kx.^1)) );
elseif regorder == 2
   %
   % Cost function based on sum of squared 2nd derivatives
   %
   H = (1*kron(kz.^2,kron(ky.^0,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^2,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^0,kx.^2)) +...
        3*kron(kz.^1,kron(ky.^1,kx.^0)) +...
        3*kron(kz.^1,kron(ky.^0,kx.^1)) +...
        3*kron(kz.^0,kron(ky.^1,kx.^1)) );
elseif regorder == 3
   %
   % Cost function based on sum of squared 3rd derivatives
   %
   H = (1*kron(kz.^3,kron(ky.^0,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^3,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^0,kx.^3)) +...
        3*kron(kz.^2,kron(ky.^1,kx.^0)) +...
        3*kron(kz.^2,kron(ky.^0,kx.^1)) +...
        3*kron(kz.^1,kron(ky.^2,kx.^0)) +...
        3*kron(kz.^0,kron(ky.^2,kx.^1)) +...
        3*kron(kz.^1,kron(ky.^0,kx.^2)) +...
        3*kron(kz.^0,kron(ky.^1,kx.^2)) +...
        6*kron(kz.^1,kron(ky.^1,kx.^1)) );
else
   error('Invalid order of regularisation');
end
return
%_______________________________________________________________________

%_______________________________________________________________________
function AtA = build_AtA(ref,dx,dy,dz,Bx,By,Bz,dBy,ds,P)
% Now lets build up the design matrix A, or rather A'*A  since
% A itself would be forbiddingly large.

if ds.jm, spm_uw_show('StartAtA',10);
else spm_uw_show('StartAtA',6); end

nof   = prod(size(ds.fot)) + size(ds.sot,1);

AxtAx = uwAtA1(dx.*dx,Bx,By,Bz); spm_uw_show('NewAtA',1);
AytAy = uwAtA1(dy.*dy,Bx,By,Bz); spm_uw_show('NewAtA',2);
AztAz = uwAtA1(dz.*dz,Bx,By,Bz); spm_uw_show('NewAtA',3);
AxtAy = uwAtA1(dx.*dy,Bx,By,Bz); spm_uw_show('NewAtA',4);
AxtAz = uwAtA1(dx.*dz,Bx,By,Bz); spm_uw_show('NewAtA',5);
AytAz = uwAtA1(dy.*dz,Bx,By,Bz); spm_uw_show('NewAtA',6);

if ds.jm
   AjtAj = uwAtA1(ref.*ref,Bx,dBy,Bz); spm_uw_show('NewAtA',7);
   AxtAj = uwAtA2( dx.*ref,Bx, By,Bz,Bx,dBy,Bz); spm_uw_show('NewAtA',8);
   AytAj = uwAtA2( dy.*ref,Bx, By,Bz,Bx,dBy,Bz); spm_uw_show('NewAtA',9);
   AztAj = uwAtA2( dz.*ref,Bx, By,Bz,Bx,dBy,Bz); spm_uw_show('NewAtA',10);
end

R = zeros(length(P),3);
for i=1:length(P)
   tmp    = inv(P(i).mat\P(1).mat);
   R(i,:) = tmp(1:3,2)';
end

tmp  = ones(1,nof);
tmp1 = ds.q.*kron(tmp,R(:,1));
tmp2 = ds.q.*kron(tmp,R(:,2));
tmp3 = ds.q.*kron(tmp,R(:,3));
AtA  =    kron(tmp1'*tmp1,AxtAx) + kron(tmp2'*tmp2,AytAy) + kron(tmp3'*tmp3,AztAz) +...
       2*(kron(tmp1'*tmp2,AxtAy) + kron(tmp1'*tmp3,AxtAz) + kron(tmp2'*tmp3,AytAz));

if ds.jm
   tmp  = [ds.q.*kron(tmp,ones(length(P),1))];
   AtA  = AtA + kron(tmp'*tmp,AjtAj) +...
             2*(kron(tmp1'*tmp,AxtAj) + kron(tmp2'*tmp,AytAj) + kron(tmp3'*tmp,AztAj));
end
spm_uw_show('EndAtA');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function AtA = uwAtA1(y,Bx,By,Bz)
% Calculating off-diagonal block of AtA.

[nx,mx] = size(Bx);
[ny,my] = size(By);
[nz,mz] = size(Bz);
AtA     = zeros(mx*my*mz);
for sl =1:nz
   tmp  = reshape(y((sl-1)*nx*ny+1:sl*nx*ny),nx,ny);
   spm_krutil(Bz(sl,:)'*Bz(sl,:),spm_krutil(tmp,Bx,By,1),AtA);
%    AtA  = AtA + kron(Bz(sl,:)'*Bz(sl,:),spm_krutil(tmp,Bx,By,1));
end
return
%_______________________________________________________________________

%_______________________________________________________________________
function AtA = uwAtA2(y,Bx1,By1,Bz1,Bx2,By2,Bz2)
% Calculating cross-term of diagonal block of AtA
% when A is a sum of the type A1+A2 and where both
% A1 and A2 are possible to express as a kronecker
% product of lesser matrices.

[nx,mx1] = size(Bx1); [nx,mx2] = size(Bx2);
[ny,my1] = size(By1); [ny,my2] = size(By2);
[nz,mz1] = size(Bz1); [nz,mz2] = size(Bz2);
AtA      = zeros(mx1*my1*mz1,mx2*my2*mz2);
for sl =1:nz
   tmp  = reshape(y((sl-1)*nx*ny+1:sl*nx*ny),nx,ny);
   spm_krutil(Bz1(sl,:)'*Bz2(sl,:),spm_krutil(tmp,Bx1,By1,Bx2,By2),AtA);
%   AtA  = AtA + kron(Bz1(sl,:)'*Bz2(sl,:),spm_krutil(tmp,Bx1,By1,Bx2,By2));
end
return
%_______________________________________________________________________

%_______________________________________________________________________
function [Aty,yty] = build_Aty(ref,dx,dy,dz,Bx,By,Bz,dBy,ds,P,xyz,beta,sf,def,ddefa,msk)
% Building Aty.

nof = prod(size(ds.fot)) + size(ds.sot,1);
Aty = zeros(nof*prod(ds.order),1);
yty = 0;

spm_uw_show('StartAty',length(P));
for scan = 1:length(P)
   spm_uw_show('NewAty',scan);
   T    = P(scan).mat\ds.M;
   txyz = xyz*T(1:3,:)';
   if ~(all(beta == 0) && isempty(ds.sfield))
      [idef,jac] = spm_get_image_def(scan,ds,def,ddefa);
      txyz(:,2)  = txyz(:,2)+idef;
   end;
   c    = spm_bsplinc(P(scan),ds.hold);
   y    = sf(scan) * spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds.hold);
   if ds.jm && ~(all(beta == 0) && isempty(ds.sfield))
      y = y .* jac;
   end;

   y_diff       = (ref - y).*msk;
   indx         = find(isnan(y));
   y_diff(indx) = 0;
   iTcol        = inv(T(1:3,1:3));
   tmpAty       = spm_get_def(Bx',By',Bz',([dx dy dz]*iTcol(:,2)).*y_diff);
   if ds.jm ~= 0
      tmpAty = tmpAty + spm_get_def(Bx',dBy',Bz',ref.*y_diff);
   end
   for i=1:nof
      rindx      = (i-1)*prod(ds.order)+1:i*prod(ds.order);
      Aty(rindx) = Aty(rindx) + ds.q(scan,i)*tmpAty;
   end
   yty = yty + y_diff'*y_diff;
end
spm_uw_show('EndAty');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [ref,dx,dy,dz,P,ds,sf,dispP] = make_ref(ds,def,ddefa,P,xyz,M,msk,beta,dispP)
% First of all get the mean of all scans given their
% present deformation fields. When using this as the
% "reference" we will explicitly be minimising variance.
%
% First scan in P is still reference in terms of
% y-direction (in the scanner framework). A single set of
% image gradients is estimated from the mean of all scans,
% transformed into the space of P(1), and used for all scans.
%
spm_uw_show('StartRef',length(P));

rem = ds.rem;
if all(beta==0), rem=0; end;

if rem
   [m_ref,D] = get_refD(ds,def,ddefa,P,xyz,msk);
end

ref = zeros(size(xyz,1),1);
for i=1:length(P)
   spm_uw_show('NewRef',i);
   T    = P(i).mat\ds.M;
   txyz = xyz*T(1:3,:)';
   if ~(all(beta == 0) && isempty(ds.sfield))
      [idef,jac] = spm_get_image_def(i,ds,def,ddefa);
      txyz(:,2)  = txyz(:,2)+idef;
   end;
   c     = spm_bsplinc(P(i),ds.hold);
   f     = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds.hold);
   if ds.jm ~= 0 && exist('jac')==1
      f = f .* jac;
   end
   indx  = find(~isnan(f));
   sf(i) = 1 / (sum(f(indx).*msk(indx)) / sum(msk(indx)));
   ref   = ref + f;

   if rem
      indx   = find(isnan(f)); f(indx) = 0;
      Dty{i} = (((m_ref - sf(i)*f).*msk)'*D)';
   end
end
ref  = ref / length(P);
ref  = reshape(ref,dispP.dim(1:3));
indx = find(~isnan(ref));

% Scale to roughly 100 mean intensity to ensure
% consistent weighting of regularisation.
gl             = spm_global(ref);
sf             = (100 * (sum(ref(indx).*msk(indx)) / sum(msk(indx))) / gl) * sf;
ref            = (100 / gl) * ref;
dispP.dat      = ref;
c              = spm_bsplinc(ref,ds.hold);
txyz           = xyz*M(1:3,:)';
[ref,dx,dy,dz] = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds.hold);
ref            = ref.*msk;
dx             =  dx.*msk;
dy             =  dy.*msk;
dz             =  dz.*msk;

% Re-estimate (or rather nudge) movement parameters.
if rem ~= 0
   iDtD = inv(D'*D);
   for i=2:length(P)
      P(i).mat = inv(spm_matrix((iDtD*Dty{i})'))*P(i).mat;
   end
   [ds.q,ds.p0] = make_q(P,ds.fot,ds.sot,ds.exp_round);
end
spm_uw_show('EndRef');
return
%_______________________________________________________________________

%_______________________________________________________________________
function [m_ref,D] = get_refD(ds,def,ddefa,P,xyz,msk)
% Get partials w.r.t. movements from first scan.

[idef,jac] = spm_get_image_def(1,ds,def,ddefa);
T    = P(1).mat\ds.M;
txyz = xyz*T';
c    = spm_bsplinc(P(1),ds.hold);
[m_ref,dx,dy,dz] = spm_bsplins(c,txyz(:,1),...
                   txyz(:,2)+idef,txyz(:,3),ds.hold);
indx    = ~isnan(m_ref);
mref_sf = 1 / (sum(m_ref(indx).*msk(indx)) / sum(msk(indx)));
m_ref(indx) = m_ref(indx) .* mref_sf; m_ref(~indx) = 0;
dx(indx)    =    dx(indx) .* mref_sf; dx(~indx)    = 0;
dy(indx)    =    dy(indx) .* mref_sf; dy(~indx)    = 0;
dz(indx)    =    dz(indx) .* mref_sf; dz(~indx)    = 0;
if ds.jm ~= 0
   m_ref = m_ref .* jac;
   dx    = dx    .* jac;
   dy    = dy    .* jac;
   dz    = dz    .* jac;
end
D = make_D(dx.*msk,dy.*msk,dz.*msk,txyz,P(1).mat);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function D = make_D(dx,dy,dz,xyz,M)
% Utility function that creates a matrix of partial
% derivatives in units of /mm and /radian.
%
% dx,dy,dz   - Partial derivatives (/pixel) wrt x-, y- and z-translation.
% xyz        - xyz-matrix (original positions).
% M          - Current voxel->world matrix.

D = zeros(length(xyz),6);
tiny = 0.0001;
for i=1:6
   p      = [0 0 0 0 0 0 1 1 1 0 0 0];
   p(i)   = p(i)+tiny;
   T      = M\spm_matrix(p)*M;
   dxyz   = xyz*T';
   dxyz   = dxyz(:,1:3)-xyz(:,1:3);
   D(:,i) = sum(dxyz.*[dx dy dz],2)/tiny;
end
return
%_______________________________________________________________________

%_______________________________________________________________________
function cleanup(P,ds)
% Delete temporary smooth files.
%
if ~isempty(ds.fwhm) && ds.fwhm > 0
   for i=1:length(P)
      spm_unlink(P(i).fname);
      [fpath,fname,ext] = fileparts(P(i).fname);
      spm_unlink(fullfile(fpath,[fname '.hdr']));
      spm_unlink(fullfile(fpath,[fname '.mat']));
   end
end
return;
%_______________________________________________________________________

%_______________________________________________________________________
function beta = refit(P,dispP,ds,def)
% We have now estimated the fields on a grid that
% does not coincide with the voxel centers. We must
% now calculate the betas that correspond to a
% a grid centered on the voxel centers.
%
spm_uw_show('StartRefit',1);
rsP     = P(1);
p       = spm_imatrix(rsP.mat);
p       = [zeros(1,6) p(7:9) 0 0 0];
p(1)    = -mean(1:rsP.dim(1))*p(7);
p(2)    = -mean(1:rsP.dim(2))*p(8);
p(3)    = -mean(1:rsP.dim(3))*p(9);
rsP.mat = spm_matrix(p); clear p;
[x,y,z] = ndgrid(1:rsP.dim(1),1:rsP.dim(2),1:rsP.dim(3));
xyz     = [x(:) y(:) z(:) ones(numel(x),1)];
txyz    = ((inv(dispP.mat)*rsP.mat)*xyz')';
Bx      = spm_dctmtx(rsP.dim(1),ds.order(1));
By      = spm_dctmtx(rsP.dim(2),ds.order(2));
Bz      = spm_dctmtx(rsP.dim(3),ds.order(3));
nx      = rsP.dim(1); mx = ds.order(1);
ny      = rsP.dim(2); my = ds.order(2);
nz      = rsP.dim(3); mz = ds.order(3);
nof     = numel(ds.fot) + size(ds.sot,1);
for i=1:nof
   dispP.dat   = reshape(def(:,i),dispP.dim(1:3));
   c           = spm_bsplinc(dispP,ds.hold);
   field       = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds.hold);
   indx        = isnan(field);
   wgt         = ones(size(field));
   wgt(indx)   = 0;
   field(indx) = 0;
   AtA         = zeros(mx*my*mz);
   Aty         = zeros(mx*my*mz,1);
   for sl = 1:nz
      indx  = (sl-1)*nx*ny+1:sl*nx*ny;
      tmp1  = reshape(field(indx).*wgt(indx),nx,ny);
      tmp2  = reshape(wgt(indx),nx,ny);
      spm_krutil(Bz(sl,:)'*Bz(sl,:),spm_krutil(tmp2,Bx,By,1),AtA);
%      AtA   = AtA + kron(Bz(sl,:)'*Bz(sl,:),spm_krutil(tmp2,Bx,By,1));
      Aty   = Aty + kron(Bz(sl,:)',spm_krutil(tmp1,Bx,By,0));
   end
   beta((i-1)*prod(ds.order)+1:i*prod(ds.order)) = AtA\Aty;
end
spm_uw_show('EndRefit');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function msk = get_mask(P,xyz,ds,dm)
%
% Create a mask to avoid regions where data doesnt exist
% for all scans. This mask is slightly less Q&D than that
% of version 1 of Unwarp. It checks where data exist for
% all scans with present movement parameters and given
% (optionally) the static field. It creates a mask with
% non-zero values only for those voxels, and then does a
% 3D erode to preempt effects of re-estimated movement
% parameters and movement-by-susceptibility effects.
%
spm_uw_show('MaskStart',length(P));
msk = true(length(xyz),1);
for i=1:length(P)
   txyz = xyz * (P(i).mat\ds.M)';
   tmsk = (txyz(:,1)>=1 & txyz(:,1)<=P(1).dim(1) &...
           txyz(:,2)>=1 & txyz(:,2)<=P(1).dim(2) &...
           txyz(:,3)>=1 & txyz(:,3)<=P(1).dim(3));
   msk  = msk & tmsk;
   spm_uw_show('MaskUpdate',i);
end
%
% Include static field in mask estimation
% if one has been supplied.
%
if isfield(ds,'sfP') && ~isempty(ds.sfP)
   txyz = xyz * (ds.sfP.mat\ds.M)';
   tmsk = (txyz(:,1)>=1 & txyz(:,1)<=ds.sfP.dim(1) &...
           txyz(:,2)>=1 & txyz(:,2)<=ds.sfP.dim(2) &...
           txyz(:,3)>=1 & txyz(:,3)<=ds.sfP.dim(3));
   msk  = msk & tmsk;
end
msk = erode_msk(msk,dm);
spm_uw_show('MaskEnd');

% maskP                 = P(1);
% [mypath,myname,myext] = fileparts(maskP.fname);
% maskP.fname           = fullfile(mypath,['mask' myext]);
% maskP.dim             = [nx ny nz];
% maskP.dt              = P(1).dt;
% maskP                 = spm_write_vol(maskP,reshape(msk,[nx ny nz]));
return;
%_______________________________________________________________________

%_______________________________________________________________________
function omsk = erode_msk(msk,dim)
omsk = zeros(dim+[4 4 4]);
omsk(3:end-2,3:end-2,3:end-2) = reshape(msk,dim);
omsk = spm_erode(omsk(:),dim+[4 4 4]);
omsk = reshape(omsk,dim+[4 4 4]);
omsk = omsk(3:end-2,3:end-2,3:end-2);
omsk = omsk(:);
return
%_______________________________________________________________________

author: spm


text:
%Author: Guosheng Lin

function do_forward_trn(obj, inputs)

%copy from  matconvnet/matlab/+dagnn/@DagNN/eval.m
% copy this function to matconvnet/matlab/+dagnn/@DagNN/

obj.computingDerivative = true;

% -------------------------------------------------------------------------
% Forward pass
% -------------------------------------------------------------------------

% set the input values
v = obj.getVarIndex(inputs(1:2:end)) ;
if any(isnan(v))
  broken = find(isnan(v)) ;
  error('No variable of name ''%s'' could be found in the DAG.', inputs{2*broken(1)-1}) ;
end
[obj.vars(v).value] = deal(inputs{2:2:end}) ;
inputs = [] ;

obj.numPendingVarRefs = [obj.vars.fanout] ;
for l = obj.executionOrder
  time = tic ;
  obj.layers(l).block.forwardAdvanced(obj.layers(l)) ;
  obj.layers(l).forwardTime = toc(time) ;
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function fx = exp_cordic ( x, n )

%*****************************************************************************80
%
%% EXP_CORDIC evaluates the exponential function using the CORDIC method.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frederick Ruckdeschel,
%    BASIC Scientific Subroutines,
%    Volume II,
%    McGraw-Hill, 1980,
%    ISBN: 0-07-054202-3,
%    LC: QA76.95.R82.
%
%    Pitts Jarvis,
%    Implementing CORDIC Algorithms,
%    Dr Dobb's Journal,
%    October 1990.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Input, integer N, the number of steps to take.
%
%    Output, real FX, the exponential of X.
%
%  Local Parameters:
%
%    Local, real A(1:25) = exp ( (1/2)^(1:25) );
%
  a_length = 25;

  a = [ ...
    1.648721270700128, ...
    1.284025416687742, ...
    1.133148453066826, ...
    1.064494458917859, ...
    1.031743407499103, ...
    1.015747708586686, ...
    1.007843097206488, ...
    1.003913889338348, ...
    1.001955033591003, ...
    1.000977039492417, ...
    1.000488400478694, ...
    1.000244170429748, ...
    1.000122077763384, ...
    1.000061037018933, ...
    1.000030518043791, ...
    1.0000152589054785, ...
    1.0000076294236351, ...
    1.0000038147045416, ...
    1.0000019073504518, ...
    1.0000009536747712, ...
    1.0000004768372719, ...
    1.0000002384186075, ...
    1.0000001192092967, ...
    1.0000000596046466, ...
    1.0000000298023228 ];
  e = 2.718281828459045;

  x_int = floor ( x );
%
%  Determine the weights.
%
  poweroftwo = 0.5;
  z = x - x_int;

  for i = 1 : n
    w(i) = 0.0;
    if ( poweroftwo < z )
      w(i) = 1.0;
      z = z - poweroftwo;
    end
    poweroftwo = poweroftwo / 2.0;
  end
%
%  Calculate products.
%
  fx = 1.0;

  for i = 1 : n

    if ( i <= a_length )
      ai = a(i);
    else
      ai = 1.0 + ( ai - 1.0 ) / 2.0;
    end

    if ( 0.0 < w(i) )
      fx = fx * ai;
    end

  end
%
%  Perform residual multiplication.
%
  fx = fx             ...
    * ( 1.0 + z       ...
    * ( 1.0 + z / 2.0 ...
    * ( 1.0 + z / 3.0 ...
    * ( 1.0 + z / 4.0 ))));
%
%  Account for factor EXP(X_INT).
%
  if ( x_int < 0 )

    for i = 1 : -x_int
      fx = fx / e;
    end

  else

    for i = 1 : x_int
      fx = fx * e;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function haar_test ( )

%*****************************************************************************80
%
%% HAAR_TEST tests the HAAR library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HAAR_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HAAR library.\n' );

  haar_test01 ( );
  haar_test02 ( );
  haar_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HAAR_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function varargout = parsef(input, format)
%parsef   parse string value using FORTRAN formatting codes
%   [val1,val2,...valn] = parsef(input, format)
%   input is string input value
%   format is cell array of format codes

global input_
input_ = input;
varargout = getvals(1, format, 1);
clear global input_
return

% this function does the work. you probably don't want to go here
function [output, idx] = getvals(idx, format, reps)
global input_
count = 1;
output = {};
for k = 1:reps
	odx = 1;
	for i = 1:length(format)
		fmt = format{i};
		switch class(fmt)
		case 'double'
			count = fmt;
		case 'char'
			type = fmt(1);
			if type == 'X'
				idx = idx+count;
			else
				[len,cnt] = sscanf(fmt,'%*c%d',1);
				if cnt ~= 1
					error(['Invalid format specifier: ''',fmt,'''']);
				end
				switch type
				case {'I','i'}
					for j = 1:count
						[val,cnt] = sscanf(input_(idx:min(idx+len-1,end)),'%d',1);
						if cnt == 1
							output{odx}(j,k) = val;
						else
							output{odx}(j,k) = NaN;
						end
						idx = idx+len;
					end
				case {'F','f'}
					for j = 1:count
						[val,cnt] = sscanf(input_(idx:min(idx+len-1,end)),'%f',1);
						if cnt == 1
							output{odx}(j,k) = val;
						else
							output{odx}(j,k) = NaN;
						end
						idx = idx+len;
					end
				case {'E','D','G'}
					for j = 1:count
						[val,cnt] = sscanf(input_(idx:min(idx+len-1,end)),'%f%*1[DdEe]%f',2);
						if cnt == 2
							output{odx}(j,k) = val(1) * 10^val(2); %#ok<AGROW>
						elseif cnt == 1
							output{odx}(j,k) = val;
						else
							output{odx}(j,k) = NaN;
						end
						idx = idx+len;
					end
				case 'A'
					for j = 1:count
						output{odx}{j,k} = input_(idx:min(idx+len-1,end));
						idx = idx+len;
					end
				otherwise
					error(['Invalid format specifier: ''',fmt,'''']);
				end
				odx = odx+1;
			end
			count = 1;
		case 'cell'
 			[val, idx] = getvals(idx, fmt, count);
			if length(val) == 1
				output(odx) = val;
			else
				output{odx} = val;
			end
			odx = odx+1;
			count = 1;
		end
	end
end
return

author: yandld


text:
%This Matlab script can be used to reproduce Figure 6.4 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Range of number of UEs per BS
Krange = [10 20];

%Compute maximum number of UEs
Kmax = max(Krange);

%Number of BS antennas
M = 100;

%Define the range of pilot reuse factors
fRange = [1 2 4];

%Select the number of setups with random UE locations
nbrOfSetups = 1;

%Select the number of channel realizations per setup
nbrOfRealizations = 200;

%Select the range of hardware qualities. The values at the same position in
%the different vectors are considered simultaneously.
kappatUE = [0.9 0.92 0.94 0.96 0.97 0.98 0.99 0.995 1];
kapparBS = kappatUE;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Define total uplink transmit power per UE (mW)
p = 100;

%Define noise figure at BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Prepare to save simulation results
sumSE_MR = zeros(length(kappatUE),length(fRange),nbrOfSetups,length(Krange));
sumSE_RZF = zeros(length(kappatUE),length(fRange),nbrOfSetups,length(Krange));
sumSE_MMMSE = zeros(length(kappatUE),length(fRange),nbrOfSetups,length(Krange));

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,Kmax,M,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Go through all pilot reuse factors
    for s = 1:length(fRange)

        %Output simulation progress
        disp([num2str(s) ' reuse factors out of ' num2str(length(fRange))]);

        %Extract pilot reuse factor
        f = fRange(s);

        %Go through all hardware qualities
        for r = 1:length(kappatUE)

            for kind = 1:length(Krange)

                %Extract number of UEs
                K = Krange(kind);

                %Generate channel realizations with estimates and estimation
                %error correlation matrices, using LMMSE estimation
                [Hhat,C,tau_p,~,H] = functionChannelEstimates_impairments(R(:,:,1:K,:,:),channelGainOverNoise(1:K,:,:),nbrOfRealizations,M,K,L,p,f,kappatUE(r),kapparBS(r));

                %Compute SEs using Monte Carlo realizations
                [SE_MR,SE_RZF,SE_MMMSE] = functionComputeSE_UL_impairments(H,Hhat,C,tau_c,tau_p,nbrOfRealizations,M,K,L,p,kappatUE(r),kapparBS(r));

                %Save results
                sumSE_MR(r,s,n,kind) = mean(sum(SE_MR,1));
                sumSE_RZF(r,s,n,kind) = mean(sum(SE_RZF,1));
                sumSE_MMMSE(r,s,n,kind) = mean(sum(SE_MMMSE,1));

                %Delete large matrices
                clear Hhat C H;

            end

        end

    end

    %Delete large matrices
    clear R;

end

%% Plot the simulation results
for kind = 1:length(Krange)

    figure;
    hold on; box on;

    plot(kappatUE,max(mean(sumSE_MMMSE(:,:,:,kind),3),[],2),'rd-','LineWidth',1);
    plot(kappatUE,max(mean(sumSE_RZF(:,:,:,kind),3),[],2),'k-.','LineWidth',1);
    plot(kappatUE,max(mean(sumSE_MR(:,:,:,kind),3),[],2),'bs-','LineWidth',1);

    xlabel('Hardware quality');
    ylabel('Average sum SE [bit/s/Hz/cell]');
    legend('M-MMSE','RZF','MR','Location','NorthWest');
    ylim([0 ceil(max(max(mean(sumSE_MMMSE(:,:,:,kind),3),[],2))/10)*10]);

end

author: emilbjornson

repo: massivemimobook


text:
function gf = gradddivHexp(Hkt_r, varargin)
% function gf = gradddivH(Hkt_r, varargin)
% gf is 1 x k*t row vector
% Hkt_r = reshape(Hkt,1,k*t) -> row vector
% Vxt = varargin{1};  %data
% Wxk = varargin{2};  %W matrix
% Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
% Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)

alphaH=1; %for now....
Vxt = varargin{1};  %data
Wxk_tmp = varargin{2};  %W matrix
Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)
peval = varargin{5}; %parameters

if ~isfield(peval, 'w_lambda') peval.w_lambda=0; end

t=size(Vxt,2);
k=length(peval.h_dovec);

Hkt_tmp = exp(reshape(Hkt_r,k,t));

Wxk = zeros(peval.numpix, peval.ncomp+1);
Hkt = zeros(peval.ncomp+1, peval.nt);

Wxk(:,peval.w_dovec)=Wxk_tmp;
Hkt(peval.h_dovec,:)=Hkt_tmp;

Wxk(:,peval.w_fixvec)=Wxk_fix;
Hkt(peval.h_fixvec,:)=Hkt_fix;

deltasum=sum(sum(Wxk_tmp,1))-k;
if and(~isempty(Wxk_tmp), abs(deltasum)>10^-6)
    error('Wxk is not correctly normalized! (sum(Wxk_tmp,1)<>1)\n sum(Wxk_tmp,1)=%f',deltasum)
end

% gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt)))*alphaH.*Hkt; %d/dh(d-divergence)
% fprintf('!!!test in gradddivHexp!!!')
denomin=(sum((sum(Hkt,2)).^2))^2;
sumt_h = repmat(sum(Hkt,2),1, peval.nt);
gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt)))*alphaH.*Hkt - 2*sumt_h/denomin; %d/dh(d-divergence)
% one is tehre because Wxt is normalized: sum(Wxt,1)=1
gf=reshape(gfkt(1:k,:),1,k*t); %making row vector

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Crystal Directions
%
%%
%

author: mtex-toolbox

repo: mtex

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


text:
classdef TestRotatedRect
    %TestRotatedRect

    methods (Static)
        function test_RotatedRect_from3points
            pts = {[10,10], [10,20], [20,20]};
            rrect = cv.RotatedRect.from3points(pts{:});
            validateattributes(rrect, {'struct'}, {'scalar'});
            assert(all(ismember({'center','size','angle'}, fieldnames(rrect))));
        end

        function test_RotatedRect_points
            rrect = struct('center',[10,20], 'size',[5,6], 'angle',30);
            pts = cv.RotatedRect.points(rrect);
            validateattributes(pts, {'numeric'}, {'size',[4 2]});
        end

        function test_RotatedRect_boundingRect
            rrect = struct('center',[10,20], 'size',[5,6], 'angle',30);

            rect = cv.RotatedRect.boundingRect(rrect);
            validateattributes(rect, {'numeric'}, {'vector', 'numel',4, 'integer'});

            rect = cv.RotatedRect.boundingRect2f(rrect);
            validateattributes(rect, {'numeric'}, {'vector', 'numel',4, 'real'});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [OLS_Bhat, OLS_betahat, OLS_sigmahat, OLS_forecast_estimates, biclag]=arbicloop(data_endo,data_endo_a,const,p,n, m, Fperiods, Fband)

% function [OLS_Bhat{ii}, OLS_betahat{ii}, OLS_sigmahat{ii}, OLS_forecast_estimates]=arbicloop(data_endo,data_endo_a,const,p,n)
% computes individual OLS estimations of autoregressive models and record their residual variances, as stated p16 of technical guide
% inputs:  - matrix 'data_endo': matrix of endogenous data used for model estimation
%          - integer 'const': 0-1 value to determine if a constant is included in the model
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'Fperiods': number of forecast periods
%          - scalar 'Fband': confidence level for forecasts
% outputs: - cell 'OLS_Bhat' :
%          - cell 'OLS_betahat' :
%          - cell 'OLS_sigmahat': residual variance of individual AR models estimated for each endogenous variable
%          - cell 'OLS_forecast_estimates :

bic=[];
nvar=1;
%loop over lags(p)
for ii=1:p
  [sigmahatlag]=bear.arloop(data_endo,const,ii,n);
   bic=[bic -2*log(sigmahatlag)+(log(length(data_endo))*nvar)];
end

%optimising VAR for AR - for each one.
[bicmin, biclag]=min(bic', [], 1);

% Estimating and forecasting for BIC optimised model
% loop over columns of data_endo
for ii=1:n
% Estimating BIC optimised VAR and record parameters and residual variance
[OLS_Bhat{ii}, OLS_betahat{ii}, OLS_sigmahat{ii},~,~,~,~,~,~,~,~,~,~,~,~]=bear.olsvar(data_endo(:,ii),[],const,biclag(ii));
% Forecasting for BIC optimised model
[OLS_forecast_estimates{ii}]=bear.olsforecast(data_endo_a(:,ii),[],Fperiods,OLS_betahat{ii},OLS_Bhat{ii},OLS_sigmahat{ii},1,m,biclag(ii),biclag(ii)+1,const,Fband);
end

author: european-central-bank

repo: BEAR-toolbox


text:
if 0
    %% Set RTLD_NOW and RTLD_DEEPBIND
    %https://nl.mathworks.com/matlabcentral/answers/327193-calling-python-module-from-matlab-causes-segmentation-fault-in-h5py
    py.sys.setdlopenflags(int32(10));
else
    %PythonEnvironment objects contain information about the settings and status of the Python® interpreter.
    %MATLAB® communicates with the interpreter when you call a py. command.
    %Python environment information is persistent across different MATLAB sessions.

    % Out-of-Process Execution of Python Functionality
    % MATLAB® can run Python scripts and functions in a separate process. Running Python in a separate process enables you to:
    % Use some third-party libraries in the Python code that are not compatible with MATLAB.
    % Isolate the MATLAB process from crashes in the Python code.
    PythonEnvironment = pyenv("ExecutionMode","OutOfProcess");
end

try
  py.list('x','y',1)
catch e
  e.message
  if(isa(e,'matlab.exception.PyException'))
    e.ExceptionObject
  end
end

try
  py.help('pandas.read_csv')
catch e
  e.message
  if(isa(e,'matlab.exception.PyException'))
    e.ExceptionObject
  end
end

%pandas.read_csv function signature
%pandas.read_csv = read_csv(filepath_or_buffer: Union[str, pathlib.Path, IO[~AnyStr]], sep=',', delimiter=None, header='infer', names=None, index_col=None, usecols=None, squeeze=False, prefix=None, mangle_dupe_cols=True, dtype=None, engine=None, converters=None, true_values=None, false_values=None, skipinitialspace=False, skiprows=None, skipfooter=0, nrows=None, na_values=None, keep_default_na=True, na_filter=True, verbose=False, skip_blank_lines=True, parse_dates=False, infer_datetime_format=False, keep_date_col=False, date_parser=None, dayfirst=False, cache_dates=True, iterator=False, chunksize=None, compression='infer', thousands=None, decimal: str = '.', lineterminator=None, quotechar='"', quoting=0, doublequote=True, escapechar=None, comment=None, encoding=None, dialect=None, error_bad_lines=True, warn_bad_lines=True, delim_whitespace=False, low_memory=True, memory_map=False, float_precision=None)

%pyargs https://nl.mathworks.com/help/matlab/ref/pyargs.html
% Python command:
% db = pd.read_csv('./data/cache_compounds_20160818.csv',index_col='compound_id')
% MATLAB command using pyargs
db = py.pandas.read_csv(pyargs('filepath_or_buffer',[pyLibraryFolder filesep 'data' filesep 'test_compounds.csv'],'index_col','compound_id'));

order = py.dict(pyargs('soup',3.57,'bread',2.29,'bacon',3.91,'salad',5.00))
% A dictionary has pairs of keys and values. Display the menu items in the variable order using the Python keys function.
keys(order)
values(order)

if count(py.sys.path,'') == 0
    insert(py.sys.path,int32(0),'');
end

mymod = py.importlib.import_module('mymod')

N = py.list({'Jones','Johnson','James'})

names = py.mymod.search(N)

author: opencobra

repo: cobratoolbox


text:
function i1mach_test ( )

%*****************************************************************************80
%
%% I1MACH_TEST reports the constants returned by I1MACH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I1MACH_TEST\n' );
  fprintf ( 1, '  I1MACH reports the value of constants associated\n' );
  fprintf ( 1, '  with integer computer arithmetic.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Numbers associated with input/output units:\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(1) = the standard input unit.\n' );
  fprintf ( 1, '%d\n', i1mach(1) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(2) = the standard output unit.\n' );
  fprintf ( 1, '%d\n', i1mach(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(3) = the standard punch unit.\n' );
  fprintf ( 1, '%d\n', i1mach(3) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(4) = the standard error message unit.\n' );
  fprintf ( 1, '%d\n', i1mach(4) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Numbers associated with words:\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(5) = the number of bits per integer.\n' );
  fprintf ( 1, '%d\n', i1mach(5) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(6) = the number of characters per integer.\n' );
  fprintf ( 1, '%d\n', i1mach(6) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Numbers associated with integer values:\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Assume integers are represented in the S digit \n' );
  fprintf ( 1, '  base A form:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Sign * (X(S-1)*A^(S-1) + ... + X(1)*A + X(0))\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where the digits X satisfy 0 <= X(1:S-1) < A.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(7) = A, the base.\n' );
  fprintf ( 1, '%d\n', i1mach(7) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(8) = S, the number of base A digits.\n' );
  fprintf ( 1, '%d\n', i1mach(8) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(9) = A^S-1, the largest integer.\n' );
  fprintf ( 1, '%d\n', i1mach(9) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Numbers associated with floating point values:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Assume floating point numbers are represented \n' );
  fprintf ( 1, '  in the T digit base B form:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Sign * (B^E) * ((X(1)/B) + ... + (X(T)/B^T) )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where \n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    0 <= X(1:T) < B,\n' );
  fprintf ( 1, '    0 < X(1) (unless the value being represented is 0),\n' );
  fprintf ( 1, '    EMIN <= E <= EMAX.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(10) = B, the base.\n' );
  fprintf ( 1, '%d\n', i1mach(10) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Numbers associated with single precision values:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(11) = T, the number of base B digits.\n' );
  fprintf ( 1, '%d\n', i1mach(11) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(12) = EMIN, the smallest exponent E.\n' );
  fprintf ( 1, '%d\n', i1mach(12) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(13) = EMAX, the largest exponent E.\n' );
  fprintf ( 1, '%d\n', i1mach(13) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Numbers associated with double precision values:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(14) = T, the number of base B digits.\n' );
  fprintf ( 1, '%d\n', i1mach(14) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(15) = EMIN, the smallest exponent E.\n' );
  fprintf ( 1, '%d\n', i1mach(15) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  I1MACH(16) = EMAX, the largest exponent E.\n' );
  fprintf ( 1, '%d\n', i1mach(16) );

  return
end


text:
function [E, T_noise_squared, d] = error_foo(Theta, X, sigma, P_inlier, parameters)

% [E T_noise_squared d] = error_foo(Theta, X, sigma, P_inlier, parameters)
%
% DESC:
% Template to estimate the error due to the foo constraint. To
% return only the error threshold the function call should be:
%
% [dummy T_noise d] = error_foo([], [], sigma, P_inlier, parameters);
%
% INPUT:
% Theta             = foo parameter vector
% X                 = samples on the manifold
% sigma             = noise std
% P_inlier          = Chi squared probability threshold for inliers
%                     If 0 then use directly sigma.
% parameters        = the parameters used by the functions
%
% OUTPUT:
% E                 = squared error
% T_noise_squared   = squared noise threshold
% d                 = degrees of freedom of the error distribution

% compute the error obtained by the orthogonal projection of
% the data points X onto the model manifold instantiated with the
% parameters Theta
E = [];
if ~isempty(Theta) && ~isempty(X)

    % error computation

end;

% compute the error threshold
if (nargout > 1)

    if (P_inlier == 0)
        % in this case the parameter sigma coincides with the noise
        % threshold
        T_noise_squared = sigma;
    else
        % otherwise we compute the error threshold given the standard
        % deviation of the noise assuming that the errors are normally
        % distributed. Hence the sum of their squares is Chi2
        % distributed with d degrees of freedom
        d = ;

        % compute the inverse probability
        T_noise_squared = sigma^2 * chi2inv_LUT(P_inlier, d);

    end;

end;

return;

author: RANSAC

repo: RANSAC-Toolbox


text:
function test_suite = test_grShortestPath
%TEST_GRSHORTESTPATH  Test case for the file grShortestPath
%
%   Test case for the file grShortestPath

%   Example
%   test_grShortestPath
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2014-01-13,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2014 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>
% propagate distance on a graph where a branch has smaller nodes but larger
% cumulated distance

% create graph structure
n0 = [10 10;20 10; 30 10;40 10;50 10;30 40];
e0 = [1 2;1 6; 2 3;3 4;4 5;5 6];
% weights of edges
l0 = [2 5 2 2 2 5]';
path = grShortestPath(n0, e0, 1, 5, l0);

expPath = [1 2 3 4 5]';
testCase.assertEqual(expPath, path, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [outputVar1, outputVar2] = cs_mri2scs(MRI, mriCoord)
% CS_MRI2SCS: Compute the transform to move from the MRI coordinate system (in mm) to the SCS
%
% USAGE:            [transfSCS] = cs_mri2scs(MRI);
%         [scsCoord, transfSCS] = cs_mri2scs(MRI, mriCoord);
%
% DEPRECATED: Replace with cs_convert()

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@

disp('BST> WARNING: Deprecated function cs_mri2scs(), use cs_convert() instead.');

if (nargin == 1)
    outputVar1 = cs_compute(MRI, 'scs');
    outputVar2 = [];
else
    outputVar1 = cs_convert(MRI, 'mri', 'scs', mriCoord' ./ 1000)' .* 1000;
    outputVar2 = MRI.SCS;
end

author: brainstorm-tools


text:
% The COBRAToolbox: testTranslateMetagenome2AGORA.m
%
% Purpose:
%     - tests the basic functionality of translateMetagenome2AGORA
%
% Authors:
%     - Almut Heinken - March 2019
%
% Ensure that the organisms in the input file were correctly translated.
% This is the case if all organisms in the output file correspond to an
% AGORA strain model or AGORA pan-model.
% Note that this test cannot tell if the biological identity of the
% translated organisms is correct as this requires manual inspection.
% Moreover, input files with abundance data from other source may fail due
% to differences in formatting and nomenclature.

if ispc && exist('AGORA_infoFile.xlsx','file')
    % Read in the info file with all AGORA strains and taxa
    [~, infoFile, ~] = xlsread('AGORA_infoFile.xlsx');

    % Download an example file with abundance data from EMBL-EBI
    system('curl -O https://www.ebi.ac.uk/metagenomics/api/v1/studies/MGYS00001248/pipelines/3.0/file/SRP065497_taxonomy_abundances_v3.0.tsv')

    % Test the output file on different taxonomical levels
    % Note: On phylum and class level, some samples sum up to zero due to too
    % few corresponding data points in the input file. Thus, these taxon levels
    % cannot be tested.
    taxLevels={
        'Species'
        'Genus'
        'Family'
        'Order'
        };

    for t=1:size(taxLevels,1)
        fprintf([' > Running translateMetagenome2AGORA with ' taxLevels{t,1} '...\n']);
        % Define the taxon level of the input file
        orgList=unique(infoFile(2:end,find(strcmp(taxLevels{t,1},infoFile(1,:)))));
        orgList(strncmp('unclassified', orgList, 12)) = [];
        [translatedAbundances,normalizedAbundances,unmappedRows]=translateMetagenome2AGORA('SRP065497_taxonomy_abundances_v3.0.tsv',taxLevels{t,1});

        % Verify that the output file is not empty
        assert(size(translatedAbundances,1)>1)

        fprintf(' > Output file is not empty.\n');

        if size(translatedAbundances,1)>1
            translatedOrgs=translatedAbundances(2:end,1);
            translatedOrgs=strrep(translatedOrgs,'pan','');
            translatedOrgs=strrep(translatedOrgs,'_',' ');

            % Find the overlap between translated organisms and AGORA organisms
            C = intersect(translatedOrgs,orgList);

            % Verify that all output organisms overlap with AGORA organisms
            assert(length(translatedOrgs) == length(C))
        end

        % Verify that the relative abundances for each sample sum up to 1
        for i=2:size(normalizedAbundances,2)
            assert(sum(str2double(normalizedAbundances(2:end,i)))-1 < 0.0001)
        end

        fprintf([' > Testing tax level ' taxLevels{t,1} ' done.\n']);
    end

    % output a success message
    fprintf('Done.\n');
else
    fprintf('Could not run testTranslateMetagenome2AGORA because AGORA_infoFile.xlsx not found or this is not a windows pc')
end


text:
function [ xmin, xmax ] = r8_gaml ( )

%*****************************************************************************80
%
%% R8_GAML evaluates bounds for an R8 argument of the gamma function.
%
%  Discussion:
%
%    This function calculates the minimum and maximum legal bounds
%    for X in the evaluation of GAMMA ( X ).
%
%    XMIN and XMAX are not the only bounds, but they are the only
%    non-trivial ones to calculate.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Output, real XMIN, XMAX, the bounds.
%
  alnsml = log ( r8_mach ( 1 ) );
  xmin = - alnsml;

  for i = 1 : 10

    xold = xmin;
    xln = log ( xmin );
    xmin = xmin - xmin * ( ( xmin + 0.5 ) * xln - xmin ...
      - 0.2258 + alnsml ) / ( xmin * xln + 0.5 );

    if ( abs ( xmin - xold ) < 0.005 )

      xmin = - xmin + 0.01;

      alnbig = log ( r8_mach ( 2 ) );
      xmax = alnbig;

      for j = 1 : 10

        xold = xmax;
        xln = log ( xmax );
        xmax = xmax - xmax * ( ( xmax - 0.5 ) * xln - xmax ...
          + 0.9189 - alnbig ) / ( xmax * xln - 0.5 );

        if ( abs ( xmax - xold ) < 0.005 )
          xmax = xmax - 0.01;
          xmin = max ( xmin, - xmax + 1.0 );
          return
        end

      end

      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8_GAML - Fatal error!\n' );
      fprintf ( 1, '  Unable to find XMAX.\n' );
      error ( 'R8_GAML - Fatal error!' )

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_GAML - Fatal error!\n' );
  fprintf ( 1, '  Unable to find XMIN.\n' );

  error ( 'R8_GAML - Fatal error!' )
end


text:
function legendre_polynomial_test10 ( p )

%*****************************************************************************80
%
%% LEGENDRE_POLYNOMIAL_TEST10 tests PN_PAIR_PRODUCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polynomial
%    factors.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_POLYNOMIAL_TEST10\n' );
  fprintf ( 1, '  Compute a pair product table for Pn(n,x):\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tij = integral ( -1 <= x <= +1 ) Pn(i,x) Pn(j,x) dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The Pn(n,x) polynomials are orthonormal,\n' );
  fprintf ( 1, '  so T should be the identity matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum degree P = %d\n', p );

  table = pn_pair_product ( p );

  r8mat_print ( p + 1, p + 1, table, '  Pair product table:' );

  return
end


text:
function [payload, interproc_data ] = aip_detector( ECG_matrix, ECG_header, ECG_start_offset, progress_handle, payload_in, interproc_data)
% aip_detector Arbitrary impulsive pseudoperiodic detector
%   Detailed explanation goes here
    payload = [];

    payload.series_quality.AnnNames = {};
    payload.series_quality.ratios = [];
    payload.series_quality.estimated_labs = {};

    %% config parsing

    if( ~isfield(payload_in, 'trgt_width') )
        % Default QRS width
        payload_in.trgt_width = 0.06; % seconds
    end

    if( ~isfield(payload_in, 'trgt_min_pattern_separation') )
        % Default QRS width
        payload_in.trgt_min_pattern_separation = 0.3; % seconds
    end

    if( ~isfield(payload_in, 'trgt_max_pattern_separation') )
        % Default QRS width
        payload_in.trgt_max_pattern_separation = 2; % seconds
    end

    if( ~isfield(payload_in, 'stable_RR_time_win') )
        % Default time window to consider an stable rhythm
        payload_in.stable_RR_time_win = 2; % seconds
    end

    if( ~isfield(payload_in, 'final_build_time_win') )
        % Default time window to build the final detection based on each
        % pattern found in the recording
        payload_in.final_build_time_win = 20; % seconds
    end

    if( ~isfield(payload_in, 'powerline_interference') )
        % Default auto
        payload_in.powerline_interference = nan; % Hz
    end

    if( ~isfield(payload_in, 'max_patterns_found') )
        % Default QRS width
        payload_in.max_patterns_found = 3; % # patterns
    end

    if( ~isfield(payload_in, 'sig_idx') )
        % Default QRS width
        payload_in.sig_idx = 1:ECG_header.nsig; % # patterns
    end

    %% preprocessing

    % lead names desambiguation
    str_aux = regexprep(cellstr(ECG_header.desc), '\W*(\w+)\W*', '$1');
    lead_names = regexprep(str_aux, '\W', '_');

    [str_aux2, ~ , aux_idx] = unique( lead_names );
    aux_val = length(str_aux2);

    if( aux_val ~= ECG_header.nsig )
        for ii = 1:aux_val
            bAux = aux_idx==ii;
            aux_matches = sum(bAux);
            if( sum(bAux) > 1 )
                lead_names(bAux) = strcat( lead_names(bAux), repmat({'v'}, aux_matches,1), cellstr(num2str((1:aux_matches)')) );
            end
        end
    end

    lead_names = regexprep(lead_names, '\W*(\w+)\W*', '$1');
    lead_names = regexprep(lead_names, '\W', '_');

    progress_handle.checkpoint('Powerline filtering');

    ECG_matrix = double(ECG_matrix);

    if( isnan(payload_in.powerline_interference) )
        f_armonic = 1:3;
        f_armonic = f_armonic*50;
        f_armonic = f_armonic( f_armonic < (ECG_header.freq / 2) );
        f50_idx = round( f_armonic * ECG_header.nsamp / ECG_header.freq ) + 1;
        pwl50_e = abs(goertzel( ECG_matrix(:,payload_in.sig_idx), f50_idx ));

        f_armonic = 1:3;
        f_armonic = f_armonic*60;
        f_armonic = f_armonic( f_armonic < (ECG_header.freq / 2) );
        f60_idx = round( f_armonic * ECG_header.nsamp /ECG_header.freq ) + 1;
        pwl60_e = abs(goertzel( ECG_matrix(:,payload_in.sig_idx), f60_idx ));

        if( sum(sum(pwl50_e)) > sum(sum(pwl60_e)) )
            payload_in.powerline_interference = 50; % Hz
        else
            payload_in.powerline_interference = 60; % Hz
        end
    end

    % Resample for the closest smaller sampling freq
    target_fs = floor(ECG_header.freq/payload_in.powerline_interference) * payload_in.powerline_interference;

    if( target_fs == ECG_header.freq )
        ECG_matrix_res = ECG_matrix(:,payload_in.sig_idx);
    else
        [ups, downs] = rat(target_fs / ECG_header.freq);
        ECG_matrix_res = resample(ECG_matrix(:,payload_in.sig_idx), ups, downs );
    end

    pwl_size = round(target_fs / payload_in.powerline_interference);

    ECG_matrix_res = filter( ones(pwl_size,1)/pwl_size, 1, flipud(ECG_matrix_res) );
    ECG_matrix_res = filter( ones(pwl_size,1)/pwl_size, 1, flipud(ECG_matrix_res) );

    if( target_fs ~= ECG_header.freq )
        ECG_matrix_res = resample(ECG_matrix_res, downs, ups );
        aux_val = size(ECG_matrix_res,1);

        % Pad or trim after resampling
        if( aux_val < ECG_header.nsamp )
            ECG_matrix_res = [ECG_matrix_res; zeros( ECG_header.nsamp - aux_val, size(ECG_matrix_res,1) )];
        elseif( aux_val > ECG_header.nsamp )
            ECG_matrix_res = ECG_matrix_res(1:ECG_header.nsamp,:);
        end
    end

    ECG_matrix(:,payload_in.sig_idx) = ECG_matrix_res;

    %% Start
    pattern_size = 2*round(payload_in.trgt_width/2*ECG_header.freq)+1; % force odd number

    pattern_coeffs = diff(gausswin(pattern_size+1)) .* gausswin(pattern_size);
    first_pattern_coeffs = pattern_coeffs;

    progress_handle.Loops2Do = length(payload_in.sig_idx);

    lp_size = round(1.2*pattern_size);

    for this_sig_idx = rowvec(payload_in.sig_idx)

        progress_handle.start_loop();

        pb_str_prefix = [ lead_names{this_sig_idx} ': ' ];
        %% Initial Pattern search

        progress_handle.checkpoint([ pb_str_prefix 'First pattern guess' ]);

        rise_detector = filter( first_pattern_coeffs, 1, flipud(ECG_matrix(:,this_sig_idx)) );
        rise_detector = filter( first_pattern_coeffs, 1, flipud(rise_detector) );
        rise_detector = filter( ones(lp_size,1)/lp_size, 1, flipud(abs(rise_detector)) );
        rise_detector = filter( ones(lp_size,1)/lp_size, 1, flipud(rise_detector) );

        % figure(4); plot(rise_detector ./ max(rise_detector)); hold on; plot(ECG_matrix(:,this_sig_idx) ./ max(ECG_matrix(:,this_sig_idx))); legend('rise detector','ECG');

        % figure(3); plot(ECG_matrix(:,1) ); xlims = xlim(); ylims = ylim(); box_h = ylims + [0.01 -0.1] * diff(ylims); aux_sc = 0.8*diff(ylims)/(max(rise_detector)-min(rise_detector));aux_off = ylims(1) + 0.6*diff(ylims); hold on; plot(rise_detector * aux_sc + aux_off); hold off; ylim(ylims);

        progress_handle.checkpoint([ pb_str_prefix 'Peak detection first guess' ]);

        actual_thr = thr_calc(rise_detector, payload_in.trgt_min_pattern_separation, ECG_header.freq);

        if isnan(actual_thr)
            cprintf('Red', '\nRecording %s, lead %s, no patterns above threshold found.(Current %f seconds). \n\n',ECG_header.recname , lead_names{this_sig_idx}, payload_in.stable_RR_time_win);
            continue
        end

        first_detection_idx = modmax(rise_detector, 1, actual_thr, 1, round(payload_in.trgt_min_pattern_separation * ECG_header.freq));

        %% Look for stable segments
        if numel(first_detection_idx) <= 1
            cprintf('Red', '\nRecording %s, lead %s, Not enough patterns found over noise floor. \n\n',ECG_header.recname , lead_names{this_sig_idx});
            continue
        end

        RRserie = RR_calculation(first_detection_idx, ECG_header.freq);

        RRserie_filt = MedianFiltSequence(first_detection_idx, RRserie, round(5 * ECG_header.freq));

        RR_scatter = abs(RRserie - RRserie_filt)./RRserie_filt;

        % Find the most stable regions
        RR_thr = prctile(RR_scatter, 50);

        str_aux = [ 'aip_guess' '_' lead_names{this_sig_idx} ];
        payload.(str_aux).time = first_detection_idx + ECG_start_offset - 1;
        payload.series_quality.AnnNames = [ payload.series_quality.AnnNames ; {str_aux} {'time'} ];
        RRserie_mean_sq_error = mean((RRserie - RRserie_filt).^2);
        payload.series_quality.ratios = [payload.series_quality.ratios; RRserie_mean_sq_error];
        payload.series_quality.estimated_labs = [ payload.series_quality.estimated_labs; {[]} ];

%         stable_rhythm_regions = get_segments_from_sequence(first_detection_idx, RR_scatter, round(payload_in.stable_RR_time_win * ECG_header.freq), RR_thr );
%
% %     figure(3); plot(first_detection_idx, [RRserie RRserie_filt] ); ylims = ylim(); ylims = ylims + [0.1 -0.1] * diff(ylims); hold on; plot( [ stable_rhythm_regions'; flipud(stable_rhythm_regions'); stable_rhythm_regions(:,1)' ],  repmat([ylims(1); ylims(1); ylims(2);ylims(2);ylims(1)],1, size(stable_rhythm_regions,1) ), 'b:' ); hold off;
% %     figure(3); plot(first_detection_idx, [RR_scatter], ':xb' ); xlims = xlim(); ylims = ylim(); ylims = ylims + [0.01 -0.1] * diff(ylims); hold on; plot(xlims, [RR_thr RR_thr], '--r'); plot( [ stable_rhythm_regions'; flipud(stable_rhythm_regions'); stable_rhythm_regions(:,1)' ],  repmat([ylims(1); ylims(1); ylims(2);ylims(2);ylims(1)],1, size(stable_rhythm_regions,1) ), 'b:' ); hold off;
% %     figure(3); plot(first_detection_idx, [RRserie RRserie_filt] ); xlims = xlim(); ylims = ylim(); ylims = ylims + [0.1 -0.1] * diff(ylims); aux_sc = 0.25*diff(ylims)/(max(RR_scatter)-min(RR_scatter));aux_off = ylims(1) + 0.65*diff(ylims); hold on; plot(first_detection_idx, (RR_scatter*aux_sc) +aux_off , ':xb' ); plot(xlims, ([RR_thr RR_thr]*aux_sc ) + aux_off, '--r'); plot( [ stable_rhythm_regions'; flipud(stable_rhythm_regions'); stable_rhythm_regions(:,1)' ],  repmat([ylims(1); ylims(1); ylims(2);ylims(2);ylims(1)],1, size(stable_rhythm_regions,1) ), 'b:' ); hold off;
%
%         if( isempty(stable_rhythm_regions) )
%             cprintf('Red', '\nRecording %s, lead %s, no stable rhythm segments found. Consider decreasing "stable_RR_time_win" (Current %f seconds). \n\n',ECG_header.recname , lead_names{this_sig_idx}, payload_in.stable_RR_time_win);
%             continue
%         end
%
%         lstable_rhythm_regions = size(stable_rhythm_regions,1);
%
%         % Get the larger or more stable segments first
%         [~, aux_idx] = sort(diff(stable_rhythm_regions,1,2), 'descend');
%
%         larger_stable_rank = arrayfun(@(a)(find(a==aux_idx)), 1:lstable_rhythm_regions);
%
%         pattern_prewin = round(payload_in.trgt_width/2*ECG_header.freq); % force odd number
%
%         % Get the segments with less change in morphology
%         pack_variance = repmat(realmax,lstable_rhythm_regions,1);
%
%         % From the longest - most rhythm stable regions
%         for ii = 1: round(lstable_rhythm_regions/4)
%
%             jj = aux_idx(ii);
%
%             % Woody method refinement
%             aux_idx2 = find(first_detection_idx >= stable_rhythm_regions(jj,1) & first_detection_idx <= stable_rhythm_regions(jj,2));
%             avg_pack = pack_signal(ECG_matrix(:,1), first_detection_idx(aux_idx2), [ pattern_prewin pattern_prewin ], true);
%
%             % Mean variance across the ensemble
%             pack_variance(jj) = nanmean(nanvar(squeeze(avg_pack), [], 2));
%
%         end
%         [~, aux_idx] = sort(pack_variance);
%
%         morphology_rank = arrayfun(@(a)(find(a==aux_idx)), 1:lstable_rhythm_regions);
%
%         % Final rank
%         [~, aux_idx] = sort(larger_stable_rank + morphology_rank);
%
%         stable_rhythm_regions = stable_rhythm_regions(aux_idx( 1:min(payload_in.max_patterns_found , size(stable_rhythm_regions,1))), : );
%
%         % figure(3); plot(first_detection_idx, [RRserie RRserie_filt] ); xlims = xlim(); ylims = ylim(); ylims = ylims + [0.1 -0.1] * diff(ylims); aux_sc = 0.25*diff(ylims)/(max(RR_scatter)-min(RR_scatter));aux_off = ylims(1) + 0.65*diff(ylims); hold on; plot(first_detection_idx, (RR_scatter*aux_sc) +aux_off , ':xb' ); plot(xlims, ([RR_thr RR_thr]*aux_sc ) + aux_off, '--r'); plot( [ stable_rhythm_regions'; flipud(stable_rhythm_regions'); stable_rhythm_regions(:,1)' ],  repmat([ylims(1); ylims(1); ylims(2);ylims(2);ylims(1)],1, size(stable_rhythm_regions,1) ), 'b:' ); arrayfun(@(a)(text( stable_rhythm_regions(a,1), ylims(2), num2str(a))), 1:min(payload_in.max_patterns_found )); hold off;
%
%         %% save first detections as a separate time serie
%
%         str_aux = [ 'aip_guess' '_' lead_names{this_sig_idx} ];
%         payload.(str_aux).time = first_detection_idx + ECG_start_offset - 1;
%         payload.series_quality.AnnNames = [ payload.series_quality.AnnNames ; {str_aux} {'time'} ];
%         RRserie_mean_sq_error = mean((RRserie - RRserie_filt).^2);
%         payload.series_quality.ratios = [payload.series_quality.ratios; RRserie_mean_sq_error];
%         payload.series_quality.estimated_labs = [ payload.series_quality.estimated_labs; {[]} ];
%
%         k_woody = 3;
%         win_limits = [ pattern_prewin (pattern_prewin+1) ];
%         extended_win_limits = [k_woody * win_limits(1), k_woody * win_limits(1)+1];
%         aux_offset = (extended_win_limits(1) - win_limits(1) + 1); % ???
%
%         aux_pat_coeff = zeros(sum(win_limits),size(stable_rhythm_regions,1));
%
%         for ii = 1:size(stable_rhythm_regions,1)
%
%             % refinamos con el metodo de woody
%             aux_idx2 = find(first_detection_idx >= stable_rhythm_regions(ii,1) & first_detection_idx <= stable_rhythm_regions(ii,2));
%
%     % figure(3); avg_pack = pack_signal(ECG_matrix(:,this_sig_idx), first_detection_idx(aux_idx2), 10*[ pattern_prewin pattern_prewin ], true); plot_ecg_mosaic(avg_pack)
%
%             if isempty(aux_idx2)
%                 % Cannot found any stable regions to apply the woody_method
%                 cprintf('Red', '\nRecording %s, lead %s, Cannot found any stable regions. \n\n',ECG_header.recname , lead_names{this_sig_idx});
%                 continue
%             end
%
%             pattern_coeffs = woody_method(ECG_matrix(:,this_sig_idx), first_detection_idx(aux_idx2), extended_win_limits, 0.95);
%             % trim the pattern
%             aux_pat_coeff(:,ii) = pattern_coeffs( aux_offset:(aux_offset+size(aux_pat_coeff,1)-1) );
%
%             % figure(5); plot(first_pattern_coeffs ./ max(abs(first_pattern_coeffs))); hold on; plot(aux_pat_coeff(:,ii) ./ max(abs(aux_pat_coeff(:,ii))))
%
%             progress_handle.checkpoint([ pb_str_prefix 'Pattern match ' num2str(ii) ]);
%
%         %         rise_detector = filter(ones(lp_size,1)/lp_size,1, abs([0; diff(double(ECG_matrix(:,ss)))]));
%         %         rise_detector = [zeros(((lp_size-1)/2)-1,1); rise_detector((lp_size-1)/2:end)];
%             rise_detector = filter( pattern_coeffs, 1, flipud(ECG_matrix(:,this_sig_idx)) );
%             rise_detector = filter( pattern_coeffs, 1, flipud(rise_detector) );
%             rise_detector = filter( ones(lp_size,1)/lp_size, 1, flipud(abs(rise_detector)) );
%             rise_detector = filter( ones(lp_size,1)/lp_size, 1, flipud(rise_detector) );
%
%             % figure(6); plot(rise_detector ./ max(rise_detector)); hold on; plot(ECG_matrix(:,this_sig_idx) ./ max(ECG_matrix(:,this_sig_idx)))
%
%             progress_handle.checkpoint([ pb_str_prefix 'Peak detection pattern ' num2str(ii) ]);
%
%             actual_thr = thr_calc(rise_detector, payload_in.trgt_min_pattern_separation, ECG_header.freq);
%
%             if isnan(actual_thr)
%                 cprintf('Red', '\nRecording %s, lead %s, no patterns above threshold found.(Current %f seconds). \n\n',ECG_header.recname , lead_names{this_sig_idx}, payload_in.stable_RR_time_win);
%                 continue
%             end
%
%             % Here the MIN time restriction is important
%             this_idx = modmax(rise_detector, 1, actual_thr, 1, round(payload_in.trgt_min_pattern_separation * ECG_header.freq));
%
%             if numel(first_detection_idx) <= 1
%                 cprintf('Red', '\nRecording %s, lead %s, Not enough patterns found over noise floor. \n\n',ECG_header.recname , lead_names{this_sig_idx});
%                 continue
%             end
%
%             RRserie = RR_calculation(this_idx, ECG_header.freq);
%
%             this_RRserie_filt = colvec(MedianFiltSequence(this_idx, RRserie, round(5 * ECG_header.freq)));
%
%             RRserie_mean_sq_error = mean((RRserie - this_RRserie_filt).^2);
%
%             % Structure of detection loading
%             str_aux = [ 'aip_patt_' num2str(ii) '_' lead_names{this_sig_idx} ];
%             payload.(str_aux).time = this_idx + ECG_start_offset - 1;
%             payload.series_quality.AnnNames = [ payload.series_quality.AnnNames ; {str_aux} {'time'} ];
%             payload.series_quality.ratios = [payload.series_quality.ratios; RRserie_mean_sq_error];
%             payload.series_quality.estimated_labs = [ payload.series_quality.estimated_labs; {[]} ];
%
%         end

        % Como el ratio se usa para rankear mediciones, hacemos un ratio
        % ficticio para cada detección por el ranking de "parsimoniosidad" que
        % sacó
        payload.series_quality.ratios = 1-(payload.series_quality.ratios * 1./max(payload.series_quality.ratios));

        progress_handle.end_loop();

    end
end

function [actual_thr] = thr_calc(rise_detector, trgt_min_pattern_separation, freq)
    initial_thr = 30;

    actual_thr = prctile(rise_detector, initial_thr);

    [~, max_values ] = modmax(rise_detector, 1, actual_thr, 1, round(trgt_min_pattern_separation * freq));

    if isempty(max_values)

        % possible flat recording.
        actual_thr = nan;

    else

        prctile_grid = prctile(max_values,1:100);

        grid_step = median(diff(prctile_grid));

        if grid_step == 0
            actual_thr = nan;
        else
            thr_grid = actual_thr:grid_step:max(max_values);

            hist_max_values = histcounts(max_values, thr_grid);

            first_bin_idx = 1;

            [thr_idx, thr_max ] = modmax( colvec( hist_max_values ) , first_bin_idx, 0, 0, [], [] );

            if( isempty(thr_max) )
                actual_thr = nan;
            else
                % Mass center of the distribution, probably the value where the
                % patterns under search are located.
                thr_idx_expected = floor(rowvec(thr_idx) * colvec(thr_max) *1/sum(thr_max));

                aux_seq = 1:length(thr_grid);

                % MPS0018 recording from Basel VIII database forced to add the
                % lesser or equal to aux_seq <= thr_idx_expected. Probably a
                % recording where no clear maximum above the noise threshold were
                % found.
                min_hist_max_values = min( hist_max_values( aux_seq >= first_bin_idx & aux_seq < thr_idx_expected) );

                % In case several indexes match the boolean condition, the mean
                % index is the center of all those indexes. Other criteria such as
                % min or max can be explored
                % thr_min_idx = round(mean(find(bAux & [hist_max_values 0] == min_hist_max_values)));
                thr_min_idx = round(mean(find(aux_seq >= first_bin_idx & aux_seq < thr_idx_expected & [hist_max_values 0] == min_hist_max_values)));

                actual_thr = thr_grid(thr_min_idx);
            end
        end

    end

end

author: marianux

repo: ecg-kit

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


text:
%
% Test the optimization algorithm that fits the parse
%  from one image of a character to another image "token" of that character
%
nsamp = 10; % total number of MCMC samples
nsamp_mcmc = 10; % number of samples to keep
auto_affine = true; % start optimization after affine fitting
ns = 1; % number of strokes

ps = defaultps;
load(ps.libname,'lib');
motor_program = generate_character(lib,ns);
M1 = motor_program();
M2 = motor_program();
while M1.ink_off_page || M2.ink_off_page
   motor_program = generate_character(lib,ns);
   M1 = motor_program();
   M2 = motor_program();
end
M1.I = M1.pimg > 0.5;
M2.I = M2.pimg > 0.5;

% visualize
figure(1);
clf
subplot(1,3,1);
plot_motor_to_image(M1.I,M1.motor_warped);
title(['image score ',num2str(scoreMP(M1,lib,'image',true,'type',false,'token',false),3)]);
xlabel('example 1');
subplot(1,3,2);
plot_motor_to_image(M2.I,M2.motor_warped);
title(['image score ',num2str(scoreMP(M2,lib,'image',true,'type',false,'token',false),3)]);
xlabel('example 2');
pause(0.1);
drawnow

% compute samples for type-level variables
all_samples = mcmc_all(M1,lib,nsamp_mcmc,'type');
indx = round(linspace(nsamp,nsamp_mcmc,nsamp));
subset_samples = all_samples(indx);

% run the optimization
Mfit2 = FitNewExemplar(M2.I,subset_samples,lib,auto_affine);

subplot(1,3,3);
plot_motor_to_image(Mfit2.I,Mfit2.motor_warped);
title(['image score ',num2str(scoreMP(Mfit2,lib,'image',true,'type',false,'token',false),3)]);
xlabel('ex 1 fit to ex 2');

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function varargout=readsph(varargin)
%
% For calling details please see v_readsph.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_readsph.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: readsph.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_readsph(varargin{:});
else
    v_readsph(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
 classdef GemText < GemUserInterfaceObject
    % GemText GEM class for a read-only, clickable text object
    %
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        FontSize
        Bold
        FontColour
        HorizontalAlignment
        SelectedColour
        SelectedFontColour
        Clickable
        BackgroundColour
    end

    properties (Access = private)
        Text
        Tag
        ToolTip
        Selected
        Highlighted
    end

    events
        TextClicked
        TextRightClicked
        TextShiftClicked
    end

    methods
        function obj = GemText(parent, text, tooltip, tag)
            obj = obj@GemUserInterfaceObject(parent);
            obj.Text = text;
            obj.Bold = false;
            obj.ToolTip = tooltip;
            obj.Tag = tag;
            obj.FontSize = 11;
            obj.HorizontalAlignment = 'left';
            obj.Selected = false;
            obj.Highlighted = false;
            obj.FontColour = obj.StyleSheet.TextPrimaryColour;
            obj.SelectedColour = obj.StyleSheet.SelectedBackgroundColour;
            obj.SelectedFontColour = obj.StyleSheet.TextContrastColour;
            obj.BackgroundColour = obj.StyleSheet.BackgroundColour;

            obj.Clickable = true;
        end

        function Select(obj, selected)
            if (selected ~= obj.Selected)
                obj.Selected = selected;
                obj.UpdateBackgroundColour;
            end
        end

        function Highlight(obj, highlighted)
            if (highlighted ~= obj.Highlighted)
                obj.Highlighted = highlighted;
                obj.UpdateBackgroundColour;
            end
        end

        function CreateGuiComponent(obj, position)
            text_size = [position(1), position(2), position(3), position(4)];
            if obj.Selected
                background_colour =  obj.SelectedColour;
                text_colour = obj.SelectedFontColour;
            else
                background_colour = obj.BackgroundColour;
                text_colour = obj.FontColour;
            end

            if obj.Bold
                weight = 'bold';
            else
                weight = 'normal';
            end

            obj.GraphicalComponentHandle = uicontrol('Style', 'text', 'Parent', obj.Parent.GetContainerHandle, 'Units', 'pixels', 'String', obj.Text, ...
                'Tag', obj.Tag, 'ToolTipString', obj.ToolTip, ...
                'BackgroundColor', background_colour, ...
                'FontAngle', 'normal', 'ForegroundColor', text_colour, 'FontName', obj.StyleSheet.Font, 'FontUnits', 'pixels', 'FontSize', obj.FontSize, ...
                'HorizontalAlignment', obj.HorizontalAlignment, 'FontWeight', weight, ...
                'Position', text_size, 'Enable', 'inactive');
        end

        function ChangeText(obj, new_text)
            if ~strcmp(obj.Text, new_text)
                obj.Text = new_text;
                if obj.ComponentHasBeenCreated && ishandle(obj.GraphicalComponentHandle)
                    set(obj.GraphicalComponentHandle, 'String', new_text);
                end
            end
        end
    end

    methods (Access = protected)

        function input_has_been_processed = MouseDown(obj, click_point, selection_type, src, eventdata)
            % This method is called when the mouse is clicked inside the control

            if ~obj.Clickable
                input_has_been_processed = false;
                return;
            end

            input_has_been_processed = true;

            if strcmp(selection_type, 'extend')
                notify(obj, 'TextShiftClicked', CoreEventData(obj.Tag));
            elseif strcmp(selection_type, 'alt')
                notify(obj, 'TextRightClicked', CoreEventData(obj.Tag));
            else
                notify(obj, 'TextClicked', CoreEventData(obj.Tag));
            end
        end

        function UpdateBackgroundColour(obj)
            if ~isempty(obj.GraphicalComponentHandle)
                if obj.Selected
                    background_colour = obj.SelectedColour;
                    text_colour = obj.SelectedFontColour;
                else
                    background_colour = obj.BackgroundColour;
                    text_colour = obj.FontColour;
                end

                if obj.Highlighted
                    background_colour = min(1, background_colour + 0.2);
                end
                set(obj.GraphicalComponentHandle, 'BackgroundColor', background_colour, 'ForegroundColor', text_colour);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit