text:
function [s, cfg] = ft_statfun_indepsamplesZcoh(cfg, dat, design)

% FT_STATFUN_INDEPSAMPLESCOHZ calculates the independent samples coherence
% Z-statistic on the biological data in dat (the dependent variable), using the
% information on the independent variable (ivar) in design.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option
%   cfg.statistic = 'ft_statfun_indepsamplesZcoh'
%
% The samples-dimension of the dat-variable must be the result of a reshaping
% operation applied to a data structure with dimord chan_(freq_time) or
% pos_(freq_time). The configuration must contain channel labels in cfg.label or
% position information in cfg.pos. This information is used to determine the number
% of channels. The dimord of the output fields is [prod(nchancmb,nfreq,ntime),1]. The
% channel combinations are the elements of the lower diagonal of the cross-spectral
% density matrix.
%
% Configuration options
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computeprob    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or
% cfg.computeprob='yes'.
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail  = -1, 0, or 1, left, two-sided, or right (default=1)
%               cfg.tail in combination with cfg.computecritval='yes'
%               determines whether the critical value is computed at
%               quantile cfg.alpha (with cfg.tail=-1), at quantiles
%               cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%               quantile (1-cfg.alpha) (with cfg.tail=1).
%
% Design specification
%   cfg.ivar  = column number of the design that contains the labels of the conditions that
%               must be compared (default=1). The labels are the numbers 1 and 2.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

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

% set the defaults
if ~isfield(cfg, 'computestat'),       cfg.computestat='yes';     end
if ~isfield(cfg, 'computecritval'),    cfg.computecritval='no';   end
if ~isfield(cfg, 'computeprob'),       cfg.computeprob='no';      end
if ~isfield(cfg, 'alpha'),             cfg.alpha=0.05;            end
if ~isfield(cfg, 'tail'),              cfg.tail=1;                end

% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
    ft_error('P-values can only be calculated if the test statistics are calculated.');
end
if isfield(cfg,'uvar') && ~isempty(cfg.uvar)
    ft_error('cfg.uvar should not exist for an independent samples statistic');
end
% if ~isfield(cfg, 'label') && ~isfield(cfg, 'pos')
%   ft_error('the configuration needs to contain either a label or a pos field');
% elseif isfield(cfg, 'label') && isfield(cfg, 'pos') && ~isempty(cfg.label) && ~isempty(cfg.pos)
%   ft_error('the configuration needs to contain either a non-empty label or a non-empty pos field');
% elseif isfield(cfg, 'label') && ~isempty(cfg.label)
%   nchan = length(cfg.label);
% elseif isfield(cfg, 'pos') && ~isempty(cfg.pos)
%   nchan = size(cfg.pos,1);
% end
nchan = cfg.dim(1);

% perform some checks on the design
selc1 = find(design(cfg.ivar,:)==1);
selc2 = find(design(cfg.ivar,:)==2);
nreplc1 = length(selc1);
nreplc2 = length(selc2);
nrepl = nreplc1 + nreplc2;
if nrepl<size(design,1)
  ft_error('Invalid specification of the independent variable in the design array.');
end
if nreplc1<2 || nreplc2<2
    ft_error('Every condition must contain at least two trials/tapers.');
end
dfc1 = nreplc1*2;
dfc2 = nreplc2*2;

if strcmp(cfg.computestat, 'yes')
  % compute the statistic
  nsamples = size(dat,1);
  %nchan = length(cfg.label); %this is computed earlier
  chancmbsel = find(tril(ones(nchan),-1));
  nfreqtim = nsamples/nchan;
  nchancmb = length(chancmbsel);
  nnewsamples = nchancmb*nfreqtim;
  s.stat = zeros(nnewsamples,1);
  for freqtimindx=1:nfreqtim
    chansel=((freqtimindx-1)*nchan + 1):(freqtimindx*nchan);
    csdc1=dat(chansel,selc1)*dat(chansel,selc1)'/nreplc1;
    powerc1=diag(csdc1);
    normmat=diag(powerc1.^(-1/2));
    csdc1=normmat*csdc1*normmat;
    csdc2=dat(chansel,selc2)*dat(chansel,selc2)'/nreplc2;
    powerc2=diag(csdc2);
    normmat=diag(powerc2.^(-1/2));
    csdc2=normmat*csdc2*normmat;
    biasc1=1./(dfc1-2);
    biasc2=1./(dfc2-2);
    denomZ=sqrt(1./(dfc1-2) + 1./(dfc2-2));
    tempstat=(atanh(abs(csdc1))-biasc1-atanh(abs(csdc2))+biasc2)./denomZ;
    s.stat(((freqtimindx-1)*nchancmb + 1):(freqtimindx*nchancmb))=tempstat(chancmbsel);
  end
  s.stat = reshape(s.stat, [nchancmb nfreqtim]);
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  if cfg.tail==-1
    s.critval = norminv(cfg.alpha);
  elseif  cfg.tail==0
    s.critval = [norminv(cfg.alpha/2),norminv(1-cfg.alpha/2)];
  elseif cfg.tail==1
    s.critval = norminv(1-cfg.alpha);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  if cfg.tail==-1
    s.prob = normcdf(s.stat);
  elseif  cfg.tail==0
    s.prob = 2*normcdf(-abs(s.stat));
  elseif cfg.tail==1
    s.prob = 1-normcdf(s.stat);
  end
  s.prob = reshape(s.prob, [nchancmb nfreqtim]);
end

% adjust the dimord
if strcmp(cfg.dimord(1:3), 'pos')
  cfg.dimord = ['poscmb',cfg.dimord(4:end)];
elseif strcmp(cfg.dimord(1:4), 'chan')
  cfg.dimord = ['chancmb',cfg.dimord(5:end)];
end

% append an indexing matrix to the cfg to be able to recover the channel
% combinations
chanindx = tril(true(nchan),-1);
cmbindx1 = repmat((1:nchan)', [1 nchan]);
cmbindx2 = repmat((1:nchan),  [nchan 1]);
cfg.chancmbindx(:,1) = cmbindx1(chanindx);
cfg.chancmbindx(:,2) = cmbindx2(chanindx);

author: spm


text:
% LoadGUI

% manually check that directories are correct in 'FishExplorer2\dir setup'

clear all;close all;clc
% path(pathdef);

%%
% scriptName = mfilename('fullpath');
% [currentpath, filename, fileextension]= fileparts(scriptName);
% code_dir = currentpath;

%% setup dir
file_masterdir = 'C:\Users\Xiu\Dropbox (Personal)';
if ~exist(file_masterdir,'dir')
    file_masterdir = 'C:\Users\xiuye\Dropbox';
end
code_dir = fullfile(file_masterdir,'FishExplorer2');

% SET CUSTOM PATH example:
% code_dir = 'C:\Users\xiuye\Dropbox\FishExplorer-master';

addpath(genpath(code_dir));

cd(code_dir)

%% setup
global VAR; %#ok<NUSED>
load_dir = GetOutputDataDir();
load(fullfile(load_dir,'VAR_new.mat'),'VAR'); % stores all clustering indices
VAR = CompileVARnames(VAR);

set(0, 'defaultUicontrolFontName', 'Arial');
set(0, 'defaultUitableFontName', 'Arial');
set(0, 'defaultUipanelFontName', 'Arial');
set(0, 'defaultAxesFontName', 'Arial');
set(0, 'defaultTextFontName', 'Arial');

%% Start GUI
% Run this, or press F5 in 'GUI_FishExplorer.m'
[hGUI,fcns] = GUI_FishExplorer;

% to push new functions to GUI, called with function 'push_cIX_gIX':
% hfig = GUI_FishExplorer(data_dir,name_CONSTs,name_MASKs,VAR,flag_script,var_script)

%%
% to get handle of local functions from GUI to use in workspace
% (call by using 'f.FunctionName')
f = [];
for i = 1:length(fcns),
    eval(['f.' char(fcns{i}) ' = fcns{i};']);
end

%%
% to recover after closing figure:
% global EXPORT_autorecover;

%% demo: load data without visualizaing in GUI (only using same infrastructure)
if false,
    %%
    h0 = figure;
    InitializeAppData(h0);
    i_fish = 3;
    [cIX,gIX,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,h0);

    figure('Position',[50,100,800,1000]);
    I = LoadCurrentFishForAnatPlot(h0);
    DrawCellsOnAnat(I);

    % filename = fullfile(saveDir, num2str(i_fish));
    % saveas(gcf, filename, 'png');
end

%%
% MASKs = getappdata(hGUI,'MASKs');
% Masks = MASKs.MaskDatabaseNames';

author: xiuyechen

repo: FishExplorer


text:
function cS = garnet(varargin)
% simple morphology
%

if nargin == 0 || ~isa(varargin{1},'symmetry')
  cs = crystalSymmetry('m3m','mineral','garnet');
else
  cs = varargin{1};
end

if check_option(varargin,'simple')

  N= Miller({1,0,0},{2,1,1},cs);
  dist = [0.45, 1];
  cS = crystalShape(N./dist);

  %N = Miller({1,1,0},{2,1,1},cs);
  %cS = crystalShape(N,1.5);

else

  N = Miller({1,0,0},{1,1,0},{4,3,1},{3,1,1},cs);
  dist = [0.92, 1.02, 3.93, 2.93];
  cS = crystalShape(N./dist);

end

author: mtex-toolbox

repo: mtex

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


text:
function [ y ] = FFT2D( x )
%FFT2D Summary of this function goes here
%   Detailed explanation goes here

y = fftshift(fft2(ifftshift(x)));

end

author: thomaskuestner


text:
function r = exprand(mu,m,n)
% EXPRAND Random matrices from exponential distribution.
%
%    R = EXPRAND(MU) returns a matrix of random numbers chosen
%    from the exponential distribution with parameter MU.
%    The size of R is the size of MU.
%    Alternatively, R = EXPRAND(MU,M,N) returns an M by N matrix.
%
% Copyright (c) 1998-2004 Aki Vehtari

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

error('No mex-file for this architecture. See Matlab help and convert.m in ./linuxCsource or ./winCsource for help.')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
classdef MOEADM2M_F3 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing MOEA/D-M2M

%------------------------------- Reference --------------------------------
% H. Liu, F. Gu, and Q. Zhang, Decomposition of a multiobjective
% optimization problem into a number of simple multiobjective subproblems,
% IEEE Transactions on Evolutionary Computation, 2014, 18(3): 450-455.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,2:end) - repmat(sin(pi/2*X(:,1)),1,size(X,2)-1);
            g = 10*sin(pi/2*X(:,1)).*sum(abs(t)./(1+exp(5*abs(t))),2);
            PopObj(:,1) = (1+g).*cos(pi/2*X(:,1));
            PopObj(:,2) = (1+g).*sin(pi/2*X(:,1));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = sqrt(1-R(:,1).^2);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% SP_INTEGRAL_PRESERVING_TRANSFORM: apply the integral-preserving transform to the functions in the parametric domain
%
%     sp = sp_integral_preserving_transform (space, msh)
%
% INPUTS:
%
%    space:   structure with the information in the parametric domain (see sp_scalar/sp_evaluate_col)
%    msh:     msh structure containing the information of the parametrization
%              in the points where basis functions have to be computed (see msh_cartesian/msh_evaluate_col)
%
%              Name     |   Default value |  Meaning
%           ------------+-----------------+----------------------------------
%            value      |      true       |  compute shape_functions
%
% OUTPUT:
%
%    sp: struct representing the discrete function space, with the following fields:
%              (see the article for a detailed description)
%
%    FIELD_NAME      (SIZE)                                 DESCRIPTION
%    ncomp           (scalar)                               number of components of the functions of the space (actually, 1)
%    ndof            (scalar)                               total number of degrees of freedom
%    ndof_dir        (ncomp x ndim matrix)                  for each component, number of degrees of freedom along each direction
%    nsh_max         (scalar)                               maximum number of shape functions per element
%    nsh             (1 x msh_col.nel vector)               actual number of shape functions per each element
%    connectivity    (nsh_max x msh_col.nel vector)         indices of basis functions that do not vanish in each element
%    shape_functions (msh_col.nqn x nsh_max x msh_col.nel)  basis functions evaluated at each quadrature node in each element
%
% Copyright (C) 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp = sp_integral_preserving_transform (sp, msh, value)

  if (nargin < 3 && isfield (sp, 'shape_functions'))
    value = true;
  end

  if (value)
    jacdet = reshape (geopdes_det__ (msh.geo_map_jac), msh.nqn, 1, msh.nel);
    sp.shape_functions = bsxfun (@rdivide, sp.shape_functions, jacdet);
    sp.shape_functions = sp.shape_functions;
    if (isfield (msh, 'side_number'))
      sp.shape_functions = sp.shape_functions * (-1)^msh.side_number;
    end
  end

  if (isfield (sp, 'shape_function_gradients'))
    sp = rmfield (sp, 'shape_function_gradients');
  end
  if (isfield (sp, 'shape_function_hessians'))
    sp = rmfield (sp, 'shape_function_hessians');
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
% Author:
% - Mehrtash Harandi (mehrtash.harandi at gmail dot com)
%
% This file is provided without any warranty of
% fitness for any purpose. You can redistribute
% this file and/or modify it under the terms of
% the GNU General Public License (GPL) as published
% by the Free Software Foundation, either version 3
% of the License or (at your option) any later version.

function [outCost,outGrad,covD_Struct] = supervised_WB_CostGrad(U,covD_Struct)
outCost = 0;
dF = zeros(size(U));

nPoints = length(covD_Struct.y);

I_r = eye(covD_Struct.r);
UXU = zeros(covD_Struct.r,covD_Struct.r,nPoints);
inv_UXU = zeros(covD_Struct.r,covD_Struct.r,nPoints);
for tmpC1 = 1:nPoints
    UXU(:,:,tmpC1) = U'*covD_Struct.X(:,:,tmpC1)*U;
    inv_UXU(:,:,tmpC1) = I_r/UXU(:,:,tmpC1);
end

for i = 1:nPoints
    X_i = covD_Struct.X(:,:,i);
    for j = 1:nPoints
        if (covD_Struct.G(i,j) == 0)
            continue;
        end

        X_j = covD_Struct.X(:,:,j);
        switch (covD_Struct.Metric_Flag)

            case 1
                %AIRM
                outCost = outCost + covD_Struct.G(i,j)*Compute_AIRM_Metric(UXU(:,:,i) , UXU(:,:,j));
                log_XY_INV = logm(UXU(:,:,i)*inv_UXU(:,:,j));

                dF = dF + 4*covD_Struct.G(i,j)*((X_i*U)*inv_UXU(:,:,i)  ...
                    -(X_j*U)*inv_UXU(:,:,j) )*log_XY_INV;
             case 2
                %Stein  metric
                outCost = outCost + covD_Struct.G(i,j)*Compute_Stein_Metric(UXU(:,:,i) , UXU(:,:,j));

                X_ij = 0.5*(X_i + X_j);
                dF = dF + covD_Struct.G(i,j)*(2*(X_ij*U)/(U'*X_ij*U)  ...
                    - (X_i*U)*inv_UXU(:,:,i) - (X_j*U)*inv_UXU(:,:,j));
            otherwise
                error('The metric is not implemented.');
        end %end switch

    end
end

outGrad = (eye(size(U,1)) - U*U')*dF;

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [x,z]=constrainedLSEq(A,b,B,d,algorithm)
%%CONSTRAINEDLSEQ Find x to minimize norm(A*x-b)^2 under the constraint
%                 that B*x=d. It is assumed that the (m+p)Xn stacked matrix
%                 [A;B] has linearly independent columns (left invertible)
%                 and B has linearly independent rows (right invertible).
%                 This means that p<=n<=m+p. This is equality constrained
%                 least squares.
%
%INPUTS: A A real mXn matrix
%        b A real mX1 vector.
%        B A real pXn matrix.
%        d A real pX1 vector.
% algorithm An optional prameter specfying the algorithm to use. Possible
%          values are:
%          0 (The default if omitted or an empty matrix is passed) Use
%            Algorithm  6.2.2 in Chapter 6.2.3 of [1].
%          1 In Chapter 14 of [2], it is shown that x solves the
%            constrained least squares problem if and only if there exists
%            a z such that [A'*A, B'; B, zeros(p,p)]*[x;z]=[A'*b;d]
%            This is a set of n+p equations and unknowns. This solution
%            just solves that system.
%          2 Use the algorithm of Chapter 6.2.5 of [1], which makes use of
%            the generalized singular value decomposition (GSVD).
%
%OUTPUTS: x The solution to the equality constrained least squares problem.
%           An mX1 vector.
%         z If algorithm=0, this is an (n-p)X1 unconstrained solution that
%           is an intermediate result. If algorithm=1, this is the pX1 set
%           of Lagrangian multipliers used in the solution. If algorithm=2,
%           this is an empty matrix.
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%[2] S. Boyd and L. Vandenberghe, Vectors, Matrices, and Least Squares,
%    Sep. 2015, draft version of book. [Online].
%    Available: http://www.seas.ucla.edu/~vandenbe/133A/133A-textbook.pdf
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(algorithm))
    algorithm=0;
end

switch(algorithm)
    case 0%Algorithm 6.2.2 in Chapter 6.2.3 of [1].
        [~,n1]=size(A);
        [m2,~]=size(B);

        [Q,R]=qr(B');
        opts.UT=false;
        opts.LT=true;
        opts.RECT=false;
        y=linsolve(R(1:m2,1:m2)',d,opts);

        A=A*Q;
        z=linearLeastSquares(A(:,(m2+1):n1),b-A(:,1:m2)*y);
        x=Q(:,1:m2)*y+Q(:,(m2+1):n1)*z;
    case 1%From Chapter 14 of [2].
        m=size(A,2);
        p=size(B,1);
        H=A'*A;
        c=A'*b;
        xz=[H,B';
            B,zeros(p,p)]\[c;d];

        x=xz(1:m);
        z=xz((m+1):end); %The Lagrangian multipliers.
    case 2%The GSVD approach of Chapter 6.2.5 of [1].
        [~,n1]=size(A);
        [m2,~]=size(B);

        [U1,U2,X,DA,DB]=gsvd(A,B);
        bTilde=U1'*b;
        dTilde=U2'*d;

        %The definition of the GSVD used in Chapter 6.1.6 of [1] is not the
        %same as that implemented by Matlab's gsvd function. Specifically,
        %we have to modify X as follows:
        X=inv(X');

        alphaVals=diag(DA);
        betaVals=diag(DB);

        x=zeros(n1,1);
        for i=1:m2
            x=x+(dTilde(i)/betaVals(i))*X(:,i);
        end

        for i=(m2+1):n1
            x=x+(bTilde(i)/alphaVals(i))*X(:,i);
        end
        z=[];
    otherwise
        error('Unknown algorithm specified.')
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function imStack=analyze2mrLoadRet3TSeries(inFileRoot,outFileRoot,nVols,firstVolIndex,doRotate,scaleFact,flipFlag)
% function imStack=analyze2mrLoadRet3TSeries(inFileRoot,outFileRoot,nVols,firstVolIndex,doRotate,scaleFact,flipFlag)
% Converts from analyze functional image data to mrLoadRet TSeries format
% Analyze functional data are stored as N individual volumes, each with S slices
% mrLoadRet has S individual files with N acquisitions in each.
% the doRotate param allows to you roate the functional data by doRotate*90 degrees
% Try to read the the first volume to see if it's there and get the dimensions of all the rest

% TODO: 1:  Rotate these so that they are in the same orientation as the anatomicals
% (Maybe make this a flag)
% 2: Why can't we see the anatomicals? Claims not to be able to read from /Raw/Anatomy/Inplane
% 3: Why can't we load in the corAnal? We do this by hand perfectly well but it won't do it from mrLoadRet
% 4: Do we want to resample the functional data? Or is this the moment to go to mrLoadRet3?
%
%
% ARW 032703 : Now saves out data in mrLoadRet3.0 format (.mat as opposed to .dat files)

if (~exist('firstVolIndex','var'))
    firstVolIndex=1;
end

if (~exist('scaleFact','var'))
    scaleFact=[1 1]; % No interpolation. (Scaling=1)
end

if (length(scaleFact)~=2) % If we just get a scalar for the scale factor, assume that it applies in both dimensions
    scaleFact=repmat(scaleFact(1),2);
end

if (~exist('doRotate','var'))
    doRotate=1; % This is on by default. Rotates 1*90 degrees
end

if (~exist('flipFlag','var'))
    flipFlag=1; % This is on by default. Flips up/down after rotation
end
suffix=sprintf('%04d',firstVolIndex);

fileName=[inFileRoot,suffix,'.hdr'];

V=spm_vol(fileName);
im=spm_read_vols(V);

[y,x,nSlices]=size(im);
fprintf('Read in a volume of size %d, %d, %d',y,x,nSlices);

% Pre-allocate memory: This is a monster and will fail on many machines.
fprintf('\nTrying to allocate an array with %d elements...\n',y*x*nSlices*nVols);

if (mod(doRotate,2)) % When we rotate by 180 degrees the x and y dimensions remain unchanged
    funcVol=zeros(y,x,nSlices,nVols);
else
    funcVol=zeros(x,y,nSlices,nVols);
end
 fprintf('Rotating by %d x 90, flipFlag=%d',doRotate,flipFlag);

for t=0:(nVols-1)
    thisImIndex=t+firstVolIndex;
    suffix=sprintf('%04d',thisImIndex);
    fileName=[inFileRoot,suffix];
    V=spm_vol(fileName);
    im=spm_read_vols(V);

    % Do the rotation and scaling

    if (mod(doRotate,2)) % When we rotate by 180 degrees the x and y dimensions remain unchanged
        im2=zeros(y,x,nSlices);
    else
        im2=zeros(x,y,nSlices);
    end
    fprintf('\nVol=%d',thisImIndex);

    for thisSlice=1:nSlices
        imSlice=squeeze(im(:,:,thisSlice));
    %s    imSlice=imresize(imSlice,[scaleFact(1)*y,scaleFact(2)*x],'nearest');

        im2(:,:,thisSlice)=rot90(imSlice,doRotate);
        if (flipFlag)

            im2(:,:,thisSlice)=flipud(im2(:,:,thisSlice));
        end

    end % next imSlice

    funcVol(:,:,:,t+1)=im2;
    %fprintf('.');

end
size(funcVol);
[y x nSlices nVols]=size(funcVol);

% Now write them out in a different format
fprintf('\nDone reading data: Writing now...\n');
for t=1:nSlices
    suffix=int2str(t);
    tSeries=squeeze(funcVol(:,:,t,:));
    [a,b,c]=size(tSeries);
    fprintf('\nSize before: %d by %d by %d vols ',a,b,c);

    tSeries=squeeze(shiftdim(tSeries,2));

    [a,b,c]=size(tSeries);
    fprintf('\nSize after:%d samples by %d by %d',a,b*scaleFact(1),c*scaleFact(2));

    % Reshape tSeries here
    tSeries=reshape(tSeries,a,(b*c));

    if ~exist(outFileRoot,'dir')
      mkdir(outFileRoot);
    end

    pathStr = fullfile(outFileRoot,['tSeries',num2str(t)]);
%disp(['Saving: ',pathStr]);
    save(pathStr,'tSeries');

    fprintf('_');
end
fprintf('\nDone\n');
imStack='pathStr';

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [futileCycleReactions,FutileCyclesTest]=identifyFutileCycles(model)
% This function serves to debug a refined reconstruction generated by the
% DEMETER pipeline that produced infeasible amounts of ATP. The output
% futileCycleReactions shows reactions that are involved in futile cycles.
% To eliminate the futile cycle, one of these reactions needs to be
% replaced by an irreversible version of the same reaction. To accomplish
% this, the replacement step needs to be added to the function
% removeFutileCycles. It is requested that identified solutions are
% submitted as a pull request.
%
% USAGE:
%
%   [gapfilledReactions,replacedReactions,revisedModel]=debugModel(model,testResultsFolder, inputDataFolder,reconVersion,microbeID,biomassReaction)
%
% INPUT
% model:                   COBRA model structure
%
% OUTPUT
% futileCycleReactions:    Reactions that are involved in futile cycles and
%                          need to be constrained
% FutileCyclesTest         Computed fluxes for ATP production
%
% .. Authors:
%       - Almut Heinken, 2016-2020

% load complex medium
constraints = readtable('ComplexMedium.txt', 'Delimiter', 'tab');
constraints=table2cell(constraints);
constraints=cellstr(string(constraints));

% apply complex medium
model = useDiet(model,constraints);
model=changeRxnBounds(model,'EX_o2(e)',-10,'l');
model=changeObjective(model,'DM_atp_c_');
FBAorg=optimizeCbModel(model,'max');

if FBAorg.f > 150

Rxns=printRxnFormula(model);
for i=1:length(model.rxns)
    if ~strcmp(model.rxns{i},'DM_atp_c_') && ~strncmp(model.rxns{i},'ATPS',4)
        modelTest=model;
        FutileCyclesTest{i,1}=model.rxns{i,1};
        modelTest.lb(i,1)=0;
        modelTest.ub(i,1)=0;
        FBA=optimizeCbModel(modelTest,'max');
        FutileCyclesTest{i,2}=Rxns{i,1};
        FutileCyclesTest{i,3}=model.lb(i,1);
        FutileCyclesTest{i,4}=model.ub(i,1);
        FutileCyclesTest{i,5}=FBA.f;
        FutileCyclesTest{i,6}=FBAorg.x(i,1);
    end
end

FutileCyclesTest(find(cellfun(@isempty, FutileCyclesTest(:,1))),:)=[];
futileCycleReactions=FutileCyclesTest(cell2mat(FutileCyclesTest(:,5))<150,1);

else
    FutileCyclesTest = {};
    futileCycleReactions = {};
end

end

author: opencobra

repo: cobratoolbox


text:
function SeqsNew = SimulationFast_ConditionalThinning_ExpHP(SeqsOld, ...
                                                        para, options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The fast simulation of Hawkes processes with exponential kernels
% conditioned on history
%
% Reference:
% Dassios, Angelos, and Hongbiao Zhao.
% "Exact simulation of Hawkes process with exponentially decaying intensity."
% Electronic Communications in Probability 18.62 (2013): 1-13.
%
% Provider:
% Hongteng Xu @ Georgia Tech
% June 13, 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SeqsNew = struct('Time', [], ...
              'Mark', [], ...
              'Start', [], ...
              'Stop', [], ...
              'Feature', []);

tic
for n = 1:length(SeqsOld)

    t=SeqsOld(n).Stop;
    History = [SeqsOld(n).Time; SeqsOld(n).Mark];

    lambdat = Intensity_HP(t, History, para);
    mt = sum(lambdat);

    while t<options.Tmax && size(History, 2)<options.Nmax

        s = random('exp', 1/mt);
        U = rand;

        lambda_ts = Intensity_Recurrent_HP(t+s, [], t, lambdat, para);
        mts = sum(lambda_ts);

        %fprintf('s=%f, v=%f\n', s, mts/mt);
        if t+s>options.Tmax || U>mts/mt
            t = t+s;
            lambdat = lambda_ts;
        else
            u = rand*mts;
            sumIs = 0;
            for d=1:length(lambda_ts)
                sumIs = sumIs + lambda_ts(d);
                if sumIs >= u
                    break;
                end
            end
            index = d;

            lambdat = Intensity_Recurrent_HP(t+s, index(1), t, lambdat, para);
            t = t+s;
            History = [History,[t;index(1)]];
        end

        mt = sum(lambdat);

    end

    SeqsNew(n).Time = History(1,:);
    SeqsNew(n).Mark = History(2,:);
    SeqsNew(n).Start = SeqsOld(n).Stop;
    SeqsNew(n).Stop = options.Tmax;
    index = find(SeqsOld(n).Stop<=SeqsNew(n).Time & ...
        SeqsNew(n).Time<=options.Tmax);
    SeqsNew(n).Time = SeqsNew(n).Time(index);
    SeqsNew(n).Mark = SeqsNew(n).Mark(index);

    if mod(n, 10)==0 || n==options.N
        fprintf('#seq=%d/%d, #event=%d, time=%.2fsec\n', ...
            n, options.N, length(SeqsNew(n).Mark), toc);
    end
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function SO3F = norm(SO3VF)
% pointwise norm of the vector field
%
% Syntax
%   SO3F = norm(SO3VF)
%
% Input
%  SO3VF - @SO3VectorField
%
% Output
%  SO3F - @SO3FunHarmonic
%

SO3F = SO3FunHarmonic.quadrature(@(rot) norm(SO3VF.eval(rot)),SO3VF.CS,SO3VF.SS);

end


text:
function varargout = pivotplot(varargin)
%PIVOTPLOT   Semilogy plot of pivot values.
%   PIVOTPLOT( F ) semilogy plot related to the pivots values taken during
%   the construction of the SPHEREFUN F.
%
%   H = PIVOTPLOT( F ) returns a handle H to the figure.
%
%   PIVOTPLOT( F, S ) allows further plotting options, such as linestyle,
%   linecolor, etc. If S contains a string 'LOGLOG', the pseudo sig will be
%   displayed on a log-log scale.

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

[varargout{1:nargout}] = pivotplot@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function A = cs_sparse2 (i,j,x)                                             %#ok
%CS_SPARSE2 same as cs_sparse, to test cs_entry function
% A = cs_sparse2 (i,j,x), removing duplicates and numerically zero entries,
% and returning A sorted (test cs_entry)
%
% Example:
%   A = cs_sparse2 (i,j,x)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_sparse2 mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function A = four_corners(m,n)
%FOUR_CORNERS Construct a m x n sparse matrix with ones in each of its four
%corners.
%
% t = four_corners(m,n);
%
% Inputs:
%  m,n  the dimension of the matrix requested
% Outputs:
%  A  the sparse matrix with ones in their four corners
%

A = sparse(m,n);
A(1,1) = 1;
A(m,1) = 1;
A(1,n) = 1;
A(m,n) = 1;

end

author: odedstein

repo: sgi-introduction-course


text:
classdef twix_map_obj < handle

% class to hold information about raw data from siemens MRI scanners
% (currently VB and VD software versions are supported and tested).
%
% Author: Philipp Ehses (philipp.ehses@tuebingen.mpg.de), Aug/19/2011
%
%
% Modified by Wolf Blecher (wolf.blecher@tuebingen.mpg.de), Apr/26/2012
% Added reorder index to indicate which lines are reflected
% Added slice position for sorting, Mai/15/2012
%
% Order of many mdh parameters are now stored (including the reflected ADC
% bit); PE, Jun/14/2012
%
% data is now 'memory mapped' and not read until demanded;
% (see mapVBVD for a description) PE, Aug/02/2012
%
% twix_obj.image.unsorted now returns the data in its acq. order
% [NCol,NCha,nsamples in acq. order], all average flags don't have an
% influence on the output, but 'flagRemoveOS' still works, PE, Sep/04/13

properties(Dependent=true)
    readerVersion
    % flags
    flagRemoveOS        % removes oversampling in read (col) during read operation
    flagDoAverage       % averages over all avg during read operation
    flagAverageReps     % averages over all repetitions
    flagAverageSets     % averages over all sets
    flagIgnoreSeg       % sum over all segments during read operation
    flagSkipToFirstLine % skips lines/partitions up to the first
                        % actually acquired line/partition
                        % (e.g. only the center k-space is acquired in
                        % refscans, we don't want all the leading zeros
                        % in our data)
                        % this is the default behaviour for everything
                        % but image scans (but can be changed manually)
    flagRampSampRegrid  % perform on-the-fly ramp sampling regridding
    flagDoRawDataCorrect     %SRY: apply raw data correction factors during read operation

    RawDataCorrectionFactors %SRY: allow the user to set/get the factors
end

properties(GetAccess='public', SetAccess='public')
    flagAverageDim      % new: flags that determines whether certain dim. should be averaged/ignored
    filename
    softwareVersion
    dataType
    rampSampTrj
end

properties(Dependent=true)
    dataSize % this is the current output size, depends on fullSize + some flags
    sqzSize
    sqzDims
end

properties(GetAccess='public', SetAccess='protected')
    dataDims

    NCol  % mdh information
    NCha  % mdh information
    NLin  % mdh information
    NPar  % mdh information
    NSli  % mdh information
    NAve  % mdh information
    NPhs  % mdh information
    NEco  % mdh information
    NRep  % mdh information
    NSet  % mdh information
    NSeg  % mdh information
    NIda  % mdh information
    NIdb  % mdh information
    NIdc  % mdh information
    NIdd  % mdh information
    NIde  % mdh information
    NAcq  % simple counter

    % mdh information
    Lin
    Par
    Sli
    Ave
    Phs
    Eco
    Rep
    Set
    Seg
    Ida
    Idb
    Idc
    Idd
    Ide

    centerCol
    centerLin
    centerPar
    cutOff
    coilSelect
    ROoffcenter
    timeSinceRF
    IsReflected
    IsRawDataCorrect %SRY: storage for MDH flag raw data correct

    slicePos
    freeParam
    iceParam
    scancounter
    timestamp
    pmutime

    % memory position in file
    memPos

    isBrokenFile % errors when parsing?
end

properties(Hidden=true, SetAccess='protected')
    arg  % arguments

    fullSize % this is the full size of the data set according to the mdhs, i.e. flags
             % like 'reduceOS' have no influence on it

    freadInfo

    skipLin
    skipPar
end

methods
    % Constructor:
    function this = twix_map_obj(arg,dataType,fname,version,rstraj)

        if ~exist('dataType','var')
            this.dataType = 'image';
        else
            this.dataType = lower(dataType);
        end

        this.filename         = fname;
        this.softwareVersion  = version;

        this.IsReflected      = logical([]);
        this.IsRawDataCorrect = logical([]); %SRY
        this.NAcq             = 0;
        this.isBrokenFile     = false;

        this.dataDims = {'Col','Cha','Lin','Par','Sli','Ave','Phs',...
            'Eco','Rep','Set','Seg','Ida','Idb','Idc','Idd','Ide'};

        this.setDefaultFlags();
        if exist('arg','var')
            % copy relevant arguments from mapVBVD argument list
            names=fieldnames(arg);
            for k=1:numel(names)
                if isfield(this.arg,names{k})
                    this.arg.(names{k}) = arg.(names{k});
                end
            end
        end

        this.flagAverageDim(ismember(this.dataDims,'Ave')) = this.arg.doAverage;
        this.flagAverageDim(ismember(this.dataDims,'Rep')) = this.arg.averageReps;
        this.flagAverageDim(ismember(this.dataDims,'Set')) = this.arg.averageSets;
        this.flagAverageDim(ismember(this.dataDims,'Seg')) = this.arg.ignoreSeg;

        switch this.softwareVersion
            case 'vb'
                % every channel has its own full mdh
                this.freadInfo.szScanHeader    =   0; % [bytes]
                this.freadInfo.szChannelHeader = 128; % [bytes]
                this.freadInfo.iceParamSz      =   4;
            case 'vd'
                if ( this.arg.doRawDataCorrect )
                    error('raw data correction for VD not supported/tested yet');
                end
                this.freadInfo.szScanHeader    = 192; % [bytes]
                this.freadInfo.szChannelHeader =  32; % [bytes]
                this.freadInfo.iceParamSz      =  24; % vd version supports up to 24 ice params
            otherwise
                error('software version not supported');
        end

        if exist('rstraj','var')
            this.rampSampTrj = rstraj;
        else
            this.rampSampTrj        = [];
            this.arg.rampSampRegrid = false;
        end
    end

    % Copy function - replacement for matlab.mixin.Copyable.copy() to create object copies
    % from http://undocumentedmatlab.com/blog/general-use-object-copy
    function newObj = copy(this)
        isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
        if isOctave || verLessThan('matlab','7.11')
            % R2010a or earlier - serialize via temp file (slower)
            fname = [tempname '.mat'];
            save(fname, 'obj');
            newObj = load(fname);
            newObj = newObj.obj;
            delete(fname);
        else
            % R2010b or newer - directly in memory (faster)
            objByteArray = getByteStreamFromArray(this);
            newObj = getArrayFromByteStream(objByteArray);
        end
    end

    function this = readMDH(this, mdh, filePos )
        % extract all values in all MDHs at once
        %
        % data types:
        % Use double for everything non-logical, both ints and floats. Seems the
        % most robust way to avoid unexpected cast-issues with very nasty side effects.
        % Examples: eps(single(16777216)) == 2
        %           uint32( 10 ) - uint32( 20 ) == 0
        %           uint16(100) + 1e5 == 65535
        %           size(1 : 10000 * uint16(1000)) ==  [1  65535]
        %
        % The 1st example always hits the timestamps.

        if ~isstruct( mdh ) || isempty( mdh )
            return
        end

        this.NAcq     = numel( filePos );
        sLC           = double( mdh.sLC ) + 1;  % +1: convert to matlab index style
        evalInfoMask1 = double( mdh.aulEvalInfoMask(:,1) ).';

        % save mdh information for each line
        this.NCol       = double( mdh.ushSamplesInScan ).';
        this.NCha       = double( mdh.ushUsedChannels ).';
        this.Lin        = sLC(:,1).' ;
        this.Ave        = sLC(:,2).' ;
        this.Sli        = sLC(:,3).' ;
        this.Par        = sLC(:,4).' ;
        this.Eco        = sLC(:,5).' ;
        this.Phs        = sLC(:,6).' ;
        this.Rep        = sLC(:,7).' ;
        this.Set        = sLC(:,8).' ;
        this.Seg        = sLC(:,9).' ;
        this.Ida        = sLC(:,10).';
        this.Idb        = sLC(:,11).';
        this.Idc        = sLC(:,12).';
        this.Idd        = sLC(:,13).';
        this.Ide        = sLC(:,14).';

        this.centerCol   = double( mdh.ushKSpaceCentreColumn ).' + 1;
        this.centerLin   = double( mdh.ushKSpaceCentreLineNo ).' + 1;
        this.centerPar   = double( mdh.ushKSpaceCentrePartitionNo ).' + 1;
        this.cutOff      = double( mdh.sCutOff ).';
        this.coilSelect  = double( mdh.ushCoilSelect ).';
        this.ROoffcenter = double( mdh.fReadOutOffcentre ).';
        this.timeSinceRF = double( mdh.ulTimeSinceLastRF ).';
        this.IsReflected = logical(min(bitand(evalInfoMask1,2^24),1));
        this.scancounter = double( mdh.ulScanCounter ).';
        this.timestamp   = double( mdh.ulTimeStamp ).';
        this.pmutime     = double( mdh.ulPMUTimeStamp ).';
        this.IsRawDataCorrect = logical(min(bitand(evalInfoMask1,2^10),1)); %SRY
        this.slicePos    = double( mdh.SlicePos ).';
        this.iceParam    = double( mdh.aushIceProgramPara ).';
        this.freeParam   = double( mdh.aushFreePara ).';

        this.memPos = filePos;

    end % of readMDH

    function this = tryAndFixLastMdh(this)
        eofWarning = [mfilename() ':UnxpctdEOF'];   % We have it inside this.readData()
        warning( 'off', eofWarning )    % silence warnings for read...
        warning( 'off', 'foo:bar'  )    % ... and a stupid placeholder

        isLastAcqGood = false;
        cnt = 0;

        while ~isLastAcqGood  &&  this.NAcq > 0  && cnt < 100
            warning( 'foo:bar', 'baz') % make sure that lastwarn() does not return eofWarning
            try
                this.clean();
                this.unsorted(this.NAcq);
                [~, warnid] = lastwarn();
                if strcmp( warnid, eofWarning )
                    error( 'Make sure to go to the catch block.')
                end
                isLastAcqGood = true;
            catch
                this.isBrokenFile = true;
                this.NAcq = this.NAcq-1;
            end
            cnt = cnt + 1;
        end

%         if this.NAcq == 0  ||  cnt > 99    % everything is garbage
%             warning(  )
%         end

        warning( 'on', eofWarning )
    end

    function this = clean(this)

        if this.NAcq == 0
            return;
        end

        % Cut mdh data to actual size. Maybe we rejected acquisitions at the end
        % due to read errors.
        fields = { 'NCol', 'NCha',                                                  ...
                   'Lin', 'Par', 'Sli', 'Ave', 'Phs', 'Eco', 'Rep',                 ...
                   'Set', 'Seg', 'Ida', 'Idb', 'Idc', 'Idd', 'Ide',                 ...
                   'centerCol'  ,   'centerLin',   'centerPar',           'cutOff', ...
                   'coilSelect' , 'ROoffcenter', 'timeSinceRF',      'IsReflected', ...
                   'scancounter',   'timestamp',     'pmutime', 'IsRawDataCorrect', ...
                   'slicePos'   ,    'iceParam',   'freeParam',           'memPos'  };

        nack = this.NAcq;
        idx = 1:nack;

        for f = fields
            f1 = f{1};
            if size(this.(f1),2) > nack   % rarely
                this.(f1) = this.(f1)(:,idx);       % 1st dim: samples,  2nd dim acquisitions
            end
        end

        this.NLin = max(this.Lin);
        this.NPar = max(this.Par);
        this.NSli = max(this.Sli);
        this.NAve = max(this.Ave);
        this.NPhs = max(this.Phs);
        this.NEco = max(this.Eco);
        this.NRep = max(this.Rep);
        this.NSet = max(this.Set);
        this.NSeg = max(this.Seg);
        this.NIda = max(this.Ida);
        this.NIdb = max(this.Idb);
        this.NIdc = max(this.Idc);
        this.NIdd = max(this.Idd);
        this.NIde = max(this.Ide);

        % ok, let us assume for now that all NCol and NCha entries are
        % the same for all mdhs:
        this.NCol = this.NCol(1);
        this.NCha = this.NCha(1);

        if strcmp(this.dataType,'refscan')
            %pehses: check for lines with 'negative' line/partition numbers
            %this can happen when the reference scan line/partition range
            %exceeds the one of the actual imaging scan
            if this.NLin>65500  %uint overflow check
                this.Lin  = mod(this.Lin + (65536 - min(this.Lin(this.Lin>65500))),65536)+1;
                this.NLin = max(this.Lin);
            end
            if this.NPar>65500  %uint overflow check
                this.Par  = mod(this.Par + (65536 - min(this.Par(this.Par>65500))),65536)+1;
                this.NPar = max(this.Par);
            end
        end

        % to reduce the matrix sizes of non-image scans, the size
        % of the refscan_obj()-matrix is reduced to the area of the
        % actually scanned acs lines (the outer part of k-space
        % that is not scanned is not filled with zeros)
        % this behaviour is controlled by flagSkipToFirstLine which is
        % set to true by default for everything but image scans
        if ~this.flagSkipToFirstLine
            % the output matrix should include all leading zeros
            this.skipLin = 0;
            this.skipPar = 0;
        else
            % otherwise, cut the matrix size to the start of the
            % first actually scanned line/partition (e.g. the acs/
            % phasecor data is only acquired in the k-space center)
            this.skipLin = min(this.Lin)-1;
            this.skipPar = min(this.Par)-1;
        end
        NLinAlloc = max(1, this.NLin - this.skipLin);
        NParAlloc = max(1, this.NPar - this.skipPar);

        this.fullSize = [ this.NCol this.NCha NLinAlloc NParAlloc...
                          this.NSli this.NAve this.NPhs this.NEco...
                          this.NRep this.NSet this.NSeg this.NIda...
                          this.NIdb this.NIdc this.NIdd this.NIde ];

        nByte = this.NCha*(this.freadInfo.szChannelHeader+8*this.NCol);

        % size for fread
        this.freadInfo.sz    = [2 nByte/8];
        % reshape size
        this.freadInfo.shape = [this.NCol+this.freadInfo.szChannelHeader/8 ...
                               , this.NCha];
        % we need to cut MDHs from fread data
        this.freadInfo.cut   = this.freadInfo.szChannelHeader/8 + (1 : this.NCol);

    end % of clean

    function varargout = subsref(this, S)
        % this is where the magic happens
        % Now seriously. Overloading of the subsref-method and working
        % with a gazillion indices got really messy really fast. At
        % some point, I should probably clean this code up a bit. But
        % good news everyone: It seems to work.
        switch S(1).type
            case '.'
                % We don't want to manage method/variable calls, so we'll
                % simply call the built-in subsref-function in this case.
                if nargout == 0
                    varargout{1} = builtin('subsref', this, S); % CTR fix.
                else
                    varargout      = cell(1, nargout);
                    [varargout{:}] = builtin('subsref', this, S);
                end
                return;
            case {'()','{}'}
            otherwise
                error('operator not supported');
        end

        [selRange,selRangeSz,outSize] = this.calcRange(S(1));

    	% calculate page table (virtual to physical addresses)
        % this is now done every time, i.e. result is no longer saved in
        % a property - slower but safer (and easier to keep track of updates)
        ixToRaw = this.calcIndices;

        tmp = reshape(1:prod(double(this.fullSize(3:end))), this.fullSize(3:end));
        tmp = tmp(selRange{3:end});
        ixToRaw = ixToRaw(tmp); clear tmp;
        ixToRaw = ixToRaw(:);
        % delete all entries that point to zero (the "NULL"-pointer)
        notAcquired = (ixToRaw == 0);
        ixToRaw (notAcquired) = []; clear notAcquired;

        % calculate ixToTarg for possibly smaller, shifted + segmented
        % target matrix:
        cIx = ones(14, numel(ixToRaw), 'single');
        if ~this.flagAverageDim(3)
            cIx( 1,:) = this.Lin(ixToRaw) - this.skipLin;
        end
        if ~this.flagAverageDim(4)
            cIx( 2,:) = this.Par(ixToRaw) - this.skipPar;
        end
        if ~this.flagAverageDim(5)
            cIx( 3,:) = this.Sli(ixToRaw);
        end
        if ~this.flagAverageDim(6)
            cIx( 4,:) = this.Ave(ixToRaw);
        end
        if ~this.flagAverageDim(7)
            cIx( 5,:) = this.Phs(ixToRaw);
        end
        if ~this.flagAverageDim(8)
            cIx( 6,:) = this.Eco(ixToRaw);
        end
        if ~this.flagAverageDim(9)
            cIx( 7,:) = this.Rep(ixToRaw);
        end
        if ~this.flagAverageDim(10)
            cIx( 8,:) = this.Set(ixToRaw);
        end
        if ~this.flagAverageDim(11)
            cIx( 9,:) = this.Seg(ixToRaw);
        end
        if ~this.flagAverageDim(12)
            cIx(10,:) = this.Ida(ixToRaw);
        end
        if ~this.flagAverageDim(13)
            cIx(11,:) = this.Idb(ixToRaw);
        end
        if ~this.flagAverageDim(14)
            cIx(12,:) = this.Idc(ixToRaw);
        end
        if ~this.flagAverageDim(15)
            cIx(13,:) = this.Idd(ixToRaw);
        end
        if ~this.flagAverageDim(16)
            cIx(14,:) = this.Ide(ixToRaw);
        end

        % make sure that indices fit inside selection range
        for k=3:numel(selRange)
            tmp = cIx(k-2,:);
            for l=1:numel(selRange{k})
                cIx(k-2,tmp==selRange{k}(l)) = l;
            end
        end

        sz = selRangeSz(3:end); % extra variable needed for octave compatibility
        ixToTarg = this.sub2ind_double(sz, cIx(1,:),cIx(2,:),cIx(3,:),...
            cIx(4,:),cIx(5,:),cIx(6,:),cIx(7,:),cIx(8,:),cIx(9,:),...
            cIx(10,:),cIx(11,:),cIx(12,:),cIx(13,:),cIx(14,:));

        mem = this.memPos(ixToRaw);
        % sort mem for quicker access, sort cIxToTarg/Raw accordingly
        [mem,ix]  = sort(mem);
        ixToTarg = ixToTarg(ix);
        ixToRaw  = ixToRaw(ix);
        clear ix;

        % For a call of type data{:,:,1:3} matlab expects more than one
        % output variable (three in this case) and will throw an error
        % otherwise. This is a lazy way (and the only one I know of) to
        % fix this.
        varargout    = cell(1, nargout);
        varargout{1} = this.readData(mem,ixToTarg,ixToRaw,selRange,selRangeSz,outSize);
    end % of subsref

    function out = unsorted(this,ival)
        % returns the unsorted data [NCol,NCha,#samples in acq. order]
        if ~exist('ival','var')
            mem = this.memPos;
        else
            mem = this.memPos(ival);
        end
        out = this.readData(mem);
    end

    function out = readData(this,mem,cIxToTarg,cIxToRaw,selRange,selRangeSz,outSize)

        if ~exist('outSize','var')
            selRange{1} = ':';
            selRange{2} = ':';
            outSize = [this.dataSize(1:2),numel(mem)];
            selRangeSz = outSize;
            cIxToTarg = 1:selRangeSz(3);
            cIxToRaw  = cIxToTarg;
        else
            if isequal( selRange{1}(:), (1:this.dataSize(1)).' )
                selRange{1} = ':';
            end
            if isequal( selRange{2}(:), (1:this.dataSize(2)).' )
                selRange{2} = ':';
            end
        end
        out = complex(zeros(outSize,'single'));
        out = reshape(out, selRangeSz(1), selRangeSz(2), []);

        if isempty( mem )
            out = reshape(out,outSize);
            return
        end

        cIxToTarg = this.cast2MinimalUint( cIxToTarg );

        % subsref overloading makes this.that-calls slow, so we need to
        % avoid them whenever possible
        szScanHeader = this.freadInfo.szScanHeader;
        readSize     = this.freadInfo.sz;
        readShape    = this.freadInfo.shape;
        readCut      = this.freadInfo.cut;
        keepOS       = [1:this.NCol/4, 1+this.NCol*3/4:this.NCol];
        bRemoveOS    = this.arg.removeOS;
        bIsReflected = this.IsReflected(cIxToRaw);
        bRegrid      = this.flagRampSampRegrid && numel(this.rampSampTrj);
        slicedata    = this.slicePos(:, cIxToRaw);
        %SRY store information about raw data correction
        bDoRawDataCorrect = this.arg.doRawDataCorrect;
        bIsRawDataCorrect = this.IsRawDataCorrect( cIxToRaw );
        isBrokenRead      = false;
        if (bDoRawDataCorrect)
            rawDataCorrect = this.arg.rawDataCorrectionFactors;
        end

        % MiVö: Raw data are read line-by-line in portions of 2xNColxNCha float32 points (2 for complex).
        % Computing and sorting(!) on these small portions is quite expensive, esp. when
        % it employs non-sequential memory paths. Examples are non-linear k-space acquisition
        % or reflected lines.
        % This can be sped up if slightly larger blocks of raw data are collected, first.
        % Whenever a block is full, we do all those operations and save it in the final "out" array.
        % What's a good block size? Depends on data size and machine (probably L2/L3/L4 cache sizes).
        % So...? Start with a small block, measure the time-per-line and double block size until
        % a minimum is found. Seems sufficiently robust to end up in a close-to-optimal size for every
        % machine and data.
        blockSz   = 2;          % size of blocks; must be 2^n; will be increased
        doLockblockSz = false;  % whether blockSZ should be left untouched
        tprev     = inf;        % previous time-per-line
        blockCtr  = 0;
        blockInit = -inf(readShape(1), readShape(2), blockSz, 'single'); %init with garbage
        blockInit = complex( blockInit );
        block     = blockInit;

        if bRegrid
            v1       = single(1:selRangeSz(2));
            v2       = single(1:blockSz);
            rsTrj    = {this.rampSampTrj,v1,v2};
            trgTrj   = linspace(min(this.rampSampTrj),max(this.rampSampTrj),this.NCol);
            trgTrj   = {trgTrj,v1,v2};
        end

        % counter for proper scaling of averages/segments
        count_ave = zeros([1 1 size(out,3)],'single');
        kMax      = numel( mem );   % max loop index

        fid = this.fileopen();

        for k = 1:kMax
            % skip scan header
            fseek(fid,mem(k) + szScanHeader,'bof');
            raw = fread(fid, readSize, 'float=>single').';

            % MiVö: With incomplete files fread() returns less than readSize points. The subsequent reshape will therefore error out.
            %       We could check if numel(raw) == prod(readSize), but people recommend exception handling for performance
            %       reasons. Do it.
            try
                raw = reshape( complex(raw(:,1), raw(:,2)), readShape);
            catch exc
                offset_bytes = mem(k) + szScanHeader;
                %remainingSz = readSize(2) - size(raw,1);
                warning( [mfilename() ':UnxpctdEOF'],  ...
                          [ '\nAn unexpected read error occurred at this byte offset: %d (%g GiB)\n'...
                            'Actual read size is [%s], desired size was: [%s]\n'                    ...
                            'Will ignore this line and stop reading.\n'                             ...
                            '=== MATLABs error message ================\n'                          ...
                            exc.message                                                             ...
                            '\n=== end of error =========================\n'                        ...
                            ], offset_bytes, offset_bytes/1024^3, num2str(size(raw)), num2str(readSize.') )

                % Reject this data fragment. To do so, init with the values of blockInit
                clear raw
                raw( 1:prod(readShape) ) = blockInit(1);
                raw = reshape( raw, readShape );
                isBrokenRead = true;   % remember it and bail out later
            end

            blockCtr = blockCtr + 1;
            block(:,:,blockCtr) = raw;  % fast serial storage in a cache array

            % Do expensive computations and reorderings on the gathered block.
            % Unfortunately, a lot of code is necessary, but that is executed much less
            % frequent, so its worthwhile for speed.
            % TODO: Do *everything* block-by-block
            if blockCtr == blockSz || k == kMax || (isBrokenRead && blockCtr > 1)
                s = tic;    % measure the time to process a block of data

                % remove MDH data from block:
                block = block(readCut,:,:);

                if bRegrid
                    % correct for readout shifts
                    % the nco frequency is always scaled to the max.
                    % gradient amp and does account for ramp-sampling
                    ro_shift = this.calcOffcenterShiftRO(slicedata(:,k));
                    deltak = max(abs(diff(rsTrj{1})));
                    phase = (0:this.NCol-1).' * deltak * ro_shift;
                    phase_factor = exp(1j*2*pi*(phase - ro_shift * rsTrj{1}));
                    block = bsxfun(@times, block, phase_factor);

                    % grid the data
                    F = griddedInterpolant(rsTrj, block);
                    block = F(trgTrj);
                end

                ix = 1 + k - blockCtr : k;
                if blockCtr ~= blockSz
                    block = block(:,:,1:blockCtr);
                end

                if sum(isnan(block(:)))>0
                    keyboard
                end

                if bRemoveOS % remove oversampling in read
                    block = ifft(block, [], 1);
                    block =  fft(block(keepOS,:,:), [], 1);
                end

                if ( bDoRawDataCorrect && bIsRawDataCorrect(k) )
                    %SRY apply raw data correction if necessary
                    block = bsxfun(@times, block, rawDataCorrect);
                end

                isRefl = bIsReflected(ix);
                block(:,:,isRefl) = block(end:-1:1,:,isRefl);

                if ~isequal(selRange{1},':') || ~isequal(selRange{2},':')
                    block = block( selRange{1}, selRange{2}, : );    % a bit slow
                end

                [sortIdx, I] = sort( cIxToTarg(ix), 'ascend' );
                block = block(:,:,I);     % reorder according to sorted target indices

                % Mark duplicate indices with 1; we'll have to treat them special for proper averaging
                % Bonus: The very first storage can be made much faster, because it's in-place.
                %        Matlab urgently needs a "+=" operater, which makes "A(:,:,idx) = A(:,:,idx) + B"
                %        in-place and more readable.
                isDupe  = [ false, diff(sortIdx) == 0 ];

                idx1 = sortIdx(~isDupe);     % acquired once in this block
                idxN = sortIdx( isDupe);     % acquired multiple times

                count_ave(idx1) = count_ave(idx1) + 1;

                if isempty( idxN )
                    % no duplicates
                    if all( count_ave(idx1) == 1 )  % first acquisition of this line
                        out(:,:,idx1) = block;                              % fast
                    else
                        out(:,:,idx1) = out(:,:,idx1) + block;              % slow
                    end
                else
                    out(:,:,idx1) = out(:,:,idx1) + block(:,:,~isDupe);     % slower

                    block = block(:,:,isDupe);
                    for n = 1:numel(idxN)
                        out(:,:,idxN(n)) = out(:,:,idxN(n)) + block(:,:,n); % snail :-)
                        count_ave(idxN(n)) = count_ave(idxN(n)) + 1;
                    end
                end

                % At the first few iterations, evaluate the spent time-per-line and decide
                % what to do with the block size.
                if ~doLockblockSz
                    t = 1e6 * toc(s)/blockSz;   % micro seconds

                    if t <= 1.1 * tprev % allow 10% inaccuracy. Usually bigger == better
                        % New block size was faster. Go a step further.
                        blockSz = blockSz * 2;
                        blockInit = cat(3, blockInit, blockInit);
                    else
                        % regression; reset size and lock it
                        blockSz = max( blockSz/2, 1 );
                        blockInit = blockInit(:,:,1:blockSz);
                        doLockblockSz = true;
                    end
                    if bRegrid
                        rsTrj{3}  = single(1:blockSz);
                        trgTrj{3} = rsTrj{3};
                    end
                    tprev = t;
                end

                blockCtr = 0;
                block = blockInit;  % reset to garbage
            end

            if isBrokenRead
                this.isBrokenFile = true;
                break
            end
        end

        fclose(fid);

        % proper scaling (we don't want to sum our data but average it)
        % For large "out" bsxfun(@rdivide,out,count_ave) is incredibly faster than
        % bsxfun(@times,out,count_ave)!
        % @rdivide is also running in parallel, while @times is not. :-/
        if any( reshape(count_ave,[],1) > 1 )
            clearvars -except  out  count_ave  outSize
            count_ave = max( 1, count_ave );
            out       = bsxfun( @rdivide, out, count_ave);
        end

        out = reshape(out,outSize);
    end % of readData

    function setDefaultFlags(this)
        % method to set flags to default values
        this.arg.removeOS            = false;
        this.arg.rampSampRegrid      = false;
        this.arg.doAverage           = false;
        this.arg.averageReps         = false;
        this.arg.averageSets         = false;
        this.arg.ignoreSeg           = false;
        this.arg.doRawDataCorrect    = false;
        this.flagAverageDim          = false(1, 16);

        if strcmp(this.dataType,'image') || strcmp(this.dataType,'phasecor') || strcmp(this.dataType,'phasestab')
            this.arg.skipToFirstLine = false;
        else
            this.arg.skipToFirstLine = true;
        end
        if ~isfield(this.arg,'rawDataCorrectionFactors')
            this.arg.rawDataCorrectionFactors = [];
        end
    end

    function dummy = resetFlags(this)
        % method to reset flags to default values
        this.flagRemoveOS            = false;
        this.flagRampSampRegrid      = false;
        this.flagDoRawDataCorrect    = false;
        this.flagAverageDim          = false(1,16);

        if strcmp(this.dataType,'image') || strcmp(this.dataType,'phasecor') || strcmp(this.dataType,'phasestab')
            this.arg.skipToFirstLine = false;
        else
            this.arg.skipToFirstLine = true;
        end
        dummy = [];
    end

    function versiontime = get.readerVersion(~)
        % returns utc-unixtime of last commit (from file precommit-unixtime)
        p = fileparts(mfilename('fullpath'));
        fid = fopen(fullfile(p, 'precommit_unixtime'));
        versiontime = uint64(str2double(fgetl(fid)));
        fclose(fid);
    end

    function out = get.dataSize(this)
        out = this.fullSize;

        if this.arg.removeOS
            ix = ismember(this.dataDims, 'Col');
            out(ix) = this.NCol/2;
        end

        if this.flagAverageDim(1) || this.flagAverageDim(2)
            warning('averaging in col and cha dim not supported, resetting flag');
            this.flagAverageDim(1:2) = false;
        end

        out(this.flagAverageDim) = 1;
    end

    function out = get.sqzDims(this)
        out = this.dataDims(this.dataSize>1);
    end

    function out = get.sqzSize(this)
        out = this.dataSize(this.dataSize>1);
    end

    function set.flagRemoveOS(this,val)
        % set method for removeOS
        this.arg.removeOS = logical(val);
    end

    function out = get.flagRemoveOS(this)
        out = this.arg.removeOS;
    end

    function set.flagDoAverage(this,val)
        ix = ismember(this.dataDims, 'Ave');
        this.flagAverageDim(ix) = val;
    end

    function out = get.flagDoAverage(this)
        ix = ismember(this.dataDims, 'Ave');
        out = this.flagAverageDim(ix);
    end

    function set.flagAverageReps(this,val)
        ix = ismember(this.dataDims, 'Rep');
        this.flagAverageDim(ix) = val;
    end

    function out = get.flagAverageReps(this)
        ix = ismember(this.dataDims, 'Rep');
        out = this.flagAverageDim(ix);
    end

    function set.flagAverageSets(this,val)
        ix = ismember(this.dataDims, 'Set');
        this.flagAverageDim(ix) = val;
    end

    function out = get.flagAverageSets(this)
        ix = ismember(this.dataDims, 'Set');
        out = this.flagAverageDim(ix);
    end

    function set.flagIgnoreSeg(this,val)
        ix = ismember(this.dataDims, 'Seg');
        this.flagAverageDim(ix) = val;
    end

    function out = get.flagIgnoreSeg(this)
        ix = ismember(this.dataDims, 'Seg');
        out = this.flagAverageDim(ix);
    end

    function set.flagSkipToFirstLine(this,val)
        val = logical(val);
        if val ~= this.arg.skipToFirstLine
            this.arg.skipToFirstLine = val;

            if this.arg.skipToFirstLine
                this.skipLin = min(this.Lin)-1;
                this.skipPar = min(this.Par)-1;
            else
                this.skipLin = 0;
                this.skipPar = 0;
            end
            NLinAlloc = max(1, this.NLin - this.skipLin);
            NParAlloc = max(1, this.NPar - this.skipPar);
            this.fullSize(3:4) = [NLinAlloc NParAlloc];
        end
    end

    function out = get.flagSkipToFirstLine(this)
        out = this.arg.skipToFirstLine;
    end

    function out = get.flagRampSampRegrid(this)
        out = this.arg.rampSampRegrid;
    end

    function set.flagRampSampRegrid(this, val)
        val = logical(val);
        if (val == true && isempty(this.rampSampTrj))
            error('No trajectory for regridding available');
        end
        this.arg.rampSampRegrid = val;
    end

    %SRY: accessor methods for raw data correction
    function out = get.flagDoRawDataCorrect(this)
        out = this.arg.doRawDataCorrect;
    end

    function set.flagDoRawDataCorrect(this, val)
        val = logical(val);
        if (val == true && strcmp(this.softwareVersion, 'vd'))
            error('raw data correction for VD not supported/tested yet');
        end

        this.arg.doRawDataCorrect = val;
    end

    function out = get.RawDataCorrectionFactors(this)
        out = this.arg.rawDataCorrectionFactors;
    end

    function set.RawDataCorrectionFactors(this, val)
        %this may not work if trying to set the factors before NCha has
        %a meaningful value (ie before calling clean)
        if (~isrow(val) || length(val) ~= this.NCha)
            error('RawDataCorrectionFactors must be a 1xNCha row vector');
        end
        this.arg.rawDataCorrectionFactors = val;
    end

end

methods (Access='protected')
    % helper functions

    function fid = fileopen(this)
        % look out for unlikely event that someone is switching between
        % windows and unix systems:
        [path,name,ext] = fileparts(this.filename);
        this.filename   = fullfile(path,[name ext]);

        % test access
        if numel(dir(this.filename))==0
            % update path when file of same name can be found in current
            % working dir. -- otherwise throw error
            [oldpath,name,ext] = fileparts(this.filename);
            newloc = fullfile(pwd,[name ext]);
            if numel(dir(newloc))==1
                fprintf('Warning: File location updated from "%s" to current working directory.\n',oldpath);
                this.filename = newloc;
            else
                error(['File "' this.filename '" not found.']);
            end
        end
        fid = fopen(this.filename);
    end

    function [selRange,selRangeSz,outSize] = calcRange(this,S)

        switch S.type
            case '()'
                bSqueeze = false;
            case '{}'
                bSqueeze = true;
        end

        selRange = num2cell(ones(1,numel(this.dataSize)));
        outSize  = ones(1,numel(this.dataSize));

        if ( isempty(S.subs) || strcmpi(S.subs(1),'') )
            % obj(): shortcut to select all data
            % unfortunately, matlab does not allow the statement
            % obj{}, so we can't use it...
            % alternative: obj{''} (obj('') also works)
            for k=1:numel(this.dataSize)
                selRange{k}   = 1:this.dataSize(k);
            end
            if ~bSqueeze
                outSize = this.dataSize;
            else
                outSize = this.sqzSize;
            end
        else
            for k=1:numel(S.subs)
                if ~bSqueeze
                    cDim = k; % nothing to do
                else
                    % we need to rearrange selRange from squeezed
                    % to original order
                    cDim = find(strcmp(this.dataDims,this.sqzDims{k}) == 1);
                end
                if strcmp(S.subs{k},':')
                    if k<numel(S.subs)
                        selRange  {cDim} = 1:this.dataSize(cDim);
                    else % all later dimensions selected and 'vectorized'!
                        for l=cDim:numel(this.dataSize)
                            selRange{l} = 1:this.dataSize(l);
                        end
                        outSize(k) = prod(double(this.dataSize(cDim:end)));
                        break; % jump out ouf for-loop
                    end
                elseif isnumeric(S.subs{k})
                    selRange{cDim} = single(S.subs{k});
                else
                    error('unknown string in brackets (e.g. 1:end does not work here)');
                end
                outSize(k) = numel(selRange{cDim});
            end
            for k=1:numel(selRange)
                if max(selRange{k}) > this.dataSize(k)
                    error('selection out of range');
                end
            end
        end

        selRangeSz = ones(1,numel(this.dataSize));
        for k=1:numel(selRange)
            selRangeSz(k) = numel(selRange{k});
        end

        % now select all indices for the dims that are averaged
        for k = find(this.flagAverageDim)
            selRange{k} = 1:this.fullSize(k);
        end
    end

    function [ixToRaw, ixToTarget] = calcIndices(this)
        % calculate indices to target & source(raw)
        LinIx     = this.Lin - this.skipLin;
        ParIx     = this.Par - this.skipPar;
        sz = this.fullSize(3:end); % extra variable needed for octave compatibility
        ixToTarget = this.sub2ind_double(sz,...
            LinIx, ParIx, this.Sli, this.Ave, this.Phs, this.Eco,...
            this.Rep, this.Set, this.Seg, this.Ida, this.Idb,...
            this.Idc, this.Idd, this.Ide);

        % now calc. inverse index (page table: virtual to physical addresses)
        % indices of lines that are not measured are zero
        ixToRaw = zeros(1,prod(this.fullSize(3:end)),'double');

        ixToRaw(ixToTarget) = 1:numel(ixToTarget);
    end
end

methods (Static)
    % helper functions, accessible from outside via <classname>.function()
    % without an instance of the class.

    function ro_shift = calcOffcenterShiftRO(slicedata)
        % calculate ro offcenter shift from mdh's slicedata

        % slice position
        pos = slicedata(1:3);

        %quaternion
        a = slicedata(5);
        b = slicedata(6);
        c = slicedata(7);
        d = slicedata(4);

        read_dir = zeros(3,1);
        read_dir(1) = 2 * (a * b - c * d);
        read_dir(2) = 1 - 2 * (a^2 + c^2);
        read_dir(3) = 2 * (b * c + a * d);

        ro_shift = dot(pos, read_dir);
    end

    function ndx = sub2ind_double(sz, varargin)
        %SUB2IND_double Linear index from multiple subscripts.
        %   Works like sub2ind but always returns double
        %   also slightly faster, but no checks
        %========================================
        sz  = double(sz);
        ndx = double(varargin{end}) - 1;
        for i = length(sz)-1:-1:1
            ix  = double(varargin{i});
            ndx = sz(i)*ndx + ix-1;
        end
        ndx = ndx + 1;
    end % sub2ind_double

    function N = cast2MinimalUint( N )
        Nmax = max( reshape(N,[],1) );
        Nmin = min( reshape(N,[],1) );
        if Nmin < 0 || Nmax > intmax('uint64')
            return
        end

        if Nmax > intmax('uint32')
            idxClass = 'uint64';
        elseif Nmax > intmax('uint16')
            idxClass = 'uint32';
        else
            idxClass = 'uint16';
        end

        N = cast( N, idxClass );
    end % cast2MinimalUint()

end % of methods (Static)

end % classdef

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function mappedMedia = mapMediumData2AGORA(strainGrowth,inputMedia)
% This function extracts experimentally determined growth of species on
% multiple media and converts the in silico minimal medium the strain can
% grow on into an input fit for DEMETER. The input data was retrieved from
% Tramontano et al., Nat Microbiol 2019 (PMID: 29556107).
%
% USAGE
%       mappedMedia = mapMediumData2AGORA(strainGrowth,inputMedia)
%
% INPUTS
% strainGrowth          Growth of strains on the different media reported
%                       by Tramontano et al.
% inputMedia            Growth media tested by Tramontano et al. converted
%                       into a computationally readable format
%
% OUTPUT
% mappedMedia           Data from Tramontano et al. converted into a format
%                       that can be added to GrowthRequirementsTable
%
% .. Author:
%       - Almut Heinken, 03/2020

% delete organisms that did not grow on any medium
rmStrain=[];
cnt=1;
for i=2:size(strainGrowth,1)
    if sum(str2double(strainGrowth(i,2:end)))==0
        rmStrain(cnt,1)=i;
        cnt=cnt+1;
    end
end
strainGrowth(rmStrain,:)=[];

mappedMedia(1,:) = inputMedia(1,:);
for i = 2:size(strainGrowth,1)
    mappedMedia{i,1} = strainGrowth{i,1};
    % find the media it can grow on
    canGrow=strainGrowth(1,find(strcmp(strainGrowth(i,:),'1')));
    inputMediaStrain=inputMedia;
    [C,IA]=setdiff(inputMediaStrain(:,1),canGrow);
    inputMediaStrain(IA(2:end),:)=[];
    for j=2:size(inputMediaStrain,2)
        if any(strcmp(inputMediaStrain(2:end,j),'-1'))
            mappedMedia{i,j} = '-1';
        else
            mappedMedia{i,j} = '0';
        end
    end
end

end


text:
%SUBSREF Subscript reference overload of mapping
%
% This routine enables constructions like DATA = W.DATA,
% which is similar to DATA = GETDATA(W).
%
% In addition V = W(I,J) is supported for affine transformations.
% It is again an affine mapping using the [I,J] block of the
% rotation matrix and the elements J of the support vector.
%
% For arbitrary mappings just V = W(:,J) is defined by output
% selection: A*V returns just the features J of A*W.

author: marianux

repo: ecg-kit

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


text:
function type = year_to_type_hebrew ( y )

%*****************************************************************************80
%
%% YEAR_TO_TYPE_HEBREW returns the type of a Hebrew year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 332.
%
%  Parameters:
%
%    Input, integer Y, the Hebrew year.
%    Nonpositive years are illegal input.
%
%    Output, integer TYPE, the year type.
%    1, Common, Deficient, 12 months, 353 days;
%    2, Common, Regular, 12 months, 354 days;
%    3, Common, Abundant, 12 months, 355 days;
%    4, Embolismic, Deficient, 13 months, 383 days;
%    5, Embolismic, Regular, 13 months, 384 days;
%    6, Embolismic, Abundant, 13 months, 385 days.
%
  if ( y <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YEAR_TO_TYPE_HEBREW - Fatal error!\n' );
    fprintf ( 1, '  Illegal input Y = 0.\n' );
    error ( 'YEAR_TO_TYPE_HEBREW - Fatal error!' );
  end

  jed = new_year_to_jed_hebrew ( y );

  jed2 = new_year_to_jed_hebrew ( y + 1 );

  year_length = round ( jed2 - jed );

       if ( year_length == 353 )
    type = 1;
  elseif ( year_length == 354 )
    type = 2;
  elseif ( year_length == 355 )
    type = 3;
  elseif ( year_length == 383 )
    type = 4;
  elseif ( year_length == 384 )
    type = 5;
  elseif ( year_length == 385 )
    type = 6;
  else
    type = 0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YEAR_TO_TYPE_HEBREW - Fatal error!\n' );
    fprintf ( 1, '  Computed an illegal type = %d\n', type );
    error ( 'YEAR_TO_TYPE_HEBREW - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Perform the mean and variance normalization (MVN)
%       output = MVN(input,mean_d,var_d,weight)
%
% Input:
%   input: a feature vector or matrix to be processed. If the input is a
%          vector, it should be a column vector; if the input is a matrix,
%          it should be MxN, where M is the number of frames and N is the
%          feature index.
%   mean_d: the desired mean of the output. If the mean_d is not provided,
%           zero is used. If the input is a vector, mean_d should be a
%           scalar; if the input is a matrix, mean_d can be either a scalar
%           or a column vector.
%   var_d: the desired variance of the output. If the var_d is not provided,
%          one is used. If the input is a vector, var_d should be a
%          scalar; if the input is a matrix, var_d can be either a scalar
%          or a column vector.
%   weight: the weight of each frame during the normalization. The
%          function supports partial ownership, i.e. a frame may be owned by a
%          class partially. This is useful in multi-class MVN. If weight is not
%          provided, all frames will have weight 1, i.e. conventional MVN.
%
% Output:
%   output: the mean-variance-normalized features
%
% Author: Xiao Xiong, School of Computer Engineering, NTU, Singapore.
% Created: Jan, 2005
% Last Modified: 21 Dec, 2010
%
function [output] = MVN(input,mean_d,var_d, weight)
useWeight = 1;
if nargin < 4
    useWeight = 0;
end
if nargin < 3
    var_d = 1;  % default variance is 1
end
if nargin < 2
    mean_d = 0; % default mean is 0
end
[M, N] = size(input);   % Perform MVN for each column individually

if useWeight
    [meanX,varX] = findMeanVarainceWeighted(input, weight, 0);      % find the mean and variance of input weighted by their weights
    output = input - repmat(meanX',M,1);
    gain = diag(sqrt(var_d./varX));
    output = output * gain;
    output = output + repmat(mean_d',M,1);
else
    output = CVN(CMN(input,mean_d),var_d);
end

author: singaxiong

repo: SignalGraph


text:
function draw_images(data,wait,lim,dispsize)
% Draw data of image array to figure. Reduce size or scale image for resolutions above 800x600
%
% * Input parameters :
%    data [uint16(,,)]       image data array (height,width,count)
%    wait                    time to wait between images (optional) default: 0.01)
%    lim  [low up]           upper and lower limits to display (optional)
%                            default: [0 m], where m is maximal value calculated from first image
%    dispsize [width height] size of top-left corner to display (optional)
%                            default: full image scaled for bigger sizes
%
% * Output parameters :

 if(nargin<1)
  error('Wrong number of input arguments. Need data array');
 end

 if(nargin<2)
  wait=0.010;
 end

 if(nargin<3)
  ulim=max(max(data(10:end-10,10:end-10,1)));
  lim=[0 ulim];
 end

 disp(['waittime set to ',num2str(wait)]);
 disp(['limits set to [0 ',num2str(lim(2)),']']);

 [h,w,count]=size(data);

 if(nargin==4)
  if(dispsize(1)<w)
   w=dispsize(1);
  end
  if(dispsize(2)<h)
   h=dispsize(2);
  end
 end

 imah=draw_image(data(1:h,1:w,1),lim);
 for ima_nr=1:count
  set(imah,'CData',data(1:h,1:w,ima_nr),'CDataMapping','scaled');
  if(wait==0)
   disp('Press any key to proceed')
   pause
  else
   pause(wait);
  end
 end
 disp('Press "Enter" to close window and proceed');
 pause();
 close();
 pause(1);
 commandwindow;
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
classdef TestDrawDetectedDiamonds
    %TestDrawDetectedDiamonds

    methods (Static)
        function test_1
            [img, diamond] = get_image_markers();
            [corners, ids] = cv.detectMarkers(img, diamond.dict);
            [diamondCorners, diamondIds] = cv.detectCharucoDiamond(img, ...
                corners, ids, diamond.slen / diamond.mlen);
            out = cv.drawDetectedDiamonds(img, diamondCorners, ...
                'IDs',diamondIds, 'BorderColor',[0 0 255]);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.drawDetectedDiamonds();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function [img, d] = get_image_markers()
    % 4 markers in a diamond
    d = struct();
    d.dict = '6x6_50';
    d.ids = randi([0 49], [1 4]);
    d.slen = 80;
    d.mlen = 50;
    img = cv.drawCharucoDiamond(d.dict, d.ids, d.slen, d.mlen);
    img = repmat(img, [1 1 3]);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function index = tileIndex3d(tile)
% Return the index of a 2-by-2-by-2 3D binary configuration tile.
%
%   INDEX = tileIndex3d(TILE)
%   Compute the index of a tile, given as a 2x2x2 binary image, by adding
%   powers of two mutliplied by values of the tile. The result is comprised
%   between 0 and 255.
%   Iterate on x, y, then z directions, corresponding to directions 2, 1,
%   and 3.
%
%   Example
%   tile = zeros([2 2 2]);
%   tile(1,1:2,1) = 1;
%   tileIndex3d(tile)
%       3
%
%   See also
%     createTile3d, tileIndex
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2009-05-25,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

tile = tile>0;
index = ...
    tile(1,1,1) + 2*tile(1,2,1) + 4*tile(2,1,1) + 8*tile(2,2,1) + ...
    16*tile(1,1,2) + 32*tile(1,2,2) + 64*tile(2,1,2) + 128*tile(2,2,2);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [inv_A, r, PR, PL] = invertProjection(A, epsilon)
% Inverts a general matrix A using the pseudoinverse
%
% USAGE:
%
%    [inv_A, r, PR, PL] = invertProjection(A, epsilon)
% INPUTS:
%    A:          general matrix
%    epsilon:    default = 1e-10
%
% OUTPUTS:
%    inv_A:      the pseudoinverse of `A`
%    r:          the rank of `A`
%    PR:        the projection matrix onto the `range(A)`
%    PL:        the projection matrix onto the `null(A')`

if nargin < 2
    epsilon = 1e-10;
end
[m, n] = size(A);

if 1
    %[U, S, V] = svd(A); % not working, uncommented line below - Lemmer
    %[U, S, V] = svds(A,min(size(A))); % Bugfix due to svd convergence problems
    [U, S, V] = svd(full(A),'econ'); %from Michael Saunders code
    r = sum(sum(abs(S) > epsilon));
    inv_S = diag(1 ./ S(abs(S) > epsilon));
    inv_A = V(:, 1:r) * inv_S(1:r, 1:r) * U(:, 1:r)';
    PR = U(:, 1:r)       * U(:, 1:r)';
    PL = U(:, (r+1):end) * U(:, (r+1):end)';

else
    %Michael Saunders code -TODO integrate this properly
    [U1,D1,V1,r] = subspaceSVD(A);
    PR=U1*U1';%projection matrix onto the range(A)
    PL=eye(m) - U1*U1';%projection matrix onto null(A')
    inv_A=pinv(A,1e-12);
end
    % Michael Saunders code
    % [U1,D1,V1,r] = subspaceSVD(A);
    % PR=U1*U1';%projection matrix onto the range(A)
    % PL=eye(m) - U1*U1';%projection matrix onto the null(A')
    % inv_A=pinv(A,1e-12);


text:
function F = uminus( F )
%- Unary minus of a CHEBFUN2V
%   -F returns the CHEBFUN2V negated componentwise.
%
%   UMINUS(F) is called by the syntax -F.

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

% Empty check:
if (isempty( F ) )
    return
end

% Take uminus of each component:
F.components = cellfun( @uminus, F.components, 'UniformOutput', false );

end