text:
function nfigures = NumTotFigures(row,col,nvar)
% =======================================================================
% Determine the number of figures for a given number of subplots, rows
% and columns
% =======================================================================
% nfigures = NumTotFigures(row,col,nvar)
% -----------------------------------------------------------------------
% INPUT
%	- row: number of rows of the subplot
%	- col: number of columns of the subplot
%	- nvar: number of variables
% -----------------------------------------------------------------------
% OUTPUT
%	- nfigures: number of charts per subplot
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

NumGraphXPage = row*col;

p = floor(nvar./NumGraphXPage);
q = (nvar./NumGraphXPage);

if NumGraphXPage>=nvar
    nfigures=1;
elseif p==q
    nfigures=p;
else
    nfigures=p+1;
end
clear p q

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function f = comp_ifwt(c,g,a,J,Ls,ext)
%COMP_IFWT Compute Inverse DWT
%   Usage:  f = comp_ifwt(c,g,J,a,Ls,ext);
%
%   Input parameters:
%         c     : Cell array of length M = J*(filtNo-1)+1. Each element is Lc(m)*W array
%         g     : Synthesis wavelet filters - cell-array of length *filtNo*.
%         J     : Number of filterbank iterations.
%         a     : Upsampling factors - array of length *filtNo*.
%         Ls    : Length of the reconstructed signal.
%         ext   : 'per','zero','odd','even', Type of the forward transform boundary handling.
%
%   Output parameters:
%         f     : Reconstructed data - Ls*W array.
%

% see comp_fwt for explanantion
assert(a(1)==a(2),'First two elements of a are not equal. Such wavelet filterbank is not suported.');

% Impulse responses to a correct format.
filtNo = numel(g);
%gCell = cellfun(@(gEl) conj(flipud(gEl.h(:))),g,'UniformOutput',0);
gCell = cellfun(@(gEl) gEl.h(:),g,'UniformOutput',0);

if strcmp(ext,'per')
   % Initial shift of the filter to compensate for it's delay.
   % "Zero" delay reconstruction is produced.
   % offset = cellfun(@(gEl) gEl.offset,g);
   %offset = cellfun(@(gEl) 1-numel(gEl.h)-gEl.offset,g);
   offset = cellfun(@(gEl) gEl.offset,g);
elseif strcmp(ext,'valid')
   offset = -cellfun(@(gEl) numel(gEl.h)-1,g);
else
   % -1 + 1 = 0 is used for better readability and to be consistent
   % with the shift in comp_fwt.
   % Here we are cheating, because we are making the filters
   % anti-causal to compensate for the delay introduced by causal
   % analysis filters.
   % Instead, we could have used causal filters here and do the
   % delay compensation at the end (cropping f).
   % offset = -cellfun(@(gEl) numel(gEl),gCell) + (a -1) +1;
   offset = -(a-1);
end

Lc = cellfun(@(cEl) size(cEl,1),c);
Lc(end+1) = Ls;
tempca = c(1);
cRunPtr = 2;
for jj=1:J
   tempca=comp_ifilterbank_td([tempca;c(cRunPtr:cRunPtr+filtNo-2)],gCell,a,Lc(cRunPtr+filtNo-1),offset,ext);
   cRunPtr = cRunPtr + filtNo -1;
end
% Save reconstructed data.
f = tempca;

% for ch=1:chans
%   tempca = c(LcStart(1):LcEnd(1),ch);
%   LcRunPtr = filtNo+1;
%   cRunPtr = 2;
%   for jj=1:J
%      tempca = comp_upconv({tempca}, Lc(LcRunPtr),{tmpg{1}},a(1),skip(1),ext,0);
%      for ff=2:filtNo
%         % tempca = tempca + comp_upconv({c{cRunPtr}(:,ch)}, Lc(LcRunPtr),{tmpg},a(ff),skip,doNoExt,0);
%         tempca = tempca + comp_upconv({c(LcStart(cRunPtr):LcEnd(cRunPtr),ch)}, Lc(LcRunPtr),{tmpg{ff}},a(ff),skip(ff),ext,0);
%         cRunPtr = cRunPtr + 1;
%      end
%      LcRunPtr = LcRunPtr + filtNo -1;
%   end
%   f(:,ch) = tempca;
% end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% Tutorial for Inscribed_Rectangle Package
% *By Jarek Tuszynski*
%
% Inscribed_Rectangle package provides 2 low level computer vision / image
% analysis functions able to locate largest square or rectangle inscribed
% inside arbitrary shape defined by a binary mask (black and white image).
% Only rectangles with vertical/horizontal edges are considered. The
% functions proved can be used as tools for larger image segmentation
% problems.
%% Change History
% * 2010-07-07 - original version

%% Licence
% The package is distributed under BSD License
format compact; % viewing preference
clear variables; close all;
type('license.txt')

%% Create a mask defining a circle
[X,Y] = meshgrid(-200:200, -200:200);
BW = (X.^2 + Y.^2)<180^2;
figure(1); imshow(BW)

%% Run FindLargestSquares and return results.
% S stores for each pixel the size of the largest all-white square
% with its upper-left corner at that pixel
S = FindLargestSquares(BW);
imagesc(S); axis off

%% Use S to find the largest Square inscribed in the circle
imshow(BW)
[~, pos] = max(S(:));
[r c] = ind2sub(size(S), pos);
rectangle('Position',[c,r,S(r,c),S(r,c)], 'EdgeColor','r', 'LineWidth',3);

%% Run FindLargestRectangles and return results.
% C stores for each pixel the area of the largest all-white rectangle
% with its upper-left corner at that pixel. W and H store width and height
% of those rectangles
[C H W] = FindLargestRectangles(BW, [0 0 1]);
subplot(2,2,1); imagesc(BW); axis off; axis equal; colormap gray
title('Circle mask');
subplot(2,2,2); imagesc(H);  axis off; axis equal; colormap gray
title('Height of the rectangles');
subplot(2,2,3); imagesc(W);  axis off; axis equal; colormap gray
title('Width of the rectangles');
subplot(2,2,4); imagesc(C);  axis off; axis equal; colormap gray
title('Area of the rectangles');

%% Find the largest Rectangle inscribed in the circle, with size measured by area
close all; imshow(BW)
[~, pos] = max(C(:));
[r c] = ind2sub(size(S), pos);
rectangle('Position',[c,r,W(r,c),H(r,c)], 'EdgeColor','r', 'LineWidth',3);

%% Find the largest Rectangle inscribed in the circle with size measured as rectangle circumference with vertical edges having twice the weight of the horizontal edges
[C H W] = FindLargestRectangles(BW, [2 1 0]);
imshow(BW)
[~, pos] = max(C(:));
[r c] = ind2sub(size(S), pos);
rectangle('Position',[c,r,W(r,c),H(r,c)], 'EdgeColor','r', 'LineWidth',3);

%% Find the largest Rectangle inscribed in the circle with size measured as rectangle circumference with horizontal edges having 3 times the weight of the vertical edges
[C H W] = FindLargestRectangles(BW, [1 3 0]);
imshow(BW)
[tmp pos] = max(C(:));
[r c] = ind2sub(size(S), pos);
rectangle('Position',[c,r,W(r,c),H(r,c)], 'EdgeColor','r', 'LineWidth',3);

%% Load an image of Jenga tower and create a boolean mask of the tower shape.
% Image by "Jason7825" was copied from http://commons.wikimedia.org/wiki/File:Jenga_arrangement.jpg and is distributed under CC-BY-SA-3.0 & GDFL licenses
I = imread('Jenga_arrangement.jpg');
imshow(I)
BW = rgb2gray(I)>115;
BW = imfill(BW,'holes');
figure(2); imshow(BW)

%% Use FindLargestRectangles to locate rectangular blocks in the image
imshow(I)
for i = 1:10
  [C H W] = FindLargestRectangles(BW, [0 10 1]);
  [tmp pos] = max(C(:));
  [r c] = ind2sub(size(C), pos);
  rectangle('Position',[c,r,W(r,c),H(r,c)], 'EdgeColor','r', 'LineWidth',3);
  BW( r:(r+H(r,c)-1), c:(c+W(r,c)-1) ) = 0;
end

author: Sable

repo: mcbench-benchmarks


text:
% std_chaninds() - look up channel indices in a EEG structure
%
% Usage:
%         >> inds = std_chaninds(EEG, channames);
% Inputs:
%         EEG   - EEG structure containing a chanlocs substructure.
%                 the chanlocs structure may also be used as input.
%     channames - [cell] channel names. May also be a string containing
%                 one or several channel names.
%
% Outputs:
%       inds - [integer array] channel indices
%
% Author: Arnaud Delorme, CERCO, 2009-

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function finalinds = eeg_chaninds(EEG, channames, errorifnotfound);

    if nargin < 2
        help eeg_chaninds;
        return;
    end;
    if nargin < 3
        errorifnotfound = 1;
    end;

    if isfield(EEG, 'chanlocs')
         chanlocs = EEG.chanlocs;
    else chanlocs = EEG;
    end;

    % decode string if necessary
    % --------------------------
    if isstr(channames)
        channames = parsetxt( channames );
    end;

    finalinds   = [];
    if isempty(chanlocs)
         tmpallchans = [];
    else tmpallchans = lower({ chanlocs.labels });
    end;
    if isempty(channames), finalinds = [1:length(chanlocs)]; return; end;
    for c = 1:length(channames)
        chanind = strmatch( lower(channames{c}), tmpallchans, 'exact');
        if isempty(chanind),
            chanind = str2double(channames{c});
            if isnan(chanind), chanind = []; end;
            if errorifnotfound && isempty(chanind)
                error(sprintf('Channel %s not found', channames{c}));
            end;
        end;
        finalinds   = [ finalinds chanind ];
    end;
    finalinds = sort(finalinds);

author: PatternRecognition

repo: OpenBMI


text:
function [outim] = ft_sliceinterp(cfg, ininterp)

% FT_SLICEINTERP plots a 2D-montage of source reconstruction and anatomical MRI
% after these have been interpolated onto the same grid.
%
% Use as
%   ft_sliceinterp(cfg, interp)
%      or
%   [rgbimage] = ft_sliceinterp(cfg, interp), rgbimage is the monatage image
%
% where interp is the output of sourceinterpolate and cfg is a structure
% with any of the following fields:
%
% cfg.funparameter  string with the functional parameter of interest (default = 'source')
% cfg.maskparameter parameter used as opacity mask (default = 'none')
% cfg.clipmin       value or 'auto' (clipping of source data)
% cfg.clipmax       value or 'auto' (clipping of source data)
% cfg.clipsym       'yes' or 'no' (default) symmetrical clipping
% cfg.colormap      colormap for source overlay (default is jet(128))
% cfg.colmin        source value mapped to the lowest color (default = 'auto')
% cfg.colmax        source value mapped to the highest color (default = 'auto')
% cfg.maskclipmin   value or 'auto' (clipping of mask data)
% cfg.maskclipmax   value or 'auto' (clipping of mask data)
% cfg.maskclipsym   'yes' or 'no' (default) symmetrical clipping
% cfg.maskmap       opacitymap for source overlay
%                   (default is linspace(0,1,128))
% cfg.maskcolmin    mask value mapped to the lowest opacity, i.e.
%                   completely transparent (default ='auto')
% cfg.maskcolmin    mask value mapped to the highest opacity, i.e.
%                   non-transparent (default = 'auto')
% cfg.alpha         value between 0 and 1 or 'adaptive' (default)
% cfg.nslices       integer value, default is 20
% cfg.dim           integer value, default is 3 (dimension to slice)
% cfg.spacemin      'auto' (default) or integer (first slice position)
% cfg.spacemax      'auto' (default) or integer (last slice position)
% cfg.resample      integer value, default is 1 (for resolution reduction)
% cfg.rotate        number of ccw 90 deg slice rotations (default = 0)
% cfg.title         optional title (default is '')
% cfg.whitebg       'yes' or 'no' (default = 'yes')
% cfg.flipdim       flip data along the sliced dimension, 'yes' or 'no'
%                   (default = 'no')
% cfg.marker        [Nx3] array defining N marker positions to display
% cfg.markersize    radius of markers (default = 5);
% cfg.markercolor   [1x3] marker color in RGB (default = [1 1 1], i.e. white)
% cfg.interactive   'yes' or 'no' (default), interactive coordinates
%                   and source values
%
% if cfg.alpha is set to 'adaptive' the opacity of the source overlay
% linearly follows the source value: maxima are opaque and minima are
% transparent.
%
% if cfg.spacemin and/or cfg.spacemax are set to 'auto' the sliced
% space is automatically restricted to the evaluated source-space
%
% if cfg.colmin and/or cfg.colmax are set to 'auto' the colormap is mapped
% to source values the following way: if source values are either all
% positive or all negative the colormap is mapped to from
% min(source) to max(source). If source values are negative and positive
% the colormap is symmetrical mapped around 0 from -max(abs(source)) to
% +max(abs(source)).
%
% If cfg.maskparameter specifies a parameter to be used as an opacity mask
% cfg.alpha is not used. Instead the mask values are maped to an opacitymap
% that can be specified using cfg.maskmap. The mapping onto that
% opacitymap is controlled as for the functional data using the
% corresponding clipping and min/max options.
%
% if cfg.whitebg is set to 'yes' the function estimates the head volume and
% displays a white background outside the head, which can save a lot of black
% printer toner.
%
% if cfg.interactive is set to 'yes' a button will be displayed for
% interactive data evaluation and coordinate reading. After clicking the
% button named 'coords' you can click on any position in the slice montage.
% After clicking these coordinates and their source value are displayed in
% a text box below the button. The coordinates correspond to indeces in the
% input data array:
%
%   f = interp.source(coord_1,coord_2,coord_3)
%
% The coordinates are not affected by any transformations used for displaying
% the data such as cfg.dim, cfg.rotate,cfg.flipdim or cfg.resample.
%
% See also FT_SOURCEANALYSIS, FT_VOLUMERESLICE

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

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

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

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

% check if the input data is valid for this function
ininterp = ft_checkdata(ininterp, 'datatype', 'volume', 'feedback', 'yes');

% set the defaults
if ~isfield(cfg, 'clipmin');      cfg.clipmin = 'auto';        end
if ~isfield(cfg, 'clipmax');      cfg.clipmax = 'auto';        end
if ~isfield(cfg, 'clipsym');      cfg.clipsym = 'no';          end
if ~isfield(cfg, 'alpha');        cfg.alpha   = 'adaptive';    end
if ~isfield(cfg, 'nslices');      cfg.nslices = 20;            end
if ~isfield(cfg, 'dim');          cfg.dim = 3;                 end
if ~isfield(cfg, 'colormap');     cfg.colormap = jet(128);     end
if ~isfield(cfg, 'spacemin');     cfg.spacemin = 'auto';       end
if ~isfield(cfg, 'spacemax');     cfg.spacemax = 'auto';       end
if ~isfield(cfg, 'colmin');       cfg.colmin = 'auto';         end
if ~isfield(cfg, 'colmax');       cfg.colmax = 'auto';         end
if ~isfield(cfg, 'resample');     cfg.resample = 1;            end
if ~isfield(cfg, 'rotate');       cfg.rotate = 0;              end
if ~isfield(cfg, 'title');        cfg.title = '';              end
if ~isfield(cfg, 'whitebg');      cfg.whitebg = 'no';          end
if ~isfield(cfg, 'flipdim');      cfg.flipdim = 'no';          end
if ~isfield(cfg, 'marker');       cfg.marker = [];             end
if ~isfield(cfg, 'markersize');   cfg.markersize = 5;          end
if ~isfield(cfg, 'markercolor');  cfg.markercolor = [1,1,1];   end
if ~isfield(cfg, 'interactive');  cfg.interactive = 'no';      end
if ~isfield(cfg, 'maskclipmin');  cfg.maskclipmin = 'auto';           end
if ~isfield(cfg, 'maskclipmax');  cfg.maskclipmax = 'auto';           end
if ~isfield(cfg, 'maskclipsym');  cfg.maskclipsym = 'no';             end
if ~isfield(cfg, 'maskmap');      cfg.maskmap = linspace(0,1,128);    end
if ~isfield(cfg, 'maskcolmin');   cfg.maskcolmin = 'auto';            end
if ~isfield(cfg, 'maskcolmax');   cfg.maskcolmax = 'auto';            end
if ~isfield(cfg, 'maskparameter');cfg.maskparameter = [];             end

% perform some checks on the configuration for backward compatibility
if ~isfield(cfg, 'funparameter') && isfield(ininterp, 'source')
  % if present, the default behavior should be to use this field for plotting
  cfg.funparameter = 'source';
end

% make the selection of functional and mask data consistent with the data
cfg.funparameter  = parameterselection(cfg.funparameter, ininterp);
cfg.maskparameter = parameterselection(cfg.maskparameter, ininterp);
% only a single parameter should be selected
try, cfg.funparameter  = cfg.funparameter{1}; end
try, cfg.maskparameter = cfg.maskparameter{1}; end

% check anatomical data
if isfield(ininterp,'anatomy');
  interp.anatomy = reshape(ininterp.anatomy, ininterp.dim);
else
  ft_error('no anatomical data supplied');
end

% check functional data
if ~isempty(cfg.funparameter)
  interp.source = double(reshape(getsubfield(ininterp, cfg.funparameter), ininterp.dim));
else
  ft_error('no functional data supplied');
end

% check mask data
if ~isempty(cfg.maskparameter)
  interp.mask = double(reshape(getsubfield(ininterp,cfg.maskparameter), ininterp.dim));
  maskdat = 1;
else
  fprintf('no opacity mask data supplied\n');
  interp.mask = [];
  maskdat = 0;
end

% only work with the copy of the relevant parameters in "interp"
clear ininterp;

% convert anatomy data type and optimize contrast
if isa(interp.anatomy, 'uint8') || isa(interp.anatomy, 'uint16')
  fprintf('converting anatomy to floating point values...');
  interp.anatomy = double(interp.anatomy);
  fprintf('done\n');
end
fprintf('optimizing contrast of anatomical data ...');
minana = min(interp.anatomy(:));
maxana = max(interp.anatomy(:));
interp.anatomy = (interp.anatomy-minana)./(maxana-minana);
fprintf('done\n');

% store original data if 'interactive' mode
if strcmp(cfg.interactive,'yes')
  data.source = interp.source;
end

% place markers
marker = zeros(size(interp.anatomy));
if ~isempty(cfg.marker)
  fprintf('placing markers ...');
  [x,y,z] = ndgrid([1:size(interp.anatomy,1)],[1:size(interp.anatomy,2)],[1:size(interp.anatomy,3)]);
  for imarker = 1:size(cfg.marker,1)
    marker(find(sqrt((x-cfg.marker(iarker,1)).^2 + (y-cfg.marker(imarker,2)).^2 + (z-cfg.marker(imarker,3)).^2)<=cfg.markersize)) = 1;
  end
  fprintf('done\n');
end

% shift dimensions
fprintf('sorting dimensions...');
interp.anatomy = shiftdim(interp.anatomy,cfg.dim-1);
interp.source = shiftdim(interp.source,cfg.dim-1);
interp.mask = shiftdim(interp.mask,cfg.dim-1);
marker = shiftdim(marker,cfg.dim-1);
fprintf('done\n');

% flip dimensions
if strcmp(cfg.flipdim,'yes')
  fprintf('flipping dimensions...');
  interp.anatomy = flipdim(interp.anatomy,1);
  interp.source = flipdim(interp.source,1);
  interp.mask = flipdim(interp.mask,1);
  marker = flipdim(marker,1);
  fprintf('done\n');
end

% set slice space
if ischar(cfg.spacemin)
  fprintf('setting first slice position...');
  spacemin = min(find(~isnan(max(max(interp.source,[],3),[],2))));
  fprintf('%d...done\n',spacemin);
else
  spacemin = cfg.spacemin;
end

if ischar(cfg.spacemax)
  fprintf('setting last slice position...');
  spacemax = max(find(~isnan(max(max(interp.source,[],3),[],2))));
  fprintf('%d...done\n',spacemax);
else
  spacemax = cfg.spacemax;
end

% clip funtional data
if ~ischar(cfg.clipmin)
  fprintf('clipping functional minimum...');
  switch cfg.clipsym
  case 'no'
    interp.source(find(interp.source<cfg.clipmin)) = nan;
  case 'yes'
    interp.source(find(abs(interp.source)<cfg.clipmin)) = nan;
  end
  fprintf('done\n');
end
if ~ischar(cfg.clipmax)
  fprintf('clipping functional maximum...');
  switch cfg.clipsym
  case 'no'
    interp.source(find(interp.source>cfg.clipmax)) = nan;
  case 'yes'
    interp.source(find(abs(interp.source)>cfg.clipmax)) = nan;
  end
  fprintf('done\n');
end

% clip mask data
if maskdat
  if ~ischar(cfg.maskclipmin)
    fprintf('clipping mask minimum...');
    switch cfg.maskclipsym
    case 'no'
      interp.mask(find(interp.mask<cfg.maskclipmin)) = nan;
    case 'yes'
      interp.mask(find(abs(interp.mask)<cfg.maskclipmin)) = nan;
    end
    fprintf('done\n');
  end
  if ~ischar(cfg.maskclipmax)
    fprintf('clipping mask maximum...');
    switch cfg.maskclipsym
    case 'no'
      interp.mask(find(interp.mask>cfg.maskclipmax)) = nan;
    case 'yes'
      interp.mask(find(abs(interp.mask)>cfg.maskclipmax)) = nan;
    end
    fprintf('done\n');
  end
end

% scale functional data
fprintf('scaling functional data...');
fmin = min(interp.source(:));
fmax = max(interp.source(:));
if ~ischar(cfg.colmin)
  fcolmin = cfg.colmin;
else
  if sign(fmin)==sign(fmax)
    fcolmin = fmin;
  else
    fcolmin = -max(abs([fmin,fmax]));
  end
end
if ~ischar(cfg.colmax)
  fcolmax = cfg.colmax;
else
  if sign(fmin)==sign(fmax)
    fcolmax = fmax;
  else
    fcolmax = max(abs([fmin,fmax]));
  end
end
interp.source = (interp.source-fcolmin)./(fcolmax-fcolmin);
if ~ischar(cfg.colmax)
  interp.source(find(interp.source>1)) = 1;
end
if ~ischar(cfg.colmin)
  interp.source(find(interp.source<0)) = 0;
end
fprintf('done\n');

% scale mask data
if maskdat
  fprintf('scaling mask data...');
  fmin = min(interp.mask(:));
  fmax = max(interp.mask(:));
  if ~ischar(cfg.maskcolmin)
    mcolmin = cfg.maskcolmin;
  else
    if sign(fmin)==sign(fmax)
      mcolmin = fmin;
    else
      mcolmin = -max(abs([fmin,fmax]));
    end
  end
  if ~ischar(cfg.maskcolmax)
    mcolmax = cfg.maskcolmax;
  else
    if sign(fmin)==sign(fmax)
      mcolmax = fmax;
    else
      mcolmax = max(abs([fmin,fmax]));
    end
  end
  interp.mask = (interp.mask-mcolmin)./(mcolmax-mcolmin);
  if ~ischar(cfg.maskcolmax)
    interp.mask(find(interp.mask>1)) = 1;
  end
  if ~ischar(cfg.maskcolmin)
    interp.mask(find(interp.mask<0)) = 0;
  end
  fprintf('done\n');
end

% merge anatomy, functional data and mask
fprintf('constructing overlay...');
if ischar(cfg.colormap)
  % replace string by colormap using standard MATLAB function
  cfg.colormap = colormap(cfg.colormap);
end
cmap = cfg.colormap;
cmaplength = size(cmap,1);
maskmap = cfg.maskmap(:);
maskmaplength = size(maskmap,1);
indslice = round(linspace(spacemin,spacemax,cfg.nslices));
nvox1 = length(1:cfg.resample:size(interp.anatomy,2));
nvox2 = length(1:cfg.resample:size(interp.anatomy,3));
if mod(cfg.rotate,2)
  dummy = nvox1;
  nvox1 = nvox2;
  nvox2 = dummy;
end
out = zeros(nvox1,nvox2,3,cfg.nslices);
for islice = 1:cfg.nslices
  sel1 = 1:cfg.resample:size(interp.anatomy,2);
  sel2 = 1:cfg.resample:size(interp.anatomy,3);

  dummy1 = reshape(interp.anatomy(indslice(islice),sel1,sel2), [numel(sel1) numel(sel2)]);
  dummy2 = reshape(interp.source(indslice(islice),sel1,sel2),  [numel(sel1) numel(sel2)]);
  indmarker = find(reshape(marker(indslice(islice),sel1,sel2), [numel(sel1) numel(sel2)]));
  indsource = find(~isnan(dummy2));
  if maskdat
    dummymask = reshape(interp.mask(indslice(islice),sel1,sel2), [numel(sel1) numel(sel2)]);
    indsource = find(~isnan(dummy2) & ~isnan(dummymask));
  end
  for icol = 1:3
    dummy3 = dummy1;
    if not(maskdat)
      if ~ischar(cfg.alpha)
        try
          dummy3(indsource) = ...
            (1-cfg.alpha) * dummy3(indsource) + ...
            cfg.alpha * cmap(round(dummy2(indsource)*(cmaplength-1))+1,icol);
        end
      else
        try
          dummy3(indsource) = ...
            (1-dummy2(indsource)) .* dummy3(indsource) + ...
            dummy2(indsource) .* cmap(round(dummy2(indsource)*(cmaplength-1))+1,icol);
        end
      end
    else
      dummy3(indsource) = ...
        (1-maskmap(round(dummymask(indsource)*(maskmaplength-1))+1)).* ...
        dummy3(indsource) + ...
        maskmap(round(dummymask(indsource)*(maskmaplength-1))+1) .* ...
        cmap(round(dummy2(indsource)*(cmaplength-1))+1,icol);
    end
    dummy3(indmarker) = cfg.markercolor(icol);
    out(:,:,icol,islice) = rot90(dummy3,cfg.rotate);
  end
  if strcmp(cfg.whitebg,'yes')
    bgmask = zeros(nvox1,nvox2);
    bgmask(find(conv2(mean(out(:,:,:,islice),3),ones(round((nvox1+nvox2)/8))/(round((nvox1+nvox2)/8).^2),'same')<0.1)) = 1;
    for icol = 1:3
      out(:,:,icol,islice) = bgmask.*ones(nvox1,nvox2) + (1-bgmask).* out(:,:,icol,islice);
    end
  end
end
fprintf('done\n');

clf;
fprintf('plotting...');
axes('position',[0.9 0.3 0.02 0.4]);
image(permute(cmap,[1 3 2]));
set(gca,'YAxisLocation','right');
set(gca,'XTick',[]);
set(gca,'YDir','normal');
set(gca,'YTick',linspace(1,cmaplength,5));
set(gca,'YTickLabel',linspace(fcolmin,fcolmax,5));
set(gca,'Box','on');
axes('position',[0.01 0.01 0.88 0.90]);
[h,nrows,ncols]=slicemon(out);
xlim=get(gca,'XLim');
ylim=get(gca,'YLim');
text(diff(xlim)/2,-diff(ylim)/100,cfg.title,'HorizontalAlignment','center','Interpreter','none');
drawnow;
fprintf('done\n');
if nargout > 0
  outim=get(h,'CData');
end

if strcmp(cfg.interactive,'yes')
  data.sin = size(interp.source);
  data.nrows = nrows;
  data.ncols = ncols;
  data.out = out;
  data.indslice = indslice;
  data.cfg = cfg;
  data.hfig = gcf;
  uicontrol('Units','norm', 'Position', [0.9 0.2 0.08 0.05], 'Style','pushbutton', 'String','coords',...
    'Callback',@getcoords,'FontSize',7);
  data.hcoords = uicontrol('Units','norm', 'Position', [0.9 0.05 0.08 0.13], 'Style','text', 'String','','HorizontalAlign','left','FontSize',7);
  guidata(data.hfig,data);
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble history ininterp
ft_postamble provenance

% ---------------- subfunctions ----------------

function getcoords(h,eventdata,handles,varargin)
data = guidata(gcf);
[xi,yi] = ginput(1);

co(2,1) = round(mod(yi,size(data.out,1)));
co(3,1) = round(mod(xi,size(data.out,2)));
switch mod(data.cfg.rotate,4)
case 1
  t1 = co(2);
  co(2) = co(3);
  co(3) = data.sin(3)-t1;
case 2
  co(2) = data.sin(2)-co(2);
  co(3) = data.sin(3)-co(3);
case 3
  t1 = co(3);
  co(3) = co(2);
  co(2) = data.sin(2)-t1;
end

try
  co(1) = data.indslice(fix(xi/size(data.out,2)) + fix(yi/size(data.out,1))*data.ncols + 1);
catch
  co(1) = NaN;
end

if strcmp(data.cfg.flipdim, 'yes')
  co(1) = data.sin(1) - co(1) + 1;
end
co = co(:);

co(2:3) = round(co(2:3)*data.cfg.resample);
for ishift = 1:data.cfg.dim-1
 co = [co(3);co(1);co(2)];
end
set(data.hcoords,'String',sprintf('1: %d\n2: %d\n3: %d\nf: %0.4f',co(1),co(2),co(3),data.source(co(1),co(2),co(3))));

function [h,nrows,ncols] = slicemon(a) % display the montage w/o image_toolbox
siz = [size(a,1) size(a,2) size(a,4)];
nn = sqrt(prod(siz))/siz(2);
mm = siz(3)/nn;
if (ceil(nn)-nn) < (ceil(mm)-mm)
  nn = ceil(nn); mm = ceil(siz(3)/nn);
else
  mm = ceil(mm); nn = ceil(siz(3)/mm);
end
b = a(1,1);
b(1,1) = 0;
b = repmat(b, [mm*siz(1), nn*siz(2), size(a,3), 1]);
rows = 1:siz(1); cols = 1:siz(2);
for i=0:mm-1
  for j=0:nn-1
    k = j+i*nn+1;
    if k<=siz(3)
      b(rows+i*siz(1),cols+j*siz(2),:) = a(:,:,:,k);
    end
  end
end
hh = image(b);
axis image;
box off;
set(gca,'XTick',[],'YTick',[],'Visible','off');
if nargout > 0
  h = hh;
  nrows = mm;
  ncols = nn;
end

author: spm


text:
function s = tifs2seq(file)
%TIFS2SEQ Create a MATLAB sequence from a multi-frame TIFF file.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Get the number of frames in the multi-frame TIFF.
frames = size(imfinfo(file), 1);

% Read the first frame, preallocate the sequence, and put the first
% in it.
i = imread(file, 1);
s = zeros([size(i) 1 frames], 'uint8');
s(:,:,:,1) = i;

% Read the remaining TIFF frames and add to the sequence.
for i = 2:frames
    s(:,:,:,i) = imread(file, i);
end

author: dipum

repo: dipum-toolbox

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


text:
function [wl,wr] = cvwindow(x)

  wr = zeros(size(x));
  wl = zeros(size(x));

  eps = 1e-16;
  sml = find(x<eps); %too small
  mid = find(x>=eps & x<=1-eps); %just right
  lrg = find(x>1-eps); %too large

  wl(sml) = 0;  wr(sml) = 1;
  wl(lrg) = 1;  wr(lrg) = 0;

  xmid = x(mid);
  a = exp(1-1./(1-exp(1-1./(1-xmid))));
  b = exp(1-1./(1-exp(1-1./xmid)));
  n = sqrt(a.^2 + b.^2);
  wl(mid) = a./n;  wr(mid) = b./n;

author: thomaskuestner


text:
function disp_(u,str)
%DISP_        Display of slope variables in "_" notation
%
%  disp_(u)
%

%Second parameter name for internal purposes: display full structure
%

% written  06/04/09     S.M. Rump
% modified 02/28/10     S.M. Rump  vector output and multi-dimensional arrays
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  loose = strcmp(get(0,'FormatSpacing'),'loose');

  numvar = size(u.s,2);
  if numvar~=INTLAB_SLOPE.NUMVAR
    warning('**** number of dependent variables in slope expansion inconsitent')
  end

  if nargin==2
    if isequal(str,'internal')
      disp([ 'slope range (internal structure) ' inputname(1) '.r = ' ])
      display(u.r,'',1);
      disp([ 'slope slope ' inputname(1) '.s = ' ])
      display(u.s,'',1);
      return
    else
      error('Invalid call of slope disp_')
    end
  end

  if loose, disp(' '); end
  disp([ 'slope intval center ' inputname(1) '.c = ' ])
  if loose, disp(' '); end
  ur = u.r(:,1);
  disp_( reshape(ur,u.size) , '' , 1 )
  if loose, disp(' '); end

  if loose, disp(' '); end
  disp([ 'slope intval range ' inputname(1) '.r = ' ])
  if loose, disp(' '); end
  ur = u.r(:,INTLAB_SLOPE.NUMVAR+1);
  disp_( reshape(ur,u.size) , '' , 1 )
  if loose, disp(' '); end

  if loose, disp(' '); end
  disp([ 'slope intval slope ' inputname(1) '.s = ' ])
  if loose, disp(' '); end
  if ( length(u.size)==2 ) & ( u.size(2)==1 )
    disp_( reshape(u.s,[u.size(1) numvar]) , '' , 1 )
  else
    disp_( reshape(u.s,[u.size numvar]) , '' , 1 )
  end
  if loose, disp(' '); end

repo: OpenMAS


text:
% This function compute the modulation spectrum of speech. The modulation
% spectrum can be derived from several types of feature trajectories: 1)
% log spectrogram; 2) log Mel filterbanks; 3) MFCC. Note that it is
% different from another definition of modulation spectrum that derived
% from the evelop of band-passed time domain signal.
% Inputs:
%   wav: the time domain signal
%   para: a structure that controls various options in the function
% Outputs:
%   modu: the modulation spectrum (complex valued) of feature blocks
%   moduAvg: the average modulation spectrum (magnitude) of the input signal
%
% Author: Xiao Xiong
% Created: 1 Dec 2014
% Last Modified: 5 Dec 2014
%
function [modu, moduAvg] = compute_modulation_stft(wav, para)
type = para.type;
fs = para.fs;
win_size = para.win_size;
win_shift = para.win_shift;

spec = wav2abs(wav, fs);

if type == 1    % compute from spectrum
    data = log(spec);
elseif type == 2    % compute from Mel filterbanks
    data = log(abs2Mel(spec, fs));
elseif type == 3    % compute from static MFCC
    data = fbank2mfcc(log(abs2Mel(spec, fs)));
end

% Perform optional normalization to the data
if isfield(para, 'norm')
    switch para.norm
        case 'CMN'
            data = CMN(data);
        case 'MVN'
            data = MVN(data);
        otherwise
            fprintf('Unknown normalization %s\n', para.norm);
    end
end

[nFr, dim] = size(data);

nBlockMax = ceil( (nFr-win_size)/win_shift ) + 1;
modu = []; moduAvg = [];
nBlock = 0;
for i=1:win_shift:nFr
    j = min(nFr, i+win_size-1);
    if j-i+1 < win_size/2   % if the remaining frames is not enough for half of the window size, discard it.
        break;
    end
    nBlock = nBlock + 1;
    modu(:,:,nBlock) = fft(data(i:j,:), win_size);
    last_block_size = (j-i)+1;
end

if size(modu,3)==0
    return;
end

modu = modu(1:win_size/2+1,:,:);

moduAvg = sum(abs(modu(:,:,1:end-1)),3) + abs(modu(:,:,end)) * last_block_size/win_size;
moduAvg = moduAvg / (size(modu,3)-1+last_block_size/win_size);

end

author: singaxiong

repo: SignalGraph


text:
function [Dt]=patchVectorTangent(F,V,D,N)

%Derive tangent contribution through dot-product with normal vectors

if isempty(N)
    [~,~,N]=patchNormal(F,V); %Get current vertex normals
end

Dn_mag=dot(D,N,2); %Allong normal displacement magnitudes
Dn=Dn_mag(:,ones(1,3)).*N; %Normal direction displacement vectors
Dt=D-Dn; %Tranverse or tangential only displacement vectors

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [vertices, faces] = triangulatePolygonPair(poly1, poly2, varargin)
%TRIANGULATEPOLYGONPAIR Compute triangulation between a pair of 3D closed curves.
%
%   [V, F] = triangulatePolygonPair(POLY1, POLY2)
%
%   [V, F] = triangulatePolygonPair(..., 'recenter', FLAG)
%   Where FLAG is a boolean, specifies whether the second curve should be
%   translated to have the same centroid as the first curve. This can
%   improve mathcing of vertices. Default is true.
%
%
%   Example
%     % triangulate a surface patch between two ellipses
%     % create two sample curves
%     poly1 = ellipseToPolygon([50 50 40 20 0], 36);
%     poly2 = ellipseToPolygon([50 50 40 20 60], 36);
%     poly1 = poly1(1:end-1,:);
%     poly2 = poly2(1:end-1,:);
%     % transform to 3D polygons / curves
%     curve1 = [poly1 10*ones(size(poly1, 1), 1)];
%     curve2 = [poly2 20*ones(size(poly2, 1), 1)];
%     % draw as 3D curves
%     figure(1); clf; hold on;
%     drawPolygon3d(curve1, 'b'); drawPoint3d(curve1, 'bo');
%     drawPolygon3d(curve2, 'g'); drawPoint3d(curve2, 'go');
%     view(3); axis equal;
%     [vertices, faces] = triangulatePolygonPair(curve1, curve2);
%     % display the resulting mesh
%     figure(2); clf; hold on;
%     drawMesh(vertices, faces);
%     drawPolygon3d(curve1, 'color', 'b', 'linewidth', 2);
%     drawPolygon3d(curve2, 'color', 'g', 'linewidth', 2);
%     view(3); axis equal;
%
%   See also
%     meshes3D, triangulatePolygonPair3d, triangulateCurvePair,
%     meshSurfaceArea
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2017-05-18, using Matlab 9.1.0.441655 (R2016b)
% Copyright 2017-2022 INRA - Cepia Software Platform

%% Settings

recenterFlag = true;
while length(varargin) > 1
    pname = varargin{1};
    if strcmpi(pname, 'recenter')
        recenterFlag = varargin{2};
    else
        error('Unknown parameter name: %s', pname);
    end
    varargin(1:2) = [];
end

%% Memory allocation

% concatenate vertex coordinates for creating mesh
vertices = [poly1 ; poly2];

% number of vertices on each polygon
n1 = size(poly1, 1);
n2 = size(poly2, 1);

% allocate the array of facets (each edge of each polygon provides a facet)
nFaces = n1 + n2;
faces = zeros(nFaces, 3);

% Translate the second polygon such that the centroids of the bounding
% boxes coincide. This is expected to improve the matching of the two
% curves.
if recenterFlag
    box1 = boundingBox3d(poly1);
    box2 = boundingBox3d(poly2);
    center1 = (box1(2:2:end) + box1(1:2:end-1)) / 2;
    center2 = (box2(2:2:end) + box2(1:2:end-1)) / 2;
    vecTrans = center1 - center2;
    trans = createTranslation3d(vecTrans);
    poly2 = transformPoint3d(poly2, trans);
end

%% Init iteration

% find the pair of points with smallest distance.
% This will be the current diagonal.
[dists, inds] = minDistancePoints(poly1, poly2);
[dummy, ind1] = min(dists); %#ok<ASGLU>
ind2 = inds(ind1);

% consider two consecutive vertices on each polygon
currentIndex1 = ind1;
currentIndex2 = ind2;

%% Main iteration
% For each diagonal, consider the two possible facets (one for each 'next'
% vertex on each polygon), each create current facet according to the
% closest one.
% Then update current diagonal for next iteration.

for iFace = 1:nFaces
    nextIndex1 = mod(currentIndex1, n1) + 1;
    nextIndex2 = mod(currentIndex2, n2) + 1;

    % compute lengths of diagonals
    dist1 = distancePoints(poly1(currentIndex1, :), poly2(nextIndex2,:));
    dist2 = distancePoints(poly1(nextIndex1, :), poly2(currentIndex2,:));

    if dist1 < dist2
        % keep current vertex of curve1, use next vertex on curve2
        face = [currentIndex1 currentIndex2+n1 nextIndex2+n1];
        currentIndex2 = nextIndex2;
    else
        % keep current vertex of curve2, use next vertex on curve1
        face = [currentIndex1 currentIndex2+n1 nextIndex1];
        currentIndex1 = nextIndex1;
    end

    % create the facet
    faces(iFace, :) = face;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ average, sd, covc ] = covariance ( dim_num, point_num, x )

%*****************************************************************************80
%
%% COVARIANCE does a covariance calculation for IHS solutions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points to be generated.
%
%    Input, integer X(DIM_NUM,POINT_NUM), the points.
%
%    Output, real AVERAGE, the average minimum distance.
%
%    Output, real SD, the standard deviation of the minimum distances.
%
%    Output, real COVC, the covariance of the minimum distances.
%

%
%  Set up the distance matrix.
%
  for i = 1 : point_num
    mindist(i) = realmax;
    for j = 1 : point_num
      if ( i ~= j )
        vec(1:dim_num) = x(1:dim_num,i) - x(1:dim_num,j);
        dist = sqrt ( vec(1:dim_num) * vec(1:dim_num)' );
        if ( dist < mindist(i) )
          mindist(i) = dist;
        end
      end
    end
  end
%
%  Find the average minimum distance.
%
  average = sum ( mindist(1:point_num) ) / point_num;
%
%  Compute the standard deviation of the distances.
%
  sd = std ( mindist(1:point_num) );
%
%  Compute the covariance.
%
  covc = sd / average;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result = transVertexToNormalAxisBase(mesh_outer, iV, coord)

% Originally written by Lucas Tamarit, tamarit@cisa.unige.ch
% Adapted to the local GI whole computation on 2007/15/11
%
% create a new referential based on the normal of the outer surface at the
% center of the outer region of interest. In order to give the distance of
% a given point (coord) the main axis of this referential.

B.origin = mesh_outer.vertices(iV,:); % set the origin of the referential
n = mesh_outer.vertexNormal(iV,:); % use the normal to the outer smoothed surface at that point

B.axe = n ./ norm(n);
B.u = getOrthogonalVector(B.axe);
B.v = cross(B.axe,B.u);
B.v = B.v ./ norm(B.v);

vt = coord - B.origin;
vt = [B.u; B.v; B.axe]*vt';

result = vt;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function outmap = read_map(filename)
% Reads a map from a file.  The file must contain two columns with
% keys in the first column and values in the second column.
% Inputs:
%   filename: The name of the text file containing the map.
% Outputs:
%   outmap: A struct with a cell array of keys and a cell array of
%   values.  The cell arrays will be the same length and will
%   contain strings.

assert(isa(filename,'char'))

fid = fopen(filename,'r');
cc = textscan(fid,'%s%s');
fclose(fid);

outmap.keySet = cc{1};
outmap.values = cc{2};

author: nesl


text:
function [beta_gibbs F_gibbs L_gibbs phi_gibbs sigma_gibbs lambda_t_gibbs sigma_t_gibbs sbar]=stvol3gibbs(Xbart,Xt,yt,B0,phi0,alpha0,delta0,f0,upsilon0,betahat,sigmahat,gamma,G,I_o,omega,T,n,k,It,Bu,pick,pickf)

% preliminary elements for the algorithm
% compute the product G'*I_gamma*G (to speed up computations of deltabar)
GIG=G'*I_o*G;
% compute alphabar
alphabar=T+alpha0;

% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices and cells
beta_gibbs=[];
F_gibbs=[];
L_gibbs=[];
phi_gibbs=[];
sigma_gibbs=[];
lambda_t_gibbs={};
sigma_t_gibbs={};

% step 1: determine initial values for the algorithm

% initial value for beta
beta=betahat;
B=reshape(beta,k,n);
% initial value for f_2,...,f_n
% obtain the triangular factorisation of sigmahat
[Fhat Lambdahat]=bear.triangf(sigmahat);
% obtain the initial value for F
F=Fhat;
% obtain the inverse of Fhat
[invFhat]=bear.invltod(Fhat,n);
% create the cell storing the different vectors of invF
Finv=cell(n,1);
% store the vectors
for ii=2:n
Finv{ii,1}=invFhat(ii,1:ii-1);
end
% initial values for L
L=zeros(T,1);
% initial values for phi
phi=1;

% step 2: determine the sbar values and Lambda
sbar=diag(Lambdahat);
Lambda=sparse(diag(sbar));
% then determine sigma^(0)
sigma=F*Lambda*F';

% step 3: recover the series of initial values for lambda_1,...,lambda_T and sigma_1,...,sigma_T
lambda_t=repmat(diag(sbar),1,1,T);
sigma_t=repmat(sigmahat,1,1,T);

hbar = bear.parfor_progressbar(It,'Progress of the Gibbs sampler');  %create the progress bar

% run the Gibbs sampler
while count<=It

   hbar.iterate(1);   % update progress by one iteration

% step 4: draw beta from its conditional posterior
% first compute the summations required for omegabar and betabar
summ1=zeros(k,k);
summ2=zeros(k,n);
   % run the summation
   for jj=1:T
   prodt=Xt{jj,1}'*exp(-L(jj,1));
   summ1=summ1+prodt*Xt{jj,1};
   summ2=summ2+prodt*yt(:,:,jj)';
   end
% then obtain the inverse of phi0
invphi0=diag(1./diag(phi0));
% obtain the inverse of phibar
invphibar=summ1+invphi0;
% recover phibar
C=chol(bear.nspd(invphibar),'Lower')';
invC=C\speye(k);
phibar=invC*invC';
% recover Bbar
Bbar=phibar*(summ2+invphi0*B0);
% draw B from its posterior
B=bear.matrixndraw(Bbar,sigma,phibar,k,n);
% finally recover beta by vectorising
beta=B(:);

% step 5: draw the series f_2,...,f_n from their conditional posteriors
   % recover first the residuals
   for jj=1:T
   epst(:,:,jj)=yt(:,:,jj)-Xbart{jj,1}*beta;
   end
   % then draw the vectors in turn
   for jj=2:n
   % first compute the summations required for upsilonbar and fbar
   summ1=zeros(jj-1,jj-1);
   summ2=zeros(jj-1,1);
      % run the summation
      for kk=1:T
      prodt=epst(1:jj-1,1,kk)*exp(-L(kk,1));
      summ1=summ1+prodt*epst(1:jj-1,1,kk)';
      summ2=summ2+prodt*epst(jj,1,kk)';
      end
   summ1=(1/sbar(jj,1))*summ1;
   summ2=(-1/sbar(jj,1))*summ2;
   % then obtain the inverse of upsilon0
   invupsilon0=diag(1./diag(upsilon0{jj,1}));
   % obtain upsilonbar
   invupsilonbar=summ1+invupsilon0;
   C=chol(bear.nspd(invupsilonbar));
   invC=C\speye(jj-1);
   upsilonbar=full(invC*invC');
   % recover fbar
   fbar=upsilonbar*(summ2+invupsilon0*f0{jj,1});
   % finally draw f_i^(-1)
   Finv{jj,1}=fbar+chol(bear.nspd(upsilonbar),'lower')*randn(jj-1,1);
   end
% recover the inverse of F
invF=eye(n);
   for jj=2:n
   invF(jj,1:jj-1)=Finv{jj,1};
   end
% eventually recover F
F=bear.invltod(invF,n);
% update sigma
sigma=F*Lambda*F';

% step 6: draw phi from its conditional posterior
% estimate deltabar
deltabar=L'*GIG*L+delta0;
% draw the value phi_i
phi=bear.igrandn(alphabar/2,deltabar/2);

% step 7: draw the series lambda_t from their conditional posteriors, t=1,...,T
   % consider periods in turn
   for kk=1:T
   % a candidate value will be drawn from N(lambdabar,phibar)
   % the definitions of lambdabar and phibar varies with the period, thus define them first
      % if the period is the first period
      if kk==1
      lambdabar=(gamma*L(2,1))/(1/omega+gamma^2);
      phibar=phi/(1/omega+gamma^2);
      % if the period is the final period
      elseif kk==T
      lambdabar=gamma*L(T-1,1);
      phibar=phi;
      % if the period is any period in-between
      else
      lambdabar=(gamma/(1+gamma^2))*(L(kk-1,1)+L(kk+1,1));
      phibar=phi/(1+gamma^2);
      end
   % now draw the candidate
   cand=lambdabar+phibar^0.5*randn;
   % compute the acceptance probability
   prob=bear.mhprob3(cand,L(kk,1),sbar,epst(:,1,kk),Finv,n);
   % draw a uniform random number
   draw=rand;
      % keep the candidate if the draw value is lower than the prob
      if draw<=prob
      L(kk,1)=cand;
      % if not, just keep the former value
      end
   end
% then recover the series of matrices lambda_t and sigma_t
for kk=1:T
lambda_t(:,:,kk)=exp(L(kk,1))*diag(sbar);
sigma_t(:,:,kk)=F*lambda_t(:,:,kk)*F';
end

% record phase
   % if the burn-in sample phase is not yet over
   if count<=Bu
   % simply add 1 to the iteration count
   count=count+1;
   % on the other hand, if the burn-in sample phase is over
   elseif count>Bu
   % adding one iteration to the count will depend on wether post-burn selection applies
      % if there is no post burn selection
      if pick==0
      % record the results
      beta_gibbs(:,count-Bu)=beta;
      F_gibbs(:,:,count-Bu)=F;
      L_gibbs(:,:,count-Bu)=L;
      phi_gibbs(count-Bu,1)=phi;
      sigma_gibbs(:,count-Bu)=sigma(:);
         for jj=1:T
         lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
         sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
         end
      % then add one to the count
      count=count+1;
      % if there is post burn selection, only one draw over 'fpick' draws will be retained
      elseif pick==1
         % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1, and do record the acceptance rate of the Metropolis-Hastings step
         if pickcount~=pickf
         pickcount=pickcount+1;
         % on the other hand, if the iteration does correspond to fpick
         elseif pickcount==pickf
         % do record the results
         beta_gibbs(:,count-Bu)=beta;
         F_gibbs(:,:,count-Bu)=F;
         L_gibbs(:,:,count-Bu)=L;
         phi_gibbs(count-Bu,1)=phi;
         sigma_gibbs(:,count-Bu)=sigma(:);
            for jj=1:T
            lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
            sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
            end
         % then increase the regular count by 1 and re-initialise pickcount
         count=count+1;
         pickcount=1;
         end
      end
   end

end

close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
function gasConst = getIdealGasConstant()
    gasConst = 8.31447; %m^3*Pa/(K*mol)
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = quiver3(F, varargin)
%QUIVER3   3-D quiver plot of a CHEBFUN3V object.
%   QUIVER3(F) plots velocity vectors as arrows with components F(1), F(2)
%   and F(3) which are CHEBFUN3 objects. QUIVER3 automatically scales the
%   arrows to fit. The arrows are plotted on a uniform grid.
%
%   QUIVER3(F, S) automatically scales the arrows to fit and then stretches
%   them by S. Use S=0 to plot the arrows with the automatic scaling.
%
%   QUIVER3(..., LINESPEC) uses the plot linestyle specified for the
%   velocity vectors.  Any marker in LINESPEC is drawn at the base instead
%   of an arrow on the tip. Use a marker of '.' to specify no marker at
%   all.
%
% See also quiver3.

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

numpts = 7;

% Empty check:
if ( isempty(F) )
    h = quiver3([], [], [], []);
    return
end

if ( isa(F, 'chebfun3v') && ...
        (isempty(varargin) || ~isa(varargin{1}, 'chebfun3v')) )

    nF = F.nComponents;
    if ( nF == 2 )
        error('CHEBFUN:CHEBFUN3V:quiver3:badInputs1',...
            'Needs three CHEBFUN3 objects.');
    else
        % Plot arrows are equally spaced locations
        dom = F.components{1}.domain;
        x = linspace(dom(1), dom(2), numpts);
        y = linspace(dom(3), dom(4), numpts);
        z = linspace(dom(5), dom(6), numpts);
        [xx, yy, zz] = ndgrid(x, y, z);
        vals1 = feval(F.components{1}, xx, yy, zz);
        vals2 = feval(F.components{2}, xx, yy, zz);
        vals3 = feval(F.components{3}, xx, yy, zz);
        h = quiver3(xx, yy, zz, vals1, vals2, vals3, varargin{:});
        axis([x(1), x(end), y(1), y(end), z(1), z(end)])
    end

else
    error('CHEBFUN:CHEBFUN3V:quiver3:badInputs3',...
        'Unrecognised input arguments.');
end

if ( nargout > 0 )
    varargout = {h};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function msm_to_mm_array_complex_general ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_ARRAY_COMPLEX_GENERAL writes a "matrix array complex general" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_ARRAY_COMPLEX_GENERAL - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_ARRAY_COMPLEX_GENERAL - Fatal error!' );
  end

  fprintf ( fid, '%%%%MatrixMarket matrix array complex general\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_ARRAY_COMPLEX_GENERAL.M\n' );
  fprintf ( fid, '  %d  %d\n', nrow, ncol );

  for j = 1 : ncol
    for i = 1 : nrow
      fprintf ( fid, '  %f  %f\n', real ( a(i,j) ), imag ( a(i,j) ) );
    end
  end

  fclose ( fid );

  return
end


text:
%% trisurfsmooth
% Below is a demonstration of the features of the |trisurfsmooth| function

%%
clear; close all; clc;

%% Syntax
% |[Vs]=trisurfsmooth(F,V,n,IND_V,par,opt);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
RGB = imread('saturn.png');
I = rgb2gray(RGB);
J = imnoise(I,'gaussian', 0, 0.005);
K = wiener2(J, [10, 10]);
H = fspecial('disk', 10);
blurred = imfilter(J, H, 'replicate');
figure, imshow(RGB);
figure, imshow(I);
figure, imshow(J);
figure, imshow(H);
figure, imshow(blurred);
figure, imshow(K);

author: UtkarshPathrabe


text:
function colorspace_demo(Cmd)
% Demo for colorspace.m - 3D visualizations of various color spaces

% Pascal Getreuer 2006

if nargin == 0
   % Create a figure with a drop-down menu
   figure('Color',[1,1,1]);
   h = uicontrol('Style','popup','Position',[15,10,90,21],...
      'BackgroundColor',[1,1,1],'Value',2,...
      'string',{'sRGB','Y''PbPr','HSV','HSL','L*a*b*','L*u*v*','L*ch'});
   set(h,'Callback',sprintf('%s(get(%.20g,''Value''))',mfilename,h));
   Cmd = 2;
end

% Plot a figure
switch Cmd
case 1
   [x,y,z,Tri] = makeshape('Cube');
   CData = [x,y,z];
   myplot((x-0.5)*0.8,(y-0.5)*0.8,(z)*0.8,Tri,CData);
   coloraxis('R''',5,0.5*0.8);
   coloraxis('G''',6,0.5*0.8);
   coloraxis('B''',3);
case 2
   [x,y,z,Tri] = makeshape('Cylinder');
   CData = colorspace('rgb<-ypbpr',[z,x/2,y/2]);
   myplot(x,y,0.8*z,Tri,CData);
   coloraxis('Y''',3);
   coloraxis('P_b',5,0.8);
   coloraxis('P_r',6,0.8);
case 3
   [x,y,z,Tri] = makeshape('Hexacone');
   CData = colorspace('rgb<-hsv',[(pi+atan2(-y,-x))*180/pi,sqrt(x.^2+y.^2),z]);
   myplot(x,y,z,Tri,CData);
   coloraxis('H',1);
   coloraxis('S',2);
   coloraxis('V',3);
case 4
   [x,y,z,Tri] = makeshape('Double Hexacone');
   CData = colorspace('rgb<-hsl',[(pi+atan2(-y,-x))*180/pi,sqrt(x.^2+y.^2),z]);
   myplot(x,y,2*z,Tri,CData);
   coloraxis('H',1);
   coloraxis('S',2);
   coloraxis('L',4);
case 5
   [x,y,z,Tri] = makeshape('Blobs');
   CData = colorspace('rgb<-lab',[z*100,x*100,y*100]);
   myplot(x,y,2*z,Tri,CData);
   coloraxis('L*',4);
   coloraxis('a*',5,2);
   coloraxis('b*',6,2);
case 6
   [x,y,z,Tri] = makeshape('Blobs');
   CData = colorspace('rgb<-luv',[z*100,x*125,y*125]);
   myplot(x,y,2*z,Tri,CData);
   coloraxis('L*',4);
   coloraxis('u*',5,2);
   coloraxis('v*',6,2);
case 7
   [x,y,z,Tri] = makeshape('Blobs');
   CData = colorspace('rgb<-lab',[z*100,x*100,y*100]);
   myplot(x,y,2*z,Tri,CData);
   coloraxis('L*',4);
   coloraxis('c',2);
   coloraxis('h',1);
end

axis equal;
axis off;
pbaspect([1,1,1]);
view(70,27);
rotate3d on;

return;

function myplot(x,y,z,Tri,CData)
% Plot a triangular mesh with color data
cla;
CData = min(max(CData,0),1);
patch('Faces',Tri,'Vertices',[x,y,z],'FaceVertexCData',CData,...
   'FaceColor','interp','EdgeColor','none');
hold on;
return;

function coloraxis(Name,Type,h)
% Draw color axes as 3D objects
FontSize = 14;

switch Type
case 1
   set(text(-0.25,-1.3,1.1,Name),'FontWeight','bold',...
      'FontSize',FontSize,'Color',[0,0,0]);
   t = linspace(0,pi*3/2,60);
   x = cos(t)*1.1;
   y = sin(t)*1.1;
   set(plot3(x,y,zeros(size(x)),'k-'),'LineWidth',2.5,'Color',[1,1,1]*0.8);
   x = [x,-0.1,0,-0.1];
   y = [y,-1.05,-1.1,-1.15];
   set(plot3(x,y,ones(size(x)),'k-'),'LineWidth',2.5);
case 2
   set(text(0,-0.6,0.15,Name),'FontWeight','bold','FontSize',FontSize,'Color',[0,0,0]);
   set(plot3([-0.05,0,0.05,0,0],[-0.9,-1,-0.9,-1,0],[0,0,0,0,0],'k-'),'LineWidth',2.5);
case 3
   set(text(0,0.15,1.3,Name),'FontWeight','bold','FontSize',FontSize,'Color',[0,0,0]);
   set(plot3([0,0,0,0,0],[-0.05,0,0.05,0,0],[1.6,1.7,1.6,1.7,0],'k-'),'LineWidth',2.5);
case 4
   set(text(0,0.15,2.4,Name),'FontWeight','bold','FontSize',FontSize,'Color',[0,0,0]);
   set(plot3([0,0,0,0,0],[-0.05,0,0.05,0,0],[2.4,2.5,2.4,2.5,0],'k-'),'LineWidth',2.5);
case 5
   set(text(0.6,0,h+0.15,Name),'FontWeight','bold','FontSize',FontSize,'Color',[0,0,0]);
   set(plot3([0.9,1,0.9,1,-1,-0.9,-1,-0.9],[-0.05,0,0.05,0,0,-0.05,0,0.05],...
      zeros(1,8)+h,'k-'),'LineWidth',2.5);
case 6
   set(text(0,0.6,h+0.15,Name),'FontWeight','bold','FontSize',FontSize,'Color',[0,0,0]);
   set(plot3([-0.05,0,0.05,0,0,-0.05,0,0.05],[0.9,1,0.9,1,-1,-0.9,-1,-0.9],...
      zeros(1,8)+h,'k-'),'LineWidth',2.5);
end
return;

function [x,y,z,Tri] = makeshape(Shape)
% Make a 3D shape: Shape = 'Cube', 'Cylinder', 'Cone', 'Double Cone', 'Blobs'

% Cube
N = 12;    % Vertices per edge
% Cylinder, Cone, Double Cone
Nth = 25;  % Vertices over angles, (Nth - 1) should be a multiple of 12
Nr = 4;    % Vertices over radius
Nz = 8;    % Vertices over z-dimension

switch lower(Shape)
case 'cube'
   [u,v] = meshgrid(linspace(0,1,N),linspace(0,1,N));
   u = u(:);
   v = v(:);
   x = [u;u;u;u;zeros(N^2,1);ones(N^2,1)];
   y = [v;v;zeros(N^2,1);ones(N^2,1);v;v];
   z = [zeros(N^2,1);ones(N^2,1);v;v;u;u];
   Tri = trigrid(N,N);
   Tri = [Tri;N^2+Tri;2*N^2+Tri;3*N^2+Tri;4*N^2+Tri;5*N^2+Tri];
case 'cylinder'
   Nth = ceil(Nth*0.75);
   [u,v] = meshgrid(linspace(0,pi*3/2,Nth),linspace(0,1,Nz));
   Tri = trigrid(Nth,Nz);
   x = cos(u(:));
   y = sin(u(:));
   z = v(:);
   [u,v] = meshgrid(linspace(0,pi*3/2,Nth),linspace(0,1,Nr));
   Tri = [Tri;Nth*Nz+trigrid(Nth,Nr);Nth*Nz+Nth*Nr+trigrid(Nth,Nr)];
   x = [x;v(:).*cos(u(:));v(:).*cos(u(:))];
   y = [y;v(:).*sin(u(:));v(:).*sin(u(:))];
   z = [z;zeros(Nth*Nr,1);ones(Nth*Nr,1)];
   [u,v] = meshgrid(linspace(0,1,Nr),linspace(0,1,Nz));
   Tri = [Tri;Nth*Nz+2*Nth*Nr+trigrid(Nr,Nz);Nth*Nz+2*Nth*Nr+Nr*Nz+trigrid(Nr,Nz)];
   x = [x;u(:);zeros(Nr*Nz,1)];
   y = [y;zeros(Nr*Nz,1);-u(:)];
   z = [z;v(:);v(:)];
case 'cone'
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nz));
   Tri = trigrid(Nth,Nz);
   x = v(:).*cos(u(:));
   y = v(:).*sin(u(:));
   z = v(:);
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nr));
   Tri = [Tri;Nth*Nz+trigrid(Nth,Nr);];
   x = [x;v(:).*cos(u(:));];
   y = [y;v(:).*sin(u(:));];
   z = [z;ones(Nth*Nr,1)];
case 'double cone'
   Nz = floor(Nz/2)*2+1;
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nz));
   Tri = trigrid(Nth,Nz);
   r = 1 - abs(2*v(:) - 1);
   x = r.*cos(u(:));
   y = r.*sin(u(:));
   z = v(:);
case 'hexacone'
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nz));
   Tri = trigrid(Nth,Nz);
   r = 0.92./max(max(abs(cos(u)),abs(cos(u - pi/3))),abs(cos(u + pi/3)));
   x = v(:).*cos(u(:)).*r(:);
   y = v(:).*sin(u(:)).*r(:);
   z = v(:);
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nr));
   Tri = [Tri;Nth*Nz+trigrid(Nth,Nr);];
   v = 0.92*v./max(max(abs(cos(u)),abs(cos(u - pi/3))),abs(cos(u + pi/3)));
   x = [x;v(:).*cos(u(:));];
   y = [y;v(:).*sin(u(:));];
   z = [z;ones(Nth*Nr,1)];
case 'double hexacone'
   Nz = floor(Nz/2)*2+1;
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nz));
   Tri = trigrid(Nth,Nz);
   r = 1 - abs(2*v - 1);
   r = 0.92*r./max(max(abs(cos(u)),abs(cos(u - pi/3))),abs(cos(u + pi/3)));
   x = r(:).*cos(u(:));
   y = r(:).*sin(u(:));
   z = v(:);
case 'blobs'
   Nz = 47;
   [u,v] = meshgrid(linspace(0,2*pi,Nth),linspace(0,1,Nz));
   Tri = trigrid(Nth,Nz);
   r = sin(v(:)*pi*3).^2.*(1 - 0.6*abs(2*v(:) - 1));
   x = r.*cos(u(:));
   y = r.*sin(u(:));
   z = v(:);
end
return;

function Tri = trigrid(Nu,Nv)
% Construct the connectivity data for a grid of triangular patches
i = (1:(Nu-1)*Nv-1).';
i(Nv:Nv:end) = [];
Tri = [i,i+1,i+Nv;i+1,i+Nv+1,i+Nv];
return;

author: taoshiqian


text:
function y = vecsym(x,K) %#ok
% y = vecsym(x,K)
%
% VECSYM  For the PSD submatrices, we let Yk = (Xk+Xk')/2
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [y]=rake_fft(x, tol)
% Discrete Fourier Transform in the QTT-Tucker format
%   function [y]=rake_fft(x, tol)
%   Computes the multidimensional DFT of the vector x in the QTT-Tucker format with
%   the approximation tolerance tol.
%   Returns the quantized dimensions in the correct order.
%
%   !!WARNING!! employs a (probably) costly rank reverse procedure.

D = size(x.tuck);
tuck = x.tuck;
crx = core2cell(x.core);
for i=1:D
    [tuck{i}, rv]=qr(tuck{i}, 'lr');
    r1 = size(crx{i}, 1); r2 = size(crx{i}, 3); rt = size(crx{i},2);
    crx{i} = permute(crx{i}, [2, 1, 3]);
    crx{i} = reshape(crx{i}, rt, r1*r2);
    crx{i} = rv*crx{i};
    rt = size(rv,1);
    crx{i} = reshape(crx{i}, rt, r1, r2);
    crx{i} = permute(crx{i}, [2, 1, 3]);
end;
% [crx,rv]=qr(crx, 'rl');
% crx = rv.'*crx;
for i=D:-1:2
    r1 = size(crx{i}, 1); rt = size(crx{i}, 2); r2 = size(crx{i}, 3);
    crx{i} = reshape(crx{i}, r1, rt*r2);
    [crx{i}, rv]=qr(crx{i}.', 0);
    r0 = size(crx{i-1}, 1); rt0 = size(crx{i-1}, 2);
    crx{i-1} = reshape(crx{i-1}, r0*rt0, r1);
    crx{i-1}=crx{i-1}*rv.';
    r1 = size(rv,1);
    crx{i-1} = reshape(crx{i-1}, r0, rt0, r1);
    crx{i} = reshape(crx{i}.', r1, rt, r2);
end;
for i=1:D
    r1 = size(crx{i}, 1); r2 = size(crx{i}, 3); rt = size(crx{i},2);
    crx{i} = permute(crx{i}, [1, 3, 2]);
    crx{i} = reshape(crx{i}, r1*r2, rt);
    [crx{i}, rv]=qr(crx{i}, 0);
    tuck{i} = tuck{i}*rv.';
    tuck{i} = qtt_fft1(tuck{i}, tol);
    % rank is reversed
    r = tuck{i}.r;
    rq1 = r(1);
    for k=1:rq1
        ek = [zeros(1,k-1), 1, zeros(1,rq1-k)];
        yk = ek*tuck{i}; % first rank is 1
        if (k==1)
            ynew = yk;
        else
            ynew = [ynew, yk];
            ynew = round(ynew, tol);
        end;
    end;
    [tuck{i}, rv]=qr(ynew, 'lr');
%     % reversing 2
%     % nonorth block is the last, but we need to SVD from the first
%     tuck{i} = core2cell(tuck{i});
%     for j=numel(tuck{i}):-1:2
%         tuck{i}{j} = reshape(tuck{i}{j}, r(j), 2*r(j+1));
%         [tuck{i}{j}, rv]=qr(tuck{i}{j}.', 0);
%         tuck{i}{j-1} = reshape(tuck{i}{j-1}, r(j-1)*2, r(j));
%         tuck{i}{j-1} = tuck{i}{j-1}*rv.';
%         r(j) = size(rv, 1);
%         tuck{i}{j} = reshape(tuck{i}{j}.', r(j), 2, r(j+1));
%         tuck{i}{j-1} = reshape(tuck{i}{j-1}, r(j-1), 2, r(j));
%     end;
%     tuck{i}{1} = permute(tuck{i}{1}, [2, 1, 3]);
%     tuck{i}{1} = reshape(tuck{i}{1}, 2, rq1*r(2));
%     [u,s,v]=svd(tuck{i}{1}, 'econ');
%     s = diag(s);
%     rnew = my_chop2(s, norm(s)*tol/sqrt(numel(tuck{i})-1));
%     u=u(:,1:rnew);
%     v=v(:,1:rnew)*diag(s(1:rnew));
%     tuck{i}{1} = reshape(u, 1, 2, rnew);
%     r(1) = 1;
%     rv = reshape(v', rnew*rq1, r(2));
%     for j=2:numel(tuck{i})
%         tuck{i}{j} = reshape(tuck{i}{j}, r(j), 2*r(j+1));
%         tuck{i}{j} = rv*tuck{i}{j};
%         r(j) = rnew;
%         tuck{i}{j} = reshape(tuck{i}{j}, r(j), rq1, 2, r(j+1));
%         tuck{i}{j} = permute(tuck{i}{j}, [1, 3, 2, 4]);
%         tuck{i}{j} = reshape(tuck{i}{j}, r(j)*2, rq1*r(j+1));
%         [u,s,v]=svd(tuck{i}{j}, 'econ');
%         s = diag(s);
%         rnew = my_chop2(s, norm(s)*tol/sqrt(numel(tuck{i})-1));
%         u=u(:,1:rnew);
%         v=v(:,1:rnew)*diag(s(1:rnew));
%         tuck{i}{j} = reshape(u, r(j), 2, rnew);
%         rv = reshape(v', rnew*rq1, r(j+1));
%     end;
%     rv = reshape(rv, rnew, rq1);
%     tuck{i} = cell2core(tt_tensor, tuck{i});

    crx{i} = crx{i}*rv.';
    rt = size(rv,1);
    crx{i} = reshape(crx{i}, r1, r2, rt);
    crx{i} = permute(crx{i}, [1, 3, 2]);
    if (i<D)
        crx{i} = reshape(crx{i}, r1*rt, r2);
        [crx{i}, rv]=qr(crx{i}, 0);
        r3 = size(crx{i+1}, 3); rt2 = size(crx{i+1},2);
        crx{i+1} = reshape(crx{i+1}, r2, rt2*r3);
        crx{i+1} = rv*crx{i+1};
        r2 = size(rv,1);
        crx{i} = reshape(crx{i}, r1, rt, r2);
        crx{i+1} = reshape(crx{i+1}, r2, rt2, r3);
    end;
end;

y = x;
y.core = cell2core(tt_tensor, crx);
y.tuck = tuck;
end

author: oseledets

repo: TT-Toolbox


text:
function x = ffwt ( n, x )

%*****************************************************************************80
%
%% FFWT performs an in-place fast Walsh transform.
%
%  Discussion:
%
%    This routine performs a fast Walsh transform on an input series X
%    leaving the transformed results in X.
%    X is dimensioned N, which must be a power of 2.
%    The results of this Walsh transform are in sequency order.
%
%    The output sequence could be normalized by dividing by N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2011
%
%  Author:
%
%    Ken Beauchamp
%
%  Reference:
%
%    Ken Beauchamp,
%    Walsh functions and their applications,
%    Academic Press, 1975,
%    ISBN: 0-12-084050-2,
%    LC: QA404.5.B33.
%
%  Parameters:
%
%    Input, integer N, the number of items in X.
%    N must be a power of 2.
%
%    Input, real X(N), the data to be transformed.
%
%    Output, real X(N), the transformed data.
%
  m = i4_log_2 ( n );

  for i = 1 : m
    two_power(i) = 2^( m - i );
  end

  for l = 1 : m

    nz = 2^( l - 1 );
    nzi = 2 * nz;
    nzn = floor ( n / nzi );
    nz2 = floor ( nz / 2 );
    if ( nz2 == 0 )
      nz2 = 1;
    end

    for i = 1 : nzn

      js = ( i - 1 ) * nzi;
      z = 1.0;

      for ii = 1 : 2

        for j = 1 : nz2
          js = js + 1;
          j2 = js + nz;
          hold = x(js) + z * x(j2);
          z = - z;
          x(j2) = x(js) + z * x(j2);
          x(js) = hold;
          z = - z;
        end

        if ( l == 1 )
          break
        end

        z = - 1.0;

      end

    end
  end
%
%  Bit reversal section.
%
  nw = 0;
  for k = 1 : n
%
%  Choose correct index and switch elements if not already switched.
%
    if ( k < nw + 1 )
      hold = x(nw+1);
      x(nw+1) = x(k);
      x(k) = hold;
    end
%
%  Bump up series by 1.
%
    for i = 1 : m

      ii = i;
      if ( nw < two_power(i) )
        break
      end
      mw = floor ( nw / two_power(i) );
      mw1 = floor ( mw / 2 );
      if ( mw <= 2 * mw1 )
        break
      end

      nw = nw - two_power(i);

    end

    nw = nw + two_power(ii);

  end

  return
end


text:
function p = gaussian_prob(x, m, C, use_log)
% GAUSSIAN_PROB Evaluate a multivariate Gaussian density.
% p = gaussian_prob(X, m, C)
% p(i) = N(X(:,i), m, C) where C = covariance matrix and each COLUMN of x is a datavector

% p = gaussian_prob(X, m, C, 1) returns log N(X(:,i), m, C) (to prevents underflow).
%
% If X has size dxN, then p has size Nx1, where N = number of examples

if nargin < 4, use_log = 0; end

if length(m)==1 % scalar
  x = x(:)';
end
[d N] = size(x);
%assert(length(m)==d); % slow
m = m(:);
M = m*ones(1,N); % replicate the mean across columns
denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
mahal = sum(((x-M)'*inv(C)).*(x-M)',2);   % Chris Bregler's trick
if any(mahal<0)
  warning('mahal < 0 => C is not psd')
end
if use_log
  p = -0.5*mahal - log(denom);
else
  p = exp(-0.5*mahal) / (denom+eps);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function j = countJumps( f )
%COUNTJUMPS Counts the jumps of a vector f
j = sum( diff(f) ~= 0);

end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function foundScans=selectInitParameters(foundScans)
% foundScans=selectInitParameters(foundScans)
% this function generates a gui to set parameters needed to preprocess fMRI
% data. The gui alows the user to change and set all the parameters needed
% to perform a preprocessing of a complete fMRI session.
%
%
%
%   I needs a foundScans(n) strcuture with n the number off scans in
%   the complete session including scout, inplanes, anaotmy or even shimming
%   founScans should contain:
%
%
%   foundScans(1).DirName='1';
%   foundScans(1).Files=0;
%   foundScans(1).Slices=0;
%   foundScans(1).Volumes=0;
%   foundScans(1).Sequence='';
%   foundScans(1).Action='';
%   foundScans(1).Filenames='';
%   foundScans(1).SkipVols='';
%   foundScans(1).Cycles='';
%   foundScans(1).Filename4d='';
%
%   a good preprocesssing scrip guesses most of these parameters before
%   calling   selectInitParameters to ajust them and add the missing
%   parameters..
%
%
%
%
% Mostly written by Mark Schira (mark@ski.org)

OptionsPerColumn=60;
if isunix
    fontSize = 10;
else
    fontSize = 9;
end
global HandlesOfGUI;
global lastedit;

ActionList={'Nothing','Inplanes','Functional','RefFunc','Anatomical','DTI','Reference'};
nDirs = size(foundScans,2);
disp(nDirs);

ncols=ceil(nDirs/OptionsPerColumn);
nOptionsPerColumn=ceil(nDirs/ncols);
disp(nOptionsPerColumn);

headerStr='Parameters';

%scale factors for x and y axis coordinates
xs = 1.8;
ys = 1.4;
colwidth=64;
colPos=[2 10 35 50 70 82 94 106 118];
%default sizes
butWidth=10;
botMargin = 0.2;
%height = nOptions+2+botMargin*2+.5;
height = nOptionsPerColumn+3.5+botMargin*2+.7;

% width = max(size(optionStr,2),length(headerStr))+2;
width = max(ncols*colwidth,length(headerStr))+2;
width = max(width,2*butWidth+2);

%open the figure
h_mainFig = figure('MenuBar','none',...
    'Units','char',...
    'Resize','off',...
    'NumberTitle','off',...
    'Position',[20, 10, width*xs,height*ys]);
bkColor = get(h_mainFig,'Color');

%___________________________________________________________
%Display title inside a frame
x = 1;
% y = nOptions+1+botMargin;
y = nOptionsPerColumn+2.5+botMargin;
uicontrol('Style','frame',...
    'Units','char',...
    'String',headerStr,...
    'Position',[x*xs,(y+.2)*ys,(width-2)*xs,ys*1.3],...
    'HorizontalAlignment','center',...
    'FontSize',fontSize);
uicontrol('Style','text',...
    'Units','char',...
    'String',headerStr,...
    'Position',[(x+.25)*xs,(y+.3)*ys,(width-2.5)*xs,ys*.9],...
    'HorizontalAlignment','center',...
    'FontSize',fontSize+2);

y=y-1.1;
putuptext('Dir',[colPos(1),(y+.3)*ys,colPos(2)-colPos(1)-1,ys*.9],fontSize);
putuptext('Protocol',[colPos(2),(y+.3)*ys,colPos(3)-colPos(2)-1,ys*.9],fontSize);
putuptext('NoFiles',[colPos(3),(y+.3)*ys,colPos(4)-colPos(3)-4,ys*.9],fontSize);
putuptext('Action',[colPos(4),(y+.3)*ys,colPos(5)-colPos(4)-1,ys*.9],fontSize);
putuptext('Volume',[colPos(5),(y+.3)*ys,colPos(6)-colPos(5)-1,ys*.9],fontSize);
putuptext('Slices',[colPos(6),(y+.3)*ys,colPos(7)-colPos(6)-1,ys*.9],fontSize);
putuptext('Skip',[colPos(7),(y+.3)*ys,colPos(8)-colPos(7)-1,ys*.9],fontSize);
putuptext('Cycles',[colPos(8),(y+.3)*ys,colPos(9)-colPos(8)-1,ys*.9],fontSize);
%putuptext('Condition # ',[colPos(9),(y+.3)*ys,colPos(10)-colPos(9)-1,ys*.9],fontSize);

y=y-0.1;

for ii=1:length(foundScans)
    y=y-1;
    %report the findings
    putuptext(int2str(ii),[colPos(1),(y+.3)*ys,colPos(2)-colPos(1)-1,ys*.9],fontSize);
    putuptext(foundScans(ii).Sequence,[colPos(2),(y+.3)*ys,colPos(3)-colPos(2)-1,ys*.9],fontSize);
    putuptext(int2str(foundScans(ii).Files),[colPos(3),(y+.3)*ys,colPos(4)-colPos(3)-1,ys*.9],fontSize);
    %the aktionMenue
    pos=[colPos(4),(y+.3)*ys,colPos(5)-colPos(4)-4,ys*.9];
    h_action(ii) = uicontrol('Units','char','HorizontalAlignment','left',...
        'Style','popupmenu','String',ActionList,'Position',pos,'BackgroundColor',[0.85 0.85 0.85]);
    if isempty(foundScans(ii).Action)
        ActionVal=1;
    else
        ActionVal=find(strcmp(foundScans(ii).Action,ActionList));
    end
    set(h_action(ii),'Value',ActionVal);

    %the parameters of the scans

    h_Volumes(ii)=uicontrol('Units','char','HorizontalAlignment','center',...
        'Style','edit','String',num2str(foundScans(ii).Volumes,4),...
        'Position',[colPos(5),(y+.3)*ys,colPos(6)-colPos(5)-1,ys*.9],'BackgroundColor',[0.85 0.85 0.85]);
    callBackStr=['global lastedit;lastedit=',char(39),'Vol_',int2str(ii),char(39),';'];
    set(h_Volumes(ii),'Callback',callBackStr);

    h_Slices(ii)=uicontrol('Units','char','HorizontalAlignment','center',...
        'Style','edit','String',int2str(foundScans(ii).Slices),...
        'Position',[colPos(6),(y+.3)*ys,colPos(7)-colPos(6)-1,ys*.9],'BackgroundColor',[0.85 0.85 0.85]);
    callBackStr=['global lastedit;lastedit=',char(39),'Sli_',int2str(ii),char(39),';'];
    set(h_Slices(ii),'Callback',callBackStr);

    h_SkipVols(ii)=uicontrol('Units','char','HorizontalAlignment','center',...
        'Style','edit','String',int2str(foundScans(ii).SkipVols),...
        'Position',[colPos(7),(y+.3)*ys,colPos(8)-colPos(7)-1,ys*.9],'BackgroundColor',[0.85 0.85 0.85]);
    callBackStr=['global lastedit;lastedit=',char(39),'Ski_',int2str(ii),char(39),';'];
    set(h_SkipVols(ii),'Callback',callBackStr);

    h_Cycles(ii)=uicontrol('Units','char','HorizontalAlignment','center',...
        'Style','edit','String',int2str(foundScans(ii).SkipVols),...
        'Position',[colPos(8),(y+.3)*ys,colPos(9)-colPos(8)-1,ys*.9],'BackgroundColor',[0.85 0.85 0.85]);
    callBackStr=['global lastedit;lastedit=',char(39),'Cyc_',int2str(ii),char(39),';'];
    set(h_Cycles(ii),'Callback',callBackStr);
end

HandlesOfGUI.h_Volumes=h_Volumes;
HandlesOfGUI.h_Slices=h_Slices;
HandlesOfGUI.h_SkipVols=h_SkipVols;
HandlesOfGUI.h_Cycles=h_Cycles;

y=y-1.4;
uicontrol('Style','pushbutton',...
    'String','Cancel',...
    'Units','char',...
    'Position',[colPos(4)-4,(y+.3)*ys,colPos(5)-colPos(4)-4,ys*.9],...
    'CallBack','uiresume',...
    'FontSize',fontSize,...
    'UserData','Cancel');

x = width-butWidth-1;
uicontrol('Style','pushbutton',...
    'String','OK',...
    'Units','char',...
    'Position',[colPos(2),(y+.3)*ys,colPos(3)-colPos(2)-4,ys*.9],...
    'CallBack','uiresume',...
    'FontSize',fontSize,...
    'UserData','OK');
global MY_GLOBAL
CopyDownString=['global lastedit; global h_Volumes; selectInitParamChange(lastedit);'];

sabutton=uicontrol('Style','pushbutton',...
    'String','Fill last Value',...
    'Units','char',...
    'Position',[colPos(6)-5,(y+.3)*ys,colPos(8)-colPos(6)-2,ys*.9],...
    'FontSize',fontSize-4,...
    'Callback',CopyDownString,...
    'UserData','OK');

%let the user select some radio buttons and
%wait for a 'uiresume' callback from OK/Cancel
uiwait

%determine which button was hit.
response = get(gco,'UserData');

%gather the status of the radio buttons if 'OK' was
%selected.  Otherwise return empty matrix.
if strcmp(response,'OK')
    for thisScan=1:length(foundScans)
        foundScans(thisScan).Action=ActionList(get(h_action(thisScan),'Value'));
        foundScans(thisScan).Volumes=getNumberFromString(get(h_Volumes(thisScan),'String'));
        foundScans(thisScan).Slices=getNumberFromString(get(h_Slices(thisScan),'String'));
        foundScans(thisScan).SkipVols=getNumberFromString(get(h_SkipVols(thisScan),'String'));
        foundScans(thisScan).Cycles=getNumberFromString(get(h_Cycles(thisScan),'String'));
    end

else
    foundScans = [];
end

close(h_mainFig);

return;
end

function putuptext(TEXT,pos,fontSize)
uicontrol('Style','text',...
    'Units','char',...
    'String',TEXT ,...
    'Position',pos,...
    'HorizontalAlignment','center',...
    'FontSize',9);
end


text:
function b = tmat_rot_vector ( a, angle, axis )

%*****************************************************************************80
%
%% TMAT_ROT_VECTOR applies an arbitrary axis rotation to the geometric transformation matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics, Principles and Practice,
%    Addison Wesley, Second Edition, 1990.
%
%  Parameters:
%
%    Input, real A(4,4), the current geometric transformation
%    matrix.
%
%    Input, real ANGLE, the angle, in degrees, of the rotation.
%
%    Input, real AXIS(3), the axis vector about which
%    rotation occurs.  AXIS may not be the zero vector.
%
%    Output, real B(4,4), the modified geometric
%    transformation matrix.
%
  dim_num = 3;

  norm = sqrt ( sum ( axis(1:dim_num).^2 ) );

  if ( norm == 0.0 )
    return
  end

  axis(1:dim_num) = axis(1:dim_num) / norm;

  angle_rad = degrees_to_radians ( angle );
  ca = cos ( angle_rad );
  sa = sin ( angle_rad );

  c = tmat_init ( );

  c(1,1) =                axis(1) * axis(1) + ca * ( 1.0 - axis(1) * axis(1) );
  c(1,2) = ( 1.0 - ca ) * axis(1) * axis(2) - sa * axis(3);
  c(1,3) = ( 1.0 - ca ) * axis(1) * axis(3) + sa * axis(2);

  c(2,1) = ( 1.0 - ca ) * axis(2) * axis(1) + sa * axis(3);
  c(2,2) =                axis(2) * axis(2) + ca * ( 1.0 - axis(2) * axis(2) );
  c(2,3) = ( 1.0 - ca ) * axis(2) * axis(3) - sa * axis(1);

  c(3,1) = ( 1.0 - ca ) * axis(3) * axis(1) - sa * axis(2);
  c(3,2) = ( 1.0 - ca ) * axis(3) * axis(2) + sa * axis(1);
  c(3,3) =                axis(3) * axis(3) + ca * ( 1.0 - axis(3) * axis(3) );

  b(1:4,1:4) = c(1:4,1:4) * a(1:4,1:4);

  return
end


text:
function module = noise_module_fixed(tau)

module.initialize = @initialize;
module.update = @update;
module.get_struct = @get_struct;

  function [Tau] = initialize(M,N)
  Tau = bsxfun(@times, tau, ones(M,N));
  end

  function S = get_struct()
  S.module = 'noise_module_fixed';
  S.tau = tau;
  end

  function [Tau, LogTau, KL_Tau] = update(iter, E2, Obs)

  I = ones(size(Obs));
  Tau = bsxfun(@times, tau, I);
  LogTau = bsxfun(@times, log(tau), I);
  KL_Tau = 0;

  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%SEQUENTIAL Sequential mapping
%
%   V = SEQUENTIAL(W1,W2)
%   B = SEQUENTIAL(A,W)
%
% INPUT
%   W,W1,W2  Mappings
%   A        Dataset
%
% OUTPUT
%   V        Sequentially combined mapping
%   B        Dataset
%
% DESCRIPTION
% The two mappings W1 and W2 are combined into a single mapping V. Note
% that SEQUENTIAL(W1,W2) is equivalent to W1*W2. If W2 is a mapping of
% a type 'combiner', it is called to make a combination attempt.
%	SEQUENTIAL(A,W) maps the dataset A by the sequential mapping W.
%
% This routine is automatically called to execute W = W1*W2 or B = A*W2 in
% case W2 is a sequential mapping. It should not be directly called by users.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function [w,varargout] = sequential(w1,w2,w3)

	% Just necessary to inform PRMAP.
	if (nargin == 0)
		w = prmapping(mfilename,'combiner');
		return;
	end
	[m1,k1] = size(w1);
	[m2,k2] = size(w2);
	if (~isa(w2,'prmapping'))
		error('Second argument should be a mapping.')
	end

	if isa(w1,'prmapping')
		% Definition
    if isempty(w1) % treat empty mapping as unity mapping
      w = w2;
    elseif (iscombiner(w2))
			% Execute the mapping W2.
      varargout = repmat({[]},[1, max((nargout-1),0)]);
			map = getmapping_file(w2);
			pars = getdata(w2);
			[w,varargout{:}] = feval(map,w1,pars{:});
		else
			% W2 is just a general mapping.
			if (k1 > 0) && (m2 > 0) && (k1 ~= m2)
				error('Inner mapping/data sizes do not agree.')
			end

			% Define the mapping type after combining W1 and W2.
			if (isuntrained(w1)) || (isuntrained(w2))
				mappingtype = 'untrained';
      elseif (isgenerator(w2))
        mappingtype = 'generator';
      %elseif (isfixed(w2)) || (isfixed_cell(w2))
      %  mappingtype = 'fixed';
			elseif (istrained(w1)) || (istrained(w2))
				mappingtype = 'trained';
			else
				mappingtype = 'fixed';
			end

			if strcmp(mappingtype,'untrained')
				labels = [];
				size_in = 0;
				size_out = 0;
			elseif (m2 == 0 || k2 == 0) && (m1 ~= 0) && (k1 ~= 0)
				% E.G. TRAINED * FIXED
				labels   = getlabels(w1);
				size_in  = getsize_in(w1);
				size_out = getsize_out(w1);

			elseif (m2 ~= 0) && (k2 ~= 0) && (m1 == 0 || k1 == 0)
				% FIXED * TRAINED
				labels   = getlabels(w2);
				size_in  = getsize_in(w2);
				size_out = getsize_out(w2);

			elseif ~istrained(w2)
				% TRAINED * FIXED
				labels   = getlabels(w1);
				size_in  = getsize_in(w1);
				size_out = getsize_out(w2);

			else
				% TRAINED * TRAINED
				labels = getlabels(w2);
				size_in = getsize_in(w1);
				size_out = getsize_out(w2);
			end
			w = prmapping(mfilename,mappingtype,{w1,w2},labels,size_in,size_out);
    end
    if ismapping(w) && (getbatch(w1) || getbatch(w2))
      [n1,b1,o1] = getbatch(w1);
      [n2,b2,o2] = getbatch(w2);
      if ~n1
        b12 = b2; o12 = o2;
      elseif ~n2
        b12 = b1; o12 = o1;
      else
        b12 = min(b1,b2); o12 = min(o1,o2);
      end
      w = setbatch(w,true,b12,o12);
    end

  elseif isempty(w2) % treat empty mapping as unity mapping
    w = w1;

  else
		% Execution. We are here, when SEQUENTIAL(A,V) is called.
		if nargin == 3 % needed as MAP breaks down sequential mappings
			w2 = w2*w3;  % restore them!
		end
		a = w1;
		if (~isa(a,'double')) && (~isa(a,'prdataset'))
			error('Just datasets or doubles can be mapped.')
		end
		% V can be a more complex mapping.
		% v = +w2; v1 = v{1}; v2 = v{2};
    [v1,v2] = getdata(w2);
		if (isuntrained(v1))
			if (isuntrained(v2))
				u = a*v1;
				w = u*(a*u*v2);
			else
				w = a*v1*v2;
			end
		else
			if (isuntrained(v2)) && (~isgenerator(v1))
				w = v1*(a*v1*v2);
				% may be v1 changed the dimensionality, reset it:
				w = setsize_in(w,size(a,2));
      else
        w1 = a*v1;
        w = w1*v2;
				%w = a*v1*v2;
				featlabels = getlabels(w2);
				if (isdataset(w)) && ~isempty(featlabels) && size(w,2) == size(featlabels,1)
					w = setfeatlab(w,featlabels);
				end
			end
		end
		if ismapping(w)
			w = setbatch(w,getbatch(w2));
		end
	end
return;

author: marianux

repo: ecg-kit

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


text:
%% Demo: Threshold-Free Cluster Enhancement (TFCE) on surface dataset
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
%    A participant made finger pressed with the index and middle finger of
%    the right hand during 4 runs in an fMRI study. Each run was divided in
%    4 blocks with presses of each finger and analyzed with the GLM,
%    resulting in 2*4*4=32 t-values
%
% This example illustrates the use of Threshold-Free Cluster Enhancement
% with a permutation test to correct for multiple comparisons.
%
% TFCE reference: Stephen M. Smith, Thomas E. Nichols, Threshold-free
% cluster enhancement: Addressing problems of smoothing, threshold
% dependence and localisation in cluster inference, NeuroImage, Volume 44,
% Issue 1, 1 January 2009, Pages 83-98.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Check externals
cosmo_check_external({'surfing','afni'});

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

% resolution parameter for input surfaces
% 64 is for high-quality results; use 16 for fast execution
ld=16;

% reset citation list
cosmo_check_external('-tic');

% load single surface
intermediate_fn=fullfile(digit_study_path,...
                            sprintf('ico%d_mh.intermediate_al.asc', ld));
[vertices,faces]=surfing_read(intermediate_fn);

%% Load functional data
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig');

targets=repmat(1:2,1,16)';    % class labels:  1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/4)+1; % half-run:      1 1 1 1 2 2 2 2 3 3 ... 8 8

vol_ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);

%% Map univariate response data to surface

% this measure averages the data near each node to get a surface dataset
radius=0;
surf_band_range=[-2 2]; % get voxel data within 2mm from surface
surf_def={vertices,faces,[-2 2]};
nbrhood=cosmo_surficial_neighborhood(vol_ds,surf_def,'radius',radius);

measure=@(x,opt) cosmo_structjoin('samples',mean(x.samples,2),'sa',x.sa);

surf_ds=cosmo_searchlight(vol_ds,nbrhood,measure);

fprintf('Univariate surface data:\n');
cosmo_disp(surf_ds);

%% Average data in each chunk
% for this example only consider the samples in the first condition
% (targets==1), and average the samples in each chunk
%
% for group analysis: set chunks to (1:nsubj)', assuming each sample is
% data from a single participant
surf_ds=cosmo_slice(surf_ds,surf_ds.sa.targets==1);
surf_ds=cosmo_average_samples(surf_ds);

fn_surf_ds=fullfile(output_path, 'digit_target1.niml.dset');

% save to disc
cosmo_map2surface(surf_ds, fn_surf_ds);
fprintf('Input data saved to %s\n', fn_surf_ds);

%% Run Threshold-Free Cluster Enhancement (TFCE)

% All data is prepared; surf_ds has 8 samples and 5124 nodes. We want to
% see if there are clusters that show a significant difference from zero in
% their response. Thus, .sa.targets is set to all ones (the same
% condition), whereas .sa.chunks is set to (1:8)', indicating that all
% samples are assumed to be independent.
%
% (While this is a within-subject analysis, exactly the same logic can be
% applied to a group-level analysis)

% define neighborhood for each feature
% (cosmo_cluster_neighborhood can be used also for meeg or volumetric
% fmri datasets)
cluster_nbrhood=cosmo_cluster_neighborhood(surf_ds,...
                                        'vertices',vertices,'faces',faces);

fprintf('Cluster neighborhood:\n');
cosmo_disp(cluster_nbrhood);

opt=struct();

% number of null iterations. for publication-quality, use >=1000;
% 10000 is even better
opt.niter=250;

% in this case we run a one-sample test against a mean of 0, and it is
% necessary to specify the mean under the null hypothesis
% (when testing classification accuracies, h0_mean should be set to chance
% level, assuming a balanced crossvalidation scheme was used)
opt.h0_mean=0;

% this example uses the data itself (with resampling) to obtain cluster
% statistcs under the null hypothesis. This is (in this case) somewhat
% conservative due to how the resampling is performed.
% Alternatively, and for better estimates (at the cost of computational
% cost), one can generate a set of (say, 50) datasets using permuted data
% e.g. using cosmo_randomize_targets), put them in a cell and provide
% them as the null argument.
opt.null=[];

fprintf('Running multiple-comparison correction with these options:\n');
cosmo_disp(opt);

% Run TFCE-based cluster correction for multiple comparisons.
% The output has z-scores for each node indicating the probablity to find
% the same, or higher, TFCE value under the null hypothesis
tfce_ds=cosmo_montecarlo_cluster_stat(surf_ds,cluster_nbrhood,opt);

%% Show results

fprintf('TFCE z-score dataset\n');
cosmo_disp(tfce_ds);

nfeatures=size(tfce_ds.samples,2);
percentiles=(1:nfeatures)/nfeatures*100;
plot(percentiles,sort(tfce_ds.samples))
title('sorted TFCE z-scores');
xlabel('feature percentile');
ylabel('z-score');

nvertices=size(vertices,1);
disp_opt=struct();
disp_opt.DataRange=[-2 2];

DispIVSurf(vertices,faces,1:nvertices,tfce_ds.samples',0,disp_opt);

% store results
fn_tfce_ds=fullfile(output_path, 'digit_target1_tfce.niml.dset');
cosmo_map2surface(tfce_ds, fn_tfce_ds);

surf_fn=fullfile(output_path, 'digit_intermediate.asc');
surfing_write(surf_fn,vertices,faces);

% show citation information
cosmo_check_external('-cite');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% mycorr2 modified version of the 2D correlation
%         for the use with im2col and col2im
%         see GETPOINT
%
%
% $Id: mycorr2.m,v 2.0 2003/06/19 12:06:52 svoboda Exp $

% Note: It written in order to gain speed. The clarity of the code suffers accordingly

function R = mycorr2(X,G,Gn,Gn2)

% Gn  = G-mean(G);
% Gn2 = sqrt(sum(Gn.^2));

mX	= repmat(mean(X),size(X,1),1);
mXn = X - mX;
smX	= sum(mXn.^2);

numerator = (mXn'*Gn)';
denominator = smX*Gn2;

R = numerator./denominator;

return

author: strawlab

repo: MultiCamSelfCal