text:
function [input] = ft_apply_montage(input, montage, varargin)

% FT_APPLY_MONTAGE changes the montage (i.e. linear combination) of a set of
% electrode or gradiometer channels. A montage can be used for EEG rereferencing, MEG
% synthetic gradients, MEG planar gradients or unmixing using ICA. This function not
% only applies the montage to the EEG or MEG data, but also applies the montage to
% the input EEG or MEG sensor array, which can subsequently be used for forward
% computation and source reconstruction of the data.
%
% Use as
%   [sens]    = ft_apply_montage(sens,     montage,  ...)
%   [data]    = ft_apply_montage(data,     montage,  ...)
%   [freq]    = ft_apply_montage(freq,     montage,  ...)
%   [montage] = ft_apply_montage(montage1, montage2, ...)
%
% A montage is specified as a structure with the fields
%   montage.tra      = MxN matrix
%   montage.labelold = Nx1 cell-array
%   montage.labelnew = Mx1 cell-array
%
% As an example, a bipolar montage could look like this
%   bipolar.labelold  = {'1',   '2',   '3',   '4'}
%   bipolar.labelnew  = {'1-2', '2-3', '3-4'}
%   bipolar.tra       = [
%     +1 -1  0  0
%      0 +1 -1  0
%      0  0 +1 -1
%   ];
%
% The montage can optionally also specify the channel type and unit of the input
% and output data with
%   montage.chantypeold = Nx1 cell-array
%   montage.chantypenew = Mx1 cell-array
%   montage.chanunitold = Nx1 cell-array
%   montage.chanunitnew = Mx1 cell-array
%
% Additional options should be specified in key-value pairs and can be
%   'keepunused'    = string, 'yes' or 'no' (default = 'no')
%   'inverse'       = string, 'yes' or 'no' (default = 'no')
%   'balancename'   = string, name of the montage (default = '')
%   'feedback'      = string, see FT_PROGRESS (default = 'text')
%   'warning'       = boolean, whether to show warnings (default = true)
%
% If the first input is a montage, then the second input montage will be
% applied to the first. In effect, the output montage will first do
% montage1, then montage2.
%
% See also FT_READ_SENS, FT_DATATYPE_SENS

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

if iscell(input)
  % this represents combined EEG, ECoG and/or MEG
  for i=1:numel(input)
    input{i} = ft_apply_montage(input{i}, montage, varargin{:});
  end
  return
end

% use "old/new" instead of "org/new"
montage = fixoldorg(montage);
input   = fixoldorg(input); % the input might be a montage or a sensor array

% get optional input arguments
keepunused    = ft_getopt(varargin, 'keepunused',  'no');
inverse       = ft_getopt(varargin, 'inverse',     'no');
feedback      = ft_getopt(varargin, 'feedback',    'text');
showwarning   = ft_getopt(varargin, 'warning',     true);
bname         = ft_getopt(varargin, 'balancename', '');

if istrue(showwarning)
  warningfun = @warning;
else
  warningfun = @nowarning;
end

% ensure that the input montage is correct, see https://github.com/fieldtrip/fieldtrip/issues/1718
assert(length(unique(montage.labelold))==length(montage.labelold), 'the montage is invalid');
assert(length(unique(montage.labelnew))==length(montage.labelnew), 'the montage is invalid');

% these are optional, at the end we will clean up the output in case they did not exist
if ~istrue(inverse)
  haschantype = (isfield(input, 'chantype') || isfield(input, 'chantypenew')) && all(isfield(montage, {'chantypeold', 'chantypenew'}));
  haschanunit = (isfield(input, 'chanunit') || isfield(input, 'chanunitnew')) && all(isfield(montage, {'chanunitold', 'chanunitnew'}));
else
  haschantype = (isfield(input, 'chantype') || isfield(input, 'chantypeold')) && all(isfield(montage, {'chantypeold', 'chantypenew'}));
  haschanunit = (isfield(input, 'chanunit') || isfield(input, 'chanunitold')) && all(isfield(montage, {'chanunitold', 'chanunitnew'}));
end

% make sure they always exist to facilitate the remainder of the code
if ~isfield(montage, 'chantypeold')
  montage.chantypeold = repmat({'unknown'}, size(montage.labelold));
  if isfield(input, 'chantype') && ~istrue(inverse)
    ft_warning('copying input chantype to montage');
    [sel1, sel2] = match_str(montage.labelold, input.label);
    montage.chantypeold(sel1) = input.chantype(sel2);
  end
end

if ~isfield(montage, 'chantypenew')
  montage.chantypenew = repmat({'unknown'}, size(montage.labelnew));
  if isfield(input, 'chantype') && istrue(inverse)
    ft_warning('copying input chantype to montage');
    [sel1, sel2] = match_str(montage.labelnew, input.label);
    montage.chantypenew(sel1) = input.chantype(sel2);
  end
end

if ~isfield(montage, 'chanunitold')
  montage.chanunitold = repmat({'unknown'}, size(montage.labelold));
  if isfield(input, 'chanunit') && ~istrue(inverse)
    ft_warning('copying input chanunit to montage');
    [sel1, sel2] = match_str(montage.labelold, input.label);
    montage.chanunitold(sel1) = input.chanunit(sel2);
  end
end

if ~isfield(montage, 'chanunitnew')
  montage.chanunitnew = repmat({'unknown'}, size(montage.labelnew));
  if isfield(input, 'chanunit') && istrue(inverse)
    ft_warning('copying input chanunit to montage');
    [sel1, sel2] = match_str(montage.labelnew, input.label);
    montage.chanunitnew(sel1) = input.chanunit(sel2);
  end
end

if ~isfield(input, 'label') && isfield(input, 'labelnew')
  % the input data structure is also a montage
  inputlabel = input.labelnew;
  if isfield(input, 'chantypenew')
    inputchantype = input.chantypenew;
  else
    inputchantype = repmat({'unknown'}, size(input.labelnew));
  end
  if isfield(input, 'chanunitnew')
    inputchanunit = input.chanunitnew;
  else
    inputchanunit = repmat({'unknown'}, size(input.labelnew));
  end
else
  % the input should describe the channel labels, and optionally the type and unit
  inputlabel = input.label;
  if isfield(input, 'chantype')
    inputchantype = input.chantype;
  else
    inputchantype = repmat({'unknown'}, size(input.label));
  end
  if isfield(input, 'chanunit')
    inputchanunit = input.chanunit;
  else
    inputchanunit = repmat({'unknown'}, size(input.label));
  end
end

% check the consistency of the montage
if ~iscell(montage.labelold) || ~iscell(montage.labelnew)
  ft_error('montage labels need to be specified in cell-arrays');
end

% check the consistency of the montage
if ~all(isfield(montage, {'tra', 'labelold', 'labelnew'}))
  ft_error('the second input argument does not correspond to a montage');
end

% check the consistency of the montage
if size(montage.tra,1)~=length(montage.labelnew)
  ft_error('the number of channels in the montage is inconsistent');
elseif size(montage.tra,2)~=length(montage.labelold)
  ft_error('the number of channels in the montage is inconsistent');
end

% use a default unit transfer from sensors to channels if not otherwise specified
if ~isfield(input, 'tra') && isfield(input, 'label')
  if     isfield(input, 'elecpos') && length(input.label)==size(input.elecpos, 1)
    nchan = length(input.label);
    input.tra = eye(nchan);
  elseif isfield(input, 'coilpos') && length(input.label)==size(input.coilpos, 1)
    nchan = length(input.label);
    input.tra = eye(nchan);
  elseif isfield(input, 'chanpos') && length(input.label)==size(input.chanpos, 1)
    nchan = length(input.label);
    input.tra = eye(nchan);
  end
end

if istrue(inverse)
  % swap the role of the old and new channels
  tmp.labelnew    = montage.labelold;
  tmp.labelold    = montage.labelnew;
  tmp.chantypenew = montage.chantypeold;
  tmp.chantypeold = montage.chantypenew;
  tmp.chanunitnew = montage.chanunitold;
  tmp.chanunitold = montage.chanunitnew;
  % apply the inverse montage, this can be used to undo a previously applied montage
  tmp.tra = full(montage.tra);
  if rank(tmp.tra) < length(tmp.tra)
    warningfun('the linear projection for the montage is not full-rank, the resulting data will have reduced dimensionality');
    tmp.tra = pinv(tmp.tra);
  else
    tmp.tra = inv(tmp.tra);
  end
  montage = tmp;
end

% keep only the columns that are not empty
selcol = ~all(montage.tra==0, 1);
if istrue(keepunused)
  for i=find(selcol==false)
    % don't remove the column if it corresponds to one of the output channels
    selcol(i) = any(strcmp(montage.labelnew, montage.labelold{i}));
  end
end
montage.tra         = montage.tra(:,selcol);
montage.labelold    = montage.labelold(selcol);
montage.chantypeold = montage.chantypeold(selcol);
montage.chanunitold = montage.chanunitold(selcol);
clear selcol

% keep only the channels that are present in the input data
remove = setdiff(montage.labelold, intersect(montage.labelold, inputlabel));
selcol = match_str(montage.labelold, remove);
% we cannot just remove the colums, all rows that depend on it should also be removed
selrow = false(length(montage.labelnew),1);
for i=1:length(selcol)
  selrow = selrow & (montage.tra(:,selcol(i))~=0);
end
% convert from indices to logical vector
selcol = indx2logical(selcol, length(montage.labelold));
% remove rows and columns
montage.labelold    = montage.labelold(~selcol);
montage.labelnew    = montage.labelnew(~selrow);
montage.chantypeold = montage.chantypeold(~selcol);
montage.chantypenew = montage.chantypenew(~selrow);
montage.chanunitold = montage.chanunitold(~selcol);
montage.chanunitnew = montage.chanunitnew(~selrow);
montage.tra         = montage.tra(~selrow, ~selcol);
clear remove selcol selrow

% add columns for channels present in the input data but that are not specified in the montage, stick to the original order in the data
[dum, ix]   = setdiff(inputlabel, montage.labelold);
addlabel    = inputlabel(sort(ix));
addchantype = inputchantype(sort(ix));
addchanunit = inputchanunit(sort(ix));
m = size(montage.tra,1);
n = size(montage.tra,2);
k = length(addlabel);

% % check for NaNs in unused channels; these will be mixed in with the rest
% % of the channels and result in NaNs in the output even when multiplied
% % with zeros or identity
% if k > 0 && isfield(input, 'trial') % check for raw data now only
%   cfg = [];
%   cfg.channel = addlabel;
%   cfg.showcallinfo = showcallinfo;
%   data_unused = ft_selectdata(cfg, input);
%   % use an anonymous function to test for the presence of NaNs in the input data
%   hasnan = @(x) any(isnan(x(:)));
%   if any(cellfun(hasnan, data_unused.trial))
%     ft_error('FieldTrip:NaNsinInputData', ['Your input data contains NaNs in channels that are unused '...
%       'in the supplied montage. This would result in undesired NaNs in the '...
%       'output data. Please remove these channels from the input data (using '...
%       'ft_selectdata) before attempting to apply the montage.']);
%   end
% end

if istrue(keepunused)
  % add the channels that are not rereferenced to the input and output of the montage
  montage.tra((m+(1:k)),(n+(1:k))) = eye(k);
  montage.labelold    = cat(1, montage.labelold(:), addlabel(:));
  montage.labelnew    = cat(1, montage.labelnew(:), addlabel(:));
  montage.chantypeold = cat(1, montage.chantypeold(:), addchantype(:));
  montage.chantypenew = cat(1, montage.chantypenew(:), addchantype(:));
  montage.chanunitold = cat(1, montage.chanunitold(:), addchanunit(:));
  montage.chanunitnew = cat(1, montage.chanunitnew(:), addchanunit(:));
else
  % add the channels that are not rereferenced to the input of the montage only
  montage.tra(:,(n+(1:k))) = zeros(m,k);
  montage.labelold    = cat(1, montage.labelold(:), addlabel(:));
  montage.chantypeold = cat(1, montage.chantypeold(:), addchantype(:));
  montage.chanunitold = cat(1, montage.chanunitold(:), addchanunit(:));
end
clear addlabel addchantype addchanunit m n k

% determine whether all channels are unique
m = size(montage.tra,1);
n = size(montage.tra,2);
if length(unique(montage.labelnew))~=m
  ft_error('not all output channels of the montage are unique');
end
if length(unique(montage.labelold))~=n
  ft_error('not all input channels of the montage are unique');
end

% determine whether all channels that have to be rereferenced are available
if length(intersect(inputlabel, montage.labelold))~=length(montage.labelold)
  ft_error('not all channels that are required in the montage are available in the data');
end

% reorder the columns of the montage matrix
[selinput, selmontage] = match_str(inputlabel, montage.labelold);
montage.tra            = montage.tra(:,selmontage);
montage.labelold       = montage.labelold(selmontage);
montage.chantypeold    = montage.chantypeold(selmontage);
montage.chanunitold    = montage.chanunitold(selmontage);

% ensure that the montage is double precision
montage.tra = double(montage.tra);

% making the tra matrix sparse will speed up subsequent multiplications, but should not result in sparse output data
% this only makes sense for matrices with a lot of zero elements; for dense matrices it is faster to keep it full
if size(montage.tra,1)>1 && nnz(montage.tra)/numel(montage.tra) < 0.3
  montage.tra = sparse(montage.tra);
else
  montage.tra = full(montage.tra);
end

% update the channel scaling if the input has different units than the montage expects
if isfield(input, 'chanunit') && ~isequal(input.chanunit, montage.chanunitold)
  scale = ft_scalingfactor(input.chanunit, montage.chanunitold);
  montage.tra = montage.tra * diag(scale);
  montage.chanunitold = input.chanunit;
elseif isfield(input, 'chanunitnew') && ~isequal(input.chanunitnew, montage.chanunitold)
  scale = ft_scalingfactor(input.chanunitnew, montage.chanunitold);
  montage.tra = montage.tra * diag(scale);
  montage.chanunitold = input.chanunitnew;
end

if isfield(input, 'chantype') && ~isequal(input.chantype, montage.chantypeold)
  ft_error('inconsistent chantype in data and montage');
elseif isfield(input, 'chantypenew') && ~isequal(input.chantypenew, montage.chantypeold)
  ft_error('inconsistent chantype in data and montage');
end

if isfield(input, 'labelold') && isfield(input, 'labelnew')
  inputtype = 'montage';
elseif isfield(input, 'tra')
  inputtype = 'sens';
elseif ft_datatype(input, 'raw')
  inputtype = 'raw';
elseif ft_datatype(input, 'timelock')
  inputtype = 'timelock';
elseif ft_datatype(input, 'freq') && isfield(input, 'fourierspctrm')
  inputtype = 'freq';
elseif ft_datatype(input, 'freq') && isfield(input, 'crsspctrm')
  inputtype = 'freq_crsspctrm';

  % attempt to convert to a chan-chan representation
  input     = ft_checkdata(input, 'cmbstyle', 'full');
else
  inputtype = 'unknown';
end

switch inputtype
  case 'montage'
    % apply the montage on top of the other montage
    if isa(input.tra, 'single')
      % sparse matrices and single precision do not match
      input.tra = full(montage.tra) * input.tra;
    else
      input.tra = montage.tra * input.tra;
    end
    input.labelnew    = montage.labelnew;
    input.chantypenew = montage.chantypenew;
    input.chanunitnew = montage.chanunitnew;

  case 'sens'
    % apply the montage to an electrode or gradiometer description
    sens = input;
    clear input

    % apply the montage to the input
    if isa(sens.tra, 'single')
      % sparse matrices and single precision do not match
      sens.tra = full(montage.tra) * sens.tra;
    else
      sens.tra = montage.tra * sens.tra;
    end

    % The montage operates on the coil weights in sens.tra, but the output channels can be different.
    % If possible, we want to keep the old channel positions and orientations.
    [sel1, sel2] = match_str(montage.labelnew, inputlabel);
    keepchans = isequal(sel1(:)', 1:numel(montage.labelnew));

    posweight = abs(montage.tra);
    posweight = diag(1./sum(posweight,2)) * posweight;

    if isfield(sens, 'chanpos')
      if keepchans
        sens.chanpos = sens.chanpos(sel2,:);
      else
        if ~isfield(sens, 'chanposold')
          % add a chanposold only if it is not there yet
          sens.chanposold  = sens.chanpos;
          %  also keep the old label, type and unit for reference
          sens.labelold    = inputlabel;
          sens.chantypeold = inputchantype;
          sens.chanunitold = inputchanunit;
        end
        % compute the new channel positions as a weighted sum of the old ones
        sens.chanpos = posweight * sens.chanpos;
      end
    end

    if isfield(sens, 'chanori')
      if keepchans
        sens.chanori = sens.chanori(sel2,:);
      else
        if ~isfield(sens, 'chanoriold')
          sens.chanoriold = sens.chanori;
        end
        % compute the new channel orientations as a weighted sum of the old ones
        sens.chanori = posweight * sens.chanori;
      end
    end

    sens.label    = montage.labelnew;
    sens.chantype = montage.chantypenew;
    sens.chanunit = montage.chanunitnew;

    % keep track of the order of the balancing and which one is the current one
    if istrue(inverse)
      if isfield(sens, 'balance')% && isfield(sens.balance, 'previous')
        if isfield(sens.balance, 'previous') && numel(sens.balance.previous)>=1
          sens.balance.current  = sens.balance.previous{1};
          sens.balance.previous = sens.balance.previous(2:end);
        elseif isfield(sens.balance, 'previous')
          sens.balance.current  = 'none';
          sens.balance          = rmfield(sens.balance, 'previous');
        else
          sens.balance.current  = 'none';
        end
      end

    elseif ~istrue(inverse) && ~isempty(bname)
      if isfield(sens, 'balance')
        % check whether a balancing montage with name bname already exist, and if so, how many
        mnt = fieldnames(sens.balance);
        sel = strmatch(bname, mnt);
        if numel(sel)==0
          % bname can stay the same
        elseif numel(sel)==1
          % the original should be renamed to 'bname1' and the new one should be 'bname2'
          sens.balance.([bname, '1']) = sens.balance.(bname);
          sens.balance                = rmfield(sens.balance, bname);
          if isfield(sens.balance, 'current') && strcmp(sens.balance.current, bname)
            sens.balance.current = [bname, '1'];
          end
          if isfield(sens.balance, 'previous')
            sel2 = strmatch(bname, sens.balance.previous);
            if ~isempty(sel2)
              sens.balance.previous{sel2} = [bname, '1'];
            end
          end
          bname = [bname, '2'];
        else
          bname = [bname, num2str(length(sel)+1)];
        end
      end

      if isfield(sens, 'balance') && isfield(sens.balance, 'current')
        if ~isfield(sens.balance, 'previous')
          sens.balance.previous = {};
        end
        sens.balance.previous = [{sens.balance.current} sens.balance.previous];
        sens.balance.current  = bname;
        sens.balance.(bname)  = montage;
      end
    end

    % rename the output variable
    input = sens;
    clear sens

  case 'raw'
    % apply the montage to the raw input data
    data = input;
    clear input

    % there are two challenges to deal with
    % 1) the input data can be single, and sparse(1) * single(0) fails
    % 2) the input data can contain nans, and 0*nan returns a nan
    % see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3035 and https://github.com/fieldtrip/fieldtrip/issues/2169

    Ntrials = numel(data.trial);
    ft_progress('init', feedback, 'processing trials');
    for i=1:Ntrials
      ft_progress(i/Ntrials, 'processing trial %d from %d\n', i, Ntrials);

      if     ~isa(data.trial{i}, 'single') && ~any(isnan(data.trial{i}(:)))
        data.trial{i} = montage.tra * data.trial{i};
      elseif ~isa(data.trial{i}, 'single') &&  any(isnan(data.trial{i}(:)))
        % do not multiply 0 in the montage with nan in the data
        tmp = zeros(size(montage.tra,1), size(data.trial{i}, 2));
        for j=1:size(montage.tra,1)
          sel = montage.tra(j,:)~=0;
          tmp(j,:) = montage.tra(j,sel) * data.trial{i}(sel,:);
        end
        data.trial{i} = tmp;
      elseif  isa(data.trial{i}, 'single') && ~any(isnan(data.trial{i}(:)))
        % sparse matrices and single precision do not match
        data.trial{i} = full(montage.tra) * data.trial{i};
      elseif  isa(data.trial{i}, 'single') &&  any(isnan(data.trial{i}(:)))
        % sparse matrices and single precision do not match
        % do not multiply 0 in the montage with nan in the data, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3035
        tmp = zeros(size(montage.tra,1), size(data.trial{i}, 2));
        for j=1:size(montage.tra,1)
          sel = montage.tra(j,:)~=0;
          tmp(j,:) = full(montage.tra(j,sel)) * data.trial{i}(sel,:);
        end
        data.trial{i} = tmp;
      end
    end % for Ntrials
    ft_progress('close');

    data.label    = montage.labelnew;
    data.chantype = montage.chantypenew;
    data.chanunit = montage.chanunitnew;

    % rename the output variable
    input = data;
    clear data

  case 'timelock'
    % apply the montage to averaged data
    timelock = input;
    clear input

    fn = {'avg', 'trial', 'individual', 'cov'};
    for i=1:numel(fn)
      if isfield(timelock, fn{i})
        switch getdimord(timelock, fn{i})
          case 'chan_time'
            timelock.(fn{i}) = montage.tra * timelock.(fn{i});
          case 'rpt_chan_time'
            siz    = getdimsiz(timelock, fn{i});
            nrpt   = siz(1);
            nchan  = siz(2);
            ntime  = siz(3);
            output = zeros(nrpt, size(montage.tra,1), ntime);
            for rptlop=1:nrpt
              output(rptlop,:,:) = montage.tra * reshape(timelock.(fn{i})(rptlop,:,:), [nchan ntime]);
            end
            timelock.(fn{i}) = output; % replace the original field
          case 'chan_chan'
            timelock.(fn{i}) = montage.tra * timelock.(fn{i}) * montage.tra';
          case 'rpt_chan_chan'
            siz    = getdimsiz(timelock, fn{i});
            nrpt   = siz(1);
            nchan  = siz(2);
            output = zeros(nrpt, size(montage.tra,1), size(montage.tra,1));
            for rptlop=1:nrpt
              output(rptlop,:,:) = montage.tra * reshape(timelock.(fn{i})(rptlop,:,:), [nchan nchan]) * montage.tra';
            end
            timelock.(fn{i}) = output; % replace the original field
          otherwise
            ft_error('unsupported dimord for %s', fn{i});
        end % switch
      end % if
    end % for

    timelock.label    = montage.labelnew;
    timelock.chantype = montage.chantypenew;
    timelock.chanunit = montage.chanunitnew;

    % rename the output variable
    input = timelock;
    clear timelock

  case 'freq'
    % apply the montage to the spectrally decomposed data
    freq = input;
    clear input

    switch getdimord(freq, 'fourierspctrm')
      case 'rpttap_chan_freq'
        siz    = [getdimsiz(freq, 'fourierspctrm') 1];
        nrpt   = siz(1);
        nchan  = siz(2);
        nfreq  = siz(3);
        output = zeros(nrpt, size(montage.tra,1), nfreq);
        for foilop=1:nfreq
          output(:,:,foilop) = freq.fourierspctrm(:,:,foilop) * montage.tra';
        end
        freq.fourierspctrm = output; % replace the original Fourier spectrum

      case 'rpttap_chan_freq_time'
        siz    = getdimsiz(freq, 'fourierspctrm');
        nrpt   = siz(1);
        nchan  = siz(2);
        nfreq  = siz(3);
        ntime  = siz(4);
        output = zeros(nrpt, size(montage.tra,1), nfreq, ntime);
        for foilop=1:nfreq
          for toilop = 1:ntime
            output(:,:,foilop,toilop) = freq.fourierspctrm(:,:,foilop,toilop) * montage.tra';
          end
        end
        freq.fourierspctrm = output; % replace the original Fourier spectrum

      otherwise
        ft_error('unsupported dimord for fourierspctrm');
    end % switch

    freq.label    = montage.labelnew;
    freq.chantype = montage.chantypenew;
    freq.chanunit = montage.chanunitnew;

    % rename the output variable
    input = freq;
    clear freq
  case 'freq_crsspctrm'
    % input freq data has a chan-chan crsspctrm, montage needs te be
    % applied to both ends

    freq = input;
    clear input
    if contains(getdimord(freq, 'crsspctrm'), 'rpt')
      % first dimension is rpt-like, so the square is dimensions 2 and 3
      siz   = [getdimsiz(freq, 'crsspctrm') 1];
      nrpt  = siz(1);
      nchan = siz(2);
      nrest = prod(siz(4:end));
      output = zeros([nrpt size(montage.tra,1).*[1 1] siz(4:end)]);
      for rptlop = 1:nrpt
        for restlop = 1:nrest
          output(rptlop,:,:,restlop) = montage.tra*reshape(freq.crsspctrm(rptlop,:,:,restlop),[nchan nchan])*montage.tra';
        end
      end
      freq.crsspctrm = output;
    else
      % the square is dimensions 1 and 2
      siz   = [getdimsiz(freq, 'crsspctrm') 1];
      nrest = prod(siz(4:end));
      output = zeros([size(montage.tra,1).*[1 1] siz(4:end)]);
      for restlop = 1:nrest
        output(:,:,restlop) = montage.tra*freq.crsspctrm(:,:,restlop)*montage.tra';
      end
      freq.crsspctrm = output;
    end

    freq.label    = montage.labelnew;
    freq.chantype = montage.chantypenew;
    freq.chanunit = montage.chanunitnew;

    % rename the output variable
    input = freq;
    clear freq
  otherwise
    ft_error('unrecognized input');
end % switch inputtype

% only retain the chantype and/or chanunit if they were present in the input
if ~haschantype
  input = removefields(input, {'chantype', 'chantypeold', 'chantypenew'});
end
if ~haschanunit
  input = removefields(input, {'chanunit', 'chanunitold', 'chanunitnew'});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HELPER FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = indx2logical(x, n)
y = false(1,n);
y(x) = true;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HELPER FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nowarning(varargin)
return


text:
function varargout=m2v(varargin)
%
% vol=m2v(node,face,Nxyz)
%  or
% vol=m2v(node,face,xi,yi,zi)
%
% shortcut for mesh2vol, rasterizing a teterahedral mesh to a volume using graphics
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input/output: please see details in the help for mesh2vol
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

[varargout{1:nargout}]=mesh2vol(varargin{:});

author: fangq


text:
function [ n_data_new, n, c ] = tau_values ( n_data )

%*****************************************************************************80
%
%% TAU_VALUES returns some values of the Tau function.
%
%  Discussion:
%
%    TAU(N) is the number of divisors of N, including 1 and N.
%
%  First values:
%
%     N   TAU(N)
%
%     1    1
%     2    2
%     3    2
%     4    3
%     5    2
%     6    4
%     7    2
%     8    4
%     9    3
%    10    4
%    11    2
%    12    6
%    13    2
%    14    4
%    15    4
%    16    5
%    17    2
%    18    6
%    19    2
%    20    6
%
%  Formula:
%
%    If the prime factorization of N is
%
%      N = P1**E1 * P2**E2 * ... * PM**EM,
%
%    then
%
%      TAU(N) = ( E1 + 1 ) * ( E2 + 1 ) * ... * ( EM + 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, integer N, the argument of the Tau function.
%
%    Output, integer C, the value of the Tau function.
%
  n_max = 20;
  c_vec = [ ...
    1,  2,  2,  3,  2,  4,  2,  4,  3,  4, ...
    2, 12, 12,  4, 18, 24,  2,  8, 14, 28 ];
  n_vec = [ ...
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, ...
     23,  72, 126, 226, 300, 480, 521, 610, 832, 960 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data_new);
    c = c_vec(n_data_new);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ALAP_demo_attenuation
% This demonstration illustrates context or state-dependent precision (i.e.
% attention), which is necessary to disambiguate between sensations
% caused exogenously and self-generated sensations. In brief, it is
% necessary to attend away from the sensory consequences of action to
% preclude sensory evidence overriding the prior beliefs that cause
% movement. This necessarily reduced the confidence in self-generated
% sensations and provides a simple (Bayes-optimal) explanation for sensory
% attenuation - in terms of the attention of sensory precision. We
% illustrate this in the setting of the force matching illusion and go on
% to show that increasing the conviction in (precision of) prior beliefs
% abolishes sensory attenuation at the expense of false (delusional)
% posterior beliefs about antagonistic external forces.
%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ALAP_demo_attenuation.m 4826 2012-08-03 16:45:09Z karl $

% process (G) and model (M)
%==========================================================================

% set dimensions for generalised coordinates
%--------------------------------------------------------------------------
G(1).E.d        = 2;                   % approximation order
G(1).E.n        = 4;                   % embedding order

M(1).E.d        = 2;                   % approximation order
M(1).E.n        = 4;                   % embedding order
M(1).E.s        = 1/2;                 % embedding order

M(1).E.method.x = 1;                   % state-dependent noise
M(1).E.method.v = 1;                   % state-dependent noise
M(1).E.method.h = 0;                   % suppress optimisation
M(1).E.method.g = 0;                   % suppress optimisation

G(1).f  = inline('tanh(a) - x/4','x','v','a','P');
G(1).g  = inline('[x; v + x]','x','v','a','P');
G(1).x  = 0;                           % hidden state
G(1).v  = [0; 0];                      % hidden cause (sensory data)
G(1).V  = exp(8);                      % precision (noise)
G(1).W  = exp(8);                      % precision (states)
G(1).U  = [exp(0) 0];                  % precision (action)

% level 2; causes
%--------------------------------------------------------------------------
G(2).v  = 0;                           % exogenous  cause
G(2).a  = 0;                           % endogenous cause (action)
G(2).V  = exp(16);

% state-dependent precision (attentional bias) in generative model (M):
%--------------------------------------------------------------------------
M(1).f  = inline('v - x/4','x','v','P');
M(1).g  = inline('[x(1); sum(x)]','x','v','P');
M(1).x  = [0; 0];                      % hidden states
M(1).v  = [0; 0];                      % hidden causes
M(1).W  = exp(4);                      % precision (states)
M(1).ph = inline('[1; 1]*(8 - h*tanh(v(1) + x(1)))','x','v','h','M');
M(1).hE = 6;

% level 2; causes
%--------------------------------------------------------------------------
M(2).v  = [0; 0];                      % hidden cause
M(2).V  = [exp(6); exp(0)];

% Demonstration of the need for sensory attenuation
%==========================================================================

% hidden cause and prior expectations
%--------------------------------------------------------------------------
N      = 32;
C      = zeros(1,N);
U(1,:) = exp(-((1:N) - N/2).^2/(4.^2))*1;
U(2,:) = zeros(1,N);

% assemble model structure
%--------------------------------------------------------------------------
DEM.M = M;
DEM.G = G;
DEM.C = C;
DEM.U = U;

hE    = (-4:1:6);
for i = 1:length(hE)

    rng('default')

    LAP         = DEM;
    LAP.M(1).hE = hE(i);
    LAP         = spm_ALAP(LAP);

    % true and perceived force exerted (endogenously)
    %----------------------------------------------------------------------
    Px(i) = max(LAP.pU.x{1}(1,:));
    Qx(i) = max(LAP.qU.x{1}(1,:));

    % plot self-generated movement
    %----------------------------------------------------------------------
    if hE(i) == 0

        spm_figure('GetWin','Figure 1: Low attenuation');
        spm_DEM_qU(LAP.qU,LAP.pU)

    elseif hE(i) == 6

        spm_figure('GetWin','Figure 2: High attenuation');
        spm_DEM_qU(LAP.qU,LAP.pU)

    end
end

% adjust axes
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2: High attenuation');
subplot(2,2,2); spm_axis tight, a = axis;
subplot(2,2,1); axis(a);
spm_figure('GetWin','Figure 1: Low attenuation');
subplot(2,2,2); axis(a);
subplot(2,2,1); axis(a);

spm_figure('GetWin','Figure 2: High attenuation');
subplot(2,2,3); spm_axis tight, a = axis;
subplot(2,2,4); axis(a);
spm_figure('GetWin','Figure 1: Low attenuation');
subplot(2,2,3); axis(a);
subplot(2,2,4); axis(a);

% plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3');

subplot(2,1,1)
plot(hE,[Px; Qx])
axis square
xlabel('attenuation of sensory precicion','FontSize',12)
ylabel('true and perceived force exerted','FontSize',12)
legend({'true','perceived'})
title('Sensory attenuation and action','FontSize',16)

% Demonstration of sensory attenuation
%==========================================================================

% replay internal force as external force
%--------------------------------------------------------------------------
rng('default')
DEM.C  = [C LAP.pU.x{1}];
DEM.U  = [U sparse(2,N)];
DEM    = spm_ALAP(DEM);

% plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4');
spm_DEM_qU(DEM.qU,DEM.pU)

% Force matching
%==========================================================================

% replay perceived (at 90% confidence) internal force
%--------------------------------------------------------------------------
for i = 1:N
    CI(i) = 1.694*sqrt(LAP.qU.S{i}(1,1));
end

DEM.C  = [C (LAP.qU.x{1}(1,:) - CI)];
DEM.U  = [U sparse(2,N)];
DEM    = spm_ALAP(DEM);

% plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 5');
spm_DEM_qU(DEM.qU,DEM.pU)
subplot(2,2,2); spm_axis tight, a = axis;
subplot(2,2,1); axis(a);
subplot(2,2,3); spm_axis tight, a = axis;
subplot(2,2,4); axis(a);

% Use a range of self-generated forces - with D = 0
%==========================================================================
F     = (1:4)/2;
for i = 1:length(F)

    DEM.C = C;
    DEM.U = U*F(i);
    DEM   = spm_ALAP(DEM);

    % self-generated and matched (inferred) force
    %----------------------------------------------------------------------
    [x j] = max(DEM.pU.x{1}(1,:));
    Sx(i) = x;
    Tx(i) = DEM.qU.x{1}(1,j) - 1.694*sqrt(DEM.qU.S{j}(1,1));

end

% repeat with (delusional) precision D = 2
%--------------------------------------------------------------------------
D           = 2;
DEM.M(1).hE = 6 - D;
DEM.M(1).W  = exp(4 + D);
DEM.M(2).V  = [exp(6 + D); exp(0)];

for i = 1:length(F)

    DEM.C = C;
    DEM.U = U*F(i);
    DEM   = spm_ALAP(DEM);

    % self-generated and matched (inferred) force
    %----------------------------------------------------------------------
    [x j] = max(DEM.pU.x{1}(1,:));
    sx(i) = x;
    tx(i) = DEM.qU.x{1}(1,j) - 1.694*sqrt(DEM.qU.S{j}(1,1));

end

% plot results of force matching
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 6');

subplot(2,1,1)
plot(Tx,Sx,'-'); hold on,
plot(Tx,Sx,'.','MarkerSize',32), hold on
plot(tx,sx,'-.'); hold on,
plot(tx,sx,'.','MarkerSize',16), hold on
plot([0 3],[0 3],':'); hold off

axis square
xlabel('external (perceptually matched) force','FontSize',12)
ylabel('self-generated force','FontSize',12)
title('Force matching illusion','FontSize',16)

% Illustrate false inference with delusional precision
%==========================================================================
D           = 4;
DEM.M(1).hE = 6 - D;
DEM.M(1).W  = exp(4 + D);
DEM.M(2).V  = [exp(6 + D); exp(0)];

% replay perceived (at 90% confidence) internal force
%--------------------------------------------------------------------------
DEM.C = C;
DEM.U = U*2;
DEM   = spm_ALAP(DEM);

for i = 1:N
    CI(i) = 1.694*sqrt(DEM.qU.S{i}(1,1));
end

DEM.C = [C (DEM.qU.x{1}(1,:) - CI)];
DEM.U = [U*2 sparse(2,N)];
DEM   = spm_ALAP(DEM);

% plot false inference under high D
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 7: False (delusional) inference');
spm_DEM_qU(DEM.qU,DEM.pU)

author: spm


text:
function p = MS_fnn(y,de,tau,th,kth);
% function nfnn = MS_fnn(y,de,tau,th,kth)
%
% determine the number of false nearest neighbours for the time
% series y embedded in dimension de with lag tau.
%
% for each pair of values (de,tau) the data y is embeded and the
% nearest neighbour to each point (excluding the immediate
% neighbourhood of n points) is determined. If the ratio of the
% distance of the next (kth) points and these points is greater than
% th then they are counted as false nearest neighbours.
%
% default:
% th=5
% kth=1
%
% p(i,j) is the proportion of false nearest neighbours for de(i)
% and tau(j).
%
% Michael Small
% michael.small@uwa.edu.au, http://school.maths.uwa.edu.au/~small/
% 3/3/2005
% For further details, please see M. Small. Applied Nonlinear Time Series
% Analysis: Applications in Physics, Physiology and Finance. Nonlinear Science
% Series A, vol. 52. World Scientific, 2005. (ISBN 981-256-117-X) and the
% references therein.
% (minor cosmetic changes made by Ben Fulcher, 2010)

if nargin < 5
  kth = 1;
  disp(['th = ',int2str(kth)]);
end
if nargin < 4
  th = 5;
  disp(['th = ',int2str(th)]);
end
if nargin < 3
  tau = MS_firstzero(y);
  disp(['tau = ',int2str(tau)]);
end
if nargin < 2
  de = [1:10];
  disp(['de = ',int2str(de(1)),':',int2str(de(end))]);
end

p=[];
for t=tau,
  px=[];
  for d=de
    %embed the data
    X=MS_embed(y,d,t);
    [dx,nx]=size(X);

    %find the nearest neighbours of each point
    ind=MS_nearest(X(:,1:(nx-kth)),tau); %whooh hooo!

    %distance between each point and its nearest neighbour
    d0=MS_rms(X(:,(1:(nx-kth)))'-X(:,ind)');
    %... and after one time step
    d1=MS_rms(X(:,(kth+1):nx)'-X(:,ind+1)');

    %exclude any coincident points
    d1(d0==0)=[];
    d0(d0==0)=[];

    %calculate the proportion fnn
    ifnn=sum((d1./d0)>th)/length(d0);

    %disp
    % disp(['tau = ', int2str(t),', de = ',int2str(d),', nfnn = ',num2str(ifnn*100),'%']);

    px=[px ifnn];
  end;

  p=[p;px];

end;

p=p';

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% The COBRAToolbox: testDetermineSignMatrix.m
%
% Purpose:
%     - test the function to determine the binary version of the stoichiometric matrix S
%
% Authors:
%     - Laurent Heirendt, November 2018
%

% save the current path and initialize the test
currentDir = fileparts(which(mfilename));

% load the model
model = createToyModelForLifting(false);

% load reference test data
load([currentDir filesep 'refData_determineBinaryMatrix.mat']);

% run the function with full output
[test_Shat, test_Shatabs, test_mconnect, test_nconnect, test_mconnectin, test_mconnectout] = determineSignMatrix(model.S);

assert(isequal(test_Shat, Shat))
assert(isequal(test_Shatabs, Shatabs))
assert(isequal(test_mconnect, mconnect))
assert(isequal(test_nconnect, nconnect))
assert(isequal(test_mconnectin, mconnectin))
assert(isequal(test_mconnectout, mconnectout))

% run the function with partial output
[test_Shat, test_Shatabs, test_mconnect, test_nconnect] = determineSignMatrix(model.S);

assert(isequal(test_Shat, Shat))
assert(isequal(test_Shatabs, Shatabs))
assert(isequal(test_mconnect, mconnect))
assert(isequal(test_nconnect, nconnect))

% run the function with minimal output
[test_Shat, test_Shatabs, test_mconnect, test_nconnect] = determineSignMatrix(model.S);

assert(isequal(test_Shat, Shat))
assert(isequal(test_Shatabs, Shatabs))

% test the sorted flag
[test_Shat, test_Shatabs, test_mconnect, test_nconnect] = determineSignMatrix(model.S, true);

assert(isequal(test_mconnect, sort(mconnect, 'descend')));
assert(isequal(test_nconnect, sort(nconnect, 'descend')));

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function test_ft_clusterplot

% MEM 12gb
% WALLTIME 00:20:00
% DEPENDENCY ft_clusterplot ft_statistics_montecarlo ft_timelockstatistics clusterstat findcluster

cd(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/eventrelatedstatistics'));

% use the tutorial on http://www.fieldtriptoolbox.org/tutorial/eventrelatedstatistics

load ERF_orig
load GA_ERF_orig

cfg = [];
cfg.method      = 'template';
cfg.template    = 'ctf151_neighb.mat';               % specify type of template
cfg.layout      = 'CTF151.lay';                      % specify layout of sensors*
cfg.feedback    = 'yes';                             % show a neighbour plot
neighbours      = ft_prepare_neighbours(cfg, GA_FC); % define neighbouring channels

cfg = [];
cfg.channel     = 'MEG';
cfg.neighbours  = neighbours; % defined as above
cfg.latency     = [-0.25 1];
cfg.avgovertime = 'no';
cfg.parameter   = 'avg';
cfg.method      = 'montecarlo';
cfg.statistic   = 'ft_statfun_depsamplesT';
cfg.alpha       = 0.05;
cfg.correctm    = 'cluster';
cfg.correcttail = 'prob';
cfg.numrandomization = 'all';
cfg.minnbchan        = 2; % minimal neighbouring channels
Nsub = 10;
cfg.design(1,1:2*Nsub)  = [ones(1,Nsub) 2*ones(1,Nsub)];
cfg.design(2,1:2*Nsub)  = [1:Nsub 1:Nsub];
cfg.ivar                = 1; % the 1st row in cfg.design contains the independent variable
cfg.uvar                = 2; % the 2nd row in cfg.design contains the subject number
stat = ft_timelockstatistics(cfg,allsubjFIC{:},allsubjFC{:});

%% chan_time, this is the original one
cfg = [];
cfg.highlightsymbolseries = ['*','*','.','.','.'];
cfg.layout = 'CTF151.lay';
cfg.contournum = 0;
cfg.markersymbol = '.';
cfg.alpha = 0.05;
cfg.parameter='stat';
cfg.zlim = [-5 5];
cfg.toi  = [-0.1:0.1:0.9];
ft_clusterplot(cfg,stat);

%% chan_freq
stat.freq = 1:376;
cfg.toi = [20:20:200];
stat = rmfield(stat, 'time');
stat.dimord = 'chan_freq';

ft_clusterplot(cfg,stat);

%% chan_freq_time with singleton time
stat.time = 1;
stat.dimord = 'chan_freq_time';
cfg.toi = [20:20:200];

ft_clusterplot(cfg,stat);

%% chan_freq_time with singleton freq
stat.time = 1:376;
stat.freq = 1;
stat.dimord = 'chan_freq_time';

% insert the singleton dimension in the data
stat.posclusterslabelmat = reshape(stat.posclusterslabelmat, [151 1 376]);
stat.negclusterslabelmat = reshape(stat.negclusterslabelmat, [151 1 376]);
stat.prob     = reshape(stat.prob,    [151 1 376]);
stat.cirange  = reshape(stat.cirange, [151 1 376]);
stat.mask     = reshape(stat.mask,    [151 1 376]);
stat.stat     = reshape(stat.stat,    [151 1 376]);
stat.ref      = reshape(stat.ref,     [151 1 376]);

cfg.toi = [20:20:200];

ft_clusterplot(cfg,stat);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function [output] = ols_reg(Y,X,options)

% this function computes the ols coefficients of Y on X
% Y TxN
% X Txk (interepct at the end)
% Y = XB + E

nindex      = find(sum(isnan([Y,X]),2)==0);
index       = find(sum(isnan([Y,X]),2)>0);
Y(index,:)  = [];
X(index,:)  = [];
[N,K]       = size(X);
robust_se_  = 0;
L           = round(N/4);

if nargin > 2
    if isfield(options,'robust_se_') == 1
        robust_se_ = options.robust_se_;
    end
    if isfield(options,'L') == 1
        L = options.L;
    end
end

ny = size(Y,2);

XX      = (X'*X);
iXX     = XX\eye(K);
Bols    = X\Y;                  % Bols   = (X'*X)\(X'*Y);
err     = Y - X*Bols;
se      = nan(K,ny);
Sols    = zeros(K*ny);

if robust_se_==2 % NW Robust SE
    %-----------------------------------------------------------------------%
    Serror     = (err'*err);
    nwWeights  = (L+1-(1:L))./(L+1);
    for j= 1 : L
        G      = (err(j+1 : N, :)'*err(1 : N-j , :));
        Serror = Serror + nwWeights(j) * (G + G');
    end
    Serror = Serror/(N-K);
    Sols   = kron(diag(diag(Serror)),iXX);
    se     = reshape(sqrt(diag(Sols)), K, ny);

elseif robust_se_ == 1 % Hamilton (1994), Ch 10 pag 282, eq (10.5.20)
    %-----------------------------------------------------------------------%
    for vv = 1: ny % equation by equation
        u= err(:,vv);
        errs=X.*u;
        V0 = [errs'*errs] / N ; %regular weighting matrix
        for ind_i = (1:L)
            S = errs(1:N-ind_i,:)'*errs(1+ind_i:N,:) / N;
            V0 = V0 + (1 - ind_i/(L+1))*(S + S');
        end
        %     D    =   inv((X'*X)/N);
        %     varb = 1/N*D*V1*D;
        Solsj =  N * iXX * V0 * iXX;
        Sols((vv-1)*K+1 : vv*K, (vv-1)*K+1 : K*vv ) = Solsj;
    end
    Serror = (err'*err)/(N-K); % not sure this is the correct Coavariance of the shocks.
    se     = reshape(sqrt(diag(Sols)), K, ny);

elseif robust_se_ == 5 % Matlab HAC function
    %-----------------------------------------------------------------------%
    if  exist('hac') ==2
        % find constant
        index = find(sum(diff(X,1,1),1)~=0);
        indexC = find(sum(diff(X,1,1),1)==0);
        for vv = 1: ny
%             [EstCoeffCov0,se0,~] = hac(X(:,index),Y(:,ny),'display','off','type','HC');
            [EstCoeffCov0,se0,~] = hac(X(:,index),Y(:,ny),'display','off','bandwidth','AR1OLS'); % remove the intercept
            se(index,vv)  = se0(2:end);
            se(indexC,vv) = se0(1); %intercept
            Sols((vv-1)*K+1 : vv*K, (vv-1)*K+1 : K*vv ) = EstCoeffCov0; % order is not correct
        end
        Serror  = 1/(N-K)*(err'*err);
%         output.TtestRobust  = coeff./se;
%         output.pvalueRobust = tpdf(output.TtestRobust,N-K);
    else
        error('Matlab Econ Toolbox missing')
    end
else
    Serror  = 1/(N-K)*(err'*err);
    Sols    = kron(diag(diag(Serror)),iXX);
    se      = reshape(sqrt(diag(Sols)), K, ny);
end

Ttest  = Bols./ reshape(sqrt(diag(Sols)), K, ny);
ESS    = diag((X*Bols - mean(Y))'*(X*Bols - mean(Y)));
RSS    = diag(err' * err);
TSS    = diag((Y - mean(Y))' * (Y - mean(Y)));
R2     = ones(length(ESS),1) - RSS ./ TSS;
adjR2  = ones(length(ESS),1) - (ones(length(ESS),1) - R2)*(N-1)/(N-K);
Ftest  = ESS/(K-1) ./ diag(Serror);

% for v = 1 : ny
% %     output.logl(v,1)   = -N/2*log(2*pi*Serror(v,v)) - RSS(v,1)/(2*Serror(v,v));
%     [output.AIC(v,1), output.SIC(v,1), output.HQIC(v,1)] = IC(output.logl(v,1), N, K);
% end

output.beta   = Bols;                   % OLS estimator
output.error  = err;                   % (TxN) matrix of Residuals
output.e_ols  = err;                   % (TxN) matrix of Residuals
output.Serror = Serror;                 % Covariance matrix of Residuals
output.Sigma_ols = Serror;                 % Covariance matrix of Residuals
output.Sols   = Sols;                   % Covariance matrix of Bols
output.Ttest  = Ttest;                  % t-statistics
output.pvalue = tpdf(Ttest,N-K);        % p-value
output.Ftest  = Ftest;                  % F-test
output.R2     = R2;                     % R2
output.adjR2  = adjR2;                  % Adjusted R2
output.yfit   = X*Bols;                 % Fitted Values
output.N      = N;                      % # of observation used
output.K      = K;                      % # of regressors
output.nindex = nindex;                 % index of missing observation
output.index  = index;                  % index of observation used
output.se     = se;
output.XX     = XX;
output.X      = X;
output.Y      = Y;

for v = 1 : ny
%     output.logl(v,1)   = -N/2*log(2*pi*Serror(v,v)) - RSS(v,1)/(2*Serror(v,v));
%     [output.AIC(v,1), output.SIC(v,1), output.HQIC(v,1)] = IC(output, N, K);
    S = Serror(v,v); E = err(:,v);
    [output.AIC(v,1), output.SIC(v,1), output.HQIC(v,1)] = IC(S, E, N, K);
end

% if nargin > 2
%     [EstCov,se,coeff] = hac(X(:,1:end-1),Y,'display','off'); % remove the intercept
%     output.serob        = se(2:end);
%     output.serob(end+1) = se(1); %intercept
%     output.TtestRobust  = coeff./se;
%     output.pvalueRobust = tpdf(output.TtestRobust,N-K);
% end


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = sum((X * theta - y) .^ 2) / (2 * m);

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
clear all

%define the global variables
global TrainSetX;
global TrainSetY;
global TrainXY;
global TestSetX;
global TestSetY;

global eta;
global gama;
global afa;
global npos;

global kneighbor_test;
global kneighbor_testp;
global kneighbor_mix;

%set the parameters
afa = 0.4;
gama = 15;
eta = 0.12;

kneighbor_testp = textread('KneighborFile.txt');
kneighbor_testp = kneighbor_testp(:,1:(size(kneighbor_testp,2)-1));
kneighbor_testp = kneighbor_testp + 1;
A = textread('data/Train1.data');
TrainSetX = spconvert(A);
A = textread('data/Test1.data');
TestSetX = spconvert(A);
clear A;
TrainSetY = textread('data/Train1.label');
TestSetY = textread('data/Test1.label');
indexi = find(TestSetY == 1);
npos = sum(TestSetY(indexi));

TrainXY = scale_cols(TrainSetX,TrainSetY);
DataX = [TrainSetX TestSetX];
pos = 0;
% read the model trained from Logistic Regression
wbest = textread('orimodel.model');
ptemp = 1./(1 + exp(-wbest'*TrainSetX));
Uconf_mix = [ptemp',1-ptemp'];
TrainA = getResult(ptemp,TrainSetY);
ptemp = 1./(1 + exp(-wbest'*TestSetX));
TestA = getResult(ptemp,TestSetY);
Uconf_mix = [Uconf_mix;[ptemp',1-ptemp']];
clear ptemp;
fprintf('The original training and testing accuracy:  %g:    %g  ...... %g\n',pos,TrainA*100,TestA*100);result = [];
[W,result] = PRP_CG(wbest,npos,gama,afa,0.1);
ptemp = 1./(1 + exp(-W'*TestSetX));
TestA = getResult(ptemp,TestSetY);
fprintf('The final testing accuracy: %g\n',TestA*100);
result(1,length(result)+1) = TestA;
xlswrite('result.xls',result');

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function pairsChart(LCO, WTI)
% PAIRSCHART: Shows the full time history for LCO and WTI

% Copyright 2011, The MathWorks, Inc.
% All rights reserved.

plot(LCO(:,1), LCO(:,4), WTI(:,1), WTI(:,4))

datetick('x')
grid on

xlabel('Date')
ylabel('Price (USD)')

title('Intraday prices for LCO and WTI')
legend('LCO', 'WTI', 'Location', 'best')

author: Sable

repo: mcbench-benchmarks


text:
function [model, L] = mlpReg(X, y, k, lambda)
% Train a multilayer perceptron neural network for regression with backpropagation
% tanh activation function is used
% Input:
%   X: d x n data matrix
%   y: 1 x n real value response vector
%   k: T x 1 vector to specify number of hidden nodes in each layer
%   lambda: regularization parameter
% Ouput:
%   model: model structure
%   L: (regularized least square) loss
% Written by Mo Chen (sth4nth@gmail.com).
if nargin < 4
    lambda = 1e-2;
end
eta = 1e-5;
tol = 1e-5;
maxiter = 50000;
L = inf(1,maxiter);

k = [size(X,1);k(:);size(y,1)];
T = numel(k)-1;
W = cell(T,1);
b = cell(T,1);
for t = 1:T
    W{t} = randn(k(t),k(t+1));
    b{t} = randn(k(t+1),1);
end
R = cell(T,1);
Z = cell(T+1,1);
Z{1} = X;
for iter = 2:maxiter
%     forward
    for t = 1:T-1
        Z{t+1} = tanh(W{t}'*Z{t}+b{t});             % 5.10 5.113
    end
    Z{T+1} = W{T}'*Z{T}+b{T};                       % 5.114

%     loss
    E = Z{T+1}-y;
    Wn = cellfun(@(x) dot(x(:),x(:)),W);            % |W|^2
    L(iter) = dot(E(:),E(:))+lambda*sum(Wn);
    if abs(L(iter)-L(iter-1)) < tol*L(iter-1); break; end

%     backward
    R{T} = E;
    for t = T-1:-1:1
        df = 1-Z{t+1}.^2;    % h'(a)
        R{t} = df.*(W{t+1}*R{t+1});    % 5.66
    end

%     gradient descent
    for t=1:T
        dW = Z{t}*R{t}'+lambda*W{t};    % 5.67
        db = sum(R{t},2);
        W{t} = W{t}-eta*dW;             % 5.43
        b{t} = b{t}-eta*db;
    end
end
L = L(2:iter);
model.W = W;
model.b = b;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
% ^   Matrix power.
%    Z = X^y is X to the y power if y is a scalar and X is square. If y
%    is an integer greater than one, the power is computed by repeated
%    squaring. For other values of y the calculation involves
%    eigenvalues and eigenvectors.
%
%    Z = x^Y is x to the Y power if Y is a square matrix and x is a scalar.
%    Computed using eigenvalues and eigenvectors.
%
%    Z = X^Y, where both X and Y are matrices, is an error.
%
%    C = MPOWER(A,B) is called for the syntax 'A ^ B' when A or B is an
%    object.
%
%    See also POWER.
%
%    Reference page in Doc Center
%       doc mpower
%
%    Other functions named mpower
%
%       codistributed/mpower    gpuArray/mpower    sym/mpower    ts/mpower
%

author: jmaih


text:
function T = qq(q,varargin)
% returns w * q' * q
%
% Input
%  q - list of quaternions
%  w - list of weights

notNaN = ~isnan(q.a);
ql = [q.a(:), q.b(:), q.c(:), q.d(:)];
ql = ql(notNaN,:);

% weigths
w = get_option(varargin,'weights');

if ~isempty(w)
  w = w(notNaN);
  w = w./sum(w);
  w = repmat(w(:).',4,1);
  T = (w .* ql.') * ql;
else
  T = ql.'*ql;
  T = T ./ size(ql,1);
end

author: mtex-toolbox

repo: mtex

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


text:
% SVMTRAIN - Trains a support vector machine incrementally
%            using the L1 soft margin approach developed by
%            Cauwenberghs for two-class problems.
%
% Syntax: [a,b,g,ind,uind,X_mer,y_mer,Rs,Q] = svmtrain(X,y,C,type,scale)
%         [a,b,g,ind,uind,X_mer,y_mer,Rs,Q] = svmtrain(X,y,C,type,scale,uind)
%         (trains a new SVM on the given examples)
%
%         [a,b,g,ind,uind,X_mer,y_mer,Rs,Q] = svmtrain(X,y,C)
%         [a,b,g,ind,uind,X_mer,y_mar,Rs,Q] = svmtrain(X,y,C,uind)
%         (trains the current SVM in memory on the given examples)
%
%      a: alpha coefficients
%      b: bias
%      g: partial derivatives of cost function w.r.t. alpha coefficients
%    ind: cell array containing indices of margin, error and reserve vectors
%         ind{1}: indices of margin vectors
%         ind{2}: indices of error vectors
%         ind{3}: indices of reserve vectors
%   uind: column vector of user-defined example indices (used for unlearning specified examples)
%  X_mer: matrix of margin, error and reserve vectors stored columnwise
%  y_mer: column vector of class labels (-1/+1) for margin, error and reserve vectors
%     Rs: inverse of extended kernel matrix for margin vectors
%      Q: extended kernel matrix for all vectors
%      X: matrix of training vectors stored columnwise
%      y: column vector of class labels (-1/+1) for training vectors
%      C: soft-margin regularization parameter(s)
%         dimensionality of C       assumption
%         1-dimensional vector      universal regularization parameter
%         2-dimensional vector      class-conditional regularization parameters (-1/+1)
%         n-dimensional vector      regularization parameter per example
%         (where n = # of examples)
%   type: kernel type
%           1: linear kernel        K(x,y) = x'*y
%         2-4: polynomial kernel    K(x,y) = (scale*x'*y + 1)^type
%           5: Gaussian kernel with variance 1/(2*scale)
%  scale: kernel scale
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function [a,b,g,ind,uind,X,y,Rs,Q] = svmtrain(X_new,y_new,C_new,varargin)

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLEARNED = 4;

% create a vector containing the regularization parameter
% for each example if necessary
if (length(C_new) == 1)              % same regularization parameter for all examples
   C_new = C_new*ones(size(y_new));
elseif (length(C_new) == 2)          % class-conditional regularization parameters
   flags = (y_new == -1);
   C_new = C_new(1)*flags + C_new(2)*(~flags);
end;

if (nargin >= 5)

   % define arguments
   type_new = varargin{1};
   scale_new = varargin{2};
   if (nargin == 6)
      uind_new = varargin(3);
   else
      uind_new = zeros(size(y_new));
   end;

   new_model = 1;
else

   % define arguments
   if (nargin == 4)
      uind_new = varargin(1);
   else
      uind_new = zeros(size(y_new));
   end;

   new_model = 0;
end;

% define global variables
global a;                     % alpha coefficients
global b;                     % bias
global C;                     % regularization parameters
global deps;                  % jitter factor in kernel matrix
global g;                     % partial derivatives of cost function w.r.t. alpha coefficients
global ind;                   % cell array containing indices of margin, error, reserve and unlearned vectors
global kernel_evals;          % kernel evaluations
global max_reserve_vectors;   % maximum number of reserve vectors stored
global perturbations;         % number of perturbations
global Q;                     % extended kernel matrix for all vectors
global Rs;                    % inverse of extended kernel matrix for margin vectors
global scale;                 % kernel scale
global type;                  % kernel type
global uind;                  % user-defined example indices
global X;                     % matrix of margin, error, reserve and unlearned vectors stored columnwise
global y;                     % column vector of class labels (-1/+1) for margin, error, reserve and unlearned vectors

% initialize variables
deps = 1e-3;
max_reserve_vectors = 3000;

if (new_model)

   num_examples = size(X_new,2);

   a = zeros(num_examples,1);
   b = 0;
   C = C_new;
   g = zeros(num_examples,1);
   ind = cell(4,1);
   ind{UNLEARNED} = 1:num_examples;
   kernel_evals = 0;
   perturbations = 0;
   Q = y_new';
   Rs = Inf;
   scale = scale_new;
   type = type_new;
   uind = uind_new;
   X = X_new;
   y = y_new;

else

   num_examples = size(X,2);
   num_new_examples = size(X_new,2);

   a = [a ; zeros(num_new_examples,1)];
   C = [C ; C_new];
   g = [g ; zeros(num_new_examples,1)];
   ind{UNLEARNED} = (1:num_new_examples) + num_examples;

   % assumes currently that there are no duplicate examples in the data - may not necessarily be true!
   Q_new = [y_new' ; (y(ind{MARGIN})*y_new').*kernel(X(:,ind{MARGIN}),X_new,type,scale)];

   Q = [Q Q_new];
   uind = [uind ; uind_new];
   X = [X X_new];
   y = [y ; y_new];

   num_examples = num_examples + num_new_examples;

end;

% begin incremental learning - enforce all constraints on each iteration
num_learned = 1;
disp('Beginning training.');
while (any(ind{UNLEARNED}))

   % randomly select example
   i = round(rand*(length(ind{UNLEARNED})-1)) + 1;
   indc = ind{UNLEARNED}(i);
%  indc = ind{UNLEARNED}(1);

	% learn example
   learn(indc,1);

   if (mod(num_learned,50) == 0)
      s = sprintf('Learned %d examples.',num_learned);
      disp(s);
   end;
   num_learned = num_learned + 1;

end;
if (mod(num_learned-1,50) ~= 0)
   s = sprintf('Learned %d examples.',num_learned-1);
   disp(s);
end;
disp('Training complete!');

% begin incremental learning - perform multiple passes through the data
% until all of the examples are learned
%while (any(ind{UNLEARNED}))
%   while (any(ind{UNLEARNED}))
%
%      % select example
%      indc = ind{UNLEARNED}(1);
%
%      % learn example
%      s = sprintf('\nLearning example %d...',indc);
%      disp(s);
%      learn(indc,0);
%
%   end;
%
%   % check to see if any reserve vectors are incorrectly classified
%   % if so, change their status to unlearned
%   ind_temp = find(g(ind{RESERVE}) < 0);
%   [ind{RESERVE},ind{UNLEARNED}] = move_ind(ind{RESERVE},ind{UNLEARNED},ind{RESERVE}(ind_temp));
%
%end;

% remove all but the closest reserve vectors from the dataset if necessary
if (length(ind{RESERVE}) == max_reserve_vectors)
   ind_keep = [ind{MARGIN} ind{ERROR} ind{RESERVE}];
   a = a(ind_keep);
   g = g(ind_keep);
   Q = Q(:,ind_keep);
   uind = uind(ind_keep);
   X = X(:,ind_keep);
   y = y(ind_keep);
   ind{MARGIN} = 1:length(ind{MARGIN});
   ind{ERROR} = length(ind{MARGIN}) + (1:length(ind{ERROR}));
   ind{RESERVE} = length(ind{MARGIN}) + length(ind{ERROR}) + (1:length(ind{RESERVE}));
end;

% summary statistics
s = sprintf('\nMargin vectors:\t\t%d',length(ind{MARGIN}));
disp(s);
s = sprintf('Error vectors:\t\t%d',length(ind{ERROR}));
disp(s);
s = sprintf('Reserve vectors:\t%d',length(ind{RESERVE}));
disp(s);
s = sprintf('Kernel evaluations:\t%d\n',kevals);
disp(s);

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [data,units] = compute_max_wing_area(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = max(trx(fly).wing_areal_mm,trx(fly).wing_arear_mm);

end
units = parseunits('mm^2');

author: kristinbranson

repo: JAABA


text:
function pass = test_div( )
% Test divergence

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% Check divergence of an empty spherefunv is an empty spherefun.
u = spherefunv;
f = div(u);
pass(1) = isempty(f) & isa(f,'spherefun');

% Check divergence of a spherefunv is a spherefun.
u = spherefunv.unormal;
f = div(u);
pass(2) = isa(f,'spherefun');

% Check divergence of the zero field is zero
f = spherefun(@(x,y,z) 0*x);
u = spherefunv(f,f,f);
pass(3) = norm(div(u), inf) < tol;

%
% Check divergence of a non-zero field gives the correct result
%

% Example 1
f = spherefun(@(x,y,z) (1+2*x.*(y-x)).*z);
g = spherefun(@(x,y,z) (-1+2*y.*(y-x)).*z);
h = spherefun(@(x,y,z) -(x-y).*(2*z.^2 - 1));
u = spherefunv(f,g,h);
divu = div(u);
% Exact divergence
exact = spherefun(@(x,y,z) -6*(x-y).*z);
pass(4) = norm(divu-exact, inf) < tol;

% Example 2
f = spherefun(@(x,y,z) 4*x.*z.*sin(4*z));
g = spherefun(@(x,y,z) 4*y.*z.*sin(4*z));
h = spherefun(@(x,y,z) -4*(1-z.^2).*sin(4*z));
u = spherefunv(f,g,h);
divu = div(u);
% Exact divergence
exact = spherefun(@(x,y,z) -8*(2*(1-z.^2).*cos(4*z) - z.*sin(4*z)));
pass(5) = norm(divu-exact, inf) < 100*tol;

% Example 3
f = spherefun(@(x,y,z) 4*(x.^2-1).*sin(4*x));
g = spherefun(@(x,y,z) 4*x.*y.*sin(4*x));
h = spherefun(@(x,y,z) 4*x.*z.*sin(4*x));
u = spherefunv(f,g,h);
divu = div(u);
% Exact divergence
exact = spherefun(@(x,y,z) -8*(2*(1-x.^2).*cos(4*x) - x.*sin(4*x)));
pass(6) = norm(divu-exact, inf) < 100*tol;

% Example 4
f = spherefun(@(x,y,z) 4*x.*y.*sin(4*y));
g = spherefun(@(x,y,z) 4*(y.^2-1).*sin(4*y));
h = spherefun(@(x,y,z) 4*y.*z.*sin(4*y));
u = spherefunv(f,g,h);
divu = div(u);
% Exact divergence
exact = spherefun(@(x,y,z) -8*(2*(1-y.^2).*cos(4*y) - y.*sin(4*y)));
pass(7) = norm(divu-exact, inf) < 200*tol;

% Check that the div and divergence give the same result.
divergenceu = divergence(u);
pass(8) = isequal(divu,divergenceu);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function b = strendswith(s, pat)
%STRENDSWITH Determines whether a string ends with a specified pattern
%
%   b = strstartswith(s, pat);
%       returns whether the string s ends with a sub-string pat.
%

%   History
%   -------
%       - Created by Dahua Lin, on Oct 9, 2008
%

%% main

sl = length(s);
pl = length(pat);

b = (sl >= pl && strcmp(s(sl-pl+1:sl), pat)) || isempty(pat);


text:
function this = finish_processing_step(this, module, varargin)
% finishes current processing step by deleting duplicate data and storing
% results of processing step
%
%   MrSeries = finish_processing_step(MrSeries, module, varargin)
%
% This is a method of class MrSeries.
%
% IN
%   module      'realign', 'smooth', ...
%
% OUT
%
% EXAMPLE
%   finish_processing_step
%
%   See also MrSeries MrSeries.init_processing_step

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-01
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% NOTE: for each new processing step added here, it has to be decided which
% files are saved additionally or which temporary files can be deleted/renamed
itemsSave = this.parameters.save.items;
doSave = ~strcmpi(itemsSave, 'none');
doSaveRaw = ismember(itemsSave, {'all'});
doSaveNifti = ismember(itemsSave, {'nii', 'all', 'processed'});
doSaveObject = ismember(itemsSave, {'object', 'all', 'processed'});

% determine where and which files have been changed from input argument
if iscell(varargin{1});
    inputImage = varargin{1}{1};
else
    inputImage = varargin{1};
end

pathSave        = inputImage.parameters.save.path;
pathRaw         = fileparts(inputImage.get_filename('prefix', 'raw'));

% delete additional, processed files...
switch module

    case 'compute_masks'
        maskImages = varargin{1};
        nImages = numel(maskImages);

        filesMask = cell(nImages,1);
        for iImage = 1:nImages
            filesMask{iImage} = ...
                maskImages{iImage}.get_filename;
            if doSaveNifti
                maskImages{iImage}.save;
            end
        end

        filesProcessed = filesMask;

    case 'compute_stat_images'
        % file names and paths already given in init_processing_step
        if doSaveNifti
            handleImageArray = this.get_all_image_objects('stats');
            for iImage = 1:numel(handleImageArray)
                handleImageArray{iImage}.save;
                filesProcessed{iImage} = ...
                    handleImageArray{iImage}.get_filename;
            end
        end

    case 'compute_tissue_probability_maps'
        createdFields = varargin{1};
        inputImage = varargin{2};
        nImages = numel(createdFields);

        filesFieldImages = cell(nImages,1);
        for iImage = 1:nImages
            filesFieldImages{iImage} = fullfile(...
                pathSave, ...
                createdFields{iImage}.parameters.save.fileName);
        end

        fileRaw     = inputImage.get_filename('prefix', 'raw');
        fileSeg8    = regexprep(fileRaw, '\.nii$', '_seg8\.mat');

        filesProcessed 	= [
            {inputImage.get_filename}
            filesFieldImages
            {fileSeg8}
            ];

    case 'coregister'
        transformedImage = varargin{1};
        equallyTransformedImages = varargin{2};

        % set files for delete
        filesProcessed = {transformedImage.get_filename};
        nImages = numel(equallyTransformedImages);
        for iImage = 1:nImages
            filesProcessed{end+1} = equallyTransformedImages{iImage}.get_filename;
            % load additionally transformed images
            handleInputImages = this.find('MrImage', 'name', ...
                equallyTransformedImages{iImage}.name);
            nameTransformed = handleInputImages{1}.name;
            handleInputImages{1}.load(equallyTransformedImages{iImage}.get_filename);
            handleInputImages{1}.name = nameTransformed;
        end

    case 'realign' % load realignment parameters into object

        fileRealignmentParameters = regexprep( ...
            tapas_uniqc_prefix_files(inputImage.get_filename('prefix', 'raw'), ...
            'rp_'), '\.nii$', '\.txt') ;
        this.glm.regressors.realign = load(fileRealignmentParameters);

        movefile(fileRealignmentParameters, pathSave);

    case 'smooth'
        filesProcessed = inputImage.get_filename();

    case 't_filter'
        if doSaveNifti
            this.data.save();
        end

    case 'specify_and_estimate_1st_level'
        % load design matrix
        spmDirectory = fullfile(this.glm.parameters.save.path, ...
            this.glm.parameters.save.spmDirectory);
        SPM = load(fullfile(spmDirectory, 'SPM.mat'));
        this.glm.designMatrix = SPM.SPM.xX.nKX;
end

% delete raw sub-folder of current processing step
if ~doSaveRaw && exist(pathRaw, 'dir')
    delete(fullfile(pathRaw, '*'));
    rmdir(pathRaw);
end

if ~doSave
    tapas_uniqc_delete_with_hdr(filesProcessed);
end

% strip object data and save ...
if doSaveObject
    fileObject = fullfile(pathSave, 'MrObject.mat');
    MrObject = this.copyobj('exclude', 'data'); % copies object without data
    save(fileObject, 'MrObject');
end

author: translationalneuromodeling

repo: tapas


text:
% Calculate the log mel filter bank from waveform
% Author: Xiao Xiong
% Created: 4 Feb 2005
% Last modified: 4 Feb 2005

function [fbank] = wav2root_fbank(is_file_name, x, frame_shift);
if is_file_name ==1
    x = readNIST(x);
end
if nargin < 3
    fbank = (abs2Mel( wav2abs(x) )).^0.1;
else
    fbank = (abs2Mel( wav2abs(x,frame_shift) )).^0.1;
end

author: singaxiong

repo: SignalGraph


text:
function save_yuv(data, new_file, height, width, h_factor, w_factor, SDR_HDR)
%%% Save YUV frame of SDR or HDR videos %%%
% 'SDR_HDR' flag should be specified with either 'SDR' or 'HDR'

% get size of data
datasize = size(data);
datasizelength = length(datasize);

% open file
fid = fopen(new_file,'a');

% subsampling of U and V
if datasizelength == 2 || h_factor == 0
    %4:0:0
    y(1:height, 1:width) = data(:, :, 1);
elseif datasizelength == 3
    y(1:height, 1:width) = double(data(:, :, 1));
    u(1:height, 1:width) = double(data(:, :, 2));
    v(1:height, 1:width) = double(data(:, :, 3));
    if w_factor == 1
        %4:1:1
        u2 = u;
        v2 = v;
    elseif h_factor == 0.5
        %4:2:0
        u2(1:height/2, 1:width/2) = u(1:2:end, 1:2:end)+u(2:2:end, 1:2:end)+u(1:2:end, 2:2:end)+u(2:2:end, 2:2:end);
        u2                         = u2/4;
        v2(1:height/2, 1:width/2) = v(1:2:end, 1:2:end)+v(2:2:end, 1:2:end)+v(1:2:end, 2:2:end)+v(2:2:end, 2:2:end);
        v2                         = v2/4;
    elseif w_factor == 0.25
        %4:1:1
        u2(1:height, 1:width/4) = u(:, 1:4:end)+u(:, 2:4:end)+u(:, 3:4:end)+u(:, 4:4:end);
        u2                       = u2/4;
        v2(1:height, 1:width/4) = v(:, 1:4:end)+v(:, 2:4:end)+v(:, 3:4:end)+v(:, 4:4:end);
        v2                       = v2/4;
    elseif w_factor == 0.5 && h_factor == 1
        %4:2:2
        u2(1:height, 1:width/2) = u(:, 1:2:end)+u(:, 2:2:end);
        u2                       = u2/2;
        v2(1:height, 1:width/2) = v(:, 1:2:end)+v(:, 2:2:end);
        v2                       = v2/2;
    end
end

if strcmp(SDR_HDR,'HDR')
    fwrite(fid,uint16(y'),'uint16'); % write Y-Data

    if h_factor ~= 0
        % write U- and V-Data if not 4:0:0 format
        fwrite(fid, uint16(u2'), 'uint16');
        fwrite(fid, uint16(v2'), 'uint16');
    end
elseif strcmp(SDR_HDR, 'SDR')
        fwrite(fid, uint8(y'), 'uchar'); % write Y-Data

    if h_factor ~= 0
        % write U- and V-Data if not 4:0:0 format
        fwrite(fid, uint8(u2'), 'uchar');
        fwrite(fid, uint8(v2'), 'uchar');
    end
end

fclose(fid);

author: sooyekim

repo: Deep-SR-ITM


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function [ucm2, ucms, elapsed_time] = contours_to_ucm(I, scales, E, O, all_parameters, param_multi, align_thr)
  % function [ucm2, ucms, elapsed_time] = contours_to_ucm(I, scales, E, O, all_parameters, param_multi, align_thr)
  % Multiscale hierarchical segmentation on RGBD images
  % Pablo Arbelaez
  % arbelaez@berkeley.edu

  if nargin<5,
    all_parameters = [5 60 4 2 0.1];
    % all_parameters = [3.3 60 4 2 0.1];
    % all_parameters = [30 60 4 2 0.09]; % NYUD2 train parameters
  end
  if nargin<6,
      param_multi.weights      = [0.7 1 1.3];
  end
  if nargin<7,
      align_thr = 0.13;
  end

  param.mult_Pb   = all_parameters(1);
  param.sat_sPb   = all_parameters(2);
  param.nvec      = all_parameters(3);
  param.dthresh   = all_parameters(4);
  param.ic_gamma  = all_parameters(5);

  [param.tx, param.ty, ~] = size(I);

  % compute ucms at multiple scales
  ucms = cell(numel(scales),1);
  elapsed_time = zeros(11,1);
  for s = 1:numel(scales),
      param.scale = scales(s);
      [ucms{s}, times] = img2ucm_scale(E{s}, O{s}, param);
      elapsed_time(3*(s-1)+1:3*s)=times;
  end

  % align ucms
  T=tic;
  ucms = project_ucms_wrap(ucms, align_thr);
  elapsed_time(10)=toc(T);

  % combine ucms
  T=tic;
  ucm2 = ucms2multi(ucms, param_multi);
  elapsed_time(11)=toc(T);

%%
function [ucm2, times] = img2ucm_scale(E, O, param)
  % compute hierarchical segmentation at each scale independently
  % I = imresize(I, param.scale, 'lanczos3');
  [ucm2, times] = img2ucm(E, O, param.mult_Pb, param.sat_sPb, param.nvec, param.dthresh, param.ic_gamma);

  % resample ucm to original image size
  T=tic;
  if ~isequal([param.tx, param.ty], size(ucm2(3:2:end, 3:2:end))),
      % TODO: Instead of resampling for each threshold, resample leaves and redo mergings
      ucm2 = resample_ucm2(ucm2, [param.tx, param.ty]);
  end
  times(3)=toc(T);

%%
function [ucm2, times] = img2ucm(E, O, mult_Pb, sat_sPb, nvec, dthresh, ic_gamma)

  % edge detection
  T=tic;
  % SAURABH: plug the rgbd contours here
  % [E,~,O] = edgesDetect( I, model );
  times(1)=toc(T);

  % continuous oriented watershed transform
  T=tic;
  [owt2, superpixels] = contours2OWT(E, O);

  % globalization
  [ sPb_thin] = spectralPb_fast(owt2 * mult_Pb, nvec, ic_gamma, dthresh) / sat_sPb;

  % ultrametric contour map with mean pb.
  ucm2 = double(ucm_mean_pb( (owt2 + sPb_thin), superpixels) );
  times(2)=toc(T);

%%
function ucm2 = ucms2multi(all_ucms, param)

  %combine ucms
  weights = param.weights;
  weights = weights ./ sum(weights);

  sz = size(all_ucms);
  W_all = repmat(repmat(weights', [1,sz(2)]),[1,1,sz(1)]); W_all = permute(W_all, [3 2 1]);
  all_ucms = all_ucms.*W_all;

  ucm2_wt = sum(all_ucms,3);

  labels2 = bwlabel(ucm2_wt == 0, 8);
  labels = labels2(2:2:end, 2:2:end) - 1; % labels begin at 0 in mex file.
  ucm2 = double(ucm_mean_pb(ucm2_wt, labels));
  % bw = (ucm2==0);
  % ucm2 = apply_sigmoid(ucm2, param.thr, param.fq);
  % ucm2(bw) = 0;

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function tests = test_imEuler2dEstimate(varargin)
% Test function for function imEuler2d
%   output = testImEuler2dEstimate(input)
%
%   Example
%   testImEuler2dEstimate
%
%   See also
%

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

tests = functiontests(localfunctions);

function testSquareInFourImages(testCase)

% create structure that does not touch borders
img = false(20, 20);
img(5:15, 5:15) = true;

% real EPC of the image
chi = imEuler2d(img);

% estimate EPC in each image using edge correction
chi1 = imEuler2dEstimate(img(1:10, 1:10));
chi2 = imEuler2dEstimate(img(1:10, 10:20));
chi3 = imEuler2dEstimate(img(10:20, 1:10));
chi4 = imEuler2dEstimate(img(10:20, 10:20));

% sum of 4 estimates
chis = chi1 + chi2 + chi3 + chi4;

assertEqual(testCase, chi, chis);

function testLabels(testCase)

% create a label image with 3 labels
img = zeros(10, 10);
img(2:3, 2:3) = 3;
img(6:8, 2:3) = 5;
img(3:5, 5:8) = 9;
img(4, 6) = 0;

[chi, labels] = imEuler2dEstimate(img);

assertElementsAlmostEqual(chi, [1 1 0]');
assertElementsAlmostEqual(labels, [3 5 9]');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [D,x0,xv,idx] = driving_function_mono_localwfs_vss(x0,xs,src,f,conf)
%DRIVING_FUNCTION_MONO_LOCALWFS_VSS driving signal for local WFS using virtual
%secondary sources
%
%   Usage: [D,xv,x0,idx] = driving_function_mono_localwfs_vss(x0,xs,src,f,conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [nx7]
%       xs          - position of virtual source or direction of plane
%                     wave / m [1x3]
%       src         - source type of the virtual source
%                         'pw' - plane wave (xs is the direction of the
%                                plane wave in this case)
%                         'ps' - point source
%                         'ls' - line source
%                         'fs' - focused source
%
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%       x0          - position, direction, and weights of the real secondary
%                     sources / m [nx7]
%       xv          - position, direction, and weights of the virtual secondary
%                     sources / m [mx7]
%       idx         - index of the selected sources from the original x0
%                     matrix [mx1]
%
%   See also: plot_sound_field, sound_field_mono_wfs
%
%   References:
%       Spors (2010) - "Local Sound Field Synthesis by Virtual Secondary
%       Sources", 40th Conference of the Audio Engineering Society, Paper 6-3,
%       http://www.aes.org/e-lib/browse.cfm?elib=15561

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargsecondarysource(x0);
isargxs(xs);
isargpositivescalar(f);
isargchar(src);
isargstruct(conf);

%% ===== Configuration ==================================================
virtualconf = conf;
virtualconf.usetapwin = conf.localwfs_vss.usetapwin;
virtualconf.tapwinlen = conf.localwfs_vss.tapwinlen;
virtualconf.secondary_sources.size = conf.localwfs_vss.size;
virtualconf.secondary_sources.center = conf.localwfs_vss.center;
virtualconf.secondary_sources.geometry = conf.localwfs_vss.geometry;
virtualconf.secondary_sources.number = conf.localwfs_vss.number;
virtualconf.wfs = conf.localwfs_vss.wfs;
virtualconf.nfchoa = conf.localwfs_vss.nfchoa;
virtualconf.driving_functions = conf.localwfs_vss.driving_functions;
method = conf.localwfs_vss.method;

%% ===== Computation ====================================================

% Determine driving functions of virtual array with different sfs methods
switch method
case 'wfs'
    % === Wave Field Synthesis ===
    % Create virtual source array
    xv = virtual_secondary_source_positions(x0,xs,src,conf);
    % Secondary_source_selection
    xv = secondary_source_selection(xv, xs, src);
    % Optional tapering
    xv = secondary_source_tapering(xv,virtualconf);
    % Driving functions for virtual source array
    Dv = driving_function_mono_wfs(xv,xs,src,f,virtualconf);
case 'nfchoa'
    % === Near-Field-Compensated Higher Order Ambisonics ===
    % Create virtual source array
    xv = secondary_source_positions(virtualconf);
    % Driving functions for virtual source array
    Dv = driving_function_mono_nfchoa(xv,xs,src,f,virtualconf);
otherwise
    error('%s: %s is not a supported method for localsfs!',upper(mfilename),method);
end

% Select secondary sources
[x0,idx] = secondary_source_selection(x0,xv(:,1:6),'vss');
% Driving functions for real source array
D = driving_function_mono_wfs_vss(x0,xv,'fs',Dv,f,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function viewGetDict = mrViewGetDictionarySetup(viewGetDict)
%
%
% USAGE:
%   This function is called upon initialization / first use to create
%   create the dictionary hash map for viewGet. This dictionary maps all of
%   the possible input combinations to the unique key that pertains to that
%   value.
%
% INPUT:
%   valueToTranslate - the orientation string, e.g. 'RAS'
%   searchString - the one character string to search for, e.g. 'S'
%   Note: the search string needs to be one of 'L','R','A','P','I','S'
%
% OUTPUT:
%   Number such that vectorString[number] = searchString
%   Note: if there are multiple occurences, then only the first one is
%   returned, along with a warning
%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function b = p01_b ( m )

%*****************************************************************************80
%
%% P01_B returns the right hand side B for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of equations.
%
%    Output, real B(M,1), the right hand side.
%
  b = [ 1.0; 2.3; 4.6; 3.1; 1.2 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ a, b ] = p05_lim ( )

%*****************************************************************************80
%
%% P05_LIM returns the integration limits for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(2), B(2), the lower and upper limits of integration.
%
  a(1:2) = 0.0;
  b(1:2) = 1.0;

  return
end


text:
function res = mig(a)
%MIG          Mignitude of point matrix (for completeness)
%
%   res = mig(a)
%

% written  11/23/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  res = abs(a);

repo: OpenMAS


text:
function [Intensity T1vector IntensityBeforeComb]=MP2RAGE_lookuptable(nimages,MPRAGE_tr,invtimesAB,flipangleABdegree,nZslices,FLASH_tr,sequence,varargin)
% first extra parameter is the inversion efficiency
% second extra parameter is the alldata
%   if ==1 all data is shown
%   if ==0 only the monotonic part is shown
alldata=0;
if nargin >=9
    if ~isempty(varargin{2})
        alldata=varargin{2};
    end
end

invtimesa=invtimesAB(1);
invtimesb=invtimesAB(2);

B1vector=1;
flipanglea=flipangleABdegree(1);
flipangleb=flipangleABdegree(2);

T1vector=0.05:0.05:5;

if length(nZslices)==2
    nZ_bef=nZslices(1);
    nZ_aft=nZslices(2);
    nZslices2=(nZslices);

    nZslices=sum(nZslices);
elseif     length(nZslices)==1
    nZ_bef=nZslices/2;
    nZ_aft=nZslices/2;
    nZslices2=(nZslices);
end;

j=0;
for T1=T1vector
    j=j+1;
    for MPRAGEtr=MPRAGE_tr
        for inversiontimesa=invtimesa
            for inversiontimesb=invtimesb
                m=0;
                for B1=B1vector
                    m=m+1;
                    %                 T1=1;%testline
                    inversiontimes2=[inversiontimesa inversiontimesb];
                    if and(and((diff(inversiontimes2))>=nZslices*FLASH_tr,inversiontimesa>=nZ_bef*FLASH_tr),inversiontimesb<=(MPRAGEtr-nZ_aft*FLASH_tr));
                        if nargin == 7
                            Signal(j,m,1:2)=1*MPRAGEfunc(nimages,MPRAGEtr,inversiontimes2,nZslices2,FLASH_tr,B1*[flipanglea flipangleb],sequence,T1);
                        else
                            if ~isempty(varargin{1})
                                Signal(j,m,1:2)=1*MPRAGEfunc(nimages,MPRAGEtr,inversiontimes2,nZslices2,FLASH_tr,B1*[flipanglea flipangleb],sequence,T1,varargin{1});
                            else
                                Signal(j,m,1:2)=1*MPRAGEfunc(nimages,MPRAGEtr,inversiontimes2,nZslices2,FLASH_tr,B1*[flipanglea flipangleb],sequence,T1);

                            end;
                        end;
                    else
                        Signal(j,m,1:2)=0;
                    end;
                end;
            end;
        end;
    end;
end;
Intensity=squeeze(real(Signal(:,:,1).*conj(Signal(:,:,2)))./(abs(Signal(:,:,1)).^2+abs(Signal(:,:,2)).^2));
T1vector=squeeze(T1vector);
if alldata==0
    [a minindex]=max(Intensity);
    [a maxindex]=min(Intensity);
    Intensity=Intensity(minindex:maxindex);
    T1vector=T1vector(minindex:maxindex);
    IntensityBeforeComb=squeeze(Signal(minindex:maxindex,1,:));
    Intensity([1 end])=[0.5 -0.5]; % pads the look up table to avoid points that fall out ot the lookuptable
else
    Intensity=squeeze(real(Signal(:,:,1).*conj(Signal(:,:,2)))./(abs(Signal(:,:,1)).^2+abs(Signal(:,:,2)).^2));
    T1vector=squeeze(T1vector);
    IntensityBeforeComb=squeeze(Signal(:,1,:));
end;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% ReadSegyTraceHeaderValue : Read a spedicifc trace header value
%
% Call:
%    % By Name
%    cdp=ReadSegyTraceHeaderValue(filename,'key','cdp');
%    SourceX=ReadSegyTraceHeaderValue(filename,'key','SourceX');
%    SourceY=ReadSegyTraceHeaderValue(filename,'key','SourceY');
%
%    % By location in Trace Header
%    SourceX=ReadSegyTraceHeaderValue(filename,'pos',72,'precision','int32');
%
%    % Call 'TraceHeaderDef(1)' to see a list of TraceHeader 'key' names
%
% See also WriteSegyTraceHeaderValue, TraceHeaderDef
%

%

function hval=ReadSegyTraceHeaderValue(filename,varargin);%pos,type)

pos=0;
precision='int32';

ninput=nargin;
% TRANSFORM VARARGING INTO PARAMETERS
cargin=1;
while (cargin<ninput)

    % ENDIAN FORMAT
    endian='ieee-be'; % Big Endian is default
    if strcmp(varargin{cargin},'endian')
       cargin=cargin+1;
       eval(['endian_tight=varargin{cargin};'])
       if endian_tight=='l',
         disp(['USING LITTLE ENDIAN TYPE'])
         endian='ieee-le';
       else
         disp(['USING BIG ENDIAN TYPE'])
       end
    end

    if strcmp(varargin{cargin},'pos')
       cargin=cargin+1;
       eval(['pos=',num2str(varargin{cargin}),';']);
       disp(['Reading at header postision : pos=',num2str(pos+1)])
    end

    if strcmp(varargin{cargin},'precision')
       cargin=cargin+1;
       eval(['precision=''',varargin{cargin},''';']);
       disp(['precision : ',precision])
    end

    if strcmp(varargin{cargin},'key')
       cargin=cargin+1;
       eval(['key=''',varargin{cargin},''';']);

       STH=TraceHeaderDef;
       try
            pos=STH.(key).pos;
            precision=STH.(key).precision;
       catch
           disp(sprintf('Trace Header Value %s not defined',key))
           pos=0; precision='int32';
           hval=[];
           return
       end
       SegymatVerbose(sprintf('key=%s, startpos=%d, precision=%s ',key,pos+1,precision))
    end

    cargin=cargin+1;
end

if nargin<2
  pos=0;
end

if nargin<3
  type='int32';
end

if exist('endian')==1,
  segyid = fopen(filename,'r',endian);
else
  segyid = fopen(filename,'r','ieee-be');  % ALL DISK FILES ARE IN BIG
end                                        % ENDIAN FORMAT, ACCORDING TO

[SegyHeader]=ReadSegyHeader(filename);

Revision=SegyHeader.SegyFormatRevisionNumber;
if Revision>0, Revision=1; end
if (SegyHeader.DataSampleFormat>length(SegyHeader.Rev(Revision+1).DataSampleFormat));
    SegymatVerbose([mfilename,' : WARNING : YOU HAVE SELECTED (OR THE FILE IS FORMATTED SUCH THAT) A DATASAMPLE FORMAT THAT IS NOT DEFINED. \nREMEBER IEEE IS NOT SPECIFIED IN THE SEGY REV0 STANDARD !'])
    if (Revision==0)
        SegymatVerbose([mfilename,' : TRYING TO USE REVISION 1 AS OPPOSED TO REVISION 0'])
        Revision=1;
        if (SegyHeader.DataSampleFormat>length(SegyHeader.Rev(Revision+1).DataSampleFormat));
            SegymatVerbose([mfilename,' : FATAL ERROR : STILL THE DATASAMPLE FORMAT IS NOT SUPPRTED - EXITING (Report error to tmh@gfy.ku.dk)'])
        else
            SegymatVerbose([mfilename,' : APPARENT SUCCES CHANING FROM Revision 0 to 1 - Continuing'])
            SegyHeader.SegyFormatRevisionNumber=1; % FORCING REVISION TO BE 1 !!!
        end
    end
end

FormatName=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).name;
Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;
BPS=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;
txt=['SegyRevision ',sprintf('%0.4g',Revision),', ',FormatName,'(',num2str(SegyHeader.DataSampleFormat),')'];

%Revision=SegyHeader.SegyFormatRevisionNumber;
%if Revision>0, Revision=1; end
%FormatName=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).name;
%Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;
%BPS=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;
%txt=['SegyRevision ',sprintf('%0.4g',Revision),',
%',FormatName,'(',num2str(SegyHeader.DataSampleFormat),')'];%

fseek(segyid,0,'eof'); DataEnd=ftell(segyid);

DataStart=3600+3200*SegyHeader.NumberOfExtTextualHeaders;
fseek(segyid,DataStart,'bof');       % Go to the beginning of the file

ntraces=(DataEnd-DataStart)./(240+(SegyHeader.ns)*(BPS/8));

hval=zeros(1,ntraces);
for itrace=1:ntraces
        if ((itrace/10000)==round(itrace/10000))
	  progress_txt(itrace,ntraces,'Trace #')
 	end

	%GOTO START OF TRACE HEADER
 	skip=DataStart+(itrace-1)*(240+(BPS/8)*SegyHeader.ns);
	fseek(segyid,skip,'bof');
	fseek(segyid,pos,'cof');
	hval(itrace)=fread(segyid,1,precision);

end
fclose(segyid);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function [sF,sDCM] = VBA_spm_split_dcm(DCM,nsplits)
% split DCM inversion using time series partition
% function [sF,sDCM] = spm_split_dcm(DCM,nsplits)
% IN:
%   - DCM: the original DCM structure
%   - nsplits: the number of splits in the partition
% OUT:
%   - sF: log-evidences for each split in the partition
%   - sDCM: cell array of DCM structures for each split

try; nsplits; catch; nsplits=4; end

% 1- split DCM structure
[sDCM] = split_DCM(DCM,nsplits);

% 2- invert split-models
sF = zeros(nsplits,1);
for i=1:nsplits
    try
        % fprintf(1,'confound removal...')
        nq = size(sDCM{i}.Y.X0,2);
        sDCM{i}.Y.X0 = sDCM{i}.Y.X0(ceil(nq./nsplits),:);
        sDCM{i}.Y.X0 = ones(sDCM{i}.v,1);
        % sDCM{i}.delays =0.*sDCM{i}.delays;
        [sDCM{i}] = VBA_spm_dcm_estimate(sDCM{i});
        sF(i) = sDCM{i}.F;
        fprintf(1,' OK.\n')
    catch
        sF(i) = NaN;
    end
end

function [splitDCM] = split_DCM(DCM,nsplits)

if nsplits==1
    splitDCM{1} = DCM;
    return
end

splitDCM = cell(nsplits,1);

nt = DCM.v;
subnt = floor(nt./nsplits);

ntu = size(DCM.U.u,1);
subntu = floor(ntu./nsplits);

for i=1:nsplits
    splitDCM{i} = DCM;
    in = subnt*(i-1)+1:subnt*i;
    splitDCM{i}.v = subnt;
    splitDCM{i}.Y.y = splitDCM{i}.Y.y(in,:);
    splitDCM{i}.Y.X0 = splitDCM{i}.Y.X0(in,:);
    nq = length(DCM.Y.Q);
    %     splitDCM{i}.Y.Q = eye(subnt);
    for j=1:nq
        tmp = zeros(nq*subnt,1);
        tmp(subnt*(j-1)+1:subnt*j) = 1;
        splitDCM{i}.Y.Q{j} = diag(tmp);
    end
    inu = subntu*(i-1)+1:subntu*i;
    splitDCM{i}.U.u = splitDCM{i}.U.u(inu,:);
end

author: MBB-team

repo: VBA-toolbox


text:
classdef RWMOP7 < PROBLEM
% <multi> <real> <constrained>
% Gear train design problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- 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
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 4;
            obj.lower    = [11.51,11.51,11.51,11.51];
            obj.upper    = [60.49,60.49,60.49,60.49];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1);
            x2 = x(:,2);
            x3 = x(:,3);
            x4 = x(:,4);
            % Objective function
            f(:,1) = abs(6.931-x3.*x4./(x1.*x2));
            f(:,2) = max(x,[],2);
            % Constraints
            g(:,1) = f(:,1)./6.931-0.5;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
         %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [3.4655000e+00   4.5417795e+01];
        end
    end
end


text:
function varargout = isempty(varargin)
%ISEMPTY   True for empty CHEBFUN2.
%   ISEMPTY(F) returns 1 if F is an empty CHEBFUN2 object and 0 otherwise.

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

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

end


text:
function plotRow(asObj, pos)

currImg = asObj.getSelectedImages(false);
figure;
cPlot(currImg(pos(1),:)');

% create plot title
title(['Row ', num2str(pos(1))]);

% create plot figure title
set(gcf,'name',['Row-Plot: ',asObj.getFigureTitle()]);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function mtrFigureDOCCSagProj(machineDir, imageDir)

subjVec = {'ss040804','mho040625','bg040719','md040714'};
paramNames = {'kLength','kSmooth','kMidSD'};

threshVec = [1000 500];
params = [-2 18 0.175];
pathsRootFile = ['paths_100k_5k_kSmooth_' num2str(params(2)) '_kLength_' num2str(params(1)) '_kMidSD_' num2str(params(3))];
pathsFile = [pathsRootFile '.dat'];

for ss = 1:length(subjVec)
    subjDir = fullfile(machineDir,subjVec{ss});

    combSTTFile = 'sttDOCCComb.nii.gz';
    fgDir = fullfile(subjDir, 'fibers');
    disp(['cd ' fgDir]);
    cd(fgDir);
    f = dir('paths_STT_LDOCC*_fd_image.nii.gz');
    if isempty(f)
        disp('Density file does not exist, creating it ...');
        mtrComputeManyFiberDensities(subjDir,'paths_STT_LDOCC.dat',3000);
    end
    sttDensityImage = dir('paths_STT_LDOCC*_fd_image.nii.gz');
    niL = niftiRead(sttDensityImage(1).name);

    f = dir('paths_STT_RDOCC*_fd_image.nii.gz');
    if isempty(f)
        disp('Density file does not exist, creating it ...');
        mtrComputeManyFiberDensities(subjDir,'paths_STT_RDOCC.dat',3000);
    end
    sttDensityImage = dir('paths_STT_RDOCC*_fd_image.nii.gz');
    niR = niftiRead(sttDensityImage(1).name);

    imgComb = double(niL.data>0) + double(niR.data>0)*2;
    disp(['Writing ' fullfile(subjDir,imageDir,combSTTFile) ' ...']);
    dtiWriteNiftiWrapper(imgComb,niL.qto_xyz,fullfile(subjDir,imageDir,combSTTFile));

    for tt = threshVec
        densityFile = [pathsRootFile '_thresh_' num2str(tt) '_fd_image.nii.gz'];
        combFile = ['ctDOCCComb_kLength_' num2str(params(1)) '_thresh_' num2str(tt) '.nii.gz'];

        fgDir = fullfile(subjDir, 'conTrack/resamp_LDOCC');
        disp(['cd ' fgDir]);
        cd(fgDir);
        % See if density file exists
        f = dir(densityFile);
        if isempty(f)
            disp('Density file does not exist, creating it ...');
            mtrComputeManyFiberDensities(subjDir,pathsFile,threshVec);
        end
        niL = niftiRead(densityFile);

        fgDir = fullfile(subjDir,'conTrack/resamp_RDOCC');
        disp(['cd ' fgDir]);
        cd(fgDir);
        % See if density file exists
        f = dir(densityFile);
        if isempty(f)
            mtrComputeManyFiberDensities(subjDir,pathsFile,threshVec);
        end
        niR = niftiRead(densityFile);

        imgComb = double(niL.data>0) + double(niR.data>0)*2;
        disp(['Writing ' fullfile(subjDir,imageDir,combFile) ' ...']);
        dtiWriteNiftiWrapper(imgComb,niL.qto_xyz,fullfile(subjDir,imageDir,combFile));
    end
end