text:
function [cfg] = ft_clusterplot(cfg, stat)

% FT_CLUSTERPLOT plots a series of topographies with highlighted clusters.
%
% Use as
%   ft_clusterplot(cfg, stat)
% where the input data is obtained from FT_TIMELOCKSTATISTICS or FT_FREQSTATISTICS.
%
% The configuration options can be
%   cfg.alpha                   = number, highest cluster p-value to be plotted max 0.3 (default = 0.05)
%   cfg.highlightseries         = 1x5 cell-array, highlight option series  with 'on', 'labels' or 'numbers' (default {'on', 'on', 'on', 'on', 'on'} for p < [0.01 0.05 0.1 0.2 0.3]
%   cfg.highlightsymbolseries   = 1x5 vector, highlight marker symbol series (default ['*', 'x', '+', 'o', '.'] for p < [0.01 0.05 0.1 0.2 0.3]
%   cfg.highlightsizeseries     = 1x5 vector, highlight marker size series   (default [6 6 6 6 6] for p < [0.01 0.05 0.1 0.2 0.3])
%   cfg.highlightcolorpos       = color of highlight marker for positive clusters (default = [0 0 0])
%   cfg.highlightcolorneg       = color of highlight marker for negative clusters (default = [0 0 0])
%   cfg.subplotsize             = layout of subplots ([h w], default [3 5])
%   cfg.saveaspng               = string, filename of the output figures (default = 'no')
%   cfg.visible                 = string, 'on' or 'off' whether figure will be visible (default = 'on')
%   cfg.position                = location and size of the figure, specified as [left bottom width height] (default is automatic)
%   cfg.renderer                = string, 'opengl', 'zbuffer', 'painters', see RENDERERINFO (default is automatic, try 'painters' when it crashes)
%   cfg.toi                     = vector, or 'all' (default) indicates which time
%                                 points (or frequency bins) are to be plotted. If specified as 'all' only the
%                                 data points with identified clusters are plotted
%
% You can also specify most configuration options that apply to FT_TOPOPLOTER or FT_TOPOPLOTTFR,
% except for cfg.xlim, any of the highlight options, cfg.comment and cfg.commentpos.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_TOPOPLOTTFR, FT_TOPOPLOTER, FT_MOVIEPLOTTFR, FT_MOVIEPLOTER

% Copyright (C) 2007-2020, F.C. Donders Centre, Ingrid Nieuwenhuis
% Copyright (C) 2021, Donders Centre for Cognitive Neuroimaging, Ingrid Nieuwenhuis and Jan-Mathijs Schoffelen
%
% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar stat
ft_preamble provenance stat

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

ws = ft_warning('off', 'FieldTrip:getdimord:warning_dimord_could_not_be_determined');

% check if the input data is valid for this function
stat = ft_checkdata(stat, 'datatype', {'timelock', 'freq'}, 'feedback', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',     {'hlmarkerseries',       'highlightsymbolseries'});
cfg = ft_checkconfig(cfg, 'renamed',     {'hlmarkersizeseries',   'highlightsizeseries'});
cfg = ft_checkconfig(cfg, 'renamed',     {'hlcolorpos',           'highlightcolorpos'});
cfg = ft_checkconfig(cfg, 'renamed',     {'hlcolorneg',           'highlightcolorneg'});
cfg = ft_checkconfig(cfg, 'renamed',     {'zparam', 'parameter'});
cfg = ft_checkconfig(cfg, 'deprecated',  {'hllinewidthseries'});
cfg = ft_checkconfig(cfg, 'deprecated',  {'xparam', 'yparam'});
cfg = ft_checkconfig(cfg, 'renamed',     {'newfigure', 'figure'});

% added several forbidden options
cfg = ft_checkconfig(cfg, 'forbidden',  {'highlight', 'highlightchannel', ...
  'highlightsymbol', 'highlightcolor', 'highlightsize', 'highlightfontsize', ...
  'xlim', 'comment', 'commentpos'});

% set the defaults
cfg.highlightseries         = ft_getopt(cfg, 'highlightseries',         {'on', 'on', 'on', 'on', 'on'});
cfg.highlightsymbolseries   = ft_getopt(cfg, 'highlightsymbolseries',   ['*', 'x', '+', 'o', '.']);
cfg.highlightsizeseries     = ft_getopt(cfg, 'highlightsizeseries',     [6 6 6 6 6]);
cfg.hllinewidthseries       = ft_getopt(cfg, 'hllinewidthseries',       [1 1 1 1 1]);
cfg.highlightfontsizeseries = ft_getopt(cfg, 'highlightfontsizeseries', [8 8 8 8 8]);
cfg.highlightcolorpos       = ft_getopt(cfg, 'highlightcolorpos',       [0 0 0]);
cfg.highlightcolorneg       = ft_getopt(cfg, 'highlightcolorneg',       [0 0 0]);
cfg.marker                  = ft_getopt(cfg, 'marker',                  'off');
cfg.alpha                   = ft_getopt(cfg, 'alpha',                   0.05);
cfg.parameter               = ft_getopt(cfg, 'parameter',               'stat');
cfg.saveaspng               = ft_getopt(cfg, 'saveaspng',               'no');
cfg.subplotsize             = ft_getopt(cfg, 'subplotsize',             [3 5]);
cfg.feedback                = ft_getopt(cfg, 'feedback',                'text');
cfg.visible                 = ft_getopt(cfg, 'visible',                 'on');
cfg.renderer                = ft_getopt(cfg, 'renderer',                []); % let MATLAB decide on the default
cfg.toi                     = ft_getopt(cfg, 'toi',                     'all');

% error if cfg.highlightseries is not a cell, for possible confusion with cfg-options
if ~iscell(cfg.highlightseries)
  ft_error('cfg.highlightseries should be a cell-array of strings')
end

% get the options that are specific for topoplotting
cfgtopo = keepfields(cfg, {'parameter', 'marker', 'markersymbol', 'markercolor', 'markersize', 'markerfontsize', 'style', 'gridscale', 'interplimits', 'interpolation', 'contournum', 'colorbar', 'shading', 'zlim'});

% prepare the layout, this only has to be done once
tmpcfg = keepfields(cfg, {'layout', 'channel', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
cfgtopo.layout = ft_prepare_layout(tmpcfg, stat);
cfgtopo.showcallinfo = 'no';
cfgtopo.feedback = 'no';

% handle with the data, it should be 1D or 2D
dimord = getdimord(stat, cfg.parameter);
dimtok = tokenize(dimord, '_');
dimsiz = getdimsiz(stat, cfg.parameter, numel(dimtok));

switch dimord
  case 'chan'
    is2D = false;
  case {'chan_time' 'chan_freq'}
    is2D = true;
  case 'chan_freq_time'
    % no more than two dimensions are supported, we can ignore singleton dimensions
    is2D = true;
    if dimsiz(2)==1
      tmpcfg = [];
      tmpcfg.avgoverfreq = 'yes';
      tmpcfg.keepfreqdim = 'no';
      tmpcfg.showcallinfo = 'no';
      stat = ft_selectdata(tmpcfg, stat);
    elseif dimsiz(3)==1
      tmpcfg = [];
      tmpcfg.avgovertime = 'yes';
      tmpcfg.keeptimedim = 'no';
      tmpcfg.showcallinfo = 'no';
      stat = ft_selectdata(tmpcfg, stat);
    else
      ft_error('if you input time-frequency data this only works if either frequency or time is a singleton dimension');
    end

  otherwise
    ft_error('unsupported dimord %s', dimord);
end % switch dimord

% these are not valid any more
clear dimord dimsiz

% this determines the labels in the figure
hastime = isfield(stat, 'time');
hasfreq = isfield(stat, 'freq');

% use the vector time, even though the 2nd dimension might be freq
if hastime
  time = stat.time;
elseif hasfreq
  time = stat.freq;
end

if issubfield(stat, 'cfg.correcttail') && ((strcmp(stat.cfg.correcttail, 'alpha') || strcmp(stat.cfg.correcttail, 'prob')) && (stat.cfg.tail == 0));
  if ~(cfg.alpha >= stat.cfg.alpha)
    ft_warning(['the pvalue you plot: cfg.alpha = ' num2str(cfg.alpha) ' is higher than the correcttail option you tested: stat.cfg.alpha = ' num2str(stat.cfg.alpha)]);
  end
end

% find significant clusters
haspos = isfield(stat, 'posclusters');
hasneg = isfield(stat, 'negclusters');

if haspos == 0 && hasneg == 0
  fprintf('%s\n', 'no clusters exceeded the nominal threshold in the data; nothing to plot')
else
  if haspos
    probpos = [stat.posclusters.prob];
    sigpos  = probpos < cfg.alpha;
  else
    sigpos  = [];
  end
  if hasneg
    probneg = [stat.negclusters.prob];
    signeg  = probneg < cfg.alpha;
  else
    signeg  = [];
  end
  sigpos  = find(sigpos == 1);
  signeg  = find(signeg == 1);
  Nsigpos = length(sigpos);
  Nsigneg = length(signeg);
  Nsigall = Nsigpos + Nsigneg;

  if Nsigall == 0
    ft_error('no clusters present with a p-value lower than the specified alpha, nothing to plot')
  end

  % make clusterslabel matrix per significant cluster
  if haspos
    posCLM = stat.posclusterslabelmat;
    sigposCLM = zeros(size(posCLM));
    for iPos = 1:Nsigpos
      sigposCLM(:,:,iPos) = (posCLM == sigpos(iPos));
      hlsignpos(iPos) = prob2hlsign(probpos(iPos), cfg.highlightsymbolseries);
    end
  else
    sigposCLM = [];
    probpos = [];
  end

  if hasneg
    negCLM = stat.negclusterslabelmat;
    signegCLM = zeros(size(negCLM));
    for iNeg = 1:Nsigneg
      signegCLM(:,:,iNeg) = (negCLM == signeg(iNeg));
      hlsignneg(iNeg) = prob2hlsign(probneg(iNeg), cfg.highlightsymbolseries);
    end
  else % no negative clusters
    signegCLM = [];
    probneg = [];
  end

  fprintf('There are %d clusters smaller than alpha (%g)\n', Nsigall, cfg.alpha);

  if is2D
    % define time or freq window per cluster
    for iPos = 1:Nsigpos
      possum_perclus = sum(sigposCLM(:,:,iPos),1); %sum over chans for each time- or freq-point
      ind_min = find(possum_perclus~=0, 1 );
      ind_max = find(possum_perclus~=0, 1, 'last' );
      time_perclus = [time(ind_min) time(ind_max)];
      if hastime
        fprintf('%s%s%s%s%s%s%s%s%s%s%s\n', 'Positive cluster: ',num2str(sigpos(iPos)), ', pvalue: ',num2str(probpos(iPos)), ' (',hlsignpos(iPos), ')', ', t = ',num2str(time_perclus(1)), ' to ',num2str(time_perclus(2)))
      elseif hasfreq
        fprintf('%s%s%s%s%s%s%s%s%s%s%s\n', 'Positive cluster: ',num2str(sigpos(iPos)), ', pvalue: ',num2str(probpos(iPos)), ' (',hlsignpos(iPos), ')', ', f = ',num2str(time_perclus(1)), ' to ',num2str(time_perclus(2)))
      end
    end
    for iNeg = 1:Nsigneg
      negsum_perclus = sum(signegCLM(:,:,iNeg),1);
      ind_min = find(negsum_perclus~=0, 1 );
      ind_max = find(negsum_perclus~=0, 1, 'last' );
      time_perclus = [time(ind_min) time(ind_max)];
      if hastime
        time_perclus = [time(ind_min) time(ind_max)];
        fprintf('%s%s%s%s%s%s%s%s%s%s%s\n', 'Negative cluster: ',num2str(signeg(iNeg)), ', pvalue: ',num2str(probneg(iNeg)), ' (',hlsignneg(iNeg), ')', ', t = ',num2str(time_perclus(1)), ' to ',num2str(time_perclus(2)))
      elseif hasfreq
        fprintf('%s%s%s%s%s%s%s%s%s%s%s\n', 'Negative cluster: ',num2str(signeg(iNeg)), ', pvalue: ',num2str(probneg(iNeg)), ' (',hlsignneg(iNeg), ')', ', f = ',num2str(time_perclus(1)), ' to ',num2str(time_perclus(2)))
      end
    end

    % define time- or freq-window containing all significant clusters
    possum = sum(sum(sigposCLM,3),1); %sum over Chans for timevector
    negsum = sum(sum(signegCLM,3),1);

    if haspos && hasneg
      allsum = possum + negsum;
    elseif haspos
      allsum = possum;
    else
      allsum = negsum;
    end

    % first and last time points of any cluster
    ind_timewin_min = find(allsum~=0, 1);
    ind_timewin_max = find(allsum~=0, 1, 'last');
    timewin = time(ind_timewin_min:ind_timewin_max);

  else
    for iPos = 1:Nsigpos
      fprintf('%s%s%s%s%s%s%s\n', 'Positive cluster: ',num2str(sigpos(iPos)), ', pvalue: ',num2str(probpos(iPos)), ' (',hlsignpos(iPos), ')')
    end
    for iNeg = 1:Nsigneg
      fprintf('%s%s%s%s%s%s%s\n', 'Negative cluster: ',num2str(signeg(iNeg)), ', pvalue: ',num2str(probneg(iNeg)), ' (',hlsignneg(iNeg), ')')
    end
  end

  % setup highlight options for all clusters and make comment for 1D data
  compos = [];
  comneg = [];
  for iPos = 1:Nsigpos
    cfgtopo.highlight(iPos)         = prob2hlsign(probpos(iPos), cfg.highlightseries);
    cfgtopo.highlightsymbol{iPos}   = prob2hlsign(probpos(iPos), cfg.highlightsymbolseries);
    cfgtopo.highlightsize{iPos}     = prob2hlsign(probpos(iPos), cfg.highlightsizeseries);
    cfgtopo.highlightfontsize{iPos} = prob2hlsign(probpos(iPos), cfg.highlightfontsizeseries);
    cfgtopo.highlightcolor{iPos}    = cfg.highlightcolorpos;
    compos = strcat(compos,cfgtopo.highlightsymbol{iPos}, 'p=',num2str(probpos(iPos)), ' '); % make comment, only used for 1D data
  end

  for iNeg = 1:Nsigneg
    cfgtopo.highlight(Nsigpos+iNeg)         = prob2hlsign(probneg(iNeg), cfg.highlightseries);
    cfgtopo.highlightsymbol{Nsigpos+iNeg}   = prob2hlsign(probneg(iNeg), cfg.highlightsymbolseries);
    cfgtopo.highlightsize{Nsigpos+iNeg}     = prob2hlsign(probneg(iNeg), cfg.highlightsizeseries);
    cfgtopo.highlightfontsize{Nsigpos+iNeg} = prob2hlsign(probneg(iNeg), cfg.highlightfontsizeseries);
    cfgtopo.highlightcolor{Nsigpos+iNeg}    = cfg.highlightcolorneg;
    comneg = strcat(comneg,cfgtopo.highlightsymbol{Nsigpos+iNeg}, 'p=',num2str(probneg(iNeg)), ' '); % make comment, only used for 1D data
  end

  if isequal(cfg.toi, 'all')
    cfg.toi = timewin;
  end

  if is2D && isequal(cfg.toi, timewin)
    % old functionality, plot every (!) time slice
    Npl = length(cfg.toi);
    begs = ind_timewin_min-1+(1:numel(cfg.toi))';
    ends = ind_timewin_min-1+(1:numel(cfg.toi))';
  elseif is2D
    % new functionality, don't plot every time slice if requested
    Npl = length(cfg.toi);

    % compute the indices in the original time axis as center points
    ix = zeros(Npl,1);
    for k = 1:Npl
      if hastime
        ix(k) = nearest(stat.time, cfg.toi(k));
      elseif hasfreq
        ix(k) = nearest(stat.freq, cfg.toi(k));
      end
    end
    win  = ceil(mean(diff(ix))/2);
    begs = max(1,ix-win);
    ends = min(numel(time), ix+win-1);
  else
    Npl  = 1;
    begs = 1;
    ends = 1;
  end

  numSubplots = prod(cfg.subplotsize);
  Nfig = ceil(Npl/numSubplots);

  % put channel indices in list
  list = cell(1,Nsigpos+Nsigneg);
  if is2D
    for iPl = 1:Npl
      for iPos = 1:Nsigpos
        list{iPl}{iPos} = find(sum(sigposCLM(:,begs(iPl):ends(iPl),iPos),2) >= 1);
      end
      for iNeg = 1:Nsigneg
        list{iPl}{Nsigpos+iNeg} = find(sum(signegCLM(:,begs(iPl):ends(iPl),iNeg),2) >= 1);
      end
    end
  else
    for iPos = 1:Nsigpos
      list{1}{iPos} = find(sigposCLM(:,iPos) == 1);
    end
    for iNeg = 1:Nsigneg
      list{1}{Nsigpos+iNeg} = find(signegCLM(:,iNeg) == 1);
    end
  end

  count = 0;
  ft_progress('init', cfg.feedback, 'making subplots...');
  ft_progress(count/Npl, 'making subplot %d from %d', count, Npl);

  % make plots
  for iPl = 1:Nfig
    % open a new figure with the specified settings, note that here it must always open a new figure
    open_figure(keepfields(cfg, {'position', 'visible', 'renderer'}));
    if is2D
      if iPl < Nfig
        for iT = 1:numSubplots
          PlN = (iPl-1)*numSubplots + iT; % plotnumber
          cfgtopo.xlim = [time(begs(PlN)) time(ends(PlN))];
          cfgtopo.highlightchannel = list{PlN};
          if hastime
            cfgtopo.comment = strcat('time: ',num2str(mean(time(begs(PlN):(ends(PlN))))), ' s');
          elseif hasfreq
            cfgtopo.comment = strcat('freq: ',num2str(mean(time(begs(PlN):(ends(PlN))))), ' Hz');
          end
          cfgtopo.commentpos = 'title';
          cfgtopo.figure = subplot(cfg.subplotsize(1), cfg.subplotsize(2), iT);
          count = count+1;
          fprintf('making subplot %d from %d\n', count, Npl);
          ft_topoplotTFR(cfgtopo, stat);
        end
      elseif iPl == Nfig
        for iT = 1:Npl-(numSubplots*(Nfig-1))
          PlN = (iPl-1)*numSubplots + iT; % plotnumber
          cfgtopo.xlim = [time(begs(PlN)) time(ends(PlN))];
          cfgtopo.highlightchannel   = list{PlN};
          if hastime
            cfgtopo.comment = strcat('time: ',num2str(mean(time(begs(PlN):(ends(PlN))))), ' s');
          elseif hasfreq
            cfgtopo.comment = strcat('freq: ',num2str(mean(time(begs(PlN):(ends(PlN))))), ' Hz');
          end
          cfgtopo.commentpos = 'title';
          cfgtopo.figure = subplot(cfg.subplotsize(1), cfg.subplotsize(2), iT);
          count = count+1;
          fprintf('making subplot %d from %d\n', count, Npl);
          ft_topoplotTFR(cfgtopo, stat);
        end
      end
    else
      cfgtopo.highlightchannel = list{1};
      cfgtopo.comment = strcat(compos, comneg);
      cfgtopo.commentpos = 'title';
      count = count+1;
      fprintf('making subplot %d from %d\n', count, Npl);
      ft_topoplotTFR(cfgtopo, stat);
    end
    if isequal(cfg.saveaspng, 'no')
      % nothing to do
    else
      % save figure
      filename = strcat(cfg.saveaspng, '_fig', num2str(iPl));
      print(gcf, '-dpng', filename);
    end
  end
end

ft_progress('close');

% return to previous warning settings
ft_warning(ws);

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% set the figure window title
if ~isempty(dataname)
  set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), mfilename, join_str(', ', dataname)));
else
  set(gcf, 'Name', sprintf('%d: %s', double(gcf), mfilename));
end
set(gcf, 'NumberTitle', 'off');

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous stat
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sign = prob2hlsign(prob, hlsign)
if prob < 0.01
  sign = hlsign(1);
elseif prob < 0.05
  sign = hlsign(2);
elseif prob < 0.1
  sign = hlsign(3);
elseif prob < 0.2
  sign = hlsign(4);
elseif prob < 0.3
  sign = hlsign(5);
end


text:
function [R, flag, timing, x0] = random_sample_precond(A, params, b)
% [R, flag, timing] = random_sample_precond(A, params)
%
% Builds a least-squares preconditioners for A based on
% random sampling.
%
% "params" - parameters governning the method.
%    params.type - type of mixing transform. Optional values: 'DCT', 'DHT', 'WHT'.
%                  Default is DHT.
%    params.gamma - gamma * n rows will be sampled (A is m-by-n).
%                   Default is 4.
%    params.preprocess_steps - number of mixing steps to do in advance.
%                               Default is 1.
%    params.maxcond - maximum condition number of the preconditioner.
%                     Default is 1 / (5 * epsilon_machine).
%
%
% Output:
%   R - the upper triangular factor of the preconditioner.
%   flag - success flag. If this function fails then probably A is rank
%          deficient.
%   timing - statistics on the time spent on various phases.
%
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.

if (nargin < 2)
    params = struct;
end

if (~isfield(params, 'type'));
    params.type = 'DHT';
end

if (~isfield(params, 'gamma'))
    params.gamma = 4;
end

if (~isfield(params, 'preprocess_steps'))
    params.preprocess_steps = 1;
end

if (~isfield(params, 'maxcond'))
    params.maxcond = 1 / (5 * eps);
end

if (~isfield(params, 'slight_coherence'))
    params.slight_coherence = 0;
end

if (~isfield(params, 'improve_start_point'))
    params.improve_start_point = false;
end

[m, n] = size(A);
mm = m;
times = 0;

timing.qr_time = 0;
timing.condest_time = 0;
timing.sample_time = 0;
timing.frut_time = 0;

if (~params.improve_start_point)
    x0 = [];
end

htimes =  params.preprocess_steps;
while (true)

    t0 = wtime;
    for i = 1:htimes
        D = sign(randn(mm, 1));
        A = fast_unitary_transform(A, D, params.type);
        mm = size(A, 1);

        if (params.improve_start_point)
            b = fast_unitary_transform(b, D, params.type);
        end
    end
    frut_time = wtime - t0;
    disp(sprintf(['\t\tRandom unit diagonal + unitary transformation time: ' ...
                  '%.2f sec'], frut_time));
    timing.frut_time = timing.frut_time + frut_time;

    t0 = wtime;
    %t = params.gamma * n / m;
    t = params.gamma * n / mm;
    s = rand(mm, 1);
    rows = find(s < t);
    R = A(rows, :);
    sample_time = wtime - t0;
    disp(sprintf('\t\tRandom sampling time: %.2f sec', sample_time));
    timing.sample_time = timing.sample_time + sample_time;

    if (params.slight_coherence > 0)
        R = [R; max(max(R))* rand(params.slight_coherence, n)];
    end

    t0 = wtime;
    if (~params.improve_start_point)
        R = mex_dgeqrf(R);
        R = triu(R(1:n, 1:n));
    else
        [R, tau] = mex_dgeqrf(R);
        Qtb = mex_dormqr(R, tau, b(rows), 'L', 'T');
        R = triu(R(1:n, 1:n));
        x0 = R \ Qtb(1:n);
    end

    qr_time = wtime - t0;
    disp(sprintf('\t\tQR on random sample time: %.2f sec', qr_time));
    timing.qr_time = timing.qr_time + qr_time;

    % Check if complete
    t0 = wtime;
    ce = mex_dtrcon(R);
    condest_time = wtime - t0;
    disp(sprintf('\t\tCondition estimation: %.2f sec', condest_time));
    timing.condest_time = timing.condest_time + condest_time;

    times = times + 1;

    if (ce > (1/ params.maxcond))
        flag = true;
        return;
    else
        if (times <= 3)
            disp(sprintf(['\t\tFailed to produced a non singular ' ...
                          'preocnditioner... applying one more FRUT...']));
            htimes = 1;
        else
            disp(sprintf('\t\tFailed to produced a non singular preocnditioner... too many FRUTs... giving up...'));
            flag = false;
            return;
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
%% EXAMPLE 9: Multitaper-Welch.
%  Multitaper-Welch estimators provide lower variance estimates at a fixed
%  frequency resolution or higher frequency resolution at similar variance
%  compared to the standard algorithm. In this example, we retain the high
%  frequency resolution of a three block Welch estimate but significantly
%  reduce the variance of the SPOD spectrum by using 10 Slepian tapers.
%
%   References:
%     [1] O. T. Schmidt, Spectral proper orthogonal decomposition using
%         multitaper estimates, Theor. Comput. Fluid Dyn., 2022, 1-14,
%         DOI 10.1007/s00162-022-00626-x, https://rdcu.be/cUtP3
%
% O. T. Schmidt (oschmidt@ucsd.edu)
% Last revision: 5-Sep-2022

clc, clear variables
addpath('utils')
disp('Loading the entire test database might take a second...')
load(fullfile('jet_data','jetLES.mat'),'p','x','r','dt');

%   trapezoidal quadrature weights for cylindrical coordinates
intWeights = trapzWeightsPolar(r(:,1),x(1,:));

%% Standard SPOD
%   SPOD with a large block size to get a high frequency resolution and
%   resolve the low-frequency regime.
nFFT    = 2048;
nOvlp   = nFFT/2;
[L,P,f] = spod(p,nFFT,intWeights,nOvlp,dt);

%   Plot the SPOD spectrum and three leading modes at the frequency of
%   interest.
f_plot  = 0.24;
[~,fi]  = min(abs(f-f_plot));
nBlk    = size(L,2);

figure
subplot(5,2,[1 3])
loglog(f,L)
xlim([f(2) f(end)]), ylim([1e-9 1e-3])
xlabel('frequency'), ylabel('SPOD mode energy')
title('Welch (standard)')
hold on
plot([f(fi) f(fi)],ylim,'k:')

count   = 5;
for mi  = 1:3
    subplot(5,2,count)
    contourf(x,r,real(squeeze(P(fi,:,:,mi))),11,'edgecolor','none'), axis equal tight, caxis(max(abs(caxis))*[-1 1])
    xlabel('x'), ylabel('r'), title(['f=' num2str(f(fi),'%.2f') ', mode ' num2str(mi)])
    xlim([0 10]); ylim([0 2])
    count = count + 2;
end
drawnow

%% Multitaper-Welch SPOD
%   SPOD using a retangular window of length 256 and 50 snaphots overlap
%   10 Slapian tapers by setting the time-halfbandwidth product to 5.5.
bw      = 5.5;
[L,P,f] = spod(p,[nFFT bw],intWeights,nOvlp,dt);

%   Plot the SPOD spectrum and modes as before. Compared to the standard
%   algorithm, the variance of the spectrum has been reduced significanlty
%   and the modes are better converged.
subplot(5,2,[2 4])
loglog(f,L(:,1:nBlk)), hold on
loglog(f,L(:,nBlk+1:end),'Color',[0.75 0.75 0.75]), hold on
xlim([f(2) f(end)]), ylim([1e-9 1e-3])
xlabel('frequency'), ylabel('SPOD mode energy')
title(['Multitaper-Welch, b_w=' num2str(bw)])

hold on
plot([f(fi) f(fi)],ylim,'k:')

count = 6;
for mi = 1:3
    subplot(5,2,count)
    contourf(x,r,real(squeeze(P(fi,:,:,mi))),11,'edgecolor','none'), axis equal tight, caxis(max(abs(caxis))*[-1 1])
    xlabel('x'), ylabel('r'), title(['f=' num2str(f(fi),'%.2f') ', mode ' num2str(mi)])
    xlim([0 10]); ylim([0 2])
    count = count + 2;
end

author: SpectralPOD


text:
%[2015]-"A human learning optimization algorithm and its application
%to multi-dimensional knapsack problems"

% (9/12/2020)

function HLO = jHumanLearningOptimization(feat,label,opts)
% Parameters
pi = 0.85;   % probability of individual learning
pr = 0.1;    % probability of exploration learning

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'pi'), pi = opts.pi; end
if isfield(opts,'pr'), pr = opts.pr; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = jInitialPopulation(N,dim);
% Fitness
fit    = zeros(1,N);
fitSKD = inf;
for i = 1:N
  fit(i) = fun(feat,label,X(i,:),opts);
  % Update SKD/gbest
  if fit(i) < fitSKD
    fitSKD = fit(i);
    SKD    = X(i,:);
  end
end
% Get IKD/pbest
fitIKD = fit;
IKD    = X;
% Pre
curve = zeros(1,max_Iter);
curve(1) = fitSKD;
t = 2;
% Generations
while t <= max_Iter
  for i = 1:N
    % Update solution (8)
    for d = 1:dim
      % Radom probability in [0,1]
      r = rand();
      if r >= 0 && r < pr
        % Random exploration learning operator (7)
        if rand() < 0.5
          X(i,d) = 0;
        else
          X(i,d) = 1;
        end
      elseif r >= pr && r < pi
        X(i,d) = IKD(i,d);
      else
        X(i,d) = SKD(d);
      end
    end
  end
  % Fitness
  for i = 1:N
    % Fitness
    fit(i) = fun(feat,label,X(i,:),opts);
    % Update IKD/pbest
    if fit(i) < fitIKD(i)
      fitIKD(i) = fit(i);
      IKD(i,:)  = X(i,:);
    end
    % Update SKD/gbest
    if fitIKD(i) < fitSKD
      fitSKD = fitIKD(i);
      SKD    = IKD(i,:);
    end
  end
  curve(t) = fitSKD;
  fprintf('\nGeneration %d Best (HLO)= %f',t,curve(t))
  t = t + 1;
end
% Select features based on selected index
Pos   = 1:dim;
Sf    = Pos(SKD == 1);
sFeat = feat(:,Sf);
% Store results
HLO.sf = Sf;
HLO.ff = sFeat;
HLO.nf = length(Sf);
HLO.c  = curve;
HLO.f  = feat;
HLO.l  = label;
end

% Binary initialization strategy
function X = jInitialPopulation(N,dim)
X = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    if rand() > 0.5
      X(i,d) = 1;
    end
  end
end
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
% ZG_rsum(X)
% row sum
%
% Machine Learning Toolbox
% Version 1.0  01-Apr-96
% Copyright (c) by Zoubin Ghahramani
% http://mlg.eng.cam.ac.uk/zoubin/software.html
%
% ------------------------------------------------------------------------------
% The MIT License (MIT)
%
% Copyright (c) 1996, Zoubin Ghahramani
%
% 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.
% ------------------------------------------------------------------------------

function Z = ZG_rsum(X)

Z = zeros(size(X(:,1)));

for i = 1:length(X(1,:))
    Z = Z + X(:,i);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function point = intersectLinePlane(line, plane, varargin)
%INTERSECTLINEPLANE Intersection point between a 3D line and a plane.
%
%   PT = intersectLinePlane(LINE, PLANE)
%   Returns the intersection point of the given line and the given plane.
%   LINE:  [x0 y0 z0 dx dy dz]
%   PLANE: [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   PT:    [xi yi zi]
%   If LINE and PLANE are parallel, return [NaN NaN NaN].
%   If LINE (or PLANE) is a matrix with 6 (or 9) columns and N rows, result
%   is an array of points with N rows and 3 columns.
%
%   PT = intersectLinePlane(LINE, PLANE, TOL)
%   Specifies the tolerance factor to test if a line is parallel to a
%   plane. Default is 1e-14.
%
%   Example
%     % define horizontal plane through origin
%     plane = [0 0 0   1 0 0   0 1 0];
%     % intersection with a vertical line
%     line = [2 3 4  0 0 1];
%     intersectLinePlane(line, plane)
%     ans =
%        2   3   0
%     % intersection with a line "parallel" to plane
%     line = [2 3 4  1 2 0];
%     intersectLinePlane(line, plane)
%     ans =
%       NaN  NaN  NaN
%
%   See also
%   lines3d, planes3d, points3d, clipLine3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-17
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

% extract tolerance if needed
tol = 1e-14;
if nargin > 2
    tol = varargin{1};
end

% unify sizes of data
nLines  = size(line, 1);
nPlanes = size(plane, 1);

% N planes and M lines not allowed
if nLines ~= nPlanes && min(nLines, nPlanes) > 1
    error('MatGeom:geom3d:intersectLinePlane', ...
        'Input must have same number of rows, or one must be 1');
end

% plane normal
n = crossProduct3d(plane(:,4:6), plane(:,7:9));

% difference between origins of plane and line
dp = bsxfun(@minus, plane(:, 1:3), line(:, 1:3));

% dot product of line direction with plane normal
denom = sum(bsxfun(@times, n, line(:,4:6)), 2);

% relative position of intersection point on line (can be inf in case of a
% line parallel to the plane)
t = sum(bsxfun(@times, n, dp),2) ./ denom;

% compute coord of intersection point
point = bsxfun(@plus, line(:,1:3),  bsxfun(@times, [t t t], line(:,4:6)));

% set indices of line and plane which are parallel to NaN
par = abs(denom) < tol;
point(par,:) = NaN;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% sets directories

dataset = 'labelme';
alg = 'occw';

switch(lower(dataset))

    case 'bsds'
        bsdsdir = '~dhoiem/data/datasets/BSDS300/';
        imdir = fullfile(bsdsdir, 'images/all/');

        basedir = '~/data/occlusion/bsds';
        gtdir = fullfile(basedir, 'gt');
        occdir = fullfile(basedir, 'occlusion');
        pbdir = fullfile(basedir, 'pb');
        pb2dir = fullfile(basedir, 'pb2');
        traindir = fullfile(basedir, 'train');

    case 'labelme'
        basedir = '~/data/occlusion/labelme';
        imdir_train = fullfile(basedir, 'train/Images');
        imdir_test = fullfile(basedir, 'test/Images');
        gtdir = fullfile(basedir, 'labels2');
        occdir = fullfile(basedir, 'occlusion');
        pbdir = fullfile(basedir, 'pb');
        pb2dir = fullfile(basedir, 'pb2');
        traindir = fullfile(basedir, 'train');

    otherwise
        error('invalid dataset')
end

author: Cloud-CV

repo: object-proposals


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function noisecgt
%NOISECGT Unit test for the function NOISECG.

%	O. Lemoine - February 1996.

N=32768;
A1=0.2; A2=0.7;
sig=noisecg(N,A1,A2);

% Mean
Mean=mean(sig);
if abs(Mean)>10/sqrt(N),
 error('noisecg test 1 failed');
end

% Variance
Var=std(sig).^2;
if abs(Var-1)>10/sqrt(N),
 error('noisecg test 2 failed');
end

% histogram
Nh=100;
[h,m]=hist(real(sig),Nh); h=h/max(h);
Nc=find(abs(m)==min(abs(m)));
pdf=amgauss(Nh,Nc,Nh/(2*sqrt(pi)))';
if any(abs(h-pdf).^2>10/sqrt(N)),
 error('noisecg test 3 failed');
end

% whiteness
sig=noisecg(N);
autocor=xcorr(sig);
Max=max(autocor);
L=length(find(abs(autocor/Max)>5e-2));
if L/N>5e-4,
 error('noisecg test 4 failed');
end

% For N=1
N=1; Np=10000;
for k=1:Np,
 sig(k)=noisecg(N);
end
Mean=mean(sig);
if abs(Mean)>10/sqrt(Np),
 error('noisecg test 5 failed');
end
Var=std(sig).^2;
if abs(Var-1)>10/sqrt(Np),
 error('noisecg test 6 failed');
end
[h,m]=hist(real(sig),Nh); h=h/max(h);
Nc=find(abs(m)==min(abs(m)));
pdf=amgauss(Nh,Nc,Nh/(2*sqrt(pi)))';
if any(abs(h-pdf).^2>10/sqrt(Np)),
 error('noisecg test 7 failed');
end

% For N=2
N=2;
for k=1:2:(Np-1),
 noise=noisecg(N);
 sig(k)=noise(1);
 sig(k+1)=noise(2);
end
Mean=mean(sig);
if abs(Mean)>10/sqrt(Np),
 error('noisecg test 8 failed');
end
Var=std(sig).^2;
if abs(Var-1)>10/sqrt(Np),
 error('noisecg test 9 failed');
end
[h,m]=hist(real(sig),Nh); h=h/max(h);
Nc=find(abs(m)==min(abs(m)));
pdf=amgauss(Nh,Nc,Nh/(2*sqrt(pi)))';
if any(abs(h-pdf).^2>10/sqrt(Np)),
 error('noisecg test 10 failed');
end

author: HeLiangHIT


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 = 1 / (2 * m) * sum( (X * theta - y) .^ 2);

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

end

author: Ayatans

repo: Machine-Learning-homework


text:
function [msh, lights] = meshVisualize(msh,id)
% Visualize a VTK mesh of the cortical surface using mrMesh
%
%  [msh] = meshVisualize(msh,[id]);
%
% This mesh structure must contain vertices and faces.
%
% A msh is input and the mesh is displayed in mrMesh window number id.  If
% no id number is sent in a new window is opened.
%
% If the input is an old style mesh, it is changed updated the best we can.
%
% Questions:
%    What are the minimal fields required for this msh to run with mrMesh?
%    Can we pass in a modified gifti mesh and convert it to the proper
%    fields and have this run?
%
% Example:
%  fName ='X:\anatomy\nakadomari\left\20050901_fixV1\left.Class';
%  fName ='/biac1/wandell/data/anatomy/dougherty/t1_class_left.nii.gz';
%  msh = meshBuildFromClass(fName);
%  msh = meshSmooth(msh);
%  meshVisualize(msh);
%
% (c) Stanford Vista Team

backColor = [1,1,1];

if notDefined('msh'), error('The mesh is required.'); end
% if ieNotDefined('mmPerVox'), mmPerVox = [1 1 1]; end
if notDefined('host'), host = 'localhost'; end
if notDefined('id'), id = -1; end

% Set initial parameters for the mesh.
if isempty(meshGet(msh,'host')), msh = meshSet(msh,'host',host); end
if isempty(meshGet(msh,'id')),   msh = meshSet(msh,'id',id); end

% If the window is already open, no harm is done.
if mrmCheckServer, mrMesh(meshGet(msh,'host'),meshGet(msh,'id'),'close'); end
msh = mrmInitHostWindow(msh);

% Initializes the mesh
[msh, lights] = mrmInitMesh(msh,backColor);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ind = subsind(grains,subs)

% dermine grains by x,y coordinates
if numel(subs)==2 && all(cellfun(@isnumeric, subs))
  ind = grains.findByLocation([subs{:}]);
  return
elseif numel(subs)==2 && ischar(subs{1}) && strcmpi(subs{1},'id')
  ind = reshape(grains.id2ind(subs{2}),size(subs{2}));
  if any(ind(:)==0)
    error('No data with the specified ids in the data set');
  end
  return
end

ind = true(length(grains),1);

for i = 1:length(subs)

  if ischar(subs{i}) && strcmpi(subs{i},'indexed')

    ind = ind & grains.isIndexed;

  elseif ischar(subs{i}) || iscellstr(subs{i})

    phases = false(length(grains.phaseMap),1);
    mineralsSubs = ensurecell(subs{i});
    phaseNumbers = cellfun(@num2str,num2cell(grains.phaseMap(:)),'Uniformoutput',false);

    for k=1:numel(mineralsSubs)
      phases = phases ...
        | strcmpi(grains.mineralList(:),mineralsSubs{k}) ...
        | strcmpi(phaseNumbers,mineralsSubs{k});
    end

    % if no complete match was found allow also for partial match
    if ~any(phases)
      for k=1:numel(mineralsSubs)
        phases = phases ...
          | strncmpi(grains.mineralList(:),mineralsSubs{k},length(mineralsSubs{k}));
      end
    end

    if ~any(phases)
      disp(' ');
      warning off backtrace
      warning(['There is no such phase "' mineralsSubs{1} '". Maybe you mispelled it?']);
      warning on backtrace
    end

    %miner = ensurecell(subs{i});
    %alt_mineral = cellfun(@num2str,num2cell(grains.phaseMap),'Uniformoutput',false);
    %for k=1:numel(miner)
    %  phases = phases | ~cellfun('isempty',regexpi(grains.mineralList(:),['^' miner{k}])) | ...
    %    strcmpi(alt_mineral(:),miner{k});
    %end
    ind = ind & phases(grains.phaseId(:));

  elseif isa(subs{i},'logical')

    sub = any(subs{i}, find(size(subs{i}')==max(size(ind)),1));

    ind = ind & reshape(sub,size(ind));

  elseif isnumeric(subs{i})

    if any(subs{i} <= 0 | subs{i} > length(grains))
      error('Out of range; index must be a positive integer or logical.')
    end

    ind = subs{i};
    return

  elseif isa(subs{i},'polygon')

    ind = ind & inpolygon(grains,subs{i})';

  elseif isa(subs{i},'crystalSymmetry')

    phaseId = grains.cs2phaseId(subs{i});
    ind = ind & grains.phaseId == phaseId;

  end
end
end

author: mtex-toolbox

repo: mtex

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


text:
% The COBRAToolbox: testROOM.m
%
% Purpose:
%     - Check the function ROOM
%
% Authors:
%     - Luis V. Valcarcel 2018-12-17
%

global CBTDIR

solversToUse = prepareTest('needsMILP', true, 'useSolversIfAvailable', {'tomlab_cplex', 'ibm_cplex', 'gurobi', 'glpk','cplex_direct'}, 'excludeSolvers', {'qpng'});
% Note: the solver QPNG cannot be used with this test

% save the current path
currentDir = pwd;

% initialize the test
testDir = fileparts(which('testROOM'));
cd(testDir);

% Create Toy example model
ReactionFormulas = {'A -> B', '2 B -> C + byp', '2 B + cof -> D + byp',...
'D -> E + cof', 'C + cof -> D', 'C -> E', ' -> A', 'E -> ', 'byp -> '};
ReactionNames = {'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'b1', 'b2', 'b3'};
lowerbounds = [0 0 0 0 0 0 0 0 0];
upperbounds = [20, 20, 20, 20, 20, 20, 20, 20, 20];
model = createModel(ReactionNames, ReactionNames, ReactionFormulas,...
'lowerBoundList', lowerbounds, 'upperBoundList', upperbounds);

% Generate reference flux
WT_flux = [10, 5, 0, 0, 0, 5, 10, 5, 5];

% Check errors when missing argument
assert(verifyCobraFunctionError('ROOM', 'inputs', {model}))
assert(verifyCobraFunctionError('ROOM', 'inputs', {model, WT_flux}))
assert(verifyCobraFunctionError('ROOM', 'inputs', {model, WT_flux, 'a'}))

% define the solver packages to be used to run this test
solverPkgs = solversToUse.MILP;

predictedROOMsol = [10 5 0 5 5 0 10 5 5]';

% Test solving for different solvers
for k = 1:length(solverPkgs)
    fprintf(' -- Running testROOM using the solver interface: %s ... ', solverPkgs{k});

    solverOK = changeCobraSolver(solverPkgs{k}, 'all', 0);

    if solverOK
        % Calculate ROOM and check solutions
        fluxROOM = ROOM(model, WT_flux, {'v6'},'delta', 1e-5, 'epsilon', 1e-5);
        assert(norm(fluxROOM-predictedROOMsol)<1e-2)

        [fluxROOM, solutionROOM, totalFluxDiff] = ROOM(model, WT_flux, {'v6'},'delta',1e-6,'epsilon',1e-6);
        assert(norm(fluxROOM-predictedROOMsol)<1e-2)
        assert(abs(totalFluxDiff^2-75)<1e-2)

        [fluxROOM2, solutionROOM2, totalFluxDiff2] = linearROOM(model, WT_flux, {'v6'},'delta',1e-6,'epsilon',1e-6);
        assert(norm(fluxROOM-predictedROOMsol)<1e-2)
        assert(abs(totalFluxDiff^2-75)<1e-2)

    else
        warning('The test testROOM cannot run using the solver interface: %s. The solver interface is not installed or not configured properly.\n', solverPkgs{k});
    end

    % output a success message
    fprintf('Done.\n');
end

% remove the file created during the test
if exist('CobraLPSolver.log','file')
    delete('CobraLPSolver.log');
end

% Set seed to default value
rng('default')
% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function output = callmaxdet(interfacedata)

[F_struc,F_blksz,G_struc,G_blksz] = sedumi2maxdet(interfacedata.F_struc,interfacedata.K);

c = interfacedata.c;
options = interfacedata.options;
x0 = interfacedata.x0;
onlyfeasible = (nnz(c)==0) & isempty(G_blksz);

if isempty(G_blksz)
    G_blksz = 1;
    G_struc = [1 zeros(1,length(c))];
else
end

if isempty(F_blksz)
    x = zeros(length(c),1)+NaN;
    problem = 11;
    infostr = yalmiperror(problems,'MAXDET failed since there are no inequality constraints in F(x)');
    solveroutput = [];
    solvertime = 0;
    D_struc = [];
end

problem = 0;
solvertimephase1=0;
D_struc = [];
if isempty(x0)
    solvertimephase1 = tic;
    showprogress('Calling MAXDET/phase1',options.showprogress);
    [x0,z0,w0,problem,infostr,solveroutput] = callmaxdetphase1(full(F_struc),full(F_blksz), full(G_struc),full(G_blksz), full(c), options);
    solvertimephase1 = toc(solvertimephase1);
end
if (problem~=0) | (onlyfeasible==1)
    if isempty(x0)
        x = repmat(nan,length(c),1);
    else
        x = x0;
    end
    solvertime = solvertimephase1;
else
    z0=zeros(size(F_struc,1),1);
    w0=zeros(size(G_struc,1),1);
    solvertime = tic;
    showprogress('Calling MAXDET',options.showprogress);
    if options.verbose==0
        evalc('[x,Z,W,ul,hist,infostr]=maxdet(full(F_struc),F_blksz,full(G_struc), G_blksz,c'',x0,z0,w0,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.gam,options.maxdet.NTiters);');
    else
        [x,Z,W,ul,hist,infostr]=maxdet(full(F_struc),F_blksz,full(G_struc), G_blksz,c',x0,z0,w0,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.gam,options.maxdet.NTiters);
    end
    solvertime = toc(solvertime)+solvertimephase1;

    D_struc = [Z;W];

    if options.savesolveroutput
        solveroutput.x = x;
        solveroutput.Z = Z;
        solveroutput.W = W;
        solveroutput.ul = ul;
        solveroutput.hist = hist;
        solveroutput.infostr = infostr;
    else
        solveroutput = [];
    end

    switch  infostr
        case {'target reached','relative accuracy reached','absolute accuracy reached','absolute tolerance reached','relative tolerance reached'}
            problem = 0;
        case 'maximum Newton iteration exceeded'
            problem = 3;
        otherwise
            problem = 1;
    end
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,interfacedata.solver.tag,[],solveroutput,solvertime);

function [x0,z0,w0,problem,infostr,solveroutput] = callmaxdetphase1(F_struc,F_blksz, G_struc,G_blksz, c, options);

try
    if options.verbose==0
        evalc('[x0,z0,w0,ul] = phase1(full(F_struc),F_blksz,full(G_struc), G_blksz,options.maxdet.gam,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.NTiters);');
    else
        [x0,z0,w0,ul] = phase1(full(F_struc),F_blksz,full(G_struc), G_blksz,options.maxdet.gam,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.NTiters);
    end
    if ul(1)<0
        problem = 0;
    elseif ul(2)>=0
        problem = 1;
    else
        problem = 8;
    end
    %  % Problems currently detected outside
%  problem = (ul(1)>0);

% In case of problems, phase1 outputs the extended vector instead
% of last primal iterate
if problem
    x0 = x0(1:length(c));
end

catch
    problem = 9;
    x0 = [];
    z0 = [];
    w0 = [];
    ul = [];
end
infostr = yalmiperror(problem,'MAXDET/phase1');
if options.savesolveroutput
    solveroutput.x = x0;
    solveroutput.Z = z0;
    solveroutput.W = w0;
    solveroutput.ul = ul;
    solveroutput.infostr = infostr;
else
    solveroutput = [];
end

% CODE TAKEN DIRECTLY FROM MAXDET/PHASE1
% THE REASON FOR HAVING A COPY HERE IS THAT
% A FILE NAMED PHASE1 EXIST FOR SP ALSO
%MAXDET, version ALPHA, April 1996.

%COPYRIGHT 1996 SHAO-PO WU, LIEVEN VANDENBERGHE AND STEPHEN BOYD
%Permission to use, copy, modify, and distribute this software for
%any purpose without fee is hereby granted, provided that this entire
%notice is included in all copies of any software which is or includes
%a copy or modification of this software and in all copies of the
%supporting documentation for such software.
%This software is being provided "as is", without any express or
%implied warranty.  In particular, the authors do not make any
%representation or warranty of any kind concerning the merchantability
%of this software or its fitness for any particular purpose.
function [x,Z,W,ul] = phase1(F,F_blkszs,G,G_blkszs,gam,abstol,reltol,NTiters);

m = size(F,2)-1;
if (m ~= size(G,2)-1)
    error('F and G must have the same number of columns.');
end
if (size(F,1) ~= sum(F_blkszs.^2))
    error('Dimensions of F do not match F_blkszs.');
end;
if (size(G,1) ~= sum(G_blkszs.^2))
    error('Dimensions of G do not match G_blkszs.');
end;

% mineigF is the smallest eigenvalue of F_0
mineigF = 0.0;
k=0; for n=F_blkszs,
    mineigF = min(mineigF, min(eig(reshape(F(k+[1:n*n],1),n,n))));
    k=k+n*n;   % k = sum n_i*n_i
end;
% mineigG is the smallest eigenvalue of G_0
mineigG = 0.0;
k=0; for n=G_blkszs,
    mineigG = min(mineigG, min(eig(reshape(G(k+[1:n*n],1),n,n))));
    k=k+n*n;   % k = sum n_i*n_i
end;

% eyeF is the identity
eyeF = zeros(size(F,1),1);
k=0; for n=F_blkszs,
    eyeF(k+[1:n*n]) = reshape(eye(n),n*n,1);   % identity
    k=k+n*n;   % k = sum n_i*n_i
end;
% eyeG is the identity
eyeG = zeros(size(G,1),1);
k=0; for n=G_blkszs,
    eyeG(k+[1:n*n]) = reshape(eye(n),n*n,1);   % identity
    k=k+n*n;   % k = sum n_i*n_i
end;

% initial x0
x0 = [zeros(m,1); max(-1.1*min(mineigF,mineigG), 1)];

% linear objective
c = [zeros(m,1); 1];

% call maxdet
[x,Z,W,ul,hist,infostr]=maxdet([F,eyeF; G,eyeG],...
    [F_blkszs(:)',G_blkszs(:)'],...
    [1 zeros(1,m+1)],1,c,x0,...
    zeros(size(F,1)+size(G,1),1),0,...
    abstol,reltol,gam,NTiters);

% prepare output
x = x(1:m);
W = Z(size(F,1)+1:size(F,1)+size(G,1));
Z = Z(1:size(F,1));
%if ul(1)<0
%  infostr = 'feasible';
%elseif ul(2)>=0
%  infostr = 'infeasible';
%else
%  infostr = 'feasibility cannot be determined';
%end

function [F_struc,blksz] = deblock(F_struc,blksz);
X = any(F_struc(end-blksz^2+1:end,:),2);
X = reshape(X,blksz,blksz);
[v,dummy,r,dummy2]=dmperm(X);
blks = diff(r);

logt = F_struc;

newlogt = [];
for i = 1:size(logt,2)
    temp = reshape(logt(:,i),blksz,blksz);
    temp = temp(v,v);
    newlogt = [newlogt temp(:)];
end
logt = newlogt;

pattern = [];
for i = 1:length(blks)
    pattern = blkdiag(pattern,ones(blks(i)));
end

F_struc = logt(find(pattern),:);
blksz = blks;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% The COBRAToolbox: testKLDis.m
%
% Purpose:
%     - test the functionality of KLDis
%
% Authors:
%     - Loic Marx
%

% save the current path
currentDir = pwd;

% initialize the test
cd(fileparts(which(mfilename)));

% set the tolerance
tol = 1e-6;

% Generate 2 matrix of the same size
P = [1 2; 2 1];
Q = [3 2; 2 3];

% run the function
KLD = KLDis(P,Q);

% compare the scaled matrix against the reference data
assert(norm(KLD - 0.12723233459191 * [1; 1]) < tol)

Q2 = [inf inf];
P2 = [inf inf];

% Test if KLDis throws an error
assert(verifyCobraFunctionError('KLDis','inputs', {P, Q(:,1)}, 'testMessage', 'the number of columns in P and Q should be the same'));
assert(verifyCobraFunctionError('KLDis','inputs', {P2, Q2}, 'testMessage', 'the inputs contain non-finite values!'));

% change the directory
cd(currentDir)


text:
function [ia,ja,sa] = find(a)
%FIND         Implements  find(a)  for Taylor
%
%   index = find(a)
%   [ia,ja] = find(a)
%   [ia,ja,sa] = find(a)
%
% functionality as Matlab function find
%

% written  05/21/09     S.M. Rump
%

  if nargout<=1
    ia = find(any(a.t));
  elseif nargout==2
    index = reshape(1:prod(a.size),a.size);
    index(~any(a.t)) = 0;
    [ia,ja] = find(index);
  else
    ia = find(any(a.t));
    %VVVV  sa = a(ia);
    s.type = '()'; s.subs = {ia}; sa = subsref(a,s);
    %AAAA  Matlab bug fix
    index = reshape(1:prod(a.size),a.size);
    index(~any(a.t)) = 0;
    [ia,ja] = find(index);
  end

repo: OpenMAS


text:
function x = r8ci_sl ( n, a, b, job )

%*****************************************************************************80
%
%% R8CI_SL solves a R8CI system.
%
%  Discussion:
%
%    The R8CI storage format is used for an N by N circulant matrix.
%    An N by N circulant matrix A has the property that the entries on
%    row I appear again on row I+1, shifted one position to the right,
%    with the final entry of row I appearing as the first of row I+1.
%    The R8CI format simply records the first row of the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2004
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N), the R8CI matrix.
%
%    Input, real B(N), the right hand side.
%
%    Input, integer JOB, specifies the system to solve.
%    0, solve A * x = b.
%    nonzero, solve A' * x = b.
%
%    Output, real X(N), the solution of the linear system.
%
  if ( job == 0 )
%
%  Solve the system with the principal minor of order 1.
%
    r1 = a(1);
    x(1) = b(1) / r1;

    r2 = 0.0;
%
%  Recurrent process for solving the system.
%
    for nsub = 2 : n
%
%  Compute multiples of the first and last columns of
%  the inverse of the principal minor of order N.
%
      r5 = a(n+2-nsub);
      r6 = a(nsub);

      if ( 2 < nsub )

        work(nsub-1) = r2;

        for i = 1 : nsub-2
          r5 = r5 + a(n+1-i) * work(nsub-i);
          r6 = r6 + a(i+1) * work(n-1+i);
        end

      end

      r2 = -r5 / r1;
      r3 = -r6 / r1;
      r1 = r1 + r5 * r3;

      if ( 2 < nsub )

        r6 = work(n);
        work(n-1+nsub-1) = 0.0;
        for i = 2 : nsub-1
          r5 = work(n-1+i);
          work(n-1+i) = work(i) * r3 + r6;
          work(i) = work(i) + r6 * r2;
          r6 = r5;
        end

      end

      work(n) = r3;
%
%  Compute the solution of the system with the principal minor of order NSUB.
%
      r5 = 0.0;
      for i = 1 : nsub-1
        r5 = r5 + a(n+1-i) * x(nsub-i);
      end

      r6 = ( b(nsub) - r5 ) / r1;
      x(1:nsub-1) = x(1:nsub-1) + work(n:n+nsub-2) * r6;
      x(nsub) = r6;

    end

  else
%
%  Solve the system with the principal minor of order 1.
%
    r1 = a(1);
    x(1) = b(1) / r1;

    r2 = 0.0;
%
%  Recurrent process for solving the system.
%
    for nsub = 2 : n
%
%  Compute multiples of the first and last columns of
%  the inverse of the principal minor of order N.
%
      r5 = a(nsub);
      r6 = a(n+2-nsub);

      if ( 2 < nsub )

        work(nsub-1) = r2;

        for i = 1 : nsub-2
          r5 = r5 + a(i+1) * work(nsub-i);
          r6 = r6 + a(n+1-i) * work(n-1+i);
        end

      end

      r2 = -r5 / r1;
      r3 = -r6 / r1;
      r1 = r1 + r5 * r3;

      if ( 2 < nsub )

        r6 = work(n);
        work(n-1+nsub-1) = 0.0;
        for i = 2 : nsub-1
          r5 = work(n-1+i);
          work(n-1+i) = work(i) * r3 + r6;
          work(i) = work(i) + r6 * r2;
          r6 = r5;
        end

      end

      work(n) = r3;
%
%  Compute the solution of the system with the principal minor of order NSUB.
%
      r5 = 0.0E+00;
      for i = 1 : nsub-1
        r5 = r5 + a(i+1) * x(nsub-i);
      end

      r6 = ( b(nsub) - r5 ) / r1;
      for i = 1 : nsub-1
        x(i) = x(i) + work(n-1+i) * r6;
      end

      x(nsub) = r6;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
  function out = full(ob, varargin)
% function out = full(ob, varargin)
%| full(ob) = ob(:,:)
%|
%| note: alternative is ob * eye(np) but this may yield the
%| wrong precision because eye() is double but ob may be single.
%|
%| optional argument(s) passed to fatrix2_subsref_colon()
%| to facilitate test_adjoint
%|
%| Copyright 2010-12-05, Jeff Fessler, University of Michigan

out = fatrix2_subsref_colon(ob, ':', varargin{:}); % ob(:,:)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function headSizeAll = getHeadSizeAll(annolist)
headSizeAll = nan(length(annolist),1);
for imgidx = 1:length(annolist)
    rect = annolist(imgidx).annorect;
    headSizeAll(imgidx) =  util_get_head_size(rect);
end
end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
%This Progra, develop by Renoald Tang
%University Teknologi Malaysia, Photogrammetry and Laser scanning group
%for academic purpose
%Email:renoald@live.com
%***For Cocone(cocone.exe) and Tight Cocone(tcocone.exe) software
%***please request download from http://www.cse.ohio-state.edu/~tamaldey/cocone.html
function varargout = Reconstruction(varargin)
% RECONSTRUCTION M-file for Reconstruction.fig
%      RECONSTRUCTION, by itself, creates a new RECONSTRUCTION or raises the existing
%      singleton*.
%
%      H = RECONSTRUCTION returns the handle to a new RECONSTRUCTION or the handle to
%      the existing singleton*.
%
%      RECONSTRUCTION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in RECONSTRUCTION.M with the given input arguments.
%
%      RECONSTRUCTION('Property','Value',...) creates a new RECONSTRUCTION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before Reconstruction_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to Reconstruction_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help Reconstruction

% Last Modified by GUIDE v2.5 10-Jan-2012 16:15:16

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Reconstruction_OpeningFcn, ...
                   'gui_OutputFcn',  @Reconstruction_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before Reconstruction is made visible.
function Reconstruction_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to Reconstruction (see VARARGIN)

% Choose default command line output for Reconstruction
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes Reconstruction wait for user response (see UIRESUME)
% uiwait(handles.figure1);
axes(handles.axes1)
%axis off
view(3)
grid on
axes(handles.axes2)
view(3)
grid on

% --- Outputs from this function are returned to the command line.
function varargout = Reconstruction_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in importbutton.
function importbutton_Callback(hObject, eventdata, handles)
% hObject    handle to importbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%Let the user select the file
[FileName,PathName] = uigetfile({'*.xyz'},'Select the Point clouds');
%combine the filename and pathname
directory=strcat(PathName ,FileName);
%.xyz check file extension
ext=findstr(FileName,'.xyz');
n=isempty(ext);
if n==0
    point=dlmread(directory);

end
dlmwrite('point', point, 'delimiter', '\t','precision', 6);
n2=length(point);
num=num2str(n2);
set(handles.text2,'String',num)
 setappdata(Reconstruction,'point1',point);
 axes(handles.axes1)
plot3(point(:,1),point(:,2),point(:,3),'r.')
grid on

% --- Executes on button press in DVbutton.
function DVbutton_Callback(hObject, eventdata, handles)
% hObject    handle to DVbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
p=getappdata(Reconstruction,'point1');
 [t,tnorm]=Crust(p);
 n=length(t);
 x=ones(n,1);

 num1=num2str(n);
axes(handles.axes2);
trisurf(t,p(:,1),p(:,2),p(:,3));
axis on
grid on
set(handles.text4,'String',num1);
 for i=1 : n
     x(i,1)=x(i,1)*3;
     t(i,1)=t(i,1)-1;
     t(i,2)=t(i,2)-1;
     t(i,3)=t(i,3)-1;
 end
t=[x t];
dlmwrite('face', t, 'delimiter', '\t','precision', 6);

% --- Executes on button press in Cbutton.
function Cbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Cbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
a1 = get(handles.C_edit,'String');
C=str2num(a1);
R1=isempty(C);
if R1==1
    C=0.393;
end
b=get(handles.S_edit,'String');
S=str2num(b);
R1=isempty(S);
if R1==1
    S=1.571;
end
d=get(handles.f_edit,'String');
f=str2num(d);
R1=isempty(f);
if R1==1
    f=1.047;
end
ee=get(handles.r_edit,'String');
r=str2num(ee);
R1=isempty(r);
if R1==1
    r=1.2;
end
command=sprintf('cocone.exe -c %2f -s %2f -f %2f -r %2f  point \n',C,S,f,r);
h = waitbar(0,'Please wait , the surface is compute......','name','surface computation');
system(command);
waitbar(1/3,h,'30 % finish');
readJV('point.jv');
waitbar(2/3,h,'60 % finish');
 readOFF('point.off');
 axes(handles.axes2);
 N=dlmread('face');
 nn=length(N);
 num1=num2str(nn);
 set(handles.text4,'String',num1);

myview(0,0,0,0)
axis on
grid on
waitbar(3/3,h,'100 % finish');
close(h);
delete('point.dl');
delete('point.jv');
delete('point.wrl');
delete('point.off');
function C_edit_Callback(hObject, eventdata, handles)
% hObject    handle to C_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of C_edit as text
%        str2double(get(hObject,'String')) returns contents of C_edit as a double
%store the contents of input1_editText as a string. if the string
%is not a number then input will be empty
input = str2num(get(hObject,'String'));

%checks to see if input is empty. if so, default input1_editText to zero
if (isempty(input))
     set(hObject,'String','0.393')
end
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function C_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to C_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function S_edit_Callback(hObject, eventdata, handles)
% hObject    handle to S_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of S_edit as text
%        str2double(get(hObject,'String')) returns contents of S_edit as a double
input = str2num(get(hObject,'String'));

%checks to see if input is empty. if so, default input1_editText to zero
if (isempty(input))
     set(hObject,'String','1.571')
end
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function S_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to S_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function f_edit_Callback(hObject, eventdata, handles)
% hObject    handle to f_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of f_edit as text
%        str2double(get(hObject,'String')) returns contents of f_edit as a double
input = str2num(get(hObject,'String'));

%checks to see if input is empty. if so, default input1_editText to zero
if (isempty(input))
     set(hObject,'String','1.047')
end
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function f_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to f_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function r_edit_Callback(hObject, eventdata, handles)
% hObject    handle to r_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of r_edit as text
%        str2double(get(hObject,'String')) returns contents of r_edit as a double
if (isempty(input))
     set(hObject,'String','1.2')
end
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function r_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to r_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function re_edit_Callback(hObject, eventdata, handles)
% hObject    handle to re_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of re_edit as text
%        str2double(get(hObject,'String')) returns contents of re_edit as a double
input = str2num(get(hObject,'String'));

%checks to see if input is empty. if so, default input1_editText to zero
if (isempty(input))
     set(hObject,'String','1.1')
end
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function re_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to re_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in Tbutton.
function Tbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Tbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
aa = get(handles.re_edit,'String');
CK=str2num(aa);
R1=isempty(CK);
if R1==1
    CK=1.1;
end
command=sprintf('tcocone.exe -m -r %2f point surface\n',CK);
h = waitbar(0,'Please wait , the surface is compute......','name','surface computation');
system(command);
waitbar(1/4,h,'25 % finish');

change('surface_surf.off')
waitbar(2/4,h,'50 % finish');
     readOFF('OFF');
     waitbar(3/4,h,'75 % finish');
     axes(handles.axes2);
     N=dlmread('face');
 nn=length(N);
 num1=num2str(nn);
 set(handles.text4,'String',num1);
      myview(0,0,0,0);
      axis on
grid on
  waitbar(4/4,h,'100 % finish');
  close(h);
delete('surface_axis.off');
delete('surface_axis.wrl');
delete('surface_surf.wrl');
delete('surface_surf.off');
delete('OFF');

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
rotate3d on

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
rotate3d off

% --- Executes on button press in expbutton.
function expbutton_Callback(hObject, eventdata, handles)
% hObject    handle to expbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[file1,path1] = uiputfile({'*.ply';'*.off';},'Save file name');
File2=strcat(path1,file1);
ext=findstr(File2,'.ply');
nn1=isempty(ext);
if nn1==0
    WritePly(File2);
end
ext1=findstr(File2,'.off');
nn2=isempty(ext1);
if nn2==0
     WriteOFF(File2);
end

% --- Executes on button press in clearbutton.
function clearbutton_Callback(hObject, eventdata, handles)
% hObject    handle to clearbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
arrayfun(@cla,findall(0,'type','axes'));
 set(handles.text2,'String','');
 set(handles.text4,'String','');
delete('face');
delete('point');


text:
matlabbatch{1}.spm.tools.physio.save_dir = {'physio_out'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'20191023_V02.puls'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'20191023_V02.resp'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'MR.27842.1230'};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.0025;
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = [];
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'first';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 24;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 1.5;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 382;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 5;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = 0.0625;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.scan_timing_log = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 120;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.respiratory.filter.passband = [0.01 2];
matlabbatch{1}.spm.tools.physio.preproc.respiratory.despike = false;
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

author: translationalneuromodeling

repo: tapas


text:
classdef FourierComponent
% This class is obsolet since MTEX 5.9. Use the class @SO3FunHarmonic instead.
% Anyway the class is preserved, so that saved @ODFs can be loaded.

methods (Static = true, Hidden=true)
  function odf = loadobj(s)
    fhat = s.f_hat;
    CS = s.CS;
    SS = s.SS;
    if isempty(fhat), fhat = 0; end
    if isempty(CS), CS = crystalSymmetry; end
    if isempty(SS), SS = specimenSymmetry; end
    odf = SO3FunHarmonic(fhat);
    % We use L2-normalized Wigner-D functions since MTEX Version 5.9.
    odf = L2normalizeFourierCoefficients(odf);
    % The ODFs has been real valued before MTEX verison 5.9.
    odf.isReal = 1;
    % do not symmetrise odf by defining the symmetries.
    odf.CS = CS;
    odf.SS = SS;
  end
end

end


text:
% SOSDEMO5 --- Upper bound for the structured singular value mu
% Section 3.5 of SOSTOOLS User's Manual
%

clear; echo on;
syms x1 x2 x3 x4 x5 x6 x7 x8;
vartable = [x1; x2; x3; x4; x5; x6; x7; x8];

% The matrix under consideration
alpha = 3 + sqrt(3);
beta = sqrt(3) - 1;
a = sqrt(2/alpha);
b = 1/sqrt(alpha);
c = b;
d = -sqrt(beta/alpha);
f = (1 + i)*sqrt(1/(alpha*beta));
U = [a 0; b b; c i*c; d f];
V = [0 a; b -b; c -i*c; -i*f -d];
M = U*V';

% Constructing A(x)'s
gam = 0.8724;

Z = monomials(vartable,1);
for i = 1:4
    H = M(i,:)'*M(i,:) - (gam^2)*sparse(i,i,1,4,4,1);
    H = [real(H) -imag(H); imag(H) real(H)];
    A{i} = (Z.')*H*Z;
end;

% =============================================
% Initialize the sum of squares program
prog = sosprogram(vartable);

% =============================================
% Define SOSP variables

% -- Q(x)'s -- : sums of squares
% Monomial vector: [x1; ... x8]
for i = 1:4
    [prog,Q{i}] = sossosvar(prog,Z);
end;

% -- r's -- : constant sum of squares
Z = monomials(vartable,0);
r = cell(4,4);
for i = 1:4
    for j = (i+1):4
        [prog,r{i,j}] = sossosvar(prog,Z,'wscoeff');
    end;
end;

% =============================================
% Next, define SOSP constraints

% Constraint : -sum(Qi(x)*Ai(x)) - sum(rij*Ai(x)*Aj(x)) + I(x) >= 0
expr = 0;
% Adding term
for i = 1:4
    expr = expr - A{i}*Q{i};
end;
for i = 1:4
    for j = (i+1):4
        expr = expr - A{i}*A{j}*r{i,j};
    end;
end;
% Constant term: I(x) = -(x1^4 + ... + x8^4)
I = -sum(vartable.^4);
expr = expr + I;

prog = sosineq(prog,expr);

% =============================================
% And call solver
prog = sossolve(prog);

% =============================================
% Program is feasible, thus 0.8724 is an upper bound for mu.

echo off

author: yu-jiang

repo: radpbook

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


text:
function reactions = similarity(formula,data,bar)
% Returns the reaction lines with similarities
%
% USAGE:
%
%    reactions = similarity(formula,data,bar)
%
% INPUTS:
%    formula:    reaction formula
%    data:       reaction list
%    bar:        `waitbar` exists then 1, else 0, default 1
%
% OUTPUT:
%    reactions:   numbers of rxn lines that share similarity
%
%
% .. Author: - Stefan G. Thorleifsson 2010
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is

if nargin < 3
    bar = 1;
end

comp = cell(1);
compd = cell(1);
reactions = [];

%Generate metabolite, compartment, coefficient and rev array
[met,coe,rev] = parseRxnFormula(formula);
S = length(met);
for n = 1:S
    k = regexpi(met{n},'\[');
    metab = met{n};
    comp{n} = metab(k+1);
    met{n} = metab(1:k-1);
end

Sdata = length(data);

if bar == 1
    showprogress(0/Sdata,'This may take a while. Checking reaction');
end

for i = 1:Sdata
    if bar == 1
        showprogress(i/Sdata);
    end

    [metd,coed,revd] = parseRxnFormula(data{i});
    %Sd = length(metd);

    if length(metd) == S % If reaction has equal number of metabolites.
        for n = 1:S %Create comparison array
            %k = regexpi(metd{n},'\[');
            %metab = metd{n};
            %compd{n} = metab(k+1);
            %metd{n} = metab(1:k-1);
            %
            compd{n} = regexprep(regexpi(metd{n},'\[.*?\]','match'),'\[*\]*','');
            metd{n} = metd{n}(1:regexp(metd{n},'\[','end')-1);

        end

        match = 0;
        for m = 1:S % Check every metabolite in met
            metab_match = 0;
            for l = 1:S %Check every metabolite in met
                if strcmp(met{m},metd{l})
                    metab_match = 1;
                end
            end
            if  metab_match == 0 % If metabolite dous not match
                break
            else
                match = match+1;
            end
        end

        if match == S
            reactions = [reactions i];
        end

    end

end


text:
function dB0=MagUser(p)
%create dB0 field from data file

global VMmg
global VObj

load(p.MagFile); % load dB0 file
Interp=p.Interp;

% Initialize display grid
xgrid=VMmg.xgrid;
ygrid=VMmg.ygrid;
zgrid=VMmg.zgrid;

Mxdims=size(VObj.Rho);
max_xgrid=((Mxdims(2)-1)/2)*VObj.XDimRes;
max_ygrid=((Mxdims(1)-1)/2)*VObj.YDimRes;
max_zgrid=((Mxdims(3)-1)/2)*VObj.ZDimRes;
min_xgrid=(-(Mxdims(2)-1)/2)*VObj.XDimRes;
min_ygrid=(-(Mxdims(1)-1)/2)*VObj.YDimRes;
min_zgrid=(-(Mxdims(3)-1)/2)*VObj.ZDimRes;

[row,col,layer] = size(dB0);

[Magx,Magy,Magz]=meshgrid(linspace(min_xgrid, max_xgrid, col),...
                          linspace(min_ygrid, max_ygrid, row),...
                          linspace(min_zgrid, max_zgrid, layer));

dB0=ba_interp3(Magx,Magy,Magz,dB0,xgrid,ygrid,zgrid,Interp);
dB0(isinf(dB0)) = 0;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [obs, hidden] = pomdp_sample(initial_prob, transmat, obsmat, act)
% SAMPLE_POMDP Generate a random sequence from a Partially Observed Markov Decision Process.
% [obs, hidden] = sample_pomdp(prior, transmat, obsmat, act)
%
% Inputs:
% prior(i) = Pr(Q(1)=i)
% transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
% obsmat(i,k) = Pr(Y(t)=k | Q(t)=i)
% act(a) = A(t), so act(1) is ignored
%
% Output:
% obs and hidden are vectors of length T=length(act)

len = length(act);
hidden = mdp_sample(initial_prob, transmat, act);
obs = zeros(1, len);
for t=1:len
  obs(t) = sample_discrete(obsmat(hidden(t),:));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pass = test_domainChck()
% Test chebfun3/domainCheck

% Default domain:
ff = @(x,y,z) cos(x+y+z);
gg = @(x,y,z) sin(x.*y.*z);
f = chebfun3(ff);
g = chebfun3(gg);
pass(1) = domainCheck(f, g);

% A different domain:
dom = [-1 2 -2 1 -3 0];
ff = @(x,y,z) cos(x+y+z);
gg = @(x,y,z) sin(x.*y.*z);
f = chebfun3(ff, dom);
g = chebfun3(gg, dom);
pass(2) = domainCheck(f, g);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% poolType = 'mean' or 'max'
% poolLayer = integer layer number after which we do the pooling
function layer = genNetworkTemporalConvDeep(para)
if isfield(para, 'LastActivation4MSE')==0
    para.LastActivation4MSE = 'linear';
end
inputDim = double(para.inputDim);

layer{1}.name = 'Input';        % this is an input layer
layer{end}.inputIdx = 1;    % specifies the index of GCC in Visible_tr
layer{end}.dim = [1 1]*inputDim;             % [input dim; output dim];

nFilter = para.nFilter;
for i=1:length(nFilter)
    if i==1 % only the first layer will use tconv, the rest uses splice and affine transform
        layer{end+1}.name = 'tconv';
        layer{end}.prev = -1;
        layer{end}.W = []; % to be initialized randomly or by pretraining
        layer{end}.b = [];
        layer{end}.dim = [nFilter(i) inputDim*para.filterLen(i)];
    else
        layer{end+1}.name = 'Splice';
        layer{end}.prev = -1;
        layer{end}.context = para.filterLen(i);
        layer{end}.dim = [layer{end}.context 1]*layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.update = 0;

        layer{end+1}.name = 'Affine';
        layer{end}.prev = -1;
        layer{end}.W = []; % to be initialized randomly or by pretraining
        layer{end}.b = [];
        layer{end}.dim = [nFilter(i) layer{length(layer)+layer{end}.prev}.dim(1)];
    end
    layer{end}.update = 1;

    layer{end+1}.name = 'tmaxpool';
    layer{end}.context = para.poolingLen(i);
    layer{end}.stride = para.poolingStride(i);
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*nFilter(i);

    if isfield(para, 'activation')
        layer{end+1}.name = para.activation;
    else
        layer{end+1}.name = 'sigmoid';
    end
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*nFilter(i);
end

layer2 = genNetworkFeedForward_v2(nFilter(end), para.hiddenLayerSizeFF, para.outputDim, para.costFn, para.LastActivation4MSE);

layer = [layer layer2(2:end)];
layer = FinishLayer(layer);

end

author: singaxiong

repo: SignalGraph


text:
function pred = COPAR_pred(Y, D, D_range_ext, opts)
% function pred = COPAR_pred(Y, D, D_range_ext, opts)
% predict label of the input Y
% INPUT:
%   opts.classify_mode = either 'GC' or 'LC'
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    if isfield(opts, 'classify_mode') == 0
        fprintf('You need to specify classification mode in opts.classify_mode: GC/LC\n');
    elseif strcmp(opts.classify_mode, 'GC')
        pred = GC(Y, D, D_range_ext, opts);
    elseif strcmp(opts.classify_mode, 'LC')
        pred = LC(Y, D, D_range_ext, opts);
    else
        fprintf('classify_mode is either GC or LC');
    end
end

function pred = GC(Y, D, D_range_ext, opts)
%     fprintf('GC mode\n');
    C = numel(D_range_ext) - 2;
    optsX = opts;
    optsX.verbose = 0;
    optsX.max_iter = 300;
    if isfield(opts, 'gamma') == 0
        fprintf('specify gamma');
    else
        X = lasso_fista(Y, D, zeros(size(D, 2), size(Y, 2)), ...
                        opts.gamma, optsX);
        DCp1_range = D_range_ext(C+1)+1: D_range_ext(C+2);
        DCp1 = D(:, DCp1_range);
        XCp1 = X(DCp1_range, :);
        Y = Y - DCp1*XCp1;
        E = zeros(C, size(Y,2));
        for i = 1:C
            Xi = get_block_row(X, i, D_range_ext);
            Di = get_block_col(D, i, D_range_ext);
            R = Y - Di*Xi;
            E(i,:) = sum(R.^2, 1);
        end
        [~,pred] = min(E);
    end

end

function pred = LC(Y, D, D_range_ext, opts)
%     fprintf('LC mode\n');
    C = numel(D_range_ext) - 2;
    opts.verbose = false;
    opts.max_iter = 300;
    if isfield(opts, 'gamma') == 0
        fprintf('specify gamma');
    else
        DCp1_range = D_range_ext(C+1)+1: D_range_ext(C+2);
        E = zeros(C, size(Y,2));
        for i = 1:C
            Dc_range = D_range_ext(i)+1: D_range_ext(i+1);
            Dchat = D(:, union(Dc_range, DCp1_range));
%             X = myLasso_fista_2(Y, Dchat, zeros(size(Dchat,2), size(Y, 2)), opts.gamma, pars);
            X = lasso_fista(Y, Dchat, [], opts.gamma, opts);

            R1 = Y - Dchat*X;
            R2 = opts.gamma*abs(X);
            E(i,:) = sum(R1.^2, 1) + sum(R2.^2, 1);
        end
        [~,pred] = min(E);
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
classdef FreeFemFilePrinter < FilePrinter

    properties (Access = private)
        printingDir
    end

    properties (Access = protected)
        linesToPrint
    end

    methods (Access = public, Static)

        function p = create(inputData)
            p = FreeFemFilePrinterFactory.create(inputData);
            p.init(inputData);
            p.createOutPutFolder();
        end

    end

    methods (Access = public)

        function print(obj)
            obj.openFile();
            obj.printingLines();
            obj.closeFile();
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.linesToPrint = d.linesToPrint;
            obj.printingDir  = d.printingDir;
            obj.fileName     = fullfile(obj.printingDir,[d.fileName,'.edp']);
        end

        function createOutPutFolder(obj)
            dir = obj.printingDir;
            if ~exist(dir,'dir')
                mkdir(dir)
                addpath(dir)
            end
        end

    end

    methods (Access = protected, Abstract)
        printingLines(obj)
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [attnInfo] = attnLayerForward(h_t, params, model, attnData, maskInfo)
%
% Attentional Layer: from lstm hidden state to softmax hidden state.
% Input:
%   attnData: require attnData.srcHidVecsOrig and attnData.srcLens
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
  assert(params.attnFunc>0); % we should have enabled attention.

  attnInfo = [];
  if params.attnGlobal % global
    srcHidVecs = attnData.srcHidVecsOrig;
    attnInfo.srcMaskedIds = find(attnData.srcMask(:, 1:params.numSrcHidVecs)'==0); % numSrcHidVecs * curBatchSize
  else % local
    % positions
    if params.attnLocalPred % predictive alignments
      [mu, attnInfo] = regressPositions(model, h_t, attnData.srcLens, params);
      srcPositions = floor(mu);
    else % monotonic alignments
      srcPositions = attnData.tgtPos*ones(1, params.curBatchSize);
      flags = srcPositions>(attnData.srcLens-1);
      srcPositions(flags) = attnData.srcLens(flags)-1;
    end

    % assert
    if params.assert
      assert(isempty(find(srcPositions<1,1)));
      assert(isempty(find(attnData.tgtLens<=1,1)));
      assert(isempty(find(srcPositions(maskInfo.unmaskedIds)>(attnData.srcLens(maskInfo.unmaskedIds)-1),1)));
    end

    % reverse
    if params.isReverse
      srcPositions = params.srcMaxLen - srcPositions;
    end

    % build context vectors
    [srcHidVecs, attnInfo] = buildSrcVecs(attnData.srcHidVecsOrig, srcPositions, maskInfo, attnData.srcLens, params.srcMaxLen, params, attnInfo);

    attnInfo.srcMaskedIds = find(attnInfo.alignMask==0);
  end % end else if attnGlobal

  % compute alignScores: numAttnPositions * curBatchSize
  % TODO: precompute for attnOpt2 and attnOpt3 (we can premultiply srcHidVecs with W_a (attnOpt2) or W_a_src (attnOpt3)
  if params.attnOpt==1 || params.attnOpt==2 % dot product or general dot product
    if params.attnOpt==1 % dot product
      [alignScores] = srcCompareLayerForward(srcHidVecs, h_t, params);
    elseif params.attnOpt==2 % general dot product
      attnInfo.transform_ht = model.W_a * h_t; % TODO: shift the multiplication to srcHidVecs
      [alignScores] = srcCompareLayerForward(srcHidVecs, attnInfo.transform_ht, params);
    end
  elseif params.attnOpt==3 % Bengio's style
    % f(H_src + W_a*h_t): lstmSize * (curBatchSize * numAttnPositions))
    attnInfo.src_ht_hid = reshape(params.nonlinear_f(bsxfun(@plus, srcHidVecs, model.W_a*h_t)), params.lstmSize, []);

    % v_a * src_ht_hid
    alignScores = linearLayerForward(model.v_a, attnInfo.src_ht_hid); % 1 * (curBatchSize * numAttnPositions)
    alignScores = reshape(alignScores, params.curBatchSize, params.numAttnPositions)'; % numAttnPositions * curBatchSize
  end

  if params.attnLocalPred && params.normLocalAttn % new approach for local attention after EMNLP'15
    [attnInfo.distWeights, attnInfo.scaleX] = distLayerForward(mu, attnInfo, params); % numAttnPositions*curBatchSize
    attnInfo.unNormAlignWeights =  alignScores .* attnInfo.distWeights; % weighted by distances
    attnInfo.alignWeights = normLayerForward(attnInfo.unNormAlignWeights, attnInfo.srcMaskedIds);
    attnInfo.alignScores = alignScores;
  else
    % normalize -> alignWeights
    attnInfo.alignWeights = normLayerForward(alignScores, attnInfo.srcMaskedIds);

    % local, regression, multiply with distWeights
    if params.attnLocalPred
      [attnInfo.distWeights, attnInfo.scaleX] = distLayerForward(mu, attnInfo, params); % numAttnPositions*curBatchSize
      attnInfo.preAlignWeights = attnInfo.alignWeights;
      attnInfo.alignWeights = attnInfo.preAlignWeights.* attnInfo.distWeights; % weighted by distances
    end
  end

  % assert
  if params.assert
    assert(computeSum(attnInfo.alignWeights(attnInfo.srcMaskedIds), params.isGPU)==0);
  end

  attnInfo.alignWeights(:, maskInfo.maskedIds) = 0;
  % alignWeights, srcHidVecs -> contextVecs
  [contextVecs] = contextLayerForward(attnInfo.alignWeights, srcHidVecs, maskInfo.unmaskedIds, params);

  % f(W_h*[context_t; h_t])
  attnInfo.input = [contextVecs; h_t];
  attnInfo.h_t = h_t;
  softmax_h = hiddenLayerForward(model.W_h, attnInfo.input, params.nonlinear_f);
  attnInfo.softmax_h = softmax_h; % attentional vectors

  % assert
  if params.assert
    assert(isequal(size(attnInfo.alignWeights), [params.numAttnPositions, params.curBatchSize]));
    assert(isequal(size(h_t), size(contextVecs))); % lstmSize * curBatchSize
  end
end

function [mu, h2sInfo] = regressPositions(model, h_t, srcLens, params)
  % h_t -> scales=sigmoid(v_pos*f(W_pos*h_t)) in [0, 1]
  [h2sInfo.scales, h2sInfo.posForwData] = scaleLayerForward(model.W_pos, model.v_pos, h_t, params);

  % scales -> srcPositions
  mu = h2sInfo.scales.*(srcLens-1) + 1;
end

author: lmthang


text:
% std_stat() - compute statistics for ERP/spectral traces or ERSP/ITC images
%              of a component or channel cluster in a STUDY.
% Usage:
%          >> std_stat( data, 'key', 'val', ...)
% Inputs:
%  data  -  [cell array] mean data for each subject group and/or data
%           condition. For example, to compute mean ERPs statistics from a
%           STUDY for epochs of 800 frames in two conditions from three
%           groups of 12 subjects:
%
%           >> data = { [800x12] [800x12] [800x12];... % 3 groups, cond 1
%                       [800x12] [800x12] [800x12] };  % 3 groups, cond 2
%           >> pcond = std_stat(data, 'condstats', 'on');
%
%           By default, parametric statistics are computed across subjects
%           in the three groups. See below and >> help statcond
%           for more information about the statistical computations.
%
% Statistics options:
%  'groupstats'  - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}
%  'condstats'   - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}
%  'statistics'  - ['param'|'perm'] Type of statistics to use: 'param' for
%                  parametric; 'perm' for permutations {default: 'param'}
%  'naccu'       - [integer] Number of surrogate averges fo accumulate when
%                  computing permutation-based statistics. For example, to
%                  test p<0.01 use naccu>=200; for p<0.001, use naccu>=2000.
%                  If a non-NaN 'threshold' is set (see below) and 'naccu'
%                  is too low, it will be automatically increased. This
%                  keyword available only from the command line {default:500}
%  'threshold'   - [NaN|p-value] threshold for computing p-value. In this
%                  function, it is only used to compute naccu above. NaN
%                  means that no threshold has been set.
%  'mcorrect'    - ['none'|'fdr'] apply correcting for multiple comparisons.
%
% Outputs:
%  pcond        - [cell] condition pvalues or mask (0 or 1). One element per group.
%  pgroup       - [cell] group pvalues or mask (0 or 1). One element per condition.
%  pinter       - [cell] three elements, condition pvalues (group pooled),
%                 group pvalues (condition pooled) and interaction pvalues.
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: statcond()

%  'statmode'    - ['subjects'|'trials'] standard statistics are
%                  'subjects' where the statistics is performed accross
%                  the mean ERSP (or ITC) of single subjects. For 'trials'
%                  statistics, the single-trial data epochs of all subjects
%                  are pooled together. This requires that they were
%                  saved on disk using option 'savetrials', 'on' at the time
%                  of computation. Note that these single-trial data
%                  may use several GB of disk space and that computation
%                  of 'trials' statistics requires a lot of RAM.

function [pcond, pgroup, pinter] = std_stat(data, varargin)

pgroup = {};
pcond  = {};
pinter = {};
if nargin < 1
    help std_stat;
    return;
end;

% decode inputs
% -------------
opt = {};
if nargin > 1
    if isstruct(varargin{1})
        g = struct(varargin{2:end});
        if isfield(g, 'condstats'),  varargin{1} = rmfield(varargin{1}, 'condstats'); end;
        if isfield(g, 'groupstats'), varargin{1} = rmfield(varargin{1}, 'groupstats'); end;
        tmpparams      = fieldnames(varargin{1}); tmpparams = tmpparams';
        tmpparams(2,:) = struct2cell(varargin{1});
        varargin       = { tmpparams{:} varargin{2:end} };
    end;
end;
opt = finputcheck( varargin, { 'threshold'   'real'    []               NaN;
                               'mcorrect'    'string'  { 'none' 'fdr' } 'none';
                               'naccu'       'integer' []               [];
                               'groupstats'  'string'  { 'on' 'off' }   'off';
                               'paired'      'cell'    { 'on' 'off' }   { 'on' 'on' };
                               'condstats'   'string'  { 'on' 'off' }   'off';
                               'statistics'  'string'  { 'param' 'perm' 'bootstrap' }       'param' }, ...
                               'std_stat', 'ignore');

if isstr(opt), error(opt); end;
if ~isnan(opt.threshold(1)) && isempty(opt.naccu), opt.naccu = 1/opt.threshold(end)*2; end;
if any(any(cellfun('size', data, 2)==1)), opt.groupstats = 'off'; opt.condstats = 'off'; end;
if strcmpi(opt.mcorrect, 'fdr'), opt.naccu = opt.naccu*20; end;
if isempty(opt.naccu), opt.naccu = 2000; end;
if strcmpi(opt.statistics, 'param') && ~isreal(data{1})
    fprintf('*** Cannot use parametric statistics for single-trial ITC significance ***\n');
    return;
end;
nc = size(data,1);
ng = size(data,2);

% compute significance mask
% --------------------------
if strcmpi(opt.condstats, 'on') && nc > 1
    for g = 1:ng
        [F df pval] = statcond(data(:,g), 'mode', opt.statistics, 'naccu', opt.naccu, 'paired', opt.paired{1});
        pcond{g} = squeeze(pval);
    end;
else
    pcond = {};
end;
if strcmpi(opt.groupstats, 'on') && ng > 1
    for c = 1:nc
        [F df pval] = statcond(data(c,:), 'mode', opt.statistics, 'naccu', opt.naccu, 'paired', opt.paired{2});
        pgroup{c} = squeeze(pval);
    end;
else
    pgroup = {};
end;
if ( strcmpi(opt.groupstats, 'on') && strcmpi(opt.condstats, 'on') ) & ng > 1 & nc > 1
    opt.paired = sort(opt.paired); % put 'off' first if present
    [F df pval] = statcond(data, 'mode', opt.statistics, 'naccu', opt.naccu, 'paired', opt.paired{1});
    for index = 1:length(pval)
        pinter{index} = squeeze(pval{index});
    end;
else
    pinter = {};
end;

if ~isempty(opt.groupstats) || ~isempty(opt.condstats)
    if strcmpi(opt.mcorrect, 'fdr'),
        disp('Applying FDR correction for multiple comparisons');
        for ind = 1:length(pcond),  pcond{ind} = fdr( pcond{ind} ); end;
        for ind = 1:length(pgroup), pgroup{ind} = fdr( pgroup{ind} ); end;
        if ~isempty(pinter),
            pinter{1} = fdr(pinter{1});
            pinter{2} = fdr(pinter{2});
            pinter{3} = fdr(pinter{3});
        end;
    end;
    if ~isnan(opt.threshold)
        for ind = 1:length(pcond),  pcond{ind}  = applythreshold(pcond{ind},  opt.threshold); end;
        for ind = 1:length(pgroup), pgroup{ind} = applythreshold(pgroup{ind}, opt.threshold); end;
        for ind = 1:length(pinter), pinter{ind} = applythreshold(pinter{ind}, opt.threshold); end;
    end;
end;

function newdata = applythreshold(data, threshold)
    threshold = sort(threshold);
    newdata = zeros(size(data));
    for index = 1:length(threshold)
        inds = data < threshold(index);
        data(inds)    = 1;
        newdata(inds) = length(threshold)-index+1;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function K = lfmjXrbfKernCompute(lfmKern, rbfKern, t1, t2)

% LFMJXRBFKERNCOMPUTE Compute cross kernel between the LFMJ and RBF kernels.
% FORMAT
% DESC computes cross kernel terms between LFMJ and RBF kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFMJ
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between LFMJ and RBF kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFMJ
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

if nargin < 4
  t2 = t1;
end
if size(t1, 2) > 1 || size(t2, 2) > 1
  error('Input can only have one column');
end

if lfmKern.inverseWidth ~= rbfKern.inverseWidth
  error('Kernels cannot be cross combined if they have different inverse widths.')
end

% Get length scale out.
sigma2 = 2/lfmKern.inverseWidth;
sigma = sqrt(sigma2);

% Parameters of the kernel
alpha = lfmKern.damper./(2*lfmKern.mass);
omega = sqrt(lfmKern.spring./lfmKern.mass - alpha.*alpha);

gamma1 = alpha + j*omega;
gamma2 = alpha - j*omega;

sK = lfmjpComputeUpsilonMatrix(gamma2,sigma2,t1, t2, 0) - ...
    lfmjpComputeUpsilonMatrix(gamma1,sigma2,t1, t2, 0);

if lfmKern.isNormalised
   K0 = lfmKern.sensitivity/(j*4*sqrt(2)*lfmKern.mass*omega);
else
   K0 = sqrt(pi)*sigma*lfmKern.sensitivity/(j*4*lfmKern.mass*omega);
end

K = K0*sK;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat