text:
% std_topomovie - make movie in the frequency domain
%
% Usage:
%  >> [STUDY] = std_movie(STUDY, ALLEEG, key1, val1, key2, val2, ...);
%
% Inputs:
%   STUDY      - STUDY structure comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - vector of EEG dataset structures for the dataset(s) in the STUDY,
%                typically created using LOAD_ALLEEG.
%   'channels' - [numeric vector]  specific channel group to plot. By
%                default, the grand mean channel spectrum is plotted (using the
%                same format as for the cluster component means described above)
%   'moviemode' - ['erp'|'spec'|'ersptime'] movie mode. Currently only
%                 'spec' is implemented.
%   'erspfreq'  - [min max] frequency range when making movie of ERSP. The
%                 ERSP values are averaged over the selected frequency
%                 range. Not implemented
%   'limitbeg'  - [min max] limits at the beginning of the movie
%   'limitend'  - [min max] limits at the end of the movie
%   'freqslim'  - [freqs] array of landmark frequencies to set color
%                 limits.
%   'movieparams' - [low inc high] lower limit, higher limit and increment. If
%                increment is omitted, movie is generate at every possible
%                increment (max freq resolution or max time resolution)
%
% Authors: Arnaud Delorme, CERCO, August, 2006
%
% See also: STD_SPECPLOT, STD_ERPPPLOT, STD_ERSPPLOT

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [STUDY, M] = std_movie(STUDY, ALLEEG, varargin);

if nargin < 2
    help std_specplot;
    return;
end

[opt addparams ] = finputcheck( varargin, ...
                             { 'erspfreq'    'real'    [] [];
                               'movieparams' 'real'    [] [];
                               'channels'    'cell'    [] {};
                               'freqslim'    'real'    [] [];
                               'limitbeg'    'real'    [] [];
                               'limitend'    'real'    [] [];
                               'moviemode'   'string'  { 'ersptime','erp','spec' } 'spec' }, 'std_movie', 'ignore');

if ischar(opt), error(opt); end
tmpchanlocs =  ALLEEG(1).chanlocs;
if isempty(opt.channels), opt.channels = { tmpchanlocs.labels }; end
if ~strcmpi(opt.moviemode, 'spec'), error('Only spec has been implemented so far'); end

% read data once
% --------------
STUDY = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', [10 11]); close;

% find first data channel with info
% ---------------------------------
for cind = 1:length(STUDY.changrp)
    if ~isempty(STUDY.changrp(cind).specdata), break; end
end

% generate movie
% --------------
if isempty(opt.movieparams),
    opt.movieparams(1) = STUDY.changrp(cind).specfreqs(1);
    opt.movieparams(2) = STUDY.changrp(cind).specfreqs(end);
end
if length(opt.movieparams) == 2
    opt.movieparams(3) = opt.movieparams(2);
    opt.movieparams(2) = STUDY.changrp(cind).specfreqs(2)-STUDY.changrp(cind).specfreqs(1);
end
if length(opt.movieparams) == 3
    opt.movieparams = [opt.movieparams(1):opt.movieparams(2):opt.movieparams(3)];
end

% find limits
% -----------
if isempty(opt.limitbeg)
    [STUDY specdata] = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(1)); close;
    opt.limitbeg = [ min(specdata{1}) max(specdata{1}) ];
    [STUDY specdata] = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(end)); close;
    opt.limitend = [ min(specdata{1}) max(specdata{1}) ];
end
lowlims  = linspace(opt.limitbeg(1), opt.limitend(1), length(opt.movieparams));
highlims = linspace(opt.limitbeg(2), opt.limitend(2), length(opt.movieparams));

% limits at specific frequencies
% ------------------------------
if ~isempty(opt.freqslim)
    if opt.freqslim(1)   ~= opt.movieparams(1)  , opt.freqslim = [ opt.movieparams(1) opt.freqslim ]; end
    if opt.freqslim(end) ~= opt.movieparams(end), opt.freqslim = [ opt.freqslim opt.movieparams(end) ]; end

    for ind = 1:length(opt.freqslim)
        [tmp indf(ind)] = min(abs(opt.freqslim(ind) - opt.movieparams));
        [STUDY specdata] = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(indf(ind))); close;
        minimum(ind) = min(specdata{1});
        maximum(ind) = max(specdata{1});
    end
    indf(1)  = 0;
    lowlims  = [ ];
    highlims = [ ];
    for ind = 2:length(opt.freqslim)
        lowlims  = [lowlims linspace(minimum(ind-1), minimum(ind), indf(ind)-indf(ind-1)) ];
        highlims = [highlims linspace(maximum(ind-1), maximum(ind), indf(ind)-indf(ind-1)) ];
    end
end

% make movie
% ----------
for ind = 1:length(opt.movieparams)
    STUDY = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(ind), 'caxis', [lowlims(ind) highlims(ind)]);
    pos = get(gcf, 'position');
    set(gcf, 'position', [pos(1) pos(2) pos(3)*2 pos(4)*2]);
    M(ind) = getframe(gcf);
    close;
end
figure; axis off; movie(M);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function hidden_probability = visible_state_to_hidden_probabilities(rbm_w, visible_state)
% <rbm_w> is a matrix of size <number of hidden units> by <number of visible units>
% <visible_state> is a binary matrix of size <number of visible units> by <number of configurations that we're handling in parallel>.
% The returned value is a matrix of size <number of hidden units> by <number of configurations that we're handling in parallel>.
% This takes in the (binary) states of the visible units, and returns the activation probabilities of the hidden units conditional on those states.
    m_ = rbm_w * visible_state;
    hidden_probability = 1 ./ (1 + e .^ (-m_));
end

repo: neural-nets


text:
function [x, info] = lewisCenter(f, x, opts)
% x = LewisCenter(f, opts, x)
% compute the analytic center for the domain {Ax=b} intersect the domain of f
%
% Input:
%    f - a ConvexProgram
%    x - a feasible initial point
%    opts - a structure for options with the following properties (optional)
%       MaxIter - maximum number of iterations
%       Output - maximum number of iterations
%       CentralityTol, FeasibilityTol - stop the following are satisfied
%           ||(A' * lambda - grad f(x)) / sqrt(hess)||_inf < CentralityTol
%           ||A x - b||_inf < FeasibilityTol
%       p - the parameter for lp Lewis weight
%       JLDim - the number of dimensions used in estimating leverage score
%
% Output:
%  x - It outputs the analytic center of f
%  info - a structure containing centrality, feasibility and iter.

defaultOpts = struct('MaxIter', 100, 'Output', @disp, 'CentralityTol', 0.1, 'FeasibilityTol', 1e-12, 'JLDim', 10, 'p', 4);
if nargin >= 3
   opts = setField(defaultOpts, opts);
else
   opts = defaultOpts;
end
x = ddouble(x);

%% prepare the printout
output = TableDisplay('Iter', '5i', 'PredStep', '13.2e', 'CorrStep', '13.2e', 'Centrality', '13.2e', 'Feasibility', '13.2e');
output.output = opts.Output;
output.header();

%% initial parameters
assert(nargin >= 2 && all(f.distance(x) > 0), 'a feasible inital point is required.');

A = f.A; At = A'; b = f.b;
y = 0*b;
lastProgress = 0; % record the last iteration making progress
bestCentrality = 1e32;
d = size(f.A, 1); n = size(f.A, 2);
w = (n-d)/n * ones(n, 1);

%% find the central path
[rs, rx, h] = updateResidual(x, y, w);

for iter = 1:opts.MaxIter
   % Compute the cholesky decomposition
   cholErr = f.solver.factorize(1./h);
   w_new = max(double(1-f.solver.leverageScore(opts.JLDim)), 0);
   w = (w + w_new)/2;

   % Compute the direction
   if (cholErr < f.solver.cholTol)
      v = f.solver.solve([A*(rs./h) rx]); %rs = A^T y - g; %rx = Ax - b;
      Atv = At * v;
      y = y - v(:,1);

      % y = y - (R\(R'\(A*(rs./hess))))
      % dx = (rs + At * (R\(R'\(rx - A*(rs./hess)))))./hess;

      dx = (rs - Atv(:,1))./h - Atv(:,2)./h;
      t = stepSize(x, dx, 0.95); % 0.95*(distance from x to the closest boundart in direction dx)
      x = x + t * dx; % move in direction dx with distance of t (measured in dx)

      [rs, rx, h] = updateResidual(x, y, w);

      centrality = max(abs(rs./sqrt(h))); %||(A^T y - grad \phi)/sqrt(h)||_Inf
      feasibility = max(abs(rx)); %||Ax-b||_Inf
      if isempty(feasibility), feasibility = 0; end % Fix the case rx is []

      % Output the error
      o = struct('Iter', iter, 'PredStep', t, 'CorrStep', t, 'Centrality', centrality, 'Feasibility', feasibility);
      output.row(o);

      if (centrality < 0.9 * bestCentrality)
         bestCentrality = centrality;
         lastProgress = iter;
      end

      % Check stop criteria
      if centrality < opts.CentralityTol && feasibility < opts.FeasibilityTol
         break;
      end
   end

   if (iter > lastProgress + 10)
      break;
   end

   if (cholErr >= f.solver.cholTol)
      opts.Output('Failed due to numerical issue.');
      f.feasible = false;
      x = []; info = [];
      return;
   end
end

info = struct('centrality', centrality, 'feasibility', feasibility, 'iter', iter, 'hess', h);

function t = stepSize(x, dx, factor)
   t = min(double(factor*min(f.barrier.distance(x, dx))),1); % min(Nan,1) = 1 for double
end

function [rs, rx, h, wp] = updateResidual(x, y, w)
   wp = w.^(1-2/opts.p);

   [gb, hb] = f.barrier.derivatives(x); % obtain gradient and hessian of log-barrier ftn \phi(x) = -\sum log(x-l) -\sum log(u-x)
   % hp = diagonals of hessian matrix of log-barrier

   gb = gb.*wp;
   hb = hb.*wp;

   gc = zeros(size(gb), class(gb));
   hc = zeros(size(hb), class(hb));

   % Below if-end is skipped while running f.normalize() in sample.m (due to c=[])
   if ~isempty(f.c)
      gc = gc + f.c;
   end

   % Below if-end is skipped while running f.normalize() in sample.m (due to df=[])
   if ~isempty(f.df)
      gc_ = f.df(f.export(x));
      gc = gc + f.scale .* gc_(f.idx);

      hc_ = f.ddf(f.export(x));
      hc = hc + f.scale2 .* hc_(f.idx);
   end

   % Thus, g=gb and h=hb in f.normalize()
   g = gb + gc; h = hb + hc;

   rs = At * y - g;
   rx = A * x - b;
end

end

author: opencobra

repo: cobratoolbox


text:
function config=createConfig()

config.path.input='/home/gneelima/work/data/datasets/PASCAL2007/testImages/';
config.path.output='/home/gneelima/work/data/output/objectProposals/randomizedPrims/PASCAL2007/proposals/';

config.opts.imageExt='.jpg';

config.params.rSeedForRun=-1;
config.params.approxFinalNBoxes=10000;
config.params.q=10;

config.params.segmentations{1}.colorspace='LAB';

config.params.segmentations{1}.superpixels.sigma=0.8;
config.params.segmentations{1}.superpixels.c=100;
config.params.segmentations{1}.superpixels.min_size=100;

config.params.segmentations{1}.simWeights.wBias=3.0017;
config.params.segmentations{1}.simWeights.wCommonBorder=-1.0029;
config.params.segmentations{1}.simWeights.wLABColorHist=-2.6864;
config.params.segmentations{1}.simWeights.wSizePer=-2.3655;

config.params.segmentations{1}.alpha=ones(1,65536);
config.params.segmentations{1}.verbose=1;

author: Cloud-CV

repo: object-proposals


text:
function prob_map_obj = atlas_get_probability_maps(obj, varargin)
% Given an atlas object, return probability maps or indicator maps based on parcels
%
% :Usage:
% ::
% prob_map_obj = atlas_get_probability_maps(obj)
%
% :Inputs:
%
%   **obj:**
%        an atlas-class object
%
% :Outputs:
%
%   **prob_map_obj:**
%        an fmri_data object containing the probability maps or indicators
%        as data
%
% Copyright 2018 Tor Wager

% Get maps if they exist
% -------------------------------------------------------------------------

if isempty(obj.probability_maps)

    obj.probability_maps = single(condf2indic(obj.dat, 'integers'));

    obj.probability_maps(isnan(obj.probability_maps)) = 0;

end

% Get maps and construct an object
% -------------------------------------------------------------------------

prob_map_obj = image_vector('dat', single(full(obj.probability_maps)), ...
    'volInfo', obj.volInfo, ...
    'removed_voxels', obj.removed_voxels, 'removed_images', obj.removed_images, ...
    'image_names', obj.image_names, 'noverbose');

prob_map_obj = fmri_data(prob_map_obj);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Map the numerical elements in the vector "v" onto the variables "s" which can
% be of any type. The number of numerical elements must match; on exit "v"
% should be empty. Non-numerical entries are just copied. See also unwrap2vec.m.

function [s v] = rewrap(s, v)

if isnumeric(s)
  if numel(v) < numel(s)
    error('The vector for conversion contains too few elements')
  end
  s = reshape(v(1:numel(s)), size(s));            % numeric values are reshaped
  v = v(numel(s)+1:end);                        % remaining arguments passed on
elseif isstruct(s)
  [s p] = orderfields(s); p(p) = 1:numel(p);      % alphabetize, store ordering
  [t v] = rewrap(struct2cell(s), v);                 % convert to cell, recurse
  s = orderfields(cell2struct(t,fieldnames(s),1),p);  % conv to struct, reorder
elseif iscell(s)
  for i = 1:numel(s)             % cell array elements are handled sequentially
    [s{i} v] = rewrap(s{i}, v);
  end
end                                             % other types are not processed

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [out_img, criterion] = TVdenoising(img, method, num_steps, lambda, clear_img, alpha, showfigs)

[H, W] = size(img);
N = H * W;

%% method aspects
switch method
    case 'ROFalg1'
        alg = 1;
        Lone = 0;
        huber = 0;
    case 'ROFalg2'
        alg = 2;
        Lone = 0;
        huber = 0;
    case 'TVL1ROFalg1'
        alg = 1;
        Lone = 1;
        huber = 0;
    case 'HuberROFalg3'
        alg = 3;
        Lone = 0;
        huber = 1;
    case 'HuberL1ROFalg1'
        alg = 1;
        Lone = 1;
        huber = 1;
    otherwise
        disp(['Unknown method: ' method]);
        return;
end

%% parameters
L = sqrt(8);

switch alg
    case 0      % Arrow-Hurwics version of Alg 2 (AHMOD)
        tau = 1/L;
        sigma = 1/L;
        gamma = 0.35 * lambda;
        theta = 0;
    case 1
        tau = 0.01;
        sigma = 1/(tau * L * L);
        theta = 1;
    case 2
        tau = 1/L;
        sigma = 1/L;
        gamma = 0.35*lambda;
    case 3
        gamma = lambda;
        delta = alpha;
        mu = 2 * sqrt(gamma * delta) / L;
        theta = 1/(1+mu);
        tau = mu / (2 * gamma);
        sigma = mu / (2 * delta);
end

%% initial solution
% primal task variables
u = img(:);
ubar = u;

% dual task variable
p = zeros(N * 2, 1);

%% precomputed
nabla = make_derivatives_mine(H, W);
divop = nabla';
% divop = make_divop(H, W);
denom = 1 + tau * lambda;

%% initial criterion value
if (Lone)
    lambda_denom = 1;
else
    lambda_denom = 2;
end

criterion = zeros(1, num_steps+1);
criterion(1) = Fval(u, img, alpha, huber) + lambda / lambda_denom * Gval(u, img, Lone);

%% plot the initial state
if showfigs
    fh1 = sfigure;

    imshow([img reshape(u, H, W)]);

    fh2 = sfigure;
    plot(0, criterion(1), 'b-');
    xlabel('step');
    ylabel('J(u)');
end

%% main loop
for step = 1:num_steps
    disp(['step: ' num2str(step)]);

    % ------ update p^n+1 ------
    p = p + sigma * nabla * ubar;

    if huber
        p = p / (1 + sigma * alpha);
    end

    % projection of p onto L2 ball
    p_len = sqrt(p(1:N).^2 + p(N+1:end).^2);
    p_len = max(1, p_len);
    p = p ./ repmat(p_len, 2, 1);

    % ----- update u^n+1 ------
    divp = divop * p;
    u_tilde = u - tau * divp;

    if Lone
        dif = u_tilde - img(:);
        idx1 = dif > tau * lambda;
        idx2 = dif < -tau * lambda;
        idx3 = abs(dif) <= tau * lambda;
        u_new = zeros(size(u));
        u_new(idx1) = u_tilde(idx1) - tau * lambda;
        u_new(idx2) = u_tilde(idx2) + tau * lambda;
        u_new(idx3) = img(idx3);
    else
        if alg == 2 || alg == 0
            denom = 1 + tau * lambda;
        end
        u_new = (u_tilde + tau * lambda * img(:)) / denom;
    end

    % ----- update tau and sigma -----
    if alg == 2
        theta = 1/(sqrt(1+2*gamma*tau));
        tau = tau * theta;
        sigma = sigma / theta;
    end
    if alg == 0
        theta_tmp = 1/(sqrt(1+2*gamma*tau));
        tau = tau * theta_tmp;
        sigma = sigma / theta_tmp;
    end

    % update ubar^n+1
    ubar = u_new + theta * (u_new - u);

    u = u_new;

    % compute the criterion function value
    criterion(step+1) = Fval(u, clear_img, alpha, huber) + lambda / lambda_denom * Gval(u, img, Lone);

    % plot current result
    if showfigs
        sfigure(fh1);
        imshow([img reshape(u, H, W)]);
        drawnow;

        sfigure(fh2);
        plot(0:step, criterion(1:step+1), 'b-');
        xlabel('step');
        ylabel('J(u)');
        drawnow;
    end
end

out_img = u;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [P, R] = mdp_example_rand (S, A, is_sparse, mask)

% mdp_example_rand    Generate a random Markov Decision Process
% Arguments -------------------------------------------------------------
%   S = number of states (> 0)
%   A = number of actions (> 0)
%   is_sparse = false to have matrices in plain format, true to have sparse matrices
%               optional (default false).
%   mask(SxS) = matrix with 0 and 1 (0 indicates a place for a zero probability),
%               optional (default, ones(S,S) )
% Evaluation -------------------------------------------------------------
%   P(SxSxA) = transition probability matrix
%   R(SxSxA) = reward matrix

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

% arguments checking
if S < 1 || A < 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: The number of states S ')
    disp('and the number of actions A must be upper than 1.')
    disp('--------------------------------------------------------')
elseif nargin >= 4 && ( size(mask,1) ~= S || size(mask,2) ~= S )
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: mask must be a SxS matrix')
    disp('--------------------------------------------------------')
else

    % initialization of optional arguments
    if nargin < 3; is_sparse = false; end;
    if nargin < 4; mask = ones(S,S); end;

    if is_sparse
        % definition of transition matrix : square stochastic matrix
        P = {};
        for a=1:A
            PP = sparse(mask .* rand(S));
            for s=1:S
                PP(s,:) = PP(s,:) / sum( PP(s,:) );
            end;
            P{a} = PP;
        end;

        % definition of reward matrix (values between -1 and +1)
        R = {};
        for a=1:A
	    R{a} = sparse(mask .* ( 2*rand(S) - ones(S,S) ));
        end;
    else
        % definition of transition matrix : square stochastic matrix
        for a=1:A
            P(:,:,a) = mask .* rand(S);
            for s=1:S
                P(s,:,a) = P(s,:,a) / sum( P(s,:,a) );
            end;
        end;

        % definition of reward matrix (values between -1 and +1)
        for a=1:A
	    R(:,:,a) = mask .* ( 2*rand(S) - ones(S,S) );
        end;
    end

end;

author: Sable

repo: mcbench-benchmarks


text:
function str = ConvertMFSK2RS(Str1)

for i = 1:length(Str1)
    if (Str1(i) == 0)
        str(i) = -Inf;
    else
        str(i) = Str1(i)-1;
    end
end


text:
function waic = gp_waic(gp, x, y, varargin)
%GP_WAIC The widely applicable information criterion (WAIC) for GP model
%
%  Description
%    WAIC = GP_WAIC(GP, X, Y) evaluates WAIC defined by
%    Watanabe(2010) given a Gaussian process model GP, training
%    inputs X and training outputs Y. Instead of Bayes loss we
%    compute the Bayes utility which is just the negative of loss
%    used by Watanabe.
%
%    WAIC is evaluated as follows when using the variance form
%
%      WAIC(n) = BUt(n) - V/n
%
%    where BUt(n) is Bayesian training utility,  V is functional variance
%    and n is the number of training inputs.
%
%      BUt = mean(log(p(yt | xt, x, y)))
%      V = sum(E[log(p(y|th))^2] - E[log(p(y|th))]^2)
%
%    When using the Gibbs training loss, WAIC is evaluated as follows
%
%          WAIC(n) = BUt(n) - 2*(BUt(n) - GUt(n))
%
%    where BUt(n) is as above and GUt is Gibbs training utility
%
%          GUt(n) = E_th[mean(log(p(y|th)))].
%
%    GP can be a Gaussian process structure, a record structure
%    from GP_MC or an array of GPs from GP_IA.
%
%   OPTIONS is optional parameter-value pair
%      method - Method to evaluate waic, 'V' = Variance method, 'G' = Gibbs
%               training utility method (default = 'V')
%      form -   Return form, 'mean' returns the mean value and 'all'
%               returns the values for all data points (default = 'mean')
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%
%   See also
%     GP_DIC, DEMO_MODELASSESMENT1, DEMO_MODELASSESMENT2
%
%   References
%
%     Watanabe(2010). Equations of states in singular statistical
%     estimation. Neural Networks 23 (2010), 20-34
%
%     Watanabe(2010). Asymptotic Equivalance of Bayes Cross Validation and
%     Widely applicable Information Criterion in Singular Learning Theory.
%     Journal of Machine Learning Research 11 (2010), 3571-3594.
%
%

% Copyright (c) 2011-2013 Ville Tolvanen

  ip=inputParser;
  ip.FunctionName = 'GP_WAIC';
  ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('method', 'V', @(x) ismember(x,{'V' 'G'}))
  ip.addParamValue('form', 'mean', @(x) ismember(x,{'mean','all'}))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.parse(gp, x, y, varargin{:});
  method=ip.Results.method;
  form=ip.Results.form;
  % pass these forward
  options=struct();
  z = ip.Results.z;
  if ~isempty(ip.Results.z)
    options.zt=ip.Results.z;
    options.z=ip.Results.z;
  end

  [tn, nin] = size(x);

  % ====================================================
  if isstruct(gp)     % Single GP or MCMC solution
    switch gp.type
      case {'FULL' 'VAR' 'DTC' 'SOR'}
        tstind = [];
      case {'FIC' 'CS+FIC'}
        tstind = 1:tn;
      case 'PIC'
        tstind = gp.tr_index;
    end

    if isfield(gp, 'etr')
      % MCMC solution
      [Ef, Varf, BUt] = gpmc_preds(gp,x,y, x, 'yt', y, 'tstind', tstind, options);
      BUt=log(mean(exp(BUt),2));
      GUt = zeros(tn,1);
      Elog = zeros(tn,1);
      Elog2 = zeros(tn,1);

      nsamples = length(gp.edata);
      if strcmp(gp.type, 'PIC')
        tr_index = gp.tr_index;
        gp = rmfield(gp, 'tr_index');
      else
        tr_index = [];
      end

      %Ef = zeros(tn, nsamples);
      %Varf = zeros(tn, nsamples);
      sigma2 = zeros(tn, nsamples);
      for j = 1:nsamples
        Gp = take_nth(gp,j);
        if  strcmp(gp.type, 'FIC') | strcmp(gp.type, 'PIC')  || strcmp(gp.type, 'CS+FIC') || strcmp(gp.type, 'VAR') || strcmp(gp.type, 'DTC') || strcmp(gp.type, 'SOR')
          Gp.X_u = reshape(Gp.X_u,length(Gp.X_u)/nin,nin);
        end
        Gp.tr_index = tr_index;

        gp_array{j} = Gp;
        %[Ef(:,j), Varf(:,j)] = gp_pred(Gp, x, y, x, 'yt', y, 'tstind', tstind, options);
        if isfield(gp.lik.fh,'trcov')
          sigma2(:,j) = repmat(Gp.lik.sigma2,1,tn);
        end
      end

      if isequal(method,'V')
        % Evaluate WAIC using the Variance method

        if isfield(gp.lik.fh,'trcov')
          % Gaussian likelihood
          for i=1:tn
%             fmin = mean(Ef(i,:) - 9*sqrt(Varf(i,:)));
%             fmax = mean(Ef(i,:) + 9*sqrt(Varf(i,:)));
%             Elog(i) = quadgk(@(f) mean(multi_npdf(f,Ef(i,:),(Varf(i,:))) ...
%                                        .*bsxfun(@minus,-bsxfun(@rdivide,(repmat((y(i)-f),nsamples,1)).^2,(2.*sigma2(i,:))'), 0.5*log(2*pi*sigma2(i,:))').^2), fmin, fmax);
%             Elog2(i) = quadgk(@(f) mean(multi_npdf(f,Ef(i,:),(Varf(i,:))) ...
%                                         .*bsxfun(@minus,-bsxfun(@rdivide,(repmat((y(i)-f),nsamples,1)).^2,(2.*sigma2(i,:))'), 0.5*log(2*pi*sigma2(i,:))')), fmin, fmax);
%
            m = Ef(i,:);
            s2 = Varf(i,:);
            m0 = 1; m1 = m; m2 = m.^2 + s2; m3 = m.*(m.^2+3*s2);
            m4 = m.^4+6.*m.^2.*s2+3*s2.^2;
            Elog2(i) = mean((-0.5.*log(2.*pi.*sigma2(i,:)) - y(i).^2./(2.*sigma2(i,:))).*m0 - 1./(2.*sigma2(i,:)) .* m2 + y(i)./sigma2(i,:) .* m1);
            Elog(i) = mean((1/4 .* m4 - y(i) .* m3 + (3.*y(i).^2./2+0.5.*log(2.*pi.*sigma2(i,:)).*sigma2(i,:)) .* m2 ...
              - (y(i).^3 + y(i).*log(2.*pi.*sigma2(i,:)).*sigma2(i,:)) .* m1 + (y(i).^4./4 + 0.5.*y(i).^2.*log(2.*pi.*sigma2(i,:)).*sigma2(i,:) ...
              + 0.25.*log(2.*pi.*sigma2(i,:)).^2.*sigma2(i,:).^2) .* m0) ./ sigma2(i,:).^2);
          end
          Elog2 = Elog2.^2;
          Vn = (Elog-Elog2);
          if strcmp(form, 'mean')
            Vn = mean(Vn);
            BUt = mean(BUt);
          end
          waic = BUt - Vn;
        else
          % non-Gaussian likelihood
          for i=1:tn
            if ~isempty(z)
              z1 = z(i);
            else
              z1 = [];
            end
            if ~isequal(gp.lik.type, 'Coxph')
              fmin = mean(Ef(i,:) - 9*sqrt(Varf(i,:)));
              fmax = mean(Ef(i,:) + 9*sqrt(Varf(i,:)));
              Elog(i) = quadgk(@(f) mean(multi_npdf(f,Ef(i,:),(Varf(i,:))) ...
                .*llvec(gp_array, y(i), f, z1).^2), fmin, fmax);
              Elog2(i) = quadgk(@(f) mean(multi_npdf(f,Ef(i,:),(Varf(i,:))) ...
                .*llvec(gp_array, y(i), f, z1)), fmin, fmax);
            else
              ntime = size(gp.lik.xtime,1);
              for i2=1:nsamples
                % Use MC to integrate over latents
                ns = 10000;
                Sigma_tmp = diag(Varf([1:ntime ntime+i],i2));
                f = mvnrnd(Ef([1:ntime ntime+i],i2), Sigma_tmp, ns);
                tmp2(i2) =  1/ns * sum(llvec(gp_array{i2}, y(i,:), f', z1));
                tmp(i2) = 1/ns * sum((llvec(gp_array{i2}, y(i,:), f', z1)).^2);
              end
              Elog2(i)=mean(tmp2);
              Elog(i)=mean(tmp);
            end
          end
          Elog2 = Elog2.^2;
          Vn = (Elog-Elog2);
          if strcmp(form, 'mean')
            Vn = mean(Vn);
            BUt = mean(BUt);
          end
          waic = BUt - Vn;
        end

      else
        % Evaluate WAIC using the expected value form via Gibbs training
        % loss

        if isfield(gp.lik.fh,'trcov')
          % Gaussian likelihood
          for i=1:tn
            fmin = mean(Ef(i,:) - 9*sqrt(Varf(i,:)));
            fmax = mean(Ef(i,:) + 9*sqrt(Varf(i,:)));
            GUt(i) = quadgk(@(f) mean(multi_npdf(f,Ef(i,:),(Varf(i,:))) ...
                                      .*bsxfun(@minus,-bsxfun(@rdivide,(repmat((y(i)-f),nsamples,1)).^2,(2.*sigma2(i,:))'), 0.5*log(2*pi*sigma2(i,:))')), fmin, fmax);
          end
          if strcmp(form, 'mean')
            GUt = mean(GUt);
            BUt = mean(BUt);
          end
          waic = BUt-2*(BUt-GUt);
        else
          % non-Gaussian likelihood
          for i=1:tn
            if ~isempty(z)
              z1 = z(i);
            else
              z1 = [];
            end
            fmin = mean(Ef(i,:) - 9*sqrt(Varf(i,:)));
            fmax = mean(Ef(i,:) + 9*sqrt(Varf(i,:)));
            GUt(i) = quadgk(@(f) mean(multi_npdf(f,Ef(i,:),(Varf(i,:))) ...
                                      .*llvec(gp_array, y(i), f, z1)), fmin, fmax);
          end
          if strcmp(form, 'mean')
            GUt = mean(GUt);
            BUt = mean(BUt);
          end
          waic = BUt-2*(BUt-GUt);
        end
      end

    else
      % A single GP solution
      [Ef, Varf, BUt] = gp_pred(gp, x, y, x, 'yt', y, 'tstind', tstind, options);

      GUt = zeros(tn,1);
      Elog = zeros(tn,1);
      Elog2 = zeros(tn,1);

      if isequal(method,'V')
        % Estimate WAIC with variance form

        if isfield(gp.lik.fh,'trcov')
          % Gaussian likelihood
          sigma2 = gp.lik.sigma2;

          for i=1:tn

            % Analytical moments for Gaussian distribution

            m0 = 1; m1 = Ef(i); m2 = Ef(i)^2 + Varf(i); m3 = Ef(i)*(Ef(i)^2+3*Varf(i));
            m4 = Ef(i)^4+6*Ef(i)^2*Varf(i)+3*Varf(i)^2;

            Elog2(i) = (-0.5*log(2*pi*sigma2) - y(i).^2./(2.*sigma2))*m0 - 1./(2.*sigma2) * m2 + y(i)./sigma2 * m1;
            Elog(i) = (1/4 * m4 - y(i) * m3 + (3*y(i).^2./2+0.5*log(2*pi*sigma2).*sigma2) * m2 ...
                       - (y(i).^3 + y(i).*log(2*pi*sigma2).*sigma2) * m1 + (y(i).^4/4 + 0.5*y(i).^2*log(2*pi*sigma2).*sigma2 ...
                                                              + 0.25*log(2*pi*sigma2).^2.*sigma2.^2) * m0) ./ sigma2.^2;

          end
          Elog2 = Elog2.^2;
          Vn = Elog-Elog2;
          if strcmp(form,'mean')
            BUt = mean(BUt);
            Vn = mean(Vn);
          end
          waic = BUt - Vn;

        else
          % Non-Gaussian likelihood
          for i=1:tn
            if ~isempty(z)
              z1 = z(i);
            else
              z1 = [];
            end
            if ~isequal(gp.lik.type, 'Coxph')
              fmin = Ef(i)-9*sqrt(Varf(i));
              fmax = Ef(i)+9*sqrt(Varf(i));
              Elog(i) = quadgk(@(f) norm_pdf(f, Ef(i), sqrt(Varf(i))).*llvec(gp, y(i), f, z1).^2 ,...
                fmin, fmax);
              Elog2(i) = quadgk(@(f) norm_pdf(f, Ef(i), sqrt(Varf(i))).*llvec(gp, y(i), f, z1) ,...
                fmin, fmax);
            else
              % Use MC to integrate over latents
              ntime = size(gp.lik.xtime,1);
              ns = 10000;
              Sigma_tmp = Varf([1:ntime ntime+i], [1:ntime ntime+i]);
              Sigma_tmp = (Sigma_tmp + Sigma_tmp') ./ 2;
              f = mvnrnd(Ef([1:ntime ntime+i]), Sigma_tmp, ns);
              Elog2(i) = 1/ns * sum(llvec(gp, y(i,:), f', z1));
              Elog(i) = 1/ns * sum((llvec(gp, y(i,:), f', z1)).^2);
            end
          end
          Elog2 = Elog2.^2;
          Vn = Elog-Elog2;
          if strcmp(form, 'mean')
            Vn = mean(Vn);
            BUt = mean(BUt);
          end
          waic = BUt - Vn;

        end

      else
        % WAIC using the expected value form via Gibbs training loss GUt

        if isfield(gp.lik.fh,'trcov')
          % Gaussian likelihood
          sigma2 = gp.lik.sigma2;
          for i=1:tn
            if Varf(i)<eps
              GUt(i)=(-0.5*log(2*pi*sigma2)- (y(i) - Ef(i)).^2/(2.*sigma2));
            else

              % GUt(i) = quadgk(@(f) norm_pdf(f,Ef(i),sqrt(Varf(i))).*(-0.5*log(2*pi*sigma2)- (y(i) - f).^2/(2.*sigma2)), fmin, fmax);

              m0 = 1; m1 = Ef(i); m2 = Ef(i)^2 + Varf(i);

              GUt(i) = (-0.5*log(2*pi*sigma2) - y(i).^2./(2.*sigma2))*m0 - 1./(2.*sigma2) * m2 + y(i)./sigma2 * m1;
            end
          end
          if strcmp(form,'mean')
            GUt = mean(GUt);
            BUt = mean(BUt);
          end
          waic = BUt-2*(BUt-GUt);
        else
          % Non-Gaussian likelihood
          for i=1:tn
            if ~isempty(z)
              z1 = z(i);
            else
              z1 = [];
            end
            if ~isequal(gp.lik.type, 'Coxph')
              fmin = Ef(i)-9*sqrt(Varf(i));
              fmax = Ef(i)+9*sqrt(Varf(i));
              GUt(i) = quadgk(@(f) norm_pdf(f, Ef(i), sqrt(Varf(i))).*llvec(gp, y(i), f, z1) ,...
                fmin, fmax);
            else
              % If likelihood coxph use mc to integrate over latents
              ntime = size(gp.lik.xtime,1);
              ns = 10000;
              Sigma_tmp = Varf([1:ntime ntime+i], [1:ntime ntime+i]);
              Sigma_tmp = (Sigma_tmp + Sigma_tmp') ./ 2;
              f = mvnrnd(Ef([1:ntime ntime+i]), Sigma_tmp, ns);
              GUt(i) = 1/ns * sum(llvec(gp, y(i), f', z1));
            end
          end
          if strcmp(form,'mean')
            GUt = mean(GUt);
            BUt = mean(BUt);
          end
          waic = BUt-2*(BUt-GUt);
        end

      end

    end

  elseif iscell(gp)

    % gp_ia solution

    switch gp{1}.type
      case {'FULL' 'VAR' 'DTC' 'SOR'}
        tstind = [];
      case {'FIC' 'CS+FIC'}
        tstind = 1:tn;
      case 'PIC'
        tstind = gp{1}.tr_index;
    end

    [tmp, tmp, BUt] = gp_pred(gp,x,y, x, 'yt', y, 'tstind', tstind, options);
    GUt = zeros(tn,1);
    Elog = zeros(tn,1);
    Elog2 = zeros(tn,1);

    nsamples = length(gp);
    for j = 1:nsamples
      Gp = gp{j};
      weight(j) = Gp.ia_weight;
      w(j,:) = gp_pak(Gp);
      [Ef(:,j), Varf(:,j)] = gp_pred(Gp, x, y, x, 'yt', y, 'tstind', tstind, options);
      if isfield(Gp.lik.fh,'trcov')
        sigma2(:,j) = repmat(Gp.lik.sigma2,1,tn);
      end
    end
    if isequal(method,'V')
      % Evaluate WAIC using the variance form

      if isfield(gp{1}.lik.fh,'trcov')
        % Gaussian likelihood
        for i=1:tn
          fmin = sum(weight.*Ef(i,:) - 9*weight.*sqrt(Varf(i,:)));
          fmax = sum(weight.*Ef(i,:) + 9*weight.*sqrt(Varf(i,:)));
          Elog(i) = quadgk(@(f) sum(bsxfun(@times, multi_npdf(f,Ef(i,:),(Varf(i,:))),weight') ...
                                    .*bsxfun(@minus,-bsxfun(@rdivide,(repmat((y(i)-f),nsamples,1)).^2,(2.*sigma2(i,:))'), 0.5*log(2*pi*sigma2(i,:))').^2), fmin, fmax);
          Elog2(i) = quadgk(@(f) sum(bsxfun(@times, multi_npdf(f,Ef(i,:),(Varf(i,:))),weight') ...
                                     .*bsxfun(@minus,-bsxfun(@rdivide,(repmat((y(i)-f),nsamples,1)).^2,(2.*sigma2(i,:))'), 0.5*log(2*pi*sigma2(i,:))')), fmin, fmax);
        end
        Elog2 = Elog2.^2;
        Vn = (Elog-Elog2);
        if strcmp(form, 'mean')
          Vn = mean(Vn);
          BUt = mean(BUt);
        end
        waic = BUt - Vn;
      else
        % non-Gaussian likelihood
        for i=1:tn
          if ~isempty(z)
            z1 = z(i);
          else
            z1 = [];
          end
          fmin = sum(weight.*Ef(i,:) - 9*weight.*sqrt(Varf(i,:)));
          fmax = sum(weight.*Ef(i,:) + 9*weight.*sqrt(Varf(i,:)));
          Elog(i) = quadgk(@(f) sum(bsxfun(@times, multi_npdf(f,Ef(i,:),(Varf(i,:))),weight') ...
                                    .*llvec(gp, y(i), f, z1).^2), fmin, fmax);
          Elog2(i) = quadgk(@(f) sum(bsxfun(@times, multi_npdf(f,Ef(i,:),(Varf(i,:))),weight') ...
                                     .*llvec(gp, y(i), f, z1)), fmin, fmax);
        end
        Elog2 = Elog2.^2;
        Vn = (Elog-Elog2);
        if strcmp(form, 'mean')
          Vn = mean(Vn);
          BUt = mean(BUt);
        end
        waic = BUt - Vn;

      end

    else
      % Evaluate WAIC using the expected value form via Gibbs training loss

      if isfield(gp{1}.lik.fh,'trcov')
        % Gaussian likelihood
        for i=1:tn
          fmin = sum(weight.*Ef(i,:) - 9*weight.*sqrt(Varf(i,:)));
          fmax = sum(weight.*Ef(i,:) + 9*weight.*sqrt(Varf(i,:)));
          GUt(i) = quadgk(@(f) sum(bsxfun(@times, multi_npdf(f,Ef(i,:),(Varf(i,:))),weight') ...
                                   .*bsxfun(@minus,-bsxfun(@rdivide,(repmat((y(i)-f),nsamples,1)).^2,(2.*sigma2(i,:))'), 0.5*log(2*pi*sigma2(i,:))')), fmin, fmax);
        end
        if strcmp(form, 'mean')
          GUt = mean(GUt);
          BUt = mean(BUt);
        end
        waic = BUt-2*(BUt-GUt);

      else
        % non-gaussian likelihood
        for i=1:tn
          if ~isempty(z)
            z1 = z(i);
          else
            z1 = [];
          end
          fmin = sum(weight.*Ef(i,:) - 9*weight.*sqrt(Varf(i,:)));
          fmax = sum(weight.*Ef(i,:) + 9*weight.*sqrt(Varf(i,:)));
          GUt(i) = quadgk(@(f) sum(bsxfun(@times, multi_npdf(f,Ef(i,:),(Varf(i,:))),weight') ...
                                   .*llvec(gp, y(i), f, z1)), fmin, fmax);
        end
        if strcmp(form, 'mean')
          GUt = mean(GUt);
          BUt = mean(BUt);
        end
        waic = BUt-2*(BUt-GUt);

      end
    end

  end

end

function lls=llvec(gp, y, fs, z)
% Compute a vector of lls for vector argument fs used by quadgk. In
% case of IA or MC, return a matrix with rows corresponding to one
% GP and columns corresponding to all of the GP's.

  if isstruct(gp)
      % single gp
      lls=zeros(1,size(fs,2));
      for i1=1:size(fs,2)
        lls(i1)=gp.lik.fh.ll(gp.lik,y,fs(:,i1),z);
      end
      %     else
      %       % mc
      %       lls=zeros(length(gp), length(fs));
      %       for i=1:numel(fs)
      %         for j=1:numel(gp.edata)
      %           Gp = take_nth(gp, j);
      %           lls(j,i) = Gp.lik.fh.ll(Gp.lik, y, fs(i), z);
      %         end
      %       end
  else
    % ia & mc
    lls=zeros(length(gp), length(fs));
    for i=1:numel(fs)
      for j=1:numel(gp)
        lls(j,i) = gp{j}.lik.fh.ll(gp{j}.lik, y, fs(i), z);
      end
    end
  end
end

function mpdf = multi_npdf(f, mean, sigma2)
% for every element in f, compute means calculated with
% norm_pdf(f(i), mean, sqrt(sigma2)). If mean and sigma2
% are vectors, returns length(mean) x length(f) matrix.

  mpdf = zeros(length(mean), length(f));
  for i=1:length(f)
    mpdf(:,i) = norm_pdf(f(i), mean, sqrt(sigma2));
  end
end

function [m_0, m_1, m_2, m_3, m_4] = moments(fun, a, b, rtol, atol, minsubs)
% QUAD_MOMENTS Calculate the 0th, 1st and 2nd moment of a given
%              (unnormalized) probability distribution
%
%   [m_0, m_1, m_2] = quad_moments(fun, a, b, varargin)
%   Inputs:
%      fun  = Function handle to the unnormalized probability distribution
%      a,b  = integration limits [a,b]
%      rtol = relative tolerance for the integration (optional, default 1e-6)
%      atol = absolute tolerance for the integration (optional, default 1e-10)
%
%   Returns the first three moments:
%      m0  = int_a^b fun(x) dx
%      m1  = int_a^b x*fun(x) dx / m0
%      m2  = int_a^b x^2*fun(x) dx / m0
%
%   The function uses an adaptive Gauss-Kronrod quadrature. The same set of
%   integration points and intervals are used for each moment. This speeds up
%   the evaluations by factor 3, since the function evaluations are done only
%   once.
%
%   The quadrature method is described by:
%   L.F. Shampine, "Vectorized Adaptive Quadrature in Matlab",
%   Journal of Computational and Applied Mathematics, 211, 2008,
%   pp. 131-140.

%   Copyright (c) 2010 Jarno Vanhatalo, Jouni Hartikainen

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

  maxsubs = 650;

  if nargin < 4
    rtol = 1.e-6;
  end
  if nargin < 5
    atol = 1.e-10;
  end
  if nargin < 6
    minsubs = 10;
  end

  rtol = max(rtol,100*eps);
  atol = max(atol,0);
  minsubs = max(minsubs,2); % At least two subintervals are needed

  % points and weights
  points15 = [0.2077849550078985; 0.4058451513773972; 0.5860872354676911; ...
              0.7415311855993944; 0.8648644233597691; 0.9491079123427585; ...
              0.9914553711208126];
  points = [-points15(end:-1:1); 0; points15];

  w15 = [0.2044329400752989, 0.1903505780647854, 0.1690047266392679, ...
         0.1406532597155259, 0.1047900103222502, 0.06309209262997855, ...
         0.02293532201052922];
  w = [w15(end:-1:1), 0.2094821410847278, w15];

  w7 = [0,0.3818300505051189,0,0.2797053914892767,0,0.1294849661688697,0];
  ew = w - [w7(end:-1:1), 0.4179591836734694, w7];

  samples = numel(w);

  % split the interval.
  if b-a <= 0
    c = a; a = b; b=c;
    warning('The start of the integration interval was less than the end of it.')
  end
  apu = a + (1:(minsubs-1))./minsubs*(b-a);
  apu = [a,apu,b];
  subs = [apu(1:end-1);apu(2:end)];

  % Initialize partial sums.
  Ifx_ok = 0;
  Ifx1_ok = 0;
  Ifx2_ok = 0;
  Ifx3_ok = 0;
  Ifx4_ok = 0;
  % The main loop
  while true
    % subintervals and their midpoints
    midpoints = sum(subs)/2;
    halfh = diff(subs)/2;
    x = bsxfun(@plus,points*halfh,midpoints);
    x = reshape(x,1,[]);

    fx = fun(x);
    fx1 = fx.*x;
    fx2 = fx.*x.^2;
    fx3 = fx.*x.^3;
    fx4 = fx.*x.^4;

    fx = reshape(fx,samples,[]);
    fx1 = reshape(fx1,samples,[]);
    fx2 = reshape(fx2,samples,[]);
    fx3 = reshape(fx3,samples,[]);
    fx4 = reshape(fx4,samples,[]);

    % Subintegrals.
    Ifxsubs = (w*fx) .* halfh;
    errsubs = (ew*fx) .* halfh;
    Ifxsubs1 = (w*fx1) .* halfh;
    Ifxsubs2 = (w*fx2) .* halfh;
    Ifxsubs3 = (w*fx3) .* halfh;
    Ifxsubs4 = (w*fx4) .* halfh;

    % Ifx and tol.
    Ifx = sum(Ifxsubs) + Ifx_ok;
    Ifx1 = sum(Ifxsubs1) + Ifx1_ok;
    Ifx2 = sum(Ifxsubs2) + Ifx2_ok;
    Ifx3 = sum(Ifxsubs3) + Ifx3_ok;
    Ifx4 = sum(Ifxsubs4) + Ifx4_ok;
    tol = max(atol,rtol*abs(Ifx));

    % determine the indices ndx of Ifxsubs for which the
    % errors are acceptable and remove those from subs
    ndx = find(abs(errsubs) <= (2/(b-a)*halfh*tol));
    subs(:,ndx) = [];
    if isempty(subs)
      break
    end

    % Update the integral.
    Ifx_ok = Ifx_ok + sum(Ifxsubs(ndx));
    Ifx1_ok = Ifx1_ok + sum(Ifxsubs1(ndx));
    Ifx2_ok = Ifx2_ok + sum(Ifxsubs2(ndx));
    Ifx3_ok = Ifx3_ok + sum(Ifxsubs3(ndx));
    Ifx4_ok = Ifx4_ok + sum(Ifxsubs4(ndx));

    % Quit if too many subintervals.
    nsubs = 2*size(subs,2);
    if nsubs > maxsubs
      warning('quad_moments: Reached the limit on the maximum number of intervals in use.');
      break
    end
    midpoints(ndx) = [];
    subs = reshape([subs(1,:); midpoints; midpoints; subs(2,:)],2,[]); % Divide the remaining subintervals in half
  end

  % Scale moments
  m_0 = Ifx;
  m_1 = Ifx1./Ifx;
  m_2 = Ifx2./Ifx;
  m_3 = Ifx3./Ifx;
  m_4 = Ifx4./Ifx;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [w, s, h] = gp_pak(gp, param)
%GP_PAK  Combine GP parameters into one vector
%
%  Description
%    W = GP_PAK(GP, PARAM) takes a Gaussian Process structure
%    GP and string PARAM defining, which parameters are packed and
%    combines the parameters into a single row vector W. If PARAM
%    is not given the function packs all parameters.
%
%    Each of the following strings in PARAM defines one group of
%    parameters to pack:
%      covariance  - pack parameters of covariance function
%      likelihood  - pack parameters of likelihood
%      inducing    - pack inducing inputs (in sparse approximations):
%                    W = gp.X_u(:)
%
%    By combining the strings one can pack more than one group of
%    parameters. For example:
%      covariance+inducing  - pack covariance function parameters
%                             and inducing inputs
%      covariance+likelih   - pack covariance function parameters
%                             and likelihood parameters
%
%    Inside each group (such as covariance functions) the
%    parameters to be packed is defined by the existence of a prior
%    structure. For example, if GP has two covariance functions but
%    only the first one has prior for its parameters then only the
%    parameters of the first one are packed. Thus, also inducing
%    inputs require prior if they are to be optimized.
%
%    GP_PAK and GP_UNPAK functions are used, e.g., when GP
%    parameters are optimized with GP_OPTIM or sampled with GP_MC.
%    See GP_SET and option 'infer_params'.
%
%    [W, WS] = GP_PAK(GP, PARAM) returns also cell array of string
%    labels for the weight vector elements, which makes diagnostics
%    easier.
%
%    [W, WS, H] = GP_PAK(GP, PARAM) returns also hierarchy level H of
%    different parameters (0 for likelihood parameters and 1 for covariance
%    function parameters)
%
%  See also
%    GP_UNPAK, GP_SET
%
% Copyright (c) 2007-2010 Jarno Vanhatalo

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

  w = []; s = {}; h = [];

  if isfield(gp,'etr') && length(gp.etr) > 1
    if strcmp(gp.type, 'PIC_BLOCK') || strcmp(gp.type, 'PIC')
      ind = gp.tr_index;           % block indeces for training points
      gp = rmfield(gp,'tr_index');
    end
    ns = length(gp.etr);
    for i1 = 1:ns
      Gp = take_nth(gp,i1);
      [w(i1,:), s, h] = gp_pak(Gp);
    end
  else

    if nargin < 2
      param = gp.infer_params;
    end

    % Pack SVI variational parameters
    if isfield(gp, 'latent_method') && isequal(gp.latent_method,'SVI') ...
        && isfield(gp, 't1')
      w=[gp.t1; gp.t2(:)]';
      if length(gp.t1)>1
        s = [s; sprintf('E[u] x %d', gp.nind); ...
          sprintf('Cov[u] x %d', gp.nind^2)];
      else
        s = [s; 'E[u]'; 'Cov[u]'];
      end
    end

    % Pack the parameters of covariance functions
    if ~isempty(strfind(param, 'covariance'))
      ncf = length(gp.cf);

      for i=1:ncf
        gpcf = gp.cf{i};
        [wi, si, hi] = gpcf.fh.pak(gpcf);
        w = [w wi];
        s = [s; si];
        h = [h hi];
      end
    end

    % Pack the parameters of likelihood function
    if ~isempty(strfind(param, 'likelihood'))
      [wi, si, hi] = gp.lik.fh.pak(gp.lik);
      if isfield(gp, 'latent_method') && isequal(gp.latent_method, 'SVI') ...
          && ~isequal(gp.lik.type, 'Gaussian') && isfield(gp, 't1') ...
          && ~isempty(gp.lik.p.sigma2)
        wi=log(gp.lik.sigma2);
        si='log(sigma2)';
        hi=0;
      end
      w = [w wi];
      s = [s; si];
      h = [h hi];
    end

    % Pack the parameters of the second likelihood function (monotonic)
    if ~isempty(strfind(param, 'likelihood')) && isfield(gp, 'lik_mono')
      [wi, si, hi] = gp.lik_mono.fh.pak(gp.lik_mono);
      if isfield(gp, 'latent_method') && isequal(gp.latent_method, 'SVI') ...
           && isfield(gp, 't1') && ~isempty(gp.lik_mono.p.sigma2)
        wi=log(gp.lik_mono.sigma2);
        si='log(sigma2_mono)';
        hi=0;
      end
      w = [w wi];
      s = [s; si];
      h = [h hi];
    end

    % Pack the inducing inputs
    if ~isempty(strfind(param, 'inducing'))
      if isfield(gp,'p') && isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
        if ~iscell(gp.p.X_u)
          % One prior for all inducing inputs
          w = [w reshape(gp.X_u', numel(gp.X_u),1)'];
%           w = [w gp.X_u'];
          s = [s; sprintf('inducing x %d',numel(gp.X_u))];
          h = [h ones(1,numel(gp.X_u))];
          [wi,si,hi]=gp.p.X_u.fh.pak(gp.p.X_u);
          w = [w wi];
          s = [s; si];
          h = [h 1+hi];
        else
          % Own prior for each inducing input
          for i=1:size(gp.X_u,1)
            w = [w gp.X_u(i,:)];
            s = [s; sprintf('inducing x %d',numel(gp.X_u(i,:)))];
            h = [h ones(1,size(gp.X_u(i,:),2))];
            [wi,si,hi]=gp.p.X_u{i}.fh.pak(gp.p.X_u{i});
            si=strcat(repmat('prior-', size(si,1),1),si);
            w = [w wi];
            s = [s; si];
            h = [h 1+hi];
          end
        end
      end
    end

    % Pack the prior weights and variances of mean functions
    if ~isempty(strfind(param, 'mean'))
      mf = length(gp.meanf);

      for i=1:mf
        gpmf = gp.meanf{i};
        [wi, si, hi] = gpmf.fh.pak(gpmf);
        w = [w wi];
        s = [s; si];
        h = [h, hi];
      end
    end

  end

end

author: gpstuff-dev

repo: gpstuff

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


text:
function i2vec_print ( n, a, title )

%*****************************************************************************80
%
%% I2VEC_PRINT prints an I2VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of components of the vector.
%
%    Input, integer A(2,N), the vector to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d: %10d  %10d\n', i, a(1:2,i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [dfdz,dfdp,sym_dfdz,sym_dfdp] = symDynJac(ode,nstates,nparam)
%SYMJAC  Returns Symbolically Differentiated Partial Derivatives of a Dynamic System
%
%   [dfdz,dfdp] = symDynJac(ode) uses the symbolic toolbox to automatically
%   generate the sensitivity partial derivatives of the function handle ode.
%   You should supply the function handle in the form @(t,z,p) z(1)^2 + p(1)
%   noting t, z and p must be the only variables used, and are indexed in the
%   equation (no vector operations).
%
%   [dfdz,dfdp] = symDynJac(ode,nstates) specifies the number of states in
%   the equation, assuming consecutive ordering. Useful if a state is not
%   specified in the original equation to pad DFDZ with zeros.
%
%   [dfdz,dfdp] = symDynJac(ode,nstates,nparam) specifies the number of
%   parameters in the equation, assuming consecutive ordering. Useful if a
%   parameter is not specified in the original equation to pad DFDP with
%   zeros.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

if(nargin < 3), nparam = 0; end
if(nargin < 2), nstates = 0; end

if(~optiCheckSymTBX())
    dfdz = []; sym_dfdz = [];
    dfdp = []; sym_dfdp = [];
    return
end

if(~isa(ode,'function_handle') && ~isa(ode,'barvec'))
    error('Fun should be a function handle!');
end
if(isa(ode,'function_handle') && nargin(ode) ~= 3)
    error('ODE should only have three input arguments (t,z,p)');
end

%Convert ODE to a symbolic expression
if(isa(ode,'function_handle'))
    [symode,ind] = func2sym(ode,{'z','p','t'});
    indz = ind{1}; indp = ind{2};
else
    symode = sym(getEq(ode)); %convert from barvec to symbolic expression
    v = char(symvar(symode));
    indz = 1:length(strfind(v,'z'));
    indp = 1:length(strfind(v,'p'));
end

%Create each symbolic partial derivative
sym_dfdz = symPartialDer(symode,'z',nstates,indz);
sym_dfdp = symPartialDer(symode,'p',nparam,indp);

%Return to function handles (note order of args important and changed)
dfdz = sym2func(sym_dfdz,{'t','z','p'});
dfdp = sym2func(sym_dfdp,{'t','z','p'});

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function v=mocov_util_platform_version()
% return the version of Matlab or Octave
%
% v=mocov_util_platform_version()
%
% Output:
%   v               Vector with two (Matlab) or three (Octave) elements,
%                   for example 8.5 (for Matlab) or 4.0.3) for Octave)
%

    version_str=version();

    parts=regexp(version_str,'\s','split');
    first_part=parts{1};

    num_parts=regexp(first_part,'\.','split');
    v=cellfun(@str2num,num_parts);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef TestFisheyeStereoRectify
    %TestFisheyeStereoRectify

    methods (Static)
        function test_1
            %TODO
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function test08 ( input_file_name )

%*****************************************************************************80
%
%% TEST08 tests STLA_READ and STLA_READ_FAST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  Compare two routines with the same functionality.\n' );
  fprintf ( 1, '  One is designed to run faster than the other,\n' );
  fprintf ( 1, '  especially on large problems with thousands of faces.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  STLA_READ reads various data from an ASCII STL file.\n' );
  fprintf ( 1, '  STLA_READ_FAST reads data from an ASCII STL file.\n' );
  fprintf ( 1, '\n' );

  [ solid_num, node_num, face_num, text_num ] = stla_size_fast ( ...
    input_file_name );

  t1 = clock;

  [ node_xyz, face_node, face_normal ] = stla_read_fast ( ...
    input_file_name, node_num, face_num );

  t2 = clock;

  diff1 = etime ( t2, t1 );

  fprintf ( 1, '  STLA_READ_FAST read the file in %f seconds.\n', diff1 );

  [ solid_num, node_num, face_num, text_num ] = stla_size ( ...
    input_file_name );

  t3 = clock;

  [ node_xyz, face_node, face_normal ] = stla_read ( ...
    input_file_name, node_num, face_num );

  t4 = clock;

  diff2 = etime ( t4, t3 );

  fprintf ( 1, '  STLA_READ      read the file in %f seconds.\n', diff2 );

  stla_size_print ( input_file_name, solid_num, node_num, face_num, ...
    text_num );

  return
end


text:
function [ score ] = compute_cosine_score( feature1, feature2 )
% compute_cosine_score Summary of this function goes here
%   This function is to compute the cosine similarity between two vectors
%
    score = feature1'*feature2/norm(feature1)/norm(feature2);

end

author: AlfredXiangWu


text:
%% Vectorization using Sparse Matrices
%
%% Vectorization

%% Sparse matrix

%% Example: assembling of discrete Laplacian operator

%% Example: node star

%% Example: eliminating hanging nodes

repo: Aero-matlab


text:
function OUTPUT = readsacfile(varargin)
   %readsacfile    Read SAC binary files.
   %    usage:  output = readsacfile('sacfile')
   %         OUTPUT is a struct with fields "header" and "amplitudes"
   %         OUTPUT.header is an array of cells, each parsed into the proper
   %         data type (ie, SINGLE, CHAR, LOGICAL, etc.).
   %         OUTPUT.amplitudes is an array of the measured data

   % VERSION: 1.1 of waveform objects
   % AUTHOR: Celso Reyes   creyes@gi.alaska.edu
   %    modified from Michael Thorne (4/2004)
   % LASTUPDATE: 9/2/2009

   try
      OUTPUT = mainFunction('little-endian',varargin{:});
   catch
      % disp('big-endian')
      OUTPUT = mainFunction('big-endian',varargin{:});
   end
end
function OUTPUT = mainFunction(endian,varargin)
   for nrecs = 1:nargin-1

      sacfile = varargin{nrecs};

      %---------------------------------------------------------------------------
      %    Default byte-order
      %    endian  = 'big-endian' byte order (e.g., UNIX)
      %            = 'little-endian' byte order (e.g., LINUX)

      %endian = 'little-endian';
      %endian = 'big-endian';

      if strcmp(endian,'big-endian')
         fid = fopen(sacfile,'r','ieee-be');
      elseif strcmp(endian,'little-endian')
         fid = fopen(sacfile,'r','ieee-le');
      end

      %preallocate
      h = cell(1,123);

      % read in single precision real header variables:
      %---------------------------------------------------------------------------
      for i=1:70
         h(i) = {fread(fid,1,'single')};
      end

      % read in single precision integer header variables:
      %---------------------------------------------------------------------------
      for i=71:105
         h(i) = {fread(fid,1,'int32')};
      end

      % Check header version = 6 and issue warning
      %---------------------------------------------------------------------------
      % If the header version is not NVHDR == 6 then the sacfile is likely of the
      % opposite byte order.  This will give h(77) some ridiculously large
      % number.  NVHDR can also be 4 or 5.  In this case it is an old SAC file
      % and rsac cannot read this file in.  To correct, read the SAC file into
      % the newest verson of SAC and w over.
      %
      if (h{77} == 4 || h{77} == 5)
         message = strcat('NVHDR = 4 or 5. File: "',sacfile,'" may be from an old version of SAC.');
         fclose(fid);
         error('Waveform:readsacfile:oldSacVersion',message)
      elseif h{77} ~= 6
         message = strcat('Current readsacfile byte order: "',endian,'". File: "',sacfile,'" may be of opposite byte-order.');
         fclose(fid);
         error('Waveform:readsacfile:wrongEndian',message)
      end

      % read in logical header variables
      %---------------------------------------------------------------------------
      for i=106:110
         h(i) = {fread(fid,1,'int32')};
      end

      % read in character header variables
      %---------------------------------------------------------------------------
      strlens = repmat(8,23,1); strlens(2) = 16;
      for i=111:133;%111:8:302
         h(i) = {(fread(fid,strlens(i-110),'char=>char'))'};
      end

      % read in amplitudes
      %---------------------------------------------------------------------------

      OUTPUT(nrecs).amplitudes = fread(fid,'single');
      fclose(fid);

      % arrange output files
      %------------------------------------------------------------------------
      %---
      OUTPUT(nrecs).header = h(:)';

   end
end

author: geoscience-community-codes

repo: GISMO


text:
function c = makesynthwaves(n);

% Add N traces of synthetic data to a correlation object. Alters properties: waves, trig,
% start, Fs.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% MAKE SINGLE SIGNAL
t = [-5:.01:14.99]';
env_orig = exp(-1*(t-5).^2/ (2*3^2) );
sprep = 0.2*sin(t+1) + .3*sin((t-3)/0.3) + 0.4*cos(t/0.8) + 0.2*cos((t-1.7)/0.2);
%sprep = sin(t);
%sprep = 0.5*sin(t) + 0.5*sin(t/.5);
s = sprep .* env_orig;

% MAKE SUITE OF SIMILAR EVENT
%n = 150; %no. of events
w = s * ones(1,n) + .01*(rand(length(s),n)-.5); % add time offset
aa = .6;    % adjust to add randomness
for i = 1:n
    w(:,i) = w(:,i) + aa*rand(1)*sin((t-rand(1))/rand(1));
end;

% VARY START TIMES
for i = 1:n
    bump = round(300*(rand(1)));
    w(:,i) = w([ length(t)-bump:length(t) 1:(length(t)-1)-bump],i);
end;

% VARY AMPLITUDES
for i = 1:n
    w(:,i) = rand * 100 * w(:,i);
end;

% DEMEAN
w = w - ones(size(w,1),1)*mean(w);

% MAKE v0 CORRELATION
d.trig = 732604 + rand(n,1);
d.start = d.trig - (5 + 0.2*rand(n,1))/86400;
d.Fs = 100;
d.w = w;

% CREATE v1 CORRELATION OBJECT
c.W = waveform;
for i = 1:length(d.trig)
    w = waveform;
    w = set(w,'Station','UNKN');
    w = set(w,'Channel','EHZ');
    w = set(w,'Start',d.start(i));
    w = set(w,'Fs',d.Fs);
    w = set(w,'Data',d.w(:,i));
    c.W(i) = w;
end;
c.W = reshape(c.W,length(c.W),1);
c.trig = d.trig;


text:
function [x,t] = tfristft(tfr,t,h,trace);
%TFRISTFT Inverse Short time Fourier transform.
%	[X,T]=TFRSTFT(tfr,T,H,TRACE) computes the inverse short-time
%	Fourier transform of a discrete-time signal X. This function
%	may be used for time-frequency synthesis of signals.
%
%	X     : signal.
%	T     : time instant(s)          (default : 1:length(X)).
%	H     : frequency smoothing window, H being normalized so as to
%	        be  of unit energy.      (default : Hamming(N/4)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR   : time-frequency decomposition (complex values). The
%	        frequency axis is graduated from -0.5 to 0.5.
%
%	Example :
%        t=200+(-128:127); sig=[fmconst(200,0.2);fmconst(200,0.4)];
%        h=hamming(57); tfr=tfrstft(sig,t,256,h,1);
%        sigsyn=tfristft(tfr,t,h,1);
%        plot(t,abs(sigsyn-sig(t)))
%
%	See also all the time-frequency representations listed in
%	the file CONTENTS (TFR*)

%	F. Auger, November 1996.
%
%  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

if (nargin<3),
 error('At least 3 parameters required');
elseif (nargin==3),
 trace=0;
end;

[N,NbPoints]=size(tfr);
[trow,tcol] =size(t);
[hrow,hcol] =size(h); Lh=(hrow-1)/2;

if (trow~=1),
 error('T must only have one row');
elseif (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
elseif (tcol~=NbPoints)
 error('tfr should have as many columns as t has rows.');
end;

Deltat=t(2:tcol)-t(1:tcol-1);
Mini=min(Deltat); Maxi=max(Deltat);
if (Mini~=1) & (Maxi~=1),
 error('The tfr must be computed at each time sample.');
end;

h=h/norm(h);

if trace, disp('Inverse Short-time Fourier transform'); end;
tfr=ifft(tfr);

x=zeros(tcol,1);

for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 valuestj=max([1,icol-N/2,icol-Lh]):min([tcol,icol+N/2,icol+Lh]);
 for tj=valuestj,
  tau=icol-tj; indices= rem(N+tau,N)+1;
  % fprintf('%g %g %g\n',tj,tau,indices);
  x(icol,1)=x(icol,1)+tfr(indices,tj)*h(Lh+1+tau);
 end;
 x(icol,1)=x(icol,1)/sum(abs(h(Lh+1+icol-valuestj)).^2);
end;

if trace, fprintf('\n'); end;

author: HeLiangHIT


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% change log:
%   040204 smw stub file
%
% file: get_ci_cs1_1s.m
%
function varargout = get_ci_cs1_1s(varargin)
  error('get_ci_cs1_1s is not part of current release.');
return;


text:
function rank = gray_code_rank ( n, t )

%*****************************************************************************80
%
%% GRAY_CODE_RANK computes the rank of a Gray code element.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of digits in each element.
%    N must be positive.
%
%    Input, integer T(N), an element of the Gray code.
%    Each entry T(I) is either 0 or 1.
%
%    Output, integer RANK, the rank of the element.
%

%
%  Check.
%
  ierror = gray_code_check ( n, t );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GRAY_CODE_RANK - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'GRAY_CODE_RANK - Fatal error!' );
  end

  rank = 0;
  b = 0;

  for i = n - 1 : -1 : 0

    if ( t(n-i) ~= 0 )
      b = 1 - b;
    end

    if ( b == 1 )
      rank = rank + 2^i;
    end

  end

  return
end


text:
function hFig = view_scouts(ResultsFiles, ScoutsArg, hFig)
% VIEW_SCOUTS: Display time series for all the scouts selected in the JList.
%
% USAGE:  view_scouts()                               : Display selected sources file time series for selected scouts
%         view_scouts(ResultsFiles, 'SelectedScouts') : Display input sources file time series for selected scouts
%         view_scouts(ResultsFiles, iScouts)          : Display input sources file time series for input scouts

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

global GlobalData;
%% ===== PARSE INPUTS =====
if (nargin < 3) || isempty(hFig)
    hFig = [];
end
% If no parameters were given
if (nargin == 0)
    % === GET SURFACES LIST ===
    hFigures = bst_figures('GetFiguresByType', '3DViz');
    % Process all figures : and keep them that have a ResultsFile defined
    ResultsFiles = {};
    for i = 1:length(hFigures)
        resFile = getappdata(hFigures(i), 'ResultsFile');
        if ~isempty(resFile)
            ResultsFiles{end+1} = resFile;
        else
            TfInfo = getappdata(hFigures(i), 'Timefreq');
            if ~isempty(TfInfo) && ~isempty(TfInfo.FileName)
                ResultsFiles{end+1} = TfInfo.FileName;
            end
        end
    end
    ResultsFiles = unique(ResultsFiles);
    ScoutsArg = 'SelectedScouts';
    clear hFigures resFile;
end
if isempty(ResultsFiles)
    % Try displaying scouts from tree selection
    tree_view_scouts();
    return
end

%% ===== GET SCOUTS LIST =====
isVolumeAtlas = 0;
% No scout
if isempty(ScoutsArg)
    return
% Use the scouts selected in the "Scouts" panel
elseif ischar(ScoutsArg) && strcmpi(ScoutsArg, 'SelectedScouts')
    % Get current atlas
    [sAtlas, iAtlas, sSurf] = panel_scout('GetAtlas');
    if isempty(sAtlas)
        return;
    end
    % Cannot display scouts from the "Source model" atlas (different regions)
    if ~isempty(strfind(sAtlas.Name, 'Source model'))
        bst_error('Cannot use the atlas "Source model" to calculate scouts time series.', 'Display scouts', 0);
        return;
    end
    % Volume scout: Get number of vertices of the atlas
    [isVolumeAtlas, nAtlasGrid] = panel_scout('ParseVolumeAtlas', sAtlas.Name);
    % Get selected scouts
    [sScouts, iScouts, sSurf] = panel_scout('GetSelectedScouts');
% Else: use directly the scout indices in argument
else
    iScouts = ScoutsArg;
    [sScouts, sSurf] = panel_scout('GetScouts', iScouts);
end
if isempty(sScouts)
    return
end
clear ScoutsArg;

%% ===== CHECK CORRESPONDANCE SCOUTS/SURFACES =====
iDroppedRes = [];
FileType = {};
% Check each results file
for i = 1:length(ResultsFiles)
    FileType{i} = file_gettype(ResultsFiles{i});
    % Load surface file from sources file
    switch (FileType{i})
        case {'link', 'results', 'presults'}
            ResMat = in_bst_results(ResultsFiles{i}, 0, 'SurfaceFile', 'HeadModelType', 'GridLoc');
            if isVolumeAtlas && (size(ResMat.GridLoc,1) ~= nAtlasGrid)
                disp(['BST> Error: The number of grid points in this atlas (' num2str(nAtlasGrid) ') does not match the loaded source file (' num2str(size(ResMat.GridLoc,1)) ').']);
                iDroppedRes(end+1) = i;
                continue;
            elseif ~isVolumeAtlas && strcmpi(ResMat.HeadModelType, 'volume')
                disp(['BST> Error: Cannot use a surface scout to extract value from a volume grid.']);
                iDroppedRes(end+1) = i;
                continue;
            end
        case {'timefreq', 'ptimefreq'}
            ResMat = in_bst_timefreq(ResultsFiles{i}, 0, 'SurfaceFile', 'HeadModelType', 'DataFile');
            if isempty(ResMat.SurfaceFile) && ~isempty(ResMat.DataFile)
                ResMat = in_bst_results(ResMat.DataFile, 0, 'SurfaceFile', 'HeadModelType');
            else
                ResMat.SurfaceFile   = sSurf.FileName;
                ResMat.HeadModelType = 'surface';
            end
    end
    % If surface is not the same as scouts' one, drop this results file
    if ~isempty(ResMat.SurfaceFile) && ~file_compare(ResMat.SurfaceFile, sSurf.FileName)    % || ismember(ResMat.HeadModelType, {'volume', 'mixed'})
        iDroppedRes = [iDroppedRes, i];
    end
end
if ~isempty(iDroppedRes)
    ResultsFiles(iDroppedRes) = [];
    FileType(iDroppedRes) = [];
end
% Check number of results files
if isempty(ResultsFiles)
    java_dialog('warning', 'No available source files.', 'Display scouts');
    return
end

%% ===== PREPARE DATA TO DISPLAY =====
% Initialize common descriptors (between all files)
SubjectFile    = '*';
StudyFile      = '*';

% Get display options
ScoutsOptions = panel_scout('GetScoutsOptions');
% if (length(iScouts) == 1)
%     ScoutsOptions.overlayScouts = 0;
% end
if (length(ResultsFiles) == 1)
    ScoutsOptions.overlayConditions = 0;
end
% Initialize data to display
scoutsActivity = cell(length(ResultsFiles), length(iScouts));
scoutsStd      = cell(length(ResultsFiles), length(iScouts));
scoutsLabels   = cell(length(ResultsFiles), length(iScouts));
scoutsColors   = cell(length(ResultsFiles), length(iScouts));
axesLabels     = cell(length(ResultsFiles), length(iScouts));
allComponents  = [];
TimeVector     = [];
issloreta      = 0;
DisplayUnits   = [];
Freqs          = [];
% Process each Results file
for iResFile = 1:length(ResultsFiles)
    % Is stat
    isStat = ismember(FileType{iResFile}, {'presults', 'ptimefreq'});
    isTimefreq = ismember(FileType{iResFile}, {'timefreq', 'ptimefreq'});
    isPsd = isTimefreq && (~isempty(strfind(ResultsFiles{iResFile}, '_psd')) || ~isempty(strfind(ResultsFiles{iResFile}, '_fft')));

    % ===== GET/CREATE RESULTS DATASET =====
    bst_progress('start', 'Display scouts time series', ['Loading results file: "' ResultsFiles{iResFile} '"...']);
    % Load results
    if ~isTimefreq
        [iDS, iResult] = bst_memory('LoadResultsFileFull', ResultsFiles{iResFile});
        if isempty(iDS)
            disp(['BST> Cannot load file : "', ResultsFiles{iResFile}, '"']);
            bst_progress('stop');
            return;
        end
        if ~isempty(strfind(lower(ResultsFiles{iResFile}), 'sloreta')) || ~isempty(strfind(lower(GlobalData.DataSet(iDS).Results(iResult).Comment), 'sloreta'))
            issloreta = 1;
        end
        fileUnits = GlobalData.DataSet(iDS).Results(iResult).DisplayUnits;
    else
        [iDS, iTimefreq, iResult] = bst_memory('LoadTimefreqFile', ResultsFiles{iResFile}, 1, 0);
        if isempty(iDS)
            disp(['BST> Cannot load file : "', ResultsFiles{iResFile}, '"']);
            bst_progress('stop');
            return;
        end
        fileUnits = GlobalData.DataSet(iDS).Timefreq(iTimefreq).DisplayUnits;
        % Check frequencies
        if isempty(Freqs)
            Freqs = GlobalData.DataSet(iDS).Timefreq(iTimefreq).Freqs;
        elseif length(Freqs) ~= length(GlobalData.DataSet(iDS).Timefreq(iTimefreq).Freqs)
            error('The frequencies must match across all files.');
        end
    end
    % Make sure that units are the same for all the files
    if ~isempty(fileUnits)
        if isempty(DisplayUnits)
            DisplayUnits = fileUnits;
        elseif ~strcmpi(fileUnits, DisplayUnits)
            DisplayUnits = 'Mixed units';
        end
    end

    % Get study
    [sStudy,   iStudy]   = bst_get('Study',   GlobalData.DataSet(iDS).StudyFile);
    if isempty(sStudy)
        error(['No study registered for file: ' strrep(ResultsFiles{iResFile},'\\','\')]);
    end
    % Get results subjectName/condition/#solInverse (FOR TITLE ONLY)
    [sSubject, iSubject] = bst_get('Subject', sStudy.BrainStormSubject);
    isInterSubject = (iStudy == -2);
    % Get identification of figure
    if ~isempty(SubjectFile)
        if (SubjectFile(1) == '*')
            SubjectFile = sStudy.BrainStormSubject;
        elseif ~file_compare(SubjectFile, sStudy.BrainStormSubject)
            SubjectFile = [];
        end
    end
    if ~isempty(StudyFile)
        if (StudyFile(1) == '*')
            StudyFile = GlobalData.DataSet(iDS).StudyFile;
        elseif ~file_compare(StudyFile, GlobalData.DataSet(iDS).StudyFile)
            StudyFile = [];
        end
    end
    DataFile = GlobalData.DataSet(iDS).DataFile;

    % ===== Prepare cell array containing time series to display =====
    TfFunction = [];
    for k = 1:length(sScouts)
        % ===== Get data to display =====
        % Get list of useful vertices in ImageGridAmp
        if ~isempty(iResult) && ~isempty(GlobalData.DataSet(iDS).Results(iResult).Atlas)
            % Atlas-based source file: use only the seed of the scout, and find it in the file atlas
            iVertices = panel_scout('GetScoutForVertex', GlobalData.DataSet(iDS).Results(iResult).Atlas, sScouts(k).Seed);
            % Error management
            if isempty(iVertices)
                disp('BST> Error: No data to display.');
                bst_progress('stop');
                return;
            end
        else
            % SORT and select unique vertices
            iVertices = sort(unique(sScouts(k).Vertices));
        end
        % Fix errors in color definition
        if isequal(size(sScouts(k).Color), [3,1])
            sScouts(k).Color = sScouts(k).Color';
        end
        % Get data (over current time window)
        if ~isTimefreq
            [DataToPlot, nComponents, DataStd] = bst_memory('GetResultsValues', iDS, iResult, iVertices, 'UserTimeWindow', 0, isVolumeAtlas);
        else
            % PSD: Get all the frequencies, display as a spectrum
            if isPsd
                iFreqs = [];
            % Other time-frequency files: Get only the selected time, display as time series
            else
                iFreqs = GlobalData.UserFrequencies.iCurrentFreq;
            end
            [DataToPlot, iTimeBands, iRow, nComponents] = bst_memory('GetTimefreqValues', iDS, iTimefreq, iVertices, iFreqs, 'UserTimeWindow', []);
            DataStd = [];
            % If the values are complex: try to get the measure from the Display tab
            if ~all(isreal(DataToPlot(:)))
                % If time-frequency function not found yet
                if isempty(TfFunction)
                    % Get current panel options
                    sOptions = panel_display('GetDisplayOptions');
                    TfFunction = sOptions.Function;
                    % Display warning
                    disp(['BST> Warning: We don''t know how to convert these complex values to real values. Trying function "' TfFunction '".']);
                end
                % Apply function
                [DataToPlot, isError] = process_tf_measure('Compute', DataToPlot, 'none', TfFunction);
                if isError
                    bst_error(['Invalid measure conversion: none => ' TfFunction], 'Display scouts', 0);
                    return;
                end
            end
        end

        % Get the time vector
        if isTimefreq && ~isempty(iTimeBands)
            TimeBands = process_tf_bands('GetBounds', GlobalData.DataSet(iDS).Timefreq(iTimefreq).TimeBands(iTimeBands,:));
            TimeVectorTmp = mean(TimeBands,2)';
        else
            [TimeVectorTmp, iTime] = bst_memory('GetTimeVector', iDS, [], 'UserTimeWindow');
            TimeVectorTmp = TimeVectorTmp(iTime);
        end
        if ~isempty(TimeVector) && (length(TimeVectorTmp) ~= length(TimeVector))
            bst_error('The time definition between the different files do not match.', 'Display scouts', 0);
            return;
        elseif isempty(TimeVector)
            TimeVector = TimeVectorTmp;
        end

        % Compute the scout values
        if ~isStat || strcmpi(sScouts(k).Function, 'All')
            ScoutFunction = sScouts(k).Function;
        else
            ScoutFunction = 'stat';
        end

        % Mixed headmodel results
        if (nComponents == 0)
            % Get atlas "Source model"
            GridAtlas = GlobalData.DataSet(iDS).Results(iResult).GridAtlas;
            % Get the vertex indices of the scout in GridLoc
            [iRows, iRegionScouts, iVertices] = bst_convert_indices(iVertices, nComponents, GridAtlas, ~isVolumeAtlas);
            % Do not accept scouts that span over multiple regions
            if isempty(iRegionScouts)
                bst_error(['Scout "' sScouts(k).Label '" is not included in the source model.' 10 'If you use this region as a volume, create a volume scout instead (menu Atlas > New atlas > Volume scouts).'], 'Display scouts', 0);
                return;
            elseif (length(iRegionScouts) > 1)
                bst_error(['Scout "' sScouts(k).Label '" spans over multiple regions of the "Source model" atlas.'], 'Display scouts', 0);
                return;
            end
            % Do not accept volume atlases with non-volume head models
            if ~isVolumeAtlas && strcmpi(GridAtlas.Scouts(iRegionScouts).Region(2), 'V')
                bst_error(['Scout "' sScouts(k).Label '" is a surface scout but region "' GridAtlas.Scouts(iRegionScouts).Label '" is a volume region.'], 'Display scouts', 0);
                return;
            elseif isVolumeAtlas && strcmpi(GridAtlas.Scouts(iRegionScouts).Region(2), 'S')
                bst_error(['Scout "' sScouts(k).Label '" is a volume scout but region "' GridAtlas.Scouts(iRegionScouts).Label '" is a surface region.'], 'Display scouts', 0);
                return;
            end
            % Set the scout computation properties based on the information in the "Source model" atlas
            if strcmpi(GridAtlas.Scouts(iRegionScouts).Region(3), 'C')
                nComponents = 1;
                VertNormals = GlobalData.DataSet(iDS).Results(iResult).GridOrient(iVertices,:);
            else
                nComponents = 3;
                VertNormals = [];
            end
        % Volume scouts
        elseif (nComponents == 3)
            VertNormals = [];
        else
            % Get vertex normals
            VertNormals = sSurf.VertNormals(iVertices,:);
        end
        % Keep track of how many components are available
        allComponents(end+1) = nComponents;

        % Only one component
        if (nComponents == 1)
            % Get meaningful tags in the results file name (without folders)
            TestResFile = file_resolve_link(ResultsFiles{iResFile});
            [tmp, TestTags] = bst_fileparts(TestResFile);
            % Do not flip sign for statistics, norms or NIRS source maps
            isFlipSign = ~isStat && ~isTimefreq && ...
                         strcmpi(GlobalData.DataSet(iDS).Results(iResult).DataType, 'results') && ...
                         isempty(strfind(TestTags, '_abs')) && ...
                         isempty(strfind(TestTags, '_norm')) && ...
                         isempty(strfind(TestTags, 'NIRS')) && ...
                         isempty(strfind(TestTags, 'Summed_sensitivities')) && ...
                         (isempty(GlobalData.DataSet(iDS).Channel) || isempty(GlobalData.DataSet(iDS).Results(iResult).GoodChannel) || ...
                          ~ismember('NIRS', {GlobalData.DataSet(iDS).Channel(GlobalData.DataSet(iDS).Results(iResult).GoodChannel).Type}));
            iTrace = k;
            scoutsActivity{iResFile,iTrace} = bst_scout_value(DataToPlot, ScoutFunction, VertNormals, nComponents, 'none', isFlipSign);
            if ~isempty(DataStd)
                scoutsStd{iResFile,iTrace} = bst_scout_value(DataStd, ScoutFunction, VertNormals, nComponents, 'none', 0);
            else
                scoutsStd{iResFile,iTrace} = [];
            end
            if ScoutsOptions.displayAbsolute
                scoutsActivity{iResFile,iTrace} = abs(scoutsActivity{iResFile,iTrace});
            end
        % More than one component & Absolute: Display the norm
        elseif ScoutsOptions.displayAbsolute
            iTrace = k;
            scoutsActivity{iResFile,iTrace} = bst_scout_value(DataToPlot, ScoutFunction, VertNormals, nComponents, 'norm');
            if ~isempty(DataStd)
                scoutsStd{iResFile,iTrace} = bst_scout_value(DataStd, ScoutFunction, VertNormals, nComponents, 'norm');
            else
                scoutsStd{iResFile,iTrace} = [];
            end
        % More than one component & Relative: Display the three components
        else
            iTrace = nComponents * (k-1) + (1:nComponents);
            tmp = bst_scout_value(DataToPlot, ScoutFunction, VertNormals, nComponents, 'none');
            for iComp = 1:nComponents
                scoutsActivity{iResFile,iTrace(iComp)} = tmp(iComp:nComponents:end,:);
            end
            % Std
            if ~isempty(DataStd)
                tmp = bst_scout_value(DataStd, ScoutFunction, VertNormals, nComponents, 'none');
                for iComp = 1:nComponents
                    scoutsStd{iResFile,iTrace(iComp)} = tmp(iComp:nComponents:end,:);
                end
            else
                scoutsStd(iResFile,iTrace) = repmat({[]}, 1, nComponents);
            end
        end

        % === AXES LABELS ===
        % === SUBJECT NAME ===
        % Format: SubjectName/CondName/(DataComment)/(Sol#iResult)/(ScoutName)
        strAxes = '';
        if ~isempty(sSubject) && (iSubject > 0)
            strAxes = sSubject.Name;
        end
        % === CONDITION ===
        % If inter-subject node
        if isInterSubject
            strAxes = [strAxes, 'Inter'];
        % Else: display conditions name
        else
            for i = 1:length(sStudy.Condition)
                strAxes = [strAxes, '/', sStudy.Condition{i}];
            end
        end
        % === DATA COMMENT ===
        % If more than one data file in this study : display current data comments
        if ~isempty(DataFile) && ((length(sStudy.Data) > 1) || isInterSubject)
            % Find DataFile in current study
            iData = find(file_compare({sStudy.Data.FileName}, DataFile));
            % Add Data comment
            if ~isempty(iData)
                strAxes = [strAxes, '/', sStudy.Data(iData).Comment];
            end
        end
        % === RESULTS COMMENT ===
        % If more than one results file in study : display indice
        if ~isempty(DataFile) && isStat && (length(sStudy.Stat) > 1)   % Stat
            % Get list of results files for current data file
            [tmp__, tmp__, iStatInStudy] = bst_get('StatForDataFile', DataFile, iStudy);
            % More than one results file for this data file
            if (length(iStatInStudy) > 1)
                strAxes = [strAxes, '/', GlobalData.DataSet(iDS).Results(iResult).Comment];
            end
        elseif ~isempty(DataFile) && isTimefreq && (length(sStudy.Timefreq) > 1)    % Time-frequency
            % Get list of results files for current data file
            [tmp__, tmp__, iTfInStudy] = bst_get('TimefreqForFile', DataFile, iStudy);
            % More than one results file for this data file
            if (length(iTfInStudy) > 1)
                strAxes = [strAxes, '/', GlobalData.DataSet(iDS).Timefreq(iTimefreq).Comment];
            end
        elseif ~isempty(DataFile) && ~isTimefreq && (length(sStudy.Result) > 1)    % Source maps
            % Get list of results files for current data file
            [tmp__, tmp__, iResInStudy] = bst_get('ResultsForDataFile', DataFile, iStudy);
            % More than one results file for this data file
            if (length(iResInStudy) > 1)
                strAxes = [strAxes, '/', GlobalData.DataSet(iDS).Results(iResult).Comment];
            end
        % Inter-subject: always display whole
        elseif (isempty(DataFile) || isInterSubject)
           strAxes = [strAxes, '/', GlobalData.DataSet(iDS).Results(iResult).Comment];
        end
        [axesLabels{iResFile,iTrace}] = deal(strAxes);

        % === SCOUTS LABELS/COLORS ===
        if (length(iTrace) == 1)
            scoutsLabels{iResFile,iTrace} = sScouts(k).Label;
            [scoutsColors{iResFile,iTrace}] = deal(sScouts(k).Color);
        else
            for iComp = 1:length(iTrace)
                scoutsLabels{iResFile,iTrace(iComp)} = sprintf('%s%d', sScouts(k).Label, iComp);
                scoutsColors{iResFile,iTrace(iComp)} = sScouts(k).Color .* (1 - .25 * (iComp-1));
            end
        end
    end
end

% ===== DISPLAY STATIC VALUES =====
% Get the number of time samples for all the scouts
nbTimes = cellfun(@(c)size(c,2), scoutsActivity(:));
% If both real scouts and static values
if ((max(nbTimes) > 2) && any(nbTimes == 2))
    iCellToExtend = find(nbTimes == 2);
    for i = 1:length(iCellToExtend)
        scoutsActivity{iCellToExtend(i)} = repmat(scoutsActivity{iCellToExtend(i)}(:,1), [1,max(nbTimes)]);
        if ~isempty(scoutsStd)
            scoutsStd{iCellToExtend(i)} = repmat(scoutsStd{iCellToExtend(i)}(:,1), [1,max(nbTimes)]);
        end
    end
end

%% ===== LEGENDS =====
% Get common beginning part in all the axesLabels
[ strAxesCommon, axesLabels ] = str_common_path( axesLabels );
% Only one time series: no overlay scout
if (size(scoutsActivity,2) == 1)
    ScoutsOptions.overlayScouts = 0;
end
% Cannot mix number of components
isMixedComponents = ~ScoutsOptions.displayAbsolute && (length(allComponents) > 1) && any(allComponents ~= allComponents(1));
if isMixedComponents
    bst_error(['Cannot display scouts with mixed numbers of components at the same time ' 10 'when the option "Values: Relative" is selected.'], 'Display scouts', 0);
    return;
end
% If at least one of the scouts functions is "All", ignore the overlay checkboxes
if ~isempty(sScouts) && any(strcmpi({sScouts.Function}, 'All'))
    ScoutsOptions.overlayScouts     = 0;
    ScoutsOptions.overlayConditions = 0;
end

% === NO OVERLAY ===
% Display all timeseries (ResultsFiles/Scouts) in different axes on the same figure
if (~ScoutsOptions.overlayScouts && ~ScoutsOptions.overlayConditions)
    % One graph for each line => Ngraph, Nscouts*Ncond
    % Scouts activity = cell-array {1, Ngraph} of doubles [1,Nt]
    scoutsActivity = scoutsActivity(:)';
    % Axes labels (subj/cond) = cell-array of strings {1, Ngraph}
    axesLabels = axesLabels(:)';
    % Scouts labels = cell-array of strings {1, Ngraph}
    scoutsLabels = scoutsLabels(:)';
    scoutsColors = scoutsColors(:)';
    % Eliminate empty entries
    iEmpty = find(cellfun(@isempty, scoutsActivity));
    if ~isempty(iEmpty)
        scoutsActivity(iEmpty) = [];
        scoutsLabels(iEmpty) = [];
        axesLabels(iEmpty) = [];
        if ~isempty(scoutsColors)
            scoutsColors(iEmpty) = [];
        end
    end
    % Std
    if ~isempty(scoutsStd)
        scoutsStd = scoutsStd(:)';
        if ~isempty(iEmpty)
            scoutsStd(iEmpty) = [];
        end
    end

% === OVERLAY SCOUTS AND CONDITIONS ===
% Only one graph with Nlines = Nscouts*Ncond
elseif (ScoutsOptions.overlayScouts && ScoutsOptions.overlayConditions)
    % Linearize entries
    scoutsActivity = scoutsActivity(:);
    scoutsLabels = scoutsLabels(:);
    axesLabels = axesLabels(:);
    % Eliminate empty entries
    iEmpty = find(cellfun(@isempty, scoutsActivity));
    if ~isempty(iEmpty)
        scoutsActivity(iEmpty) = [];
        scoutsLabels(iEmpty) = [];
        axesLabels(iEmpty) = [];
    end
    % Scouts activity = double [Nlines, Nt]
    scoutsActivity = cat(1, scoutsActivity{:});
    % Scouts labels = cell-array {Nlines}
    % => "scoutName @ subject/condition"
    for iTrace = 1:size(scoutsActivity,1)
        if ~isempty(axesLabels{iTrace})
            scoutsLabels{iTrace} = [scoutsLabels{iTrace} ' @ ' axesLabels{iTrace}];
        end
    end
    scoutsColors = [];
    % Only one graph => legend is common scouts parts
    if ~isempty(strAxesCommon)
        axesLabels = {strAxesCommon};
    else
        axesLabels = {'Mixed subjects'};
    end
    % Std
    if ~isempty(scoutsStd)
        scoutsStd = scoutsStd(:);
        if ~isempty(iEmpty)
            scoutsStd(iEmpty) = [];
        end
        scoutsStd = cat(1, scoutsStd{:});
    end

% === OVERLAY SCOUTS ONLY ===
elseif ScoutsOptions.overlayScouts
    % One graph per condition => Ngraph = Ncond
    % Scouts activity = cell-array {1, Ncond} of doubles [Nscout, Nt]
    for i = 1:size(scoutsActivity,1)
        scoutsActivityTmp(i) = {cat(1, scoutsActivity{i,:})};
    end
    scoutsActivity = scoutsActivityTmp;
    % Axes labels (subj/cond) = cell-array of strings {1, Ncond}
    axesLabels   = axesLabels(:,1)';
    % Scouts labels = cell-array of strings {Nscout, Ncond}
    scoutsLabels = scoutsLabels';
    scoutsColors = scoutsColors';
    % Std
    if ~isempty(scoutsStd)
        for i = 1:size(scoutsStd,1)
            scoutsStdTmp(i) = {cat(1, scoutsStd{i,:})};
        end
        scoutsStd = scoutsStdTmp;
    end

% === OVERLAY CONDITIONS ONLY ===
elseif ScoutsOptions.overlayConditions
    % One graph per scout => Ngraph = Nscout
    % Scouts activity = cell-array {1, Nscout} of doubles [Ncond, Nt]
    for j = 1:size(scoutsActivity,2)
        scoutsActivityTmp(j) = {cat(1, scoutsActivity{:,j})};
    end
    scoutsActivity = scoutsActivityTmp;
    % Axes labels (scout names @ common_subject/cond part) = cell-array of strings {1, Nscout}
    tmpAxesLabels = axesLabels;
    axesLabels = scoutsLabels(1,:);
    % Lines labels (subject/condition) = cell-array of strings {Ncond, Nscout}
    scoutsLabels = tmpAxesLabels;
    scoutsColors = [];
    % Std
    if ~isempty(scoutsStd)
        for j = 1:size(scoutsStd,2)
            scoutsStdTmp(j) = {cat(1, scoutsStd{:,j})};
        end
        scoutsStd = scoutsStdTmp;
    end
end
% Close progress bar
bst_progress('stop');

%% ===== CALL DISPLAY FUNCTION ====
% Get type of the first file
switch (file_gettype(ResultsFiles{1}))
    case {'results', 'link'}
        if issloreta
            Modality = 'sloreta';
        else
            Modality = 'results';
        end
    case 'timefreq'
        Modality = 'timefreq';
    case 'presults'
        Modality = 'stat';
    otherwise
        Modality = 'none';
end
% Plot spectrum
if isPsd
    if (length(scoutsActivity) > 1)
        error(['Scouts on PSD plots can only show one graph at a time.' 10 ...
               'Open only one file at a time and enable the "overlay" options.']);
    end
    hFig = view_spectrum_matrix(ResultsFiles, 'Spectrum', scoutsActivity{1}, scoutsLabels, TimeVector, Freqs, 'Scout', DisplayUnits);
% Plot time series
else
    hFig = view_timeseries_matrix(ResultsFiles, scoutsActivity, TimeVector, Modality, axesLabels, scoutsLabels, scoutsColors, hFig, scoutsStd, DisplayUnits);
end
% Associate the file with one specific result file
setappdata(hFig, 'ResultsFiles', ResultsFiles);
if (length(ResultsFiles) == 1)
    setappdata(hFig, 'ResultsFile', ResultsFiles);
else
    setappdata(hFig, 'ResultsFile', []);
end
% Update figure name
bst_figures('UpdateFigureName', hFig);
% Set the time label visible
figure_timeseries('SetTimeVisible', hFig, 1);

author: brainstorm-tools


text:
function x = IWT2_PO(wc,L,qmf)
% IWT2_PO -- Inverse 2-d MRA wavelet transform (periodized, orthogonal)
%  Usage
%    x = IWT2_PO(wc,L,qmf)
%  Inputs
%    wc    2-d wavelet transform [n by n array, n dyadic]
%    L     coarse level
%    qmf   quadrature mirror filter
%  Outputs
%    x     2-d signal reconstructed from wc
%
%  Description
%    If wc is the result of a forward 2d wavelet transform, with
%    wc = FWT2_PO(x,L,qmf), then x = IWT2_PO(wc,L,qmf) reconstructs x
%    exactly if qmf is a nice qmf, e.g. one made by MakeONFilter.
%
%  See Also
%    FWT2_PO, MakeONFilter
%
	[n,J] = quadlength(wc);
	x = wc;
	nc = 2^(L+1);
	for jscal=L:J-1, % from coarse to fine
		top = (nc/2+1):nc; bot = 1:(nc/2); all = 1:nc;
		for iy=1:nc,
			x(all,iy) =  UpDyadLo(x(bot,iy)',qmf)'  ...
					   + UpDyadHi(x(top,iy)',qmf)';
		end
		for ix=1:nc,
			x(ix,all) = UpDyadLo(x(ix,bot),qmf)  ...
					  + UpDyadHi(x(ix,top),qmf);
		end
		nc = 2*nc;
	end

%
% Copyright (c) 1993. David L. Donoho
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: thomaskuestner