text:
function y = ft_preproc_slidingrange(dat, width, normalize, varargin)

% FT_PREPROC_SLIDINGRANGE computes the range of the data in a sliding time
% window of the width specified.
%
% Use as
%   [dat] = ft_preproc_slidingrange(dat, width, normalize)
% where
%   dat        data matrix (Nchans x Ntime)
%   width      width of the smoothing kernel, this should be an odd number since the window needs to be centered on an individual sample
%   normalize  boolean, whether to normalize the range of the data with the square root of the window size (default = false)
%
% If the data contains NaNs, these are ignored for the computation, but retained in
% the output.
%
% See also PREPROC

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

if nargin>3
  % for backward compatibility, this function was first implemented as ft_preproc_slidingrange(dat, width, ...)
  % with 'normalize' as an optional key-value pair, but all other PREPROC functions take fixed input arguments
  normalize = varargin{1};
end

if nargin<3 || isempty(normalize)
  normalize = false;
end

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

if mod(width+1, 2)
  ft_error('width should be an odd number');
end

% compute half width
h = (width-1)/2;

n = size(dat,2);
minval = zeros(size(dat));
maxval = zeros(size(dat));

for i=1:n
  begsample = i-h;
  endsample = i+h;
  if begsample<1
    begsample = 1;
  end
  if endsample>n
    endsample=n;
  end
  minval(:,i) = min(dat(:,begsample:endsample),[],2);
  maxval(:,i) = max(dat(:,begsample:endsample),[],2);
end

y = maxval - minval;

if istrue(normalize)
  y = y ./ sqrt(width);
end


text:
function [Estimation_X0,FirstLandmarks] = FEKF_update(Estimation_X0, CameraMeasurementThis, Sigma_OB , FirstLandmarks)
NumberOfLandmarksObInThisStep = size(CameraMeasurementThis,1)/3;

% initialise the IndexOfFeature if possible
if size(Estimation_X0.landmarks,2) > 0
    IndexOfFeature = Estimation_X0.landmarks(4,:)';
else
    IndexOfFeature = [];
end

IndexObservedAlreadyThis = [];
IndexObservedNew = [];
for i = 1:NumberOfLandmarksObInThisStep
    % check whether the feature is observed before or not
    M = find( IndexOfFeature== CameraMeasurementThis(3*i,2) );
    if isempty(M)
        IndexObservedNew = [IndexObservedNew;CameraMeasurementThis(3*i,2) ];
    else
        IndexObservedAlreadyThis = [IndexObservedAlreadyThis;CameraMeasurementThis(3*i,2)];
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IndexObservedNew=[78; 97; 18] indicates that the
% robot firstly observes landmarks 78 97 18 in this step
% IndexObservedAlreadyThis=[19; 20; 53] indicates
% that the robot observes again landmarks 19 20 53 in this step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

orientation = Estimation_X0.orientation;
position    = Estimation_X0.position;
cov         = Estimation_X0.cov;

NumberOfFeature = size( IndexOfFeature,1);
NumberOfOldFeatureInThisStep = size(IndexObservedAlreadyThis,1);
NumberOfNewFeatureInThisStep = size(IndexObservedNew,1);

% update state vector and covariance by considering
% new feature into state and covariance
if ~isempty(IndexObservedNew)

    orientation = Estimation_X0.orientation;
    position    = Estimation_X0.position;

   % copy previous covariance
    temp    = repmat({eye(3)}, NumberOfNewFeatureInThisStep, 1 );
    tempKK  = blkdiag(temp{:});
    Sigma   = blkdiag(Estimation_X0.cov,tempKK);
    KK      = eye(6+3*(NumberOfFeature+NumberOfNewFeatureInThisStep));

 %   add new features
    for i = 1:NumberOfNewFeatureInThisStep
        indNewf = IndexObservedNew(i);
        Estimation_X0.landmarks(4,NumberOfFeature+i) = indNewf;
        m2 = find( CameraMeasurementThis(:,2) == indNewf );
        nf = CameraMeasurementThis( m2, 1 );

        Estimation_X0.landmarks(1:3,NumberOfFeature+i) = orientation*nf+position;
        Newlandmark=[ orientation*nf+position;  indNewf  ];

        FirstLandmarks=[FirstLandmarks  Newlandmark];

%         KK( 6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i,1:6 ) = [-skew(Estimation_X0.orientation*(nf)) eye(3)];
%         KK (6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i, 6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i )=Estimation_X0.orientation;

         KK( 6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i,1:6 ) = [-skew(orientation*(nf)) eye(3)];
        KK (6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i, 6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i )=orientation;
                tempKK(3*i-2:3*i,3*i-2:3*i)=diag(nf.^2)*Sigma_OB^2;
    end
    Sigma   = blkdiag(Estimation_X0.cov,tempKK);
    Estimation_X0.cov = KK*Sigma*KK';

    NumberOfFeature=NumberOfFeature+NumberOfNewFeatureInThisStep;
end

orientation = Estimation_X0.orientation;
position    = Estimation_X0.position;
cov         = Estimation_X0.cov;

% update state and covariance
if ~isempty(IndexObservedAlreadyThis)
    Z = zeros( NumberOfOldFeatureInThisStep*3 , 1);
    Y = zeros( NumberOfOldFeatureInThisStep*3 , 1);
    H = zeros(3*NumberOfOldFeatureInThisStep, 6+3*NumberOfFeature);

    temp = repmat({eye(3)}, NumberOfOldFeatureInThisStep,1 );
    R = blkdiag(temp{:});

    % update old features
    for i = 1:NumberOfOldFeatureInThisStep
        ind = find(IndexOfFeature == IndexObservedAlreadyThis(i));

        fi  = Estimation_X0.landmarks(1:3,ind);
        Y(3*i-2:3*i,1) = ObservationModel( orientation, position, fi );

        ind2 = find(CameraMeasurementThis(:,2) == IndexObservedAlreadyThis(i));
        Z(3*i-2:3*i,1 ) = CameraMeasurementThis(ind2,1);

        FirstLandmark=FirstLandmarks(1:3,ind);

        H(3*i-2:3*i, 1:6) = [-orientation'*skew(FirstLandmark-position) orientation'];
        H(3*i-2:3*i, 6+3*ind-2:6+3*ind) = -orientation';
        R(3*i-2:3*i,3*i-2:3*i) = diag(CameraMeasurementThis(ind2,1).^2)*Sigma_OB^2;
    end

    % question @RomaTeng, different computaton scheme
    z = Z-Y;
    S = H*cov*H'+R;
    K = cov*H'*inv(S);
    s = K*z;

    Estimation_X0 = SpecialAdd(Estimation_X0,-s);
    cov = ( eye(6+3*NumberOfFeature) -K*H )*cov;
    Estimation_X0.cov = cov;

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function vw = nifti2ROI(vw, mappth)
%   Convert a nifti label file into a mrVista format ROI
%
% vw = nifti2ROI([vw], [mappth])
%
% vw:       mrVista view structure (must be gray view)
%               [default = current Gray view]
% mappth:   path to nifti file with map to convert;
%               [default = dialog]
%
% Notes: the nifti format map is a 3D matrix. The mrVista parameter map is
% a vector, indexed to the gray coords. The main steps of the routine are:
%   1) Convert the nifti format 3D map to a vector and 3xN coordinate
%           matrix
%   2) Transform the coordinates to the space of the t1-weighted anatomy
%       (using header information in the two files)
%   3) Transform the coordinates to mrAnat conventions
%
% Example:  Convert a class file into ROIs
% vw = nifti2ROI(vw, viewGet(vw, 'class file', 'right'))
%
%   JW, 4/27/2015
%

mrGlobals;

% Variable check
if notDefined('vw'), vw = getSelectedGray; end

if ~exist('mappth', 'var') || ~exist(mappth, 'file')
    mappth = getPathStrDialog(dataDir(vw),'Choose nifti parameter map','*.nii.gz');
end

% read a nifti file with map to be converted to ROIs
ni = niftiRead(mappth);

% apply our canonical transform to ensure orientation is matched to t1;
ni   = niftiApplyCannonicalXform(ni);

% This is the map data
data = niftiGet(ni, 'data');
data = nifti2mrVistaAnat(data);
data = round(data);
% ensure data are integers
assert(isequal(data, round(data)));

% define ROIs
labels = setdiff(unique(data(:)),0);

fprintf('[%s]: Creating %d ROIs from file %s\n', mfilename, length(labels), mappth);

for ii = 1:length(labels)
   [x, y, z] = ind2sub(size(data), find(data == labels(ii)));
   coords = [x y z]';
   comments = sprintf('ROI defined by label %d in map %s', labels(ii), mappth);
   name = sprintf('ROI_%03d', labels(ii));
   fprintf('.');drawnow();
   vw = newROI(vw,name,1,'k',coords, comments);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = r8_lgic ( a, x, alx )

%*****************************************************************************80
%
%% R8_LGIC evaluates the log complementary incomplete gamma function for large X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, the parameter.
%
%    Input, real X, the argument.
%
%    Input, real ALX, the logarithm of X.
%
%    Output, real VALUE, the log complementary incomplete
%    gamma function.
%
  persistent eps

  if ( isempty ( eps ) )
    eps = 0.5 * r8_mach ( 3 );
  end

  xpa = x + 1.0 - a;
  xma = x - 1.0 - a;

  r = 0.0;
  p = 1.0;
  s = p;
  for k = 1 : 300
    fk = k;
    t = fk * ( a - fk ) * ( 1.0 + r );
    r = - t / ( ( xma + 2.0 * fk ) * ( xpa + 2.0 * fk ) + t );
    p = r * p;
    s = s + p;
    if ( abs ( p ) < eps * s )
      value = a * alx - x + log ( s / xpa );
      return
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_LGIC - Fatal error!\n' );
  fprintf ( 1, '  No convergence in 300 iterations.\n' );

  error ( 'R8_LGIC - Fatal error!' )
end

author: johannesgerer

repo: jburkardt-m


text:
function [pos, neg, impos] = pascal_data(cls, year)
% Get training data from the PASCAL dataset.
%   [pos, neg, impos] = pascal_data(cls, year)
%
% Return values
%   pos     Each positive example on its own
%   neg     Each negative image on its own
%   impos   Each positive image with a list of foreground boxes
%
% Arguments
%   cls     Object class to get examples for
%   year    PASCAL dataset year

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf       = voc_config('pascal.year', year);
dataset_fg = conf.training.train_set_fg;
dataset_bg = conf.training.train_set_bg;
cachedir   = conf.paths.model_dir;
VOCopts    = conf.pascal.VOCopts;

try
  load([cachedir cls '_' dataset_fg '_' year]);
catch
  % Positive examples from the foreground dataset
  ids      = textread(sprintf(VOCopts.imgsetpath, dataset_fg), '%s');
  pos      = [];
  impos    = [];
  numpos   = 0;
  numimpos = 0;
  dataid   = 0;
  for i = 1:length(ids);
    tic_toc_print('%s: parsing positives (%s %s): %d/%d\n', ...
                  cls, dataset_fg, year, i, length(ids));
    % Parse record and exclude difficult examples
    rec           = PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    clsinds       = strmatch(cls, {rec.objects(:).class}, 'exact');
    diff          = [rec.objects(clsinds).difficult];
    clsinds(diff) = [];
    count         = length(clsinds(:));
    % Skip if there are no objects in this image
    if count == 0
      continue;
    end

    % Create one entry per bounding box in the pos array
    for j = clsinds(:)'
      numpos = numpos + 1;
      dataid = dataid + 1;
      bbox   = rec.objects(j).bbox;

      pos(numpos).im      = [VOCopts.datadir rec.imgname];
      pos(numpos).x1      = bbox(1);
      pos(numpos).y1      = bbox(2);
      pos(numpos).x2      = bbox(3);
      pos(numpos).y2      = bbox(4);
      pos(numpos).boxes   = bbox;
      pos(numpos).flip    = false;
      pos(numpos).trunc   = rec.objects(j).truncated;
      pos(numpos).dataids = dataid;
      pos(numpos).sizes   = (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1);

      % Create flipped example
      numpos  = numpos + 1;
      dataid  = dataid + 1;
      oldx1   = bbox(1);
      oldx2   = bbox(3);
      bbox(1) = rec.imgsize(1) - oldx2 + 1;
      bbox(3) = rec.imgsize(1) - oldx1 + 1;

      pos(numpos).im      = [VOCopts.datadir rec.imgname];
      pos(numpos).x1      = bbox(1);
      pos(numpos).y1      = bbox(2);
      pos(numpos).x2      = bbox(3);
      pos(numpos).y2      = bbox(4);
      pos(numpos).boxes   = bbox;
      pos(numpos).flip    = true;
      pos(numpos).trunc   = rec.objects(j).truncated;
      pos(numpos).dataids = dataid;
      pos(numpos).sizes   = (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1);
    end

    % Create one entry per foreground image in the impos array
    numimpos                = numimpos + 1;
    impos(numimpos).im      = [VOCopts.datadir rec.imgname];
    impos(numimpos).boxes   = zeros(count, 4);
    impos(numimpos).dataids = zeros(count, 1);
    impos(numimpos).sizes   = zeros(count, 1);
    impos(numimpos).flip    = false;

    for j = 1:count
      dataid = dataid + 1;
      bbox   = rec.objects(clsinds(j)).bbox;

      impos(numimpos).boxes(j,:) = bbox;
      impos(numimpos).dataids(j) = dataid;
      impos(numimpos).sizes(j)   = (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1);
    end

    % Create flipped example
    numimpos                = numimpos + 1;
    impos(numimpos).im      = [VOCopts.datadir rec.imgname];
    impos(numimpos).boxes   = zeros(count, 4);
    impos(numimpos).dataids = zeros(count, 1);
    impos(numimpos).sizes   = zeros(count, 1);
    impos(numimpos).flip    = true;
    unflipped_boxes         = impos(numimpos-1).boxes;

    for j = 1:count
      dataid  = dataid + 1;
      bbox    = unflipped_boxes(j,:);
      oldx1   = bbox(1);
      oldx2   = bbox(3);
      bbox(1) = rec.imgsize(1) - oldx2 + 1;
      bbox(3) = rec.imgsize(1) - oldx1 + 1;

      impos(numimpos).boxes(j,:) = bbox;
      impos(numimpos).dataids(j) = dataid;
      impos(numimpos).sizes(j)   = (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1);
    end
  end

  % Negative examples from the background dataset
  ids    = textread(sprintf(VOCopts.imgsetpath, dataset_bg), '%s');
  neg    = [];
  numneg = 0;
  for i = 1:length(ids);
    tic_toc_print('%s: parsing negatives (%s %s): %d/%d\n', ...
                  cls, dataset_bg, year, i, length(ids));
    rec = PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    clsinds = strmatch(cls, {rec.objects(:).class}, 'exact');
    if length(clsinds) == 0
      dataid             = dataid + 1;
      numneg             = numneg+1;
      neg(numneg).im     = [VOCopts.datadir rec.imgname];
      neg(numneg).flip   = false;
      neg(numneg).dataid = dataid;
    end
  end

  save([cachedir cls '_' dataset_fg '_' year], 'pos', 'neg', 'impos');
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
%kvaddm 'Add Image1 to Image2 Pixel by Pixel'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vaddm.pane file
%
% Parameters:
% InputFile: i1 'Input Image #1', required: 'first input image'
% OutputFile: o 'Output Image ', required: 'resulting output image'
% InputFile: i2 'Input Image #2', optional: 'second input image'
% Integer: c 'Constant Image', default: 0: 'add constant'
%
% Example: o = kvaddm({i1, i2}, {'i1','';'o','';'i2','';'c',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vaddm - Add Image1 to Image2 Pixel by Pixel
%
%  DESCRIPTION
% .I vaddm
% Adds the first input image to the second input image (result = image1 + image2), with an optional operation gating image. The second input image can be entered as a constant, optionally.
% If the sum of two pixel values is more than the maximum value of the image data type, the resultant pixel value will be set to this maximum value.
%
% The two input images must be of the same size and have the same
% number of data bands.  The input
% images may be of the same or different data types.  If the input
% images are of different data types, then all input images will be
% upcast to the highest input image data type.  The output
% image data type will be the same as that of the highest data type
% of the input images.
%
%
%
%  EXAMPLES
% vaddm -i1 input.image1 -i2 input.image2 -o output.image
%
% This subtracts input.image2 from input.image1 stores the results in output.image.
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vaddm
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvaddm(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvaddm(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input';'c', 0};
maxval={0,0,1,255};
minval={0,0,1,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vaddm"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [net, stats] = cnn_toy_data(varargin)
% CNN_TOY_DATA
% Minimal demonstration of MatConNet training of a CNN on toy data.
%
% It also serves as a short tutorial on creating and using a custom imdb
% (image database).
%
% The task is to distinguish between images of triangles, squares and
% circles.

% Copyright (C) 2017 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

run([fileparts(mfilename('fullpath')) '/../../matlab/vl_setupnn.m']) ;

% Parameter defaults. You can add any custom parameters here (e.g.
% opts.alpha = 1), and change them when calling: cnn_toy_data('alpha', 2).
opts.train.batchSize = 200 ;
opts.train.numEpochs = 10 ;
opts.train.continue = true ;
opts.train.gpus = [] ;
opts.train.learningRate = 0.01 ;
opts.train.expDir = [vl_rootnn '/data/toy'] ;
opts.dataDir = [vl_rootnn '/data/toy-dataset'] ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.imdbPath = [opts.train.expDir '/imdb.mat'] ;
opts = vl_argparse(opts, varargin) ;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

% Generate images if they don't exist (this would be skipped for real data)
if ~exist(opts.dataDir, 'dir')
  mkdir(opts.dataDir) ;
  cnn_toy_data_generator(opts.dataDir) ;
end

% Create image database (imdb struct). It can be cached to a file for speed
if exist(opts.imdbPath, 'file')
  disp('Reloading image database...')
  imdb = load(opts.imdbPath) ;
else
  disp('Creating image database...')
  imdb = getImdb(opts.dataDir) ;
  mkdir(fileparts(opts.imdbPath)) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% Create network (see HELP VL_SIMPLENN)
f = 1/100 ;
net.layers = {} ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,1,5, 'single'), zeros(1, 5, 'single')}}) ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2) ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,5,10, 'single'),zeros(1,10,'single')}}) ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2) ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,10,3, 'single'),  zeros(1,3,'single')}}) ;
net.layers{end+1} = struct('type', 'softmaxloss') ;

% Fill in any values we didn't specify explicitly
net = vl_simplenn_tidy(net) ;

% --------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------

use_gpu = ~isempty(opts.train.gpus) ;

% Start training
[net, stats] = cnn_train(net, imdb, @(imdb, batch) getBatch(imdb, batch, use_gpu), ...
  'train', find(imdb.set == 1), 'val', find(imdb.set == 2), opts.train) ;

% Visualize the learned filters
figure(3) ; vl_tshow(net.layers{1}.weights{1}) ; title('Conv1 filters') ;
figure(4) ; vl_tshow(net.layers{3}.weights{1}) ; title('Conv2 filters') ;
figure(5) ; vl_tshow(net.layers{5}.weights{1}) ; title('Conv3 filters') ;

% --------------------------------------------------------------------
function [images, labels] = getBatch(imdb, batch, use_gpu)
% --------------------------------------------------------------------
% This is where we return a given set of images (and their labels) from
% our imdb structure.
% If the dataset was too large to fit in memory, getBatch could load images
% from disk instead (with indexes given in 'batch').

images = imdb.images(:,:,:,batch) ;
labels = imdb.labels(batch) ;

if use_gpu
  images = gpuArray(images) ;
end

% --------------------------------------------------------------------
function imdb = getImdb(dataDir)
% --------------------------------------------------------------------
% Initialize the imdb structure (image database).
% Note the fields are arbitrary: only your getBatch needs to understand it.
% The field imdb.set is used to distinguish between the training and
% validation sets, and is only used in the above call to cnn_train.

% The sets, and number of samples per label in each set
sets = {'train', 'val'} ;
numSamples = [1500, 150] ;

% Preallocate memory
totalSamples = 4950 ;  % 3 * 1500 + 3 * 150
images = zeros(32, 32, 1, totalSamples, 'single') ;
labels = zeros(totalSamples, 1) ;
set = ones(totalSamples, 1) ;

% Read all samples
sample = 1 ;
for s = 1:2  % Iterate sets
  for label = 1:3  % Iterate labels
    for i = 1:numSamples(s)  % Iterate samples
      % Read image
      im = imread(sprintf('%s/%s/%i/%04i.png', dataDir, sets{s}, label, i)) ;

      % Store it, along with label and train/val set information
      images(:,:,:,sample) = single(im) ;
      labels(sample) = label ;
      set(sample) = s ;
      sample = sample + 1 ;
    end
  end
end

% Show some random example images
figure(2) ;
montage(images(:,:,:,randperm(totalSamples, 100))) ;
title('Example images') ;

% Remove mean over whole dataset
images = bsxfun(@minus, images, mean(images, 4)) ;

% Store results in the imdb struct
imdb.images = images ;
imdb.labels = labels ;
imdb.set = set ;

repo: LapSRN


text:
function pbsNonSharedSimpleSubmitFcn(scheduler, job, props, ...
                                        clusterHost, remoteDataLocation)
%PBSNONSHAREDSIMPLESUBMITFCN Submits a MATLAB job to a PBS scheduler in the
% absence of a shared file system between the MATLAB client and the
% PBS scheduler.
%
% See also pbsNonSharedSimpleDecodeFcn.

% Copyright 2006-2008 The MathWorks, Inc.

if ~ischar(clusterHost)
    error('distcomp:genericscheduler:SubmitFcnError', ...
        'Hostname must be a string');
end
if ~ischar(remoteDataLocation)
    error('distcomp:genericscheduler:SubmitFcnError', ...
        'Remote Data Location must be a string');
end
scheduler.UserData = { clusterHost ; remoteDataLocation };
localDataLocation = scheduler.DataLocation;

% Set the name of the decode function which will be executed by
% the worker. The decode function must be on the path of the MATLAB
% worker when it starts up. This is typically done by placing the decode
% function in MATLABROOT/toolbox/local on the cluster nodes, or by
% prefixing commandToRun (created below) with a command to cd to the
% directory where the decode function's M file exists.
decodeFcn = 'pbsNonSharedSimpleDecodeFcn';

% Read the number of tasks which are to be created. This property
% cannot be changed.
numberOfTasks = props.NumberOfTasks;

% A unique file and directory name for the job. This is used to create
% files and a directory under scheduler.DataLocation
jobLocation = props.JobLocation;

% A cell array of unique file names for tasks. These are used to create
% files under jobLocation
taskLocations = props.TaskLocations;

% Since PBS jobs will be submitted from a UNIX host on the cluster,
% a single quoted string will protect the MATLAB command.
quote = '''';

% The MATLAB command to be run on a cluster node to execute a task.
commandToRun = [ quote props.MatlabExecutable quote ' ' ...
                 props.MatlabArguments ];

% Copy the matlab_metadata.mat file to the remote host.
localMetaDataFile = [ localDataLocation '/matlab_metadata.mat' ];
copyDataToCluster(localMetaDataFile, remoteDataLocation, clusterHost);

% Copy the local job directory to the remote host.
localJobDirectory = [ localDataLocation '/' jobLocation ];
copyDataToCluster(localJobDirectory, remoteDataLocation, clusterHost);

% Copy the local job files to the remote host.
localJobFiles = [ localDataLocation '/' jobLocation '.*' ];
copyDataToCluster(localJobFiles, remoteDataLocation, clusterHost);

% Submit tasks which the scheduler will execute by starting MATLAB workers.
for i = 1:numberOfTasks
    taskLocation = taskLocations{i};
    remoteJobDir = [ remoteDataLocation '/' jobLocation ];

    % Create a script to submit a PBS job.
    localScript = createPBSSubmitScript(commandToRun, decodeFcn, ...
        props.StorageConstructor, remoteDataLocation, ...
        jobLocation, taskLocation, quote);
    [Path, scriptName] = fileparts(localScript);

    % Copy the script to the remote host.
    copyDataToCluster(localScript, remoteJobDir, clusterHost);

    % Create the command to run on the remote host.
    remoteScriptLocation = [ remoteJobDir '/' scriptName ];
    remoteCommand = sprintf('qsub %s', remoteScriptLocation);

    % Execute the submit command on the remote host.
    runCmdOnCluster(remoteCommand, clusterHost);

    % Delete the local copy of the script
    delete(localScript);
end

function filename = createPBSSubmitScript(commandToRun, decodeFunction, ...
         storageConstructor, remoteDataLocation, jobLocation, ...
         taskLocation, quote)
%Create a PBS submit script that forwards the required environment
%variables and runs a MATLAB worker.
logFileLocation = [quote remoteDataLocation '/' taskLocation '.log' quote];

% Create the commands to set the environment variables.
% "-N Job#"  specifies the job name
% "-j oe" joins together output and error streams
% "-o ..." specifies where standard output goes to
setEnv = sprintf([ ...
    '#!/bin/sh', '\n', ...
    '#PBS -N MATLAB_', taskLocation,  '\n', ...
    '#PBS -j oe', '\n', ...
    '#PBS -o ', logFileLocation, '\n', ...
    '#PBS -Wgroup_list=matlab\n', ...   %snow
     '#PBS -q pmatlab_q\n', ... %snow
    '#PBS -lpartition=ITHACA\n', ... %snow
    '#PBS -A matlab0001\n', ... %snow
    '#PBS -v MDCE_DECODE_FUNCTION=', decodeFunction, ',', ...
    'MDCE_STORAGE_LOCATION=', remoteDataLocation, ',', ...
    'MDCE_STORAGE_CONSTRUCTOR=', storageConstructor, ',', ...
    'MDCE_JOB_LOCATION=', jobLocation, ',', ...
    'MDCE_TASK_LOCATION=', taskLocation, ',', ...
    'MDCE_DEBUG=', 'true', ...
    '\n', ...
    ]);

% Content of script.
scriptContent = sprintf('%s%s', setEnv, commandToRun);

% Create script.
filename = tempname;
% Open file in binary mode to make it cross-platform.
fid = fopen(filename, 'w');
fprintf(fid, scriptContent);
fclose(fid);


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [chi, psi] = coeff_b(k, x1, x2, a, b)
% compute chi and psi for cosine method

arg2 = k .* pi * diag((x2 - a) ./ (b - a));     % arg trig func
arg1 = k .* pi * diag((x1 - a) ./ (b - a));     % arg trig func

term1 = cos( arg2 ) * diag(exp(x2));
term2 = cos( arg1 ) * diag(exp(x1));

term3 = pi * k .* sin( arg2 ) * diag(exp(x2)./ (b-a));
term4 = pi * k .* sin( arg1 ) * diag(exp(x1)./ (b-a));

chi = 1 ./ ( 1 + ((k .* pi) *diag(1./ (b - a))).^2 ) ...
    .* ( term1 - term2 + term3 - term4 );   % modify init

chi(1,:) = (exp(x2)-exp(x1)); % init chi

psi = ((sin(arg2) - sin(arg1)) ./ (k .* pi)) *diag(b-a);    % modify init

psi(1,:) = (x2-x1);           % init psi

end

author: Sable

repo: mcbench-benchmarks


text:
function tutorial_visual_group(ProtocolName, reports_dir)
% TUTORIAL_VISUAL_GROUP: Runs the Brainstorm/SPM group analysis pipeline (group analysis, BIDS VERSION)
%
% ONLINE TUTORIALS:
%    - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle
%    - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualGroup
%
% INPUTS:
%    - ProtocolName : Name of the protocol in which the recordings for the 19 subjects have been imported (TutorialVisual or TutorialGroup)
%                     This folder must include: the run-level averages (recordings and time-frequency) and the source kernels for each run
%    - reports_dir  : If defined, exports all the reports as HTML to this folder

% @=============================================================================
% 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.
% =============================================================================@
%
% Author: Francois Tadel, Elizabeth Bock, 2016-2018

% ===== CHECK PROTOCOL =====
% Start Brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin < 1) || isempty(ProtocolName)
    ProtocolName = 'TutorialGroup';
end
% Select current protocol
iProtocol = bst_get('Protocol', ProtocolName);
if isempty(iProtocol)
    error(['Unknown protocol: ' ProtocolName]);
end
if (iProtocol ~= bst_get('iProtocol'))
    gui_brainstorm('SetCurrentProtocol', iProtocol);
end
% Process: Select results files in: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch
sFiles = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname', 'sub-01', ...
    'condition',   'sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch');
if isempty(sFiles)
    error(['No source files available in folder: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch.' 10 ...
           'You should run tutorial_visual_single first, or download the ' 10 ...
           'protocol TutorialGroup.zip from the Brainstorm website.']);
end
% Process: Select file comments with tag: Unconstr
sFiles = bst_process('CallProcess', 'process_select_tag', sFiles, [], 'tag', 'Unconstr');
if ~isempty(sFiles)
    isUnconstrained = 1;
else
    isUnconstrained = 0;
end
% Start a new report
bst_report('Start');

% ===== DELETE ALL THE EXISTING FILES =====
GroupSubjectName = bst_get('NormalizedSubjectName');
DirIntra = bst_get('DirAnalysisIntra');
% Delete group analysis
[sSubject, iSubject] = bst_get('Subject', GroupSubjectName);
if ~isempty(sSubject)
    db_delete_subjects(iSubject);
end
% Get all the files in the intra subjects
sAvgSrc = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% If there are existing files: delete the contents of the folders
if ~isempty(sAvgSrc)
    % Get the base folder of the protocol
    ProtocolInfo = bst_get('ProtocolInfo', iProtocol);
    % Get all the folders
    AllIntraDir = unique(cellfun(@bst_fileparts, {sAvgSrc.FileName}, 'UniformOutput', 0));
    % Loop over the intra folders
    for iSubj = 1:length(AllIntraDir)
        % Get all the files
        subjIntraDir = dir(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, '*.mat'));
        subjIntraFiles = setdiff({subjIntraDir.name}, 'brainstormstudy.mat');
        % Delete all the files
        for iFile = 1:length(subjIntraFiles)
            delete(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, subjIntraFiles{iFile}));
        end
    end
    % Reload protocol
    db_reload_database(iProtocol);
end

%% ===== SUBJECT AVERAGES: FAMOUS, UNFAMILIAR, SCRAMBLED ==============================
%  ====================================================================================

% ===== SUBJECT AVERAGE: MEG+EEG =====
% Process: Select data files in: */*/Avg
sAvgRun = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname', 'All', ...
    'tag',         'Avg');
% Process: Weighted Average: By trial group (subject average)
sAvgSubj = bst_process('CallProcess', 'process_average', sAvgRun, [], ...
    'avgtype',    6, ...  % By trial group (subject average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   1, ...
    'keepevents', 1);
% Process: Average: By trial group (grand average)
sAvgGroup = bst_process('CallProcess', 'process_average', sAvgSubj, [], ...
    'avgtype',    7, ...  % By trial group (grand average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== SUBJECT AVERAGES: SOURCES (EEG) =====
% Process: Select source files in: */*/EEG
sAvgRunSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'EEG');
% Process: Weighted Average: By trial group (subject average) - EEG
sAvgSubjSrcEeg = bst_process('CallProcess', 'process_average', sAvgRunSrcEeg, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 0);
% Process: Add tag: EEG
sAvgSubjSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcEeg, [], ...
    'tag',    'EEG', ...
    'output', 1);  % Add to comment

% ===== SUBJECT AVERAGES: SOURCES (MEG) =====
% Process: Select source files in: */*/MEG
sAvgRunSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'MEG');
% Process: Weighted Average: By trial group (subject average) - MEG
sAvgSubjSrcMeg = bst_process('CallProcess', 'process_average', sAvgRunSrcMeg, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: MEG
sAvgSubjSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcMeg, [], ...
    'tag',    'MEG', ...
    'output', 1);  % Add to comment

% ===== SUBJECT AVERAGES: TIMEFREQ =====
% Process: Select time-frequency files in: */*
sAvgRunTf = bst_process('CallProcess', 'process_select_files_timefreq', [], []);
% Process: Weighted Average: By trial group (subject average)
sAvgSubjTf = bst_process('CallProcess', 'process_average', sAvgRunTf, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 1);

%% ===== SUBJECT AVERAGES: FACES ======================================================
%  ====================================================================================

% ===== FACES AVERAGE: MEG+EEG ======
% Process: Select data files in: */*
sAllData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous (6 files)"
sAvgSubjFamous = bst_process('CallProcess', 'process_select_tag', sAllData, [], 'tag', 'WAvg: Avg: Famous (6 files)');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar (6 files)"
sAvgSubjUnfamiliar = bst_process('CallProcess', 'process_select_tag', sAllData, [],  'tag', 'WAvg: Avg: Unfamiliar (6 files)');
% Process: Weighted Average A&B
sAvgSubjFaces = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamous, sAvgSubjUnfamiliar, ...
    'weighted', 1);
% Process: Add tag: Faces
sAvgSubjFaces = bst_process('CallProcess', 'process_add_tag', sAvgSubjFaces, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFaces = bst_process('CallProcess', 'process_set_comment', sAvgSubjFaces, [], ...
    'tag',     'WAvg: Avg: Faces', ...
    'isindex', 0);
% Process: Average: By trial group (grand average)
sAvgGroupFaces = bst_process('CallProcess', 'process_average', sAvgSubjFaces, [], ...
    'avgtype',    7, ...  % By trial group (grand average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== FACES AVERAGE: SOURCES (EEG) ======
% Process: Select source files in: */*/EEG
sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'EEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ...
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: Faces
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcEeg, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcEeg, [], ...
    'tag',     'WAvg: Avg: Faces | EEG', ...
    'isindex', 0);

% ===== FACES AVERAGE: SOURCES (MEG) ======
% Process: Select source files in: */*
sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'MEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: Famous
sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Famous');
% Process: Select file comments with tag: Unfamiliar
sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ...
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: Faces
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcMeg, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcMeg, [], ...
    'tag',     'WAvg: Avg: Faces | MEG', ...
    'isindex', 0);

% ===== FACES AVERAGE: TIMEFREQ ======
% Process: Select time-frequency files in: */*
sAllTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesTf = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousTf, sAvgSubjUnfamiliarTf, ...
    'weighted', 1);
% Process: Add tag: Faces
sAvgSubjFacesTf = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesTf, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesTf = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesTf, [], ...
    'tag',     'WAvg: Avg: Faces', ...
    'isindex', 0);

%% ===== SUBJECT AVERAGES: WITHIN-SUBJECT DIFFERENCES =================================
%  ====================================================================================

% ===== FACES - SCRAMBLED: SOURCES EEG =====
% Process: Select source files in: */*/EEG
sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'EEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Faces"
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Faces');
% Process: Select file comments with tag: "WAvg: Avg: Scrambled"
sAvgSubjScrambledSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Scrambled');
% Process: Difference: A-B
sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcEeg, sAvgSubjScrambledSrcEeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcEeg, [], ...
    'tag',     'Faces - Scrambled | EEG', ...
    'isindex', 0);

% ===== FACES - SCRAMBLED: SOURCES MEG =====
% Process: Select source files in: */*/MEG
sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'MEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Faces"
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Faces');
% Process: Select file comments with tag: "WAvg: Avg: Scrambled"
sAvgSubjScrambledSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Scrambled');
% Process: Difference: A-B
sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcMeg, sAvgSubjScrambledSrcMeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcMeg, [], ...
    'tag',     'Faces - Scrambled | MEG', ...
    'isindex', 0);

% ===== FAMOUS - UNFAMILIAR: SOURCE EEG =====
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Difference: A-B
sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcEeg, [], ...
    'tag',     'Famous - Unfamiliar | EEG', ...
    'isindex', 0);

% ===== FAMOUS - UNFAMILIAR: SOURCE MEG =====
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Difference: A-B
sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcMeg, [], ...
    'tag',     'Famous - Unfamiliar | MEG', ...
    'isindex', 0);

%% ===== SUBJECT AVERAGES: FILTER AND NORMALIZE =======================================
%  ====================================================================================

% ===== FILTER: MEG+EEG =====
% Process: Select data files in: */Intra
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Low-pass:32Hz
sAvgData = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'sensortypes', 'MEG, EEG', ...
    'highpass',    0, ...
    'lowpass',     32, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Extract time: [-202ms,900ms]
sAvgData = bst_process('CallProcess', 'process_extract_time', sAvgData, [], ...
    'timewindow', [-0.2018, 0.9], ...
    'overwrite',  1);

% ===== FILTER AND NORMALIZE: SOURCES =====
% Process: Select source files in: */Intra
sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Low-pass:32Hz
sAvgSrc = bst_process('CallProcess', 'process_bandpass', sAvgSrc, [], ...
    'highpass',    0, ...
    'lowpass',     32, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Extract time: [-202ms,900ms]
sAvgSrc = bst_process('CallProcess', 'process_extract_time', sAvgSrc, [], ...
    'timewindow', [-0.2018, 0.9], ...
    'overwrite',  1);
% Process: Z-score transformation: [-202ms,-5ms]
sAvgSrc = bst_process('CallProcess', 'process_baseline_norm', sAvgSrc, [], ...
    'baseline',   [-0.2018, -0.005], ...
    'source_abs', 0, ...
    'method',     'zscore', ...  % Z-score transformation:    x_std = (x - μ) / σ
    'overwrite',  1);

% ===== NORMALIZE: TIMEFREQ =====
% Process: Select time-frequency files in: */Intra
sAvgTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Event-related perturbation (ERS/ERD): [-200ms,-4ms]
sAvgTf = bst_process('CallProcess', 'process_baseline_norm', sAvgTf, [], ...
    'baseline',  [-0.2, -0.004], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 1);

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_1prepare.html']));
end

%% ===== SUBJECT AVERAGES: SCREEN CAPTURES ============================================
%  ====================================================================================
% Start a new report
bst_report('Start');
% Set colormap: local color scale
bst_colormaps('SetMaxMode', 'stat2', 'local');

% ===== FACES: RECORDINGS ======
% Process: Select data files in: */*/WAvg: Avg: Faces |
sAvgFaces = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'tag',           'WAvg: Avg: Faces |', ...
    'includeintra',  1);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ...
    'target',         1, ...  % Sensor/MRI registration
    'modality',       1);     % MEG
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ...
    'target',         5, ...  % Recordings time series
    'modality',       4, ...  % EEG
    'time',           0.1055);

% ===== FACES: SOURCES MEG ======
% Process: Select source files in: */*/WAvg: Avg: Faces | MEG
sAvgFacesSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Faces | MEG', ...
    'includeintra',  1);
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sAvgFacesSrcMeg, [], ...
    'target',         8, ...  % Sources (one time)
    'orient',         6, ...  % back
    'time',           0.1055, ...
    'threshold',      10);

% ===== FACES: TIME-FREQ EEG ======
% Process: Select time-frequency files in: */*/WAvg: Avg: Faces
sAvgFacesTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'tag',           'WAvg: Avg: Faces', ...
    'includeintra',  1);
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sAvgFacesTf, [], ...
    'target',         14, ...  % Time-frequency maps
    'time',           0.1055, ...
    'rowname',        'EEG070');

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_2snapshots.html']));
end

%% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG =====================================
%  ====================================================================================
% Start a new report
bst_report('Start');
% Set colormap: global color scale
bst_colormaps('SetMaxMode', 'meg', 'global');
bst_colormaps('SetMaxMode', 'eeg', 'global');
% Set display properties for statistics: p<0.05, FDR-corrected
StatThreshOptions = bst_get('StatThreshOptions');
StatThreshOptions.pThreshold = 0.05;
StatThreshOptions.Correction = 'fdr';
StatThreshOptions.Control    = [1 2 3];
bst_set('StatThreshOptions', StatThreshOptions);

% ===== GRAND AVERAGES =====
% Process: Select data files in: Group_analysis/Intra
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   GroupSubjectName, ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Take screen captures
DataScreenCapture(sAvgData(1), 1, 1, '');
DataScreenCapture(sAvgData(2), 1, 1, '');
DataScreenCapture(sAvgData(3), 1, 1, '');
DataScreenCapture(sAvgData(4), 1, 1, '');
% Define function for screen captures
function DataScreenCapture(sFiles, isMEG, isEEG, Comment)
    if isMEG
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         5, ...  % Recordings time series
            'modality',       3, ...  % MEG (Magnetometers)
            'time',           0.1055, ...
            'Comment',        ['MEG MAG: ' Comment]);
        % Process: Snapshot: Recordings topography (contact sheet)
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         7, ...  % Recordings topography (contact sheet)
            'modality',       3, ...  % MEG (Magnetometers)
            'contact_time',   [0.050, 0.300], ...
            'contact_nimage', 6, ...
            'Comment',        ['MEG MAG: ' Comment]);
    end
    if isEEG
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         5, ...  % Recordings time series
            'modality',       4, ...  % EEG
            'time',           0.1055, ...
            'Comment',        ['EEG: ' Comment]);
        % Process: Snapshot: Recordings topography (contact sheet)
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         7, ...  % Recordings topography (contact sheet)
            'modality',       4, ...  % EEG
            'contact_time',   [0.050, 0.300], ...
            'contact_nimage', 6, ...
            'Comment',        ['EEG: ' Comment]);
    end
end

% ===== FACES-SCRAMBLED: MEAN DIFF =====
% Process: Select data files in: */DirIntra/"Avg: Faces | "
sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Faces | ', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)"
sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Scrambled (6 files)', ...
    'includeintra',  1);
% Process: Difference of means [mean]
sDiffFacesData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Set comment
sDiffFacesData = bst_process('CallProcess', 'process_set_comment', sDiffFacesData, [], ...
    'tag',     'Faces - Scrambled', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sDiffFacesData, 1, 1, sDiffFacesData.Comment);

% ===== FACES-SCRAMBLED: PARAMETRIC T-TEST =====
% Process: Select data files in: */DirIntra/"Avg: Faces | "
sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Faces | ', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)"
sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Scrambled (6 files)', ...
    'includeintra',  1);
% Process: t-test paired [ALL]          H0:(A=B), H1:(A<>B)
sStatParamFacesData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'timewindow',    [], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
    'tail',          'two');  % Two-tailed
% Process: Set comment
sStatParamFacesData = bst_process('CallProcess', 'process_set_comment', sStatParamFacesData, [], ...
    'tag',     'Faces - Scrambled: Parametric t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatParamFacesData, 1, 1, sStatParamFacesData.Comment);

% ===== FACES-SCRAMBLED: NON-PARAMETRIC T-TEST =====
% Process: Perm t-test paired [All]          H0:(A=B), H1:(A<>B)
sStatPermFacesData = bst_process('CallProcess', 'process_test_permutation2p', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'timewindow',     [], ...
    'sensortypes',    '', ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgrow',         0, ...
    'iszerobad',      0, ...
    'Comment',        '', ...
    'test_type',      'ttest_paired', ...  % Paired Student's t-test T = mean(A-B) / std(A-B) * sqrt(n)
    'randomizations', 1000, ...
    'tail',           'two');  % Two-tailed
% Process: Set comment
sStatPermFacesData = bst_process('CallProcess', 'process_set_comment', sStatPermFacesData, [], ...
    'tag',     'Faces - Scrambled: Permutation t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatPermFacesData, 1, 1, sStatPermFacesData.Comment);

% ===== FACES-SCRAMBLED: CLUSTER T-TEST =====
% Process: FT t-test paired cluster [all EEG]          H0:(A=B), H1:(A<>B)
sStatClustFacesData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'sensortypes',    'EEG', ...
    'timewindow',     [], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  2, ...  % Paired t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...  % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% If FieldTrip is available and results were returned
if ~isempty(sStatClustFacesData)
    % Process: Set comment
    sStatClustFacesData = bst_process('CallProcess', 'process_set_comment', sStatClustFacesData, [], ...
        'tag',     'Faces - Scrambled: Cluster t-test EEG', ...
        'isindex', 0);
    % Take screen captures
    DataScreenCapture(sStatClustFacesData, 0, 1, sStatClustFacesData.Comment);
end

%% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG =====================================
%  ====================================================================================
% ===== FAMOUS-UNFAMILIAR: MEAN DIFF =====
% Process: Select data files in: */DirIntra/"Avg: Famous | "
sSubjAvgFamousData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Famous (6 files)', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Unfamiliar (6 files)"
sSubjAvgUnfamiliarData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Unfamiliar (6 files)', ...
    'includeintra',  1);
% Process: Difference of means [mean]
sDiffFamousData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Set comment
sDiffFamousData = bst_process('CallProcess', 'process_set_comment', sDiffFamousData, [], ...
    'tag',     'Famous - Unfamiliar', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sDiffFamousData, 1, 1, sDiffFamousData.Comment);

% ===== FAMOUS-UNFAMILIAR: PARAMETRIC T-TEST =====
% Process: t-test paired [All]          H0:(A=B), H1:(A<>B)
sStatParamFamousData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'timewindow',    [], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
    'tail',          'two');  % Two-tailed
% Process: Set comment
sStatParamFamousData = bst_process('CallProcess', 'process_set_comment', sStatParamFamousData, [], ...
    'tag',     'Famous - Unfamiliar: Parametric t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatParamFamousData, 1, 1, sStatParamFamousData.Comment);

% ===== FAMOUS-UNFAMILIAR: CLUSTER T-TEST =====
% Process: FT t-test paired cluster [all EEG]          H0:(A=B), H1:(A<>B)
sStatClustFamousData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'sensortypes',    'EEG', ...
    'timewindow',     [], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  2, ...  % Paired t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...  % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% If FieldTrip is available and results were returned
if ~isempty(sStatClustFacesData)
    % Process: Set comment
    sStatClustFamousData = bst_process('CallProcess', 'process_set_comment', sStatClustFamousData, [], ...
        'tag',     'Famous - Unfamiliar: Cluster t-test EEG', ...
        'isindex', 0);
    % Take screen captures
    DataScreenCapture(sStatClustFamousData, 0, 1, sStatClustFamousData.Comment);
end

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_3meeg.html']));
end

%% ===== FACES-SCRAMBLED: SOURCES (MEG) ===============================================
%  ====================================================================================
% Start a new report
bst_report('Start');

% ===== PROJECT SOURCES =====
% Process: Select source files in: */Intra
sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Absolute values / Norm for unconstrained
if isUnconstrained
    % Process: Unconstrained to flat map
    sAvgSrc = bst_process('CallProcess', 'process_source_flat', sAvgSrc, [], ...
        'method', 1);  % Norm: sqrt(x^2+y^2+z^2)
else
    sAvgSrc = bst_process('CallProcess', 'process_absolute', sAvgSrc, [], ...
        'overwrite', 1);
end
% Process: Project on default anatomy
sProjSrc = bst_process('CallProcess', 'process_project_sources', sAvgSrc, [], ...
    'headmodeltype', 'surface');  % Cortex surface
% Process: Spatial smoothing (3.00,abs)
sProjSrc = bst_process('CallProcess', 'process_ssmooth_surfstat', sProjSrc, [], ...
    'fwhm',       3, ...
    'overwrite',  1, ...
    'source_abs', 1);

% ===== MOVE FILES =====
% File tags / destination folder
ListTags = {'Avg: Faces | EEG',                 'Faces_EEG'; ...
            'Avg: Faces | MEG',                 'Faces_MEG'; ...
            'Avg: Famous (6 files) | EEG',      'Famous_EEG'; ...
            'Avg: Famous (6 files) | MEG',      'Famous_MEG'; ...
            'Avg: Unfamiliar (6 files) | EEG',  'Unfamiliar_EEG'; ...
            'Avg: Unfamiliar (6 files) | MEG',  'Unfamiliar_MEG'; ...
            'Avg: Scrambled (6 files) | EEG',   'Scrambled_EEG'; ...
            'Avg: Scrambled (6 files) | MEG',   'Scrambled_MEG'; ...
            'Faces - Scrambled | EEG',          'Faces-Scrambled_EEG'; ...
            'Faces - Scrambled | MEG',          'Faces-Scrambled_MEG'; ...
            'Famous - Unfamiliar | EEG',        'Famous-Unfamiliar_EEG'; ...
            'Famous - Unfamiliar | MEG',        'Famous-Unfamiliar_MEG'};
% Move group by group
for i = 1:size(ListTags,1)
    % Find names of files to move
    iTag = find(~cellfun(@(c)isempty(strfind(c, ListTags{i,1})), {sProjSrc.Comment}));
    if isempty(iTag)
        continue;
    end
    % Process: Move files
    bst_process('CallProcess', 'process_movefile', {sProjSrc(iTag).FileName}, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      ListTags{i,2});
end

% ===== LOOP ON MODALITY =====
AllModalities = {'MEG','EEG'};
for iMod = 1:length(AllModalities)
    % Selected modality for this loop
    Mod = AllModalities{iMod};

    % ===== |FACES-SCRAMBLED|: MEAN DIFF ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_', Mod]);
    % Process: Weighted Average: By folder (grand average)
    sAbsDiffFaces = bst_process('CallProcess', 'process_average', sDiffFaces, [], ...
        'avgtype',         4, ...  % By folder (grand average)
        'avg_func',        1, ...  % Arithmetic average:  mean(x)
        'weighted',        0, ...
        'scalenormalized', 0);
    % Process: Set comment
    sAbsDiffFaces = bst_process('CallProcess', 'process_set_comment', sAbsDiffFaces, [], ...
        'tag',     ['mean(|Faces-Scrambled|) | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sAbsDiffFaces = bst_process('CallProcess', 'process_movefile', sAbsDiffFaces, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Colormap: Custom max: [-10,+10] Z
    bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces-Scrambled|)']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces-Scrambled|)']);

    % ===== |FACES-SCRAMBLED|: PARAMETRIC CHI2-TEST ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_' Mod]);
    % Process: Chi2-test [all]          H0:(|Zi| = 0)
    sChiParamFaces = bst_process('CallProcess', 'process_test_parametric1', sDiffFaces, [], ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'chi2_onesample', ...  % One-sample Chi2 test     Zi~N(0,1), i=1..nQ = sum(|Zi|^2)     Q~Chi2(n)
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sChiParamFaces = bst_process('CallProcess', 'process_set_comment', sChiParamFaces, [], ...
        'tag',     ['|Faces-Scrambled|=0: Parametric Chi2 test | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sChiParamFaces = bst_process('CallProcess', 'process_movefile', sChiParamFaces, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']);

    % ===== LOG(|FACES-SCRAMBLED|): PARAMETRIC CHI2-TEST ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_' Mod]);
    % Process: Duplicate folders: Add tag "_log"
    sDiffFacesLog = bst_process('CallProcess', 'process_duplicate', sDiffFaces, [], ...
        'target', 2, ...  % Duplicate folders
        'tag',    '_log');
    % Process: Run Matlab command
    sDiffFacesLog = bst_process('CallProcess', 'process_matlab_eval', sDiffFacesLog, [], ...
        'matlab',    'Data = log(Data);', ...
        'overwrite', 1);
    % Process: Chi2-test [all]          H0:(|Zi| = 0)
    sChiParamFacesLog = bst_process('CallProcess', 'process_test_parametric1', sDiffFacesLog, [], ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'chi2_onesample', ...  % One-sample Chi2 test     Zi~N(0,1), i=1..nQ = sum(|Zi|^2)     Q~Chi2(n)
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sChiParamFacesLog = bst_process('CallProcess', 'process_set_comment', sChiParamFacesLog, [], ...
        'tag',     ['log(|Faces-Scrambled|)=0: Parametric Chi2 test | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sChiParamFacesLog = bst_process('CallProcess', 'process_movefile', sChiParamFacesLog, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']);

    % ===== |FACES|-|SCRAMBLED|: MEAN DIFF ======
    % Process: Select source files in: Group_analysis/Faces_MOD
    sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces_' Mod]);
    % Process: Select source files in: Group_analysis/Scrambled_MOD
    sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Scrambled_' Mod]);
    % Process: Difference of means [abs(mean)]
    sDiffAbsFaces = bst_process('CallProcess', 'process_diff_mean', sAvgFaces, sAvgScrambled, ...
        'avg_func',   2, ...  % Absolute value of average abs(mean(A)) - abs(mean(B))
        'weighted',   0);
    % Process: Set comment
    sDiffAbsFaces = bst_process('CallProcess', 'process_set_comment', sDiffAbsFaces, [], ...
        'tag',     ['mean(|Faces|)-mean(|Scrambled|) | ' Mod], ...
        'isindex', 0);
    % Colormap: Custom max: [-10,+10] Z
    bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces|)-mean(|Scrambled|)']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces|)-mean(|Scrambled|)']);

    % ===== |FACES|-|SCRAMBLED|: PARAMETRIC T-TEST ======
    % Process: Select source files in: Group_analysis/Faces_MOD
    sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces_' Mod]);
    % Process: Select source files in: Group_analysis/Scrambled_MOD
    sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Scrambled_' Mod]);
    % Process: t-test paired [-202ms,900ms]          H0:(A=B), H1:(A<>B)
    sTtestParamFaces = bst_process('CallProcess', 'process_test_parametric2p', sAvgFaces, sAvgScrambled, ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sTtestParamFaces = bst_process('CallProcess', 'process_set_comment', sTtestParamFaces, [], ...
        'tag',     ['|Faces|=|Scrambled|: Parametric t-test | ' Mod], ...
        'isindex', 0);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces|-|Scrambled|: Parametric t-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces|-|Scrambled|: Parametric t-test']);
end

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_4sources.html']));
end

end

author: brainstorm-tools


text:
%INTLAB long toolbox (rudimentary, slow but correct)
%
%Long constructors
%  long         - Long constructor
%  long2dble    - Long to double
%  long2intval  - Long to intval
%  vertcat      - Vertical concatenation            [ ; ]
%  subsasgn     - Subscripted assignment A(2:4) = 1
%  subsref      - Subscripted reference r = A(3)
%
%Display of long numbers
%  display      - Command window display for long
%  disp         - Display function for pop-up windows in debugger
%
%Long arithmetic operations
%  plus         - Plus                              +
%  uplus        - Unary plus                        +
%  minus        - Minus                             -
%  uminus       - Unary minus                       -
%  mtimes       - Elementwise multiply              *
%  times        - Elementwise multiply              .*
%  mrdivide     - Elementwise right division        /
%  mldivide     - Elementwise left division         \
%  rdivide      - Elementwise right division        ./
%  ldivide      - Elementwise left division         .\
%  mpower       - Elementwise power                 ^
%  power        - Elementwise power                 .^
%  longshift    - Shift by r bits
%
%Other long operations
%  min          - minimum
%  max          - maximum
%  abs          - absolute value
%  mid          - midpoint
%  rad          - radius
%  diam         - diameter
%  inf          - infimum
%  sup          - supremum
%
%Long comparisons
%  eq           - equal                             ==
%  ne           - not equal                         ~=
%  le           - less equal                        <=
%  lt           - less than                         <
%  ge           - greater equal                     >=
%  gt           - greater than                      >
%
%Utility routines
%
%  isempty      - Long is empty in Matlab sense, i.e. []
%  isnan        - True for Not a Number
%  end          - determine last index
%
%Structural operations
%  length       - Length
%  size         - Size
%
%Other long operations
%  longprecision- Sets/gets current precision
%  addlongerror - Add error to long number
%
%Some sample functions
%  exp          - Exponential
%  longpi       - Long computation of Pi
%
%Initialization of INTLAB long package and system variables
%  longinit     - Initialization and definition of defaults
%
%Demonstration, samples
%  demolong     - Some examples for using INTLAB long package
%

% written  12/30/98     S.M. Rump
%
% Copyright (c) Siegfried M. Rump, head of the Institute for Reliable Computing,
%               Hamburg University of Technology

repo: OpenMAS


function [  ] = convertPLY( openDir, saveDir )
%CONVERTPLY Summary of this function goes here
%   Detailed explanation goes here

    % For KITTI Dataset

    fileExt = '*.bin';
    files = dir(fullfile(openDir,fileExt));

    for i=0:length(files)-1
        i

        fileID = sprintf('%06d', i);
        fileName = [openDir, fileID, '.bin' ];

        xyz = readVelodyne(fileName);
        xyz = xyz';
        xyz = xyz(:,1:3);

        pt = pointCloud(xyz);

        saveName = [saveDir, num2str(i), '.ply'];
        pcwrite(pt, saveName);

    end

end


text:
function F = metoffice_experiment_gpfa_gpkron(datanum, anomalies, validation)

%
% Load the data
%

date = datestr(now,'yyyymmdd');

if nargin < 3
  validation = true;
end

[data,dataset,folder,maskfile] = metoffice_getdata(datanum, anomalies, validation);

if anomalies
  disp('Model anomalies')
  comps = [0 1 0];
  comps_spatial = [0 1 0];
else
  disp('Don''t remove climatological averages')
  comps = [5 1 0];
  comps_spatial = [5 1 0];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GPFA inference
%

debug = false;

% Number of components
D = sum(comps);
% DEBUGGING!!
if debug
  maxiter = 1;
  N_samples = 2;
else
  maxiter = 1000;
  N_samples = 2000;
end

%
% Model for temporal X
%

covfunc_x = cell(D,1);
theta_x = cell(D,1);
is_pseudos_x = false(D,1);

% Inputs (assume uniformly spaced time instances which is not exactly correct)
in_x = linspace(data.time(1), data.time(end), length(data.time));
pseudo_x = in_x(1:10:end);

% Squared distances for covariance functions
D2_pp = sq_dist(pseudo_x);
D2_px = sq_dist(pseudo_x, in_x);
d2_x = zeros(size(in_x,2),1); %diag(D2_xx);
D_pp = sqrt(D2_pp);
D_px = sqrt(D2_px);
d_x = sqrt(d2_x);

% Distance matrices for covariance functions
d_xx = sqrt(sq_dist(in_x(:,1),in_x));

ind = 0;

% Periodic components (1 year period) with decay (rational quadratic)
% $$$ ind = 1:min(D,comps(1));
% $$$ fprintf('%d periodic components for X\n', length(ind));
% $$$ covfunc = @(D,D2) gp_cov_product(gp_cov_periodic(D, 'wavelength', 365), ...
% $$$                                  gp_cov_rq(D2));
% $$$ covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D_pp,D2_pp), 1e-6), ...
% $$$                                 covfunc(D_px,D2_px), ...
% $$$                                 covfunc(d_x,d2_x))};
% $$$ theta_x(ind) = columns_to_cells(...
% $$$     [linspace(1,1,length(ind))      % smoothness of the period
% $$$      365*linspace(20,1,length(ind)) % lengthscale of the decay (RQ)
% $$$      ones(1,length(ind))]);         % alpha (RQ)
% $$$ is_pseudos_x(ind) = true;
% Periodic components (1 year period) WITHOUT decay
if comps(1) > 0
  ind = 1:min(D,comps(1));
  fprintf('%d periodic components for X\n', length(ind));
  covfunc = @(D) gp_cov_periodic(D, 'wavelength', 365.26);
  covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D_pp), 1e-6), ...
                                  covfunc(D_px), ...
                                  covfunc(d_x))};
  theta_x(ind) = columns_to_cells(...
      [linspace(1,1,length(ind))]);      % smoothness of the period
  is_pseudos_x(ind) = true;
end

% Slow components (1-20 years): rational quadratic using pseudo inputs
if comps(2) > 0
  ind = ind(end) + (1:min(D-ind(end),comps(2)));
  fprintf('%d slow components for X\n', length(ind));
% $$$ covfunc = @(D2) gp_cov_se(D2); % SE
  covfunc = @(D2) gp_cov_rq(D2); % RQ
  covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D2_pp), 1e-3), ...
                                  covfunc(D2_px), ...
                                  covfunc(d2_x))};
% $$$ theta_x(ind) = columns_to_cells(...
% $$$     [365*linspace(20,1,length(ind))]); % lengthscale for SE
  theta_x(ind) = columns_to_cells(...
      [365*linspace(20,1,length(ind)) % lengthscale for RQ
       ones(1,length(ind))]);         % alpha for RQ
  is_pseudos_x(ind) = true;
end

% Fast components (4-18 months): piecewise polynomial in 1-D
% Take advantage of the Toeplitz structure of the covariance matrix
if comps(3) > 0
  ind = (ind(end)+1):D;
  fprintf('%d fast components for X\n', length(ind));
  covfunc_x(ind) = {gp_cov_jitter(gp_cov_toeplitz(gp_cov_pp(d_xx,1)))};
  theta_x(ind) = columns_to_cells(...
      [30*linspace(18,4,length(ind))]); % lengthscale or cut-off
  is_pseudos_x(ind) = false;
end

%
% Model for spatial W
%

covfunc_w = cell(D,1);
theta_w = cell(D,1);
is_pseudos_w = false(D,1);

% Remove land area grid points
in_w = data.coordinates;

%% Smooth components (using pseudo inputs)

% Pseudo inputs (uniformly with respect to area size)
pseudo_w = points_on_sphere(18); % uniform points by number of latitudes
% Remove pseudo inputs that are on land (the nearest grid point is land)
ind_pseudo_w = mohsst5_points_to_grid_index(pseudo_w);
mask = metoffice_get_mask(maskfile);
pseudo_w(:,~mask(ind_pseudo_w)) = [];
% $$$ % This code shows the pseudo inputs on the map
% $$$ figure
% $$$ map_projection('global-ellipse');
% $$$ map_plot(pseudo_w,'r+');
% $$$ map_coast()
% $$$ map_grid()
% $$$ return

% Transform inputs to 3-D Euclidean coordinates
in_w = geographic_to_euclidean(in_w);
pseudo_w = geographic_to_euclidean(pseudo_w);

% Squared distance matrices for the covariance functions
D2_ww = sq_dist(in_w);
D2_pp = sq_dist(pseudo_w);
D2_pw = sq_dist(pseudo_w, in_w);
d2_w = diag(D2_ww);

ind = 0;

if comps_spatial(1) > 0
  ind = ind(end) + (1:comps_spatial(1));
  fprintf('%d iid components for W\n', length(ind));

  covfunc_w(ind) = {gp_cov_scale(gp_cov_delta(size(in_w,2)))};
  theta_w(ind) = columns_to_cells(...
      [linspace(1,0.1,length(ind))]);       % magnitudes

end

if comps_spatial(2) > 0
  ind = ind(end) + (1:comps_spatial(2));
  fprintf('%d slow components for W (using %d pseudo inputs)\n', length(ind), ...
          size(pseudo_w,2));

  % Covariance function (scaled squared exponential) with pseudo inputs
  covfunc = @(D2) gp_cov_se(D2);
  covfunc_w(ind) = {gp_cov_pseudo(...
      gp_cov_scale(gp_cov_jitter(covfunc(D2_pp), 1e-3)), ...
      gp_cov_scale(covfunc(D2_pw)), ...
      gp_cov_scale(covfunc(d2_w)))};

  % Hyperparameters for the covariance functions
  theta_w(ind) = columns_to_cells(...
      [linspace(1,0.1,length(ind));       % magnitudes
       linspace(4000,1000,length(ind))]); % lengthscales
  is_pseudos_w(ind) = true;
end
%
% Process data
%

% Form the data matrix
Y = data.data;
[M,N] = size(Y);
Obs = ~isnan(Y);

% Filename for saving the results
folder = [folder '/gpfa_gpkron'];
mkdir(folder);
filename = sprintf('%s/results_rectest_%s_gpfa_D=%d_anomalies=%d_remval=%d_%s', ...
                   folder, ...
                   dataset, ...
                   D, ...
                   anomalies, ...
                   validation, ...
                   date);

% Component-wise factorization for X
X_module = factor_module_gp_factorized(N, covfunc_x, theta_x, ...
                                       'update_hyperparameters', [5 10:10:100 100:100:2000], ...
                                       'maxiter_hyperparameters', 5, ...
                                       'is_pseudo', is_pseudos_x, ...
                                       'init', zeros(D,N));

% Component-wise factorization for W
W_module = factor_module_gp_factorized(M, covfunc_w, theta_w, ...
                                       'update_hyperparameters', [5 10:10:100 100:100:2000], ...
                                       'maxiter_hyperparameters', 5, ...
                                       'is_pseudo', is_pseudos_w);

% Isotropic noise (precisions weighted proportionally to grid size)
weights = repmat(data.gridsize, [1, N]);
% $$$ figure
% $$$ mohsst5_mapplot(metoffice_add_land(weights(:,1)));
% $$$ return
noise_module = noise_module_isotropic(M, N, 1e-3, 1e-3, ...
                                      'init', 10, ...
                                      'weights', weights);

% Run GPFA
Q = gpfa(D, Y, W_module, X_module, noise_module, ...
         'maxiter', maxiter, ...
         'rotate', 1:100, ... %[1:50 60:10:2000], ...
         'autosavefile', filename, ...
         'autosave', [10:100:2000]);

Yh = Q.W'*Q.X;

% Some performance measures
fprintf('Weighted training RMSE of the reconstruction: %f\n',  ...
        rmsew(Y(Obs)-Yh(Obs),weights(Obs)));

% Save the results
save(filename, '-struct', 'Q');
fprintf('Saved GPFA results to %s\n', filename);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Short-scale GP inference for the residuals
%

% Reconstruct
if debug
  Yh = zeros(size(Y));
else
  Yh = Q.W'*Q.X;
end
Yres = Y - Yh;
clear Q;

%
% Temporal covariance function (assume uniformly spaced time instances)
%

d = abs(1-(1:length(data.time)));
covfunc1 = gp_cov_toeplitz(gp_cov_pp(d,1));
theta_temporal = [7];
% $$$ covfunc1 = gp_cov_toeplitz(gp_cov_sum(gp_cov_pp(d,1), ...
% $$$                                       gp_cov_scale(gp_cov_pp(d,1))));
% $$$ theta_temporal = [7;   % length scale 1
% $$$                   1.0; % magnitude 2
% $$$                   3];  % length scale 2

%
% Spatial covariance function
%

[LON,LAT] = meshgrid(data.lon,data.lat);

X = geographic_to_euclidean([LON(:)';LAT(:)']);

% Use block-Toeplitz structure for the covariance function
[lat,lon0] = meshgrid(data.lat,data.lon(1));
X0 = geographic_to_euclidean([lon0(:)';lat(:)']);
D = sqrt(sq_dist(X0,X));
covfunc2 = gp_cov_toeplitz_block(gp_cov_pp(D,3));
if datanum <= 4
  % 5x5 SST
  theta_spatial = [3000];  % length scale
elseif datanum == 5
  % 1x1 sea ice
  theta_spatial = [200];
end

% Select sea areas
sea = metoffice_get_mask(maskfile);
covfunc2 = gp_cov_select(covfunc2, sea);

%
% Inference
%

burnin = floor(N_samples/2);
folder_samples = sprintf('%s/samples_rectest_%s_gpkron_%s', ...
                         folder, ...
                         dataset, ...
                         date);
mkdir(folder_samples);
filename = sprintf('%s/results_rectest_%s_gpkron_D=%d_anomalies=%d_remval=%d_%s', ...
                   folder, ...
                   dataset, ...
                   sum(comps), ...
                   anomalies, ...
                   validation, ...
                   date);
filename_samples = sprintf('%s/samples_rectest_%s_gpkron_D=%d_anomalies=%d_remval=%d_%s', ...
                           folder_samples, ...
                           dataset, ...
                           sum(comps), ...
                           anomalies, ...
                           validation, ...
                           date);

a = 1e-3;
b = 1e-3;
logprior_theta = @(theta) sum(gamma_logpdf(theta, a, b));
dlogprior_theta = @(theta) gamma_dlogpdf(theta, a, b);

% Initial guess for covariance parameters
theta_init = [0.5; ...               % total magnitude
              theta_temporal(:); ... % temporal parameters
              theta_spatial(:); ...  % spatial parameters
              0.5]';                 % noise magnitude

samplefunc = get_sample_function2(numel(theta_init), N_samples, burnin, ...
                                                filename, filename_samples);

% Weights for the noise levels using the respective grid size
w = 1./sqrt(cosd(LAT));
W = repmat(w(sea), [1,size(Y,2)]);

[get_logpdf, get_dlogpdf, get_rand_model, func_theta] = ...
    gp_init_kron4(covfunc1, ...
                  solver_ldlchol(), ...
                  covfunc2, ...
                  solver_ldlchol(), ...
                  logprior_theta, ...
                  dlogprior_theta, ...
                  'samplefunc', samplefunc, ...
                  'rand_y', 'pcg', ...
                  'noise_scale', W, ...
                  'likelihood', 'whitened_prior');

% Transform to log-scale
func_theta = @(logtheta, varargin) func_theta_transformed(logtheta, ...
                                                  exp(logtheta), ...
                                                  func_theta, ...
                                                  varargin{:});
get_logpdf = @(f_theta) get_logpdf_transformed(f_theta, ...
                                               get_logpdf, ...
                                               sum(f_theta.theta_transformed));
get_dlogpdf = @(df_theta) get_dlogpdf_transformed(df_theta, ...
                                                  get_dlogpdf, ...
                                                  diag(exp(df_theta.theta_transformed)), ...
                                                  ones(size(df_theta.theta_transformed)));
theta_init = log(theta_init);

[rand_theta, f_theta] = mcmc_init_slicesampling(theta_init, ...
                                                get_logpdf, ...
                                                'fx', func_theta);

rand_model = get_rand_model(rand_theta, f_theta);

% Gibbs sampling
tic
gibbs(Yres,~Obs, N_samples, rand_model);
toc

% Check the results
res = samplefunc();

save(filename, '-struct', 'res');
disp(['Saved short-scale GP results to ', filename]);

% Mean reconstruction
F = Yh + res.F;

filename = sprintf('%s/results_rectest_%s_gpfa_gpkron_D=%d_anomalies=%d_remval=%d_%s', ...
                   folder, ...
                   dataset, ...
                   sum(comps), ...
                   anomalies, ...
                   validation, ...
                   date);
save(filename, 'F');
disp(['Saved total reconstruction to ', filename]);

fprintf('Weighted training RMSE of the reconstruction: %f\n',  ...
        rmsew(Y(Obs)-F(Obs),weights(Obs)));

if nargout < 1
  clear F;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function S = solver_ldlchol()

S.decompose = @decompose;
S.linsolve = @linsolve;
S.logdet = @logdet;
S.inv = @inv;
S.squareroot = @squareroot;

  function LD = decompose(K)
  [LD,p] = ldlchol(K);
  end

  function x = linsolve(LD, y)
  x = linsolve_ldlchol(LD,y);
  end

  function ldet = logdet(LD)
  ldet = logdet_ldlchol(LD);
  end

  function A = inv(LD)
  A = spinv_ldlchol(LD);
  end

  function L = squareroot(LD)
  L = ldlchol2lchol(LD);
  end

end

function [f_theta, df_theta] = func_theta_transformed(theta_transformed, ...
                                                  theta, func_theta, varargin)
if nargout <= 1
  f_theta = func_theta(theta, varargin{:});
  f_theta.theta_transformed = theta_transformed;
else
  [f_theta, df_theta] = func_theta(theta, varargin{:});
  f_theta.theta_transformed = theta_transformed;
  df_theta.theta_transformed = theta_transformed;
end
end

function logpdf_y = get_logpdf_transformed(fy, get_logpdf, logjacobian)
logpdf = get_logpdf(fy);
logpdf_y = @logpdf_transformed;
  function lpdf = logpdf_transformed(varargin)
  lpdf = logpdf(varargin{:}) + logjacobian;
  end
end

function dlogpdf_y = get_dlogpdf_transformed(dfy, get_dlogpdf, Jacobian, ...
                                                  dlogjacobian)
dlogpdf = get_dlogpdf(dfy);
dlogpdf_y = @dlogpdf_transformed;
  function dlpdf = dlogpdf_transformed(varargin)
  dlpdf = dlogpdf(varargin{:});
  dlpdf = Jacobian*dlpdf + dlogjacobian;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function samplefunc = get_sample_function2(D_theta, N, burnin, filename, filename_samples)
results.Y = [];
results.F = 0;
results.FF = 0;
results.theta = zeros(D_theta, N);
samplefunc = @process_sample;
n = 1;
  function res = process_sample(Y, F, theta)
  if nargin >= 1
    % Store results
    results.Y = Y;
    if true && n > burnin
      results.F = (F + (n-burnin-1)*results.F) / (n-burnin);
    end
    if true && n > burnin
      results.FF = (F.*F + (n-burnin-1)*results.FF) / (n-burnin);
    end
    results.theta(:,n) = theta(:);
    % Save results
    fprintf('Saving results to %s..', filename)
    save(filename, '-struct', 'results');
    save(sprintf('%s_F%d',filename_samples,n), 'F');
    fprintf(' done.\n')
    n = n + 1;
  end
  if nargout >= 1
    % Return results
    res = results;
  end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gibbs(Y, Imv, N, rand_model, samplefunc)

Y(Imv) = 0;

for n=1:N

  t = cputime();
  Y = rand_model(Y,Imv);
  dt = cputime() - t;
  fprintf('Iteration step %d done. (%f seconds)\n', n, dt)

end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab