text:
function ret = spm_ov_roi(varargin)
% ROI tool - plugin for spm_orthviews
%
% With ROI tool it is possible to create new or modify existing mask images
% interactively. ROI tool can be launched via the spm_orthviews image
% context menu.
% While ROI tool is active, mouse buttons have the following functions:
% left       Reposition crosshairs
% middle     Perform ROI tool box selection according to selected edit mode at
%            crosshair position
% right      context menu
%
% Menu options and prompts explained:
% Launch     Initialise ROI tool in current image
%            'Load existing ROI image? (yes/no)'
%              If you want to modify an existing mask image (e.g. mask.img from
%              a fMRI analysis), press 'yes'. You will then be prompted to
%            'Select ROI image'
%              This is the image that will be loaded as initial ROI.
%              If you want to create a new ROI image, you will first be
%              prompted to
%            'Select image defining ROI space'
%              The image dimensions, voxel sizes and slice orientation will
%              be read from this image. Thus you can edit a ROI based on a
%              image with a resolution and slice orientation different from
%              the underlying displayed image.
%
% Once ROI tool is active, the menu consists of three parts: settings,
% edit operations and load/save operations.
% Settings
% --------
% Selection  Operation performed when pressing the middle mouse button or
% mode       by clustering operations.
%            'Set selection'
%              The selection made with the following commands will
%              be included in your ROI.
%            'Clear selection'
%              The selection made with the following commands will
%              be excluded from your ROI.
% Box size   Set size of box to be (de)selected when pressing the
%            middle mouse button.
% Polygon    Set number of adjacent slices selected by one polygon
% slices     drawing.
% Cluster    Set minimum cluster size for "Cleanup clusters" and
% size       "Connected cluster" operations.
% Erosion/   During erosion/dilation operations, the binary mask will be
% dilation   smoothed. At boundaries, this will result in mask values
% threshold  that are not exactly zero or one, but somewhere in
%            between. Whether a mask will be eroded (i.e. be smaller than
%            the original) or dilated (i.e. grow) depends on this
%            threshold. A threshold below 0.5 dilates, above 0.5 erodes a
%            mask.
% Edit actions
% ------------
% Polygon    Draw an outline on one of the 3 section images. Voxels
%            within the outline will be added to the ROI. The same
%            outline can be applied to a user-defined number of
%            consecutive slices around the current crosshair position.
% Threshold  You will be prompted to enter a [min max] threshold. Only
%            those voxels in the ROI image where the intensities of the
%            underlying image are within the [min max] range will survive
%            this operation.
% Connected  Select only voxels that are connected to the voxel at
% cluster    current crosshair position through the ROI.
% Cleanup    Keep only clusters that are larger than a specified cluster
% clusters   size.
% Erode/     Erode or dilate a mask, using the current erosion/dilation
% Dilate     threshold.
% Invert     Invert currently defined ROI
% Clear      Clear ROI, but keep ROI space information
% Add ROI from file(s)
%            Add ROIs from file(s) into current ROI set. According to the
%            current edit mode voxels unequal zero will be set or
%            cleared. The image files will be resampled and thus do not
%            need to have the same orientation or voxel size as the
%            original ROI.
% Save actions
% ------------
% Save       Save ROI image
% Save As    Save ROI image under a new file name
%            The images will be rescaled to 0 (out of mask) and 1 (in
%            mask).
% Quit       Quit ROI tool
%
% This routine is a plugin to spm_orthviews. For general help about
% spm_orthviews and plugins type
%             help spm_orthviews
% at the MATLAB prompt.
%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Volkmar Glauche
% $Id: spm_ov_roi.m 6991 2017-01-19 13:09:51Z guillaume $

% Note: This plugin depends on the blobs set by spm_orthviews('addblobs',...)
% They should not be removed while ROI tool is active and no other blobs be
% added. This restriction may be removed when using the 'alpha' property
% to overlay blobs onto images.

global st;
if isempty(st)
    error('roi: This routine can only be called as a plugin for spm_orthviews!');
end

if nargin < 2
    error('roi: Wrong number of arguments. Usage: spm_orthviews(''roi'', cmd, volhandle, varargin)');
end

cmd = lower(varargin{1});
volhandle = varargin{2};

toset = [];
toclear = [];
tochange = [];
update_roi = false;

switch cmd
    case 'init'
        % spm_ov_roi('init', volhandle, Vroi, loadasroi, xyz)
        spm('pointer','watch');
        Vroi = spm_vol(varargin{3});
        switch varargin{4} % loadasroi
            case 1,
                roi = spm_read_vols(Vroi)>0;
                [x,y,z] = ndgrid(1:Vroi.dim(1),1:Vroi.dim(2),1:Vroi.dim(3));
                xyz = [x(roi(:))'; y(roi(:))'; z(roi(:))'];
            case {0,2} % ROI space image or SPM mat
                Vroi = rmfield(Vroi,'private');
                roi = false(Vroi.dim(1:3));
                Vroi.fname = fileparts(Vroi.fname); % save path
                xyz = varargin{5};
                if ~isempty(xyz)
                    ind = sub2ind(Vroi.dim(1:3),xyz(1,:),xyz(2,:),xyz(3,:));
                    roi(ind) = true;
                end;
        end;
        % reset data type to save disk space
        Vroi.dt(1) = spm_type('uint8');
        % reset scaling factor of Vroi handle
        Vroi.pinfo(1:2) = Inf;
        clear x y z

        % draw a frame only if ROI volume different from underlying GM volume
        if any(Vroi.dim(1:3)-st.vols{volhandle}.dim(1:3))|| ...
                any(Vroi.mat(:)-st.vols{volhandle}.mat(:))
            [xx1,yx1,zx1] = ndgrid(1            , 1:Vroi.dim(2), 1:Vroi.dim(3));
            [xx2,yx2,zx2] = ndgrid(Vroi.dim(1)  , 1:Vroi.dim(2), 1:Vroi.dim(3));
            [xy1,yy1,zy1] = ndgrid(1:Vroi.dim(1), 1            , 1:Vroi.dim(3));
            [xy2,yy2,zy2] = ndgrid(1:Vroi.dim(1), Vroi.dim(2)  , 1:Vroi.dim(3));
            [xz1,yz1,zz1] = ndgrid(1:Vroi.dim(1), 1:Vroi.dim(2), 1);
            [xz2,yz2,zz2] = ndgrid(1:Vroi.dim(1), 1:Vroi.dim(2), Vroi.dim(3));

            fxyz = [xx1(:)' xx2(:)' xy1(:)' xy2(:)' xz1(:)' xz2(:)'; ...
                    yx1(:)' yx2(:)' yy1(:)' yy2(:)' yz1(:)' yz2(:)'; ...
                    zx1(:)' zx2(:)' zy1(:)' zy2(:)' zz1(:)' zz2(:)'];
            clear xx1 yx1 zx1 xx2 yx2 zx2 xy1 yy1 zy1 xy2 yy2 zy2 xz1 yz1 zz1 xz2 yz2 zz2
            hframe = 1;
        else
            hframe = [];
            fxyz  = [];
        end;

        cb = cell(1,3);
        for k=1:3
            cb{k}=get(st.vols{volhandle}.ax{k}.ax,'ButtonDownFcn');
            set(st.vols{volhandle}.ax{k}.ax,...
                'ButtonDownFcn',...
                @(ob,ev)spm_ov_roi('bdfcn',volhandle,ob,ev));
        end;

        st.vols{volhandle}.roi = struct('Vroi',Vroi, 'xyz',xyz, 'roi',roi,...
                                        'hroi',1, 'fxyz',fxyz,...
                                        'hframe',hframe, 'mode','set',...
                                        'tool', 'box', ...
                                        'thresh',[60 140], 'box',[4 4 4],...
                                        'cb',[], 'polyslices',1, 'csize',5,...
                                        'erothresh',.5);
        st.vols{volhandle}.roi.cb = cb;

        if ~isempty(st.vols{volhandle}.roi.fxyz)
            if isfield(st.vols{volhandle}, 'blobs')
                st.vols{volhandle}.roi.hframe = numel(st.vols{volhandle}.blobs)+1;
            end;
            spm_orthviews('addcolouredblobs',volhandle, ...
                          st.vols{volhandle}.roi.fxyz,...
                          ones(size(st.vols{volhandle}.roi.fxyz,2),1), ...
                          st.vols{volhandle}.roi.Vroi.mat,[1 .5 .5]);
            st.vols{volhandle}.blobs{st.vols{volhandle}.roi.hframe}.max=1.3;
        end;
        update_roi=1;
        obj = findobj(0, 'Tag',  sprintf('ROI_1_%d', volhandle));
        set(obj, 'Visible', 'on');
        obj = findobj(0, 'Tag',  sprintf('ROI_0_%d', volhandle));
        set(obj, 'Visible', 'off');

    case 'edit'
        switch st.vols{volhandle}.roi.tool
            case 'box'
                spm('pointer','watch');
                pos = round(st.vols{volhandle}.roi.Vroi.mat\ ...
                            [spm_orthviews('pos'); 1]);
                tmp = round((st.vols{volhandle}.roi.box-1)/2);
                [sx,sy,sz] = meshgrid(-tmp(1):tmp(1), -tmp(2):tmp(2), -tmp(3):tmp(3));
                sel = [sx(:)';sy(:)';sz(:)']+repmat(pos(1:3), 1,prod(2*tmp+1));
                tochange = sel(:, (all(sel>0) &...
                                   sel(1,:)<=st.vols{volhandle}.roi.Vroi.dim(1) & ...
                                   sel(2,:)<=st.vols{volhandle}.roi.Vroi.dim(2) & ...
                                   sel(3,:)<=st.vols{volhandle}.roi.Vroi.dim(3)));
                update_roi = 1;

            case 'connect'
                spm_ov_roi('connect', volhandle)

            case 'poly'

                % @COPYRIGHT  :
      %             Copyright 1993,1994 Mark Wolforth and Greg Ward, McConnell
      %             Brain Imaging Centre, Montreal Neurological Institute, McGill
      %             University.
      %             Permission to use, copy, modify, and distribute this software
      %             and its documentation for any purpose and without fee is
      %             hereby granted, provided that the above copyright notice
      %             appear in all copies.  The authors and McGill University make
      %             no representations about the suitability of this software for
      %             any purpose.  It is provided "as is" without express or
      %             implied warranty.
                for k = 1:3
                    if st.vols{volhandle}.ax{k}.ax == gca
                        axhandle = k;
                        break;
                    end;
                end;
                line_color = [1 1 0];
                axes(st.vols{volhandle}.ax{axhandle}.ax);

                hold on;
                Xlimits = get (st.vols{volhandle}.ax{axhandle}.ax,'XLim');
                Ylimits = get (st.vols{volhandle}.ax{axhandle}.ax,'YLim');
                XLimMode = get(st.vols{volhandle}.ax{axhandle}.ax,'XLimMode');
                set(st.vols{volhandle}.ax{axhandle}.ax,'XLimMode','manual');
                YLimMode = get(st.vols{volhandle}.ax{axhandle}.ax,'YLimMode');
                set(st.vols{volhandle}.ax{axhandle}.ax,'YLimMode','manual');
                ButtonDownFcn = get(st.vols{volhandle}.ax{axhandle}.ax,'ButtonDownFcn');
                set(st.vols{volhandle}.ax{axhandle}.ax,'ButtonDownFcn','');
                UIContextMenu = get(st.vols{volhandle}.ax{axhandle}.ax,'UIContextMenu');
                set(st.vols{volhandle}.ax{axhandle}.ax,'UIContextMenu',[]);
                set(st.vols{volhandle}.ax{axhandle}.ax,'Selected','on');
                disp (['Please mark the ROI outline in the highlighted image' ...
                       ' display.']);
                disp ('Points outside the ROI image area will be clipped to');
                disp ('the image boundaries.');
                disp ('Left-Click on the vertices of the ROI...');
                disp ('Middle-Click to finish ROI selection...');
                disp ('Right-Click to cancel...');
                x=Xlimits(1);
                y=Ylimits(1);
                i=1;
                lineHandle = [];
                xc = 0; yc = 0; bc = 0;
                while ~isempty(bc)
                    [xc,yc,bc] = ginput(1);
                    if isempty(xc) || bc > 1
                        if bc == 3
                            x = []; y=[];
                        end;
                        if bc == 2 || bc == 3
                            bc = [];
                            break;
                        end;
                    else
                        if xc > Xlimits(2)
                            xc = Xlimits(2);
                        elseif xc < Xlimits(1)
                            xc = Xlimits(1);
                        end;
                        if yc > Ylimits(2)
                            yc = Ylimits(2);
                        elseif yc < Ylimits(1)
                            yc = Ylimits(1);
                        end;
                        x(i) = xc;
                        y(i) = yc;
                        i=i+1;
                        if ishandle(lineHandle)
                            delete(lineHandle);
                        end;
                        lineHandle = line (x,y,ones(1,length(x)), ...
                                           'Color',line_color,...
                                           'parent',st.vols{volhandle}.ax{axhandle}.ax,...
                                           'HitTest','off');
                    end;
                end

                if ishandle(lineHandle)
                    delete(lineHandle);
                end;
                if ~isempty(x)
                    spm('pointer','watch');
                    x(i)=x(1);
                    y(i)=y(1);
                    prms=spm_imatrix(st.vols{volhandle}.roi.Vroi.mat);
                    % Code from spm_orthviews('redraw') for determining image
                    % positions
                    is   = inv(st.Space);
                    cent = is(1:3,1:3)*st.centre(:) + is(1:3,4);
                    polyoff = [0 0 0];
                    switch axhandle
                        case 1,
                            M0 = [ 1 0 0 -st.bb(1,1)+1
                                   0 1 0 -st.bb(1,2)+1
                                   0 0 1 -cent(3)
                                   0 0 0 1];
                            polyoff(3) = st.vols{volhandle}.roi.polyslices/2;
                            polythick = prms(9);
                        case 2,
                            M0 = [ 1 0 0 -st.bb(1,1)+1
                                   0 0 1 -st.bb(1,3)+1
                                   0 1 0 -cent(2)
                                   0 0 0 1];
                            polyoff(2) = st.vols{volhandle}.roi.polyslices/2;
                            polythick = prms(8);
                        case 3,
                            if st.mode ==0,
                                M0 = [ 0 0 1 -st.bb(1,3)+1
                                       0 1 0 -st.bb(1,2)+1
                                       1 0 0 -cent(1)
                                       0 0 0 1];
                            else
                                M0 = [ 0 -1 0 +st.bb(2,2)+1
                                       0  0 1 -st.bb(1,3)+1
                                       1  0 0 -cent(1)
                                       0  0 0 1];
                            end;
                            polyoff(1) = st.vols{volhandle}.roi.polyslices/2;
                            polythick = abs(prms(7));
                    end;
                    polvx = st.vols{volhandle}.roi.Vroi.mat\(st.Space*(M0\...
                            [x(:)';y(:)'; zeros(size(x(:)')); ones(size(x(:)'))]));
                    % Bounding volume for polygon in ROI voxel space
                    [xbox,ybox,zbox] = ndgrid(max(min(floor(polvx(1,:)-polyoff(1))),1):...
                                              min(max(ceil(polvx(1,:)+polyoff(1))),...
                                                  st.vols{volhandle}.roi.Vroi.dim(1)),...
                                              max(min(floor(polvx(2,:)-polyoff(2))),1):...
                                              min(max(ceil(polvx(2,:)+polyoff(2))),...
                                                  st.vols{volhandle}.roi.Vroi.dim(2)),...
                                              max(min(floor(polvx(3,:)-polyoff(3))),1):...
                                              min(max(ceil(polvx(3,:)+polyoff(3))),...
                                                  st.vols{volhandle}.roi.Vroi.dim(3)));
                    % re-transform in polygon plane
                    xyzbox = M0*(st.Space\st.vols{volhandle}.roi.Vroi.mat)*[xbox(:)';ybox(:)';zbox(:)';...
                                        ones(size(xbox(:)'))];
                    xyzbox = xyzbox(:,abs(xyzbox(3,:))<=.6*polythick*...
                                    st.vols{volhandle}.roi.polyslices); % nearest neighbour to polygon
                    sel = logical(inpolygon(xyzbox(1,:),xyzbox(2,:),x,y));
                    xyz = st.vols{volhandle}.roi.Vroi.mat\(st.Space*(M0\xyzbox(:,sel)));
                    if ~isempty(xyz)
                        tochange = round(xyz(1:3,:));
                        update_roi = 1;
                    end;
                end;
                set(st.vols{volhandle}.ax{axhandle}.ax,...
                    'Selected','off', 'XLimMode',XLimMode, 'YLimMode',YLimMode,...
                    'ButtonDownFcn',ButtonDownFcn, 'UIContextMenu',UIContextMenu);
        end;
    case 'thresh'
        spm('pointer','watch');
        rind = find(st.vols{volhandle}.roi.roi);
        [x,y,z]=ind2sub(st.vols{volhandle}.roi.Vroi.dim(1:3),rind);
        tmp = round(st.vols{volhandle}.mat \ ...
                    st.vols{volhandle}.roi.Vroi.mat*[x'; y'; z'; ones(size(x'))]);
        dat = spm_sample_vol(st.vols{volhandle}, ...
                             tmp(1,:), tmp(2,:), tmp(3,:), 0);
        sel = ~((st.vols{volhandle}.roi.thresh(1) < dat) & ...
                (dat < st.vols{volhandle}.roi.thresh(2)));
        if strcmp(st.vols{volhandle}.roi.mode,'set')
            toclear = [x(sel)'; y(sel)'; z(sel)'];
        else
            toset   = [x(sel)'; y(sel)'; z(sel)'];
            toclear = st.vols{volhandle}.roi.xyz;
        end;
        update_roi = 1;

    case 'erodilate'
        spm('pointer','watch');
        V = zeros(size(st.vols{volhandle}.roi.roi));
        spm_smooth(double(st.vols{volhandle}.roi.roi), V, 2);
        [ero(1,:),ero(2,:),ero(3,:)] = ind2sub(st.vols{volhandle}.roi.Vroi.dim(1:3),...
                                               find(V(:)>st.vols{volhandle}.roi.erothresh));
        if strcmp(st.vols{volhandle}.roi.mode,'set')
            toset   = ero;
            toclear = st.vols{volhandle}.roi.xyz;
        else
            toclear = ero;
        end;
        update_roi = 1;

    case {'connect', 'cleanup'}
        spm('pointer','watch');
        [V,L] = spm_bwlabel(double(st.vols{volhandle}.roi.roi),6);
        sel = [];
        switch cmd
            case 'connect'
                pos = round(st.vols{volhandle}.roi.Vroi.mat\ ...
                            [spm_orthviews('pos'); 1]);
                sel = V(pos(1),pos(2),pos(3));
                if sel == 0
                    sel = [];
                end;
            case 'cleanup'
                numV = zeros(1,L);
                for k = 1:L
                    numV(k) = sum(V(:)==k);
                end;
                sel = find(numV>st.vols{volhandle}.roi.csize);
        end;
        if ~isempty(sel)
            ind1 = cell(1,numel(sel));
            for k=1:numel(sel)
                ind1{k} = find(V(:) == sel(k));
            end;
            ind = cat(1,ind1{:});
            conn = zeros(3,numel(ind));
            [conn(1,:),conn(2,:),conn(3,:)] = ...
                ind2sub(st.vols{volhandle}.roi.Vroi.dim(1:3),ind);

            if strcmp(st.vols{volhandle}.roi.mode,'set')
                toset   = conn;
                toclear = st.vols{volhandle}.roi.xyz;
            else
                toclear = conn;
            end;
        end;
        update_roi = 1;

    case 'invert'
        spm('pointer','watch');
        st.vols{volhandle}.roi.roi = ~st.vols{volhandle}.roi.roi;
        ind = find(st.vols{volhandle}.roi.roi);
        st.vols{volhandle}.roi.xyz = zeros(3,numel(ind));
        [st.vols{volhandle}.roi.xyz(1,:), ...
            st.vols{volhandle}.roi.xyz(2,:), ...
            st.vols{volhandle}.roi.xyz(3,:)] = ind2sub(st.vols{volhandle}.roi.Vroi.dim(1:3),ind);
        update_roi = 1;

    case 'clear'
        spm('pointer','watch');
        st.vols{volhandle}.roi.roi = false(size(st.vols{volhandle}.roi.roi));
        st.vols{volhandle}.roi.xyz=[];
        update_roi = 1;

    case 'addfile'
        V = spm_vol(spm_select([1 Inf],'image','Image(s) to add'));
        [x,y,z] = ndgrid(1:st.vols{volhandle}.roi.Vroi.dim(1),...
                         1:st.vols{volhandle}.roi.Vroi.dim(2),...
                         1:st.vols{volhandle}.roi.Vroi.dim(3));
        xyzmm = st.vols{volhandle}.roi.Vroi.mat*[x(:)';y(:)';z(:)'; ...
                            ones(1, prod(st.vols{volhandle}.roi.Vroi.dim(1:3)))];
        msk = false(1,prod(st.vols{volhandle}.roi.Vroi.dim(1:3)));
        for k = 1:numel(V)
            xyzvx = V(k).mat\xyzmm;
            dat = spm_sample_vol(V(k), xyzvx(1,:), xyzvx(2,:), xyzvx(3,:), 0);
            dat(~isfinite(dat)) = 0;
            msk = msk | logical(dat);
        end;
        [tochange(1,:),tochange(2,:),tochange(3,:)] = ind2sub(st.vols{volhandle}.roi.Vroi.dim(1:3),find(msk));
        clear xyzmm xyzvx msk
        update_roi = 1;

    case {'save','saveas'}
        if strcmp(cmd,'saveas') || ...
                exist(st.vols{volhandle}.roi.Vroi.fname, 'dir')
            flt = {'*.nii','NIfTI (1 file)';'*.img','NIfTI (2 files)'};
            [name,pth,idx] = uiputfile(flt, 'Output image');
            if ~ischar(pth)
                warning('spm:spm_ov_roi','Save cancelled');
                return;
            end;
            [p,n,e,v] = spm_fileparts(fullfile(pth,name));
            if isempty(e)
                e = flt{idx,1}(2:end);
            end;
            st.vols{volhandle}.roi.Vroi.fname = fullfile(p, [n e v]);
        end;
        spm('pointer','watch');
        spm_write_vol(st.vols{volhandle}.roi.Vroi, ...
                      st.vols{volhandle}.roi.roi);
        spm('pointer','arrow');
        return;

    case 'redraw'
        % do nothing
        return;

 %-------------------------------------------------------------------------
 % Context menu and callbacks
    case 'context_menu'
        item0 = uimenu(varargin{3}, 'Label', 'ROI tool');
        item1 = uimenu(item0, 'Label', 'Launch', 'Callback', ...
                       sprintf('%s(''context_init'', %d);', mfilename, volhandle), ...
                       'Tag',sprintf('ROI_0_%d', volhandle));
        item2 = uimenu(item0, 'Label', 'Selection mode', ...
                       'Visible', 'off', 'Tag', ['ROI_1_', num2str(volhandle)]);
        item2_1a = uimenu(item2, 'Label', 'Set selection', 'Callback', ...
                          sprintf('%s(''context_selection'', %d, ''set'');', ...
                                  mfilename, volhandle), ...
                          'Tag',sprintf('ROI_SELECTION_%d', volhandle), ...
                          'Checked','on');
        item2_1b = uimenu(item2, 'Label', 'Clear selection', 'Callback', ...
                          sprintf('%s(''context_selection'', %d,''clear'');', ...
                                  mfilename, volhandle), ...
                          'Tag', sprintf('ROI_SELECTION_%d', volhandle));
        item3 = uimenu(item0, 'Label', 'Box size', 'Callback', ...
                       sprintf('%s(''context_box'', %d);', mfilename, volhandle), ...
                       'Visible', 'off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item4 = uimenu(item0, 'Label', 'Polygon slices', 'Callback', ...
                       sprintf('%s(''context_polyslices'', %d);', mfilename, volhandle), ...
                       'Visible','off', 'Tag', sprintf('ROI_1_%d', volhandle));
        item5 = uimenu(item0, 'Label', 'Cluster size', 'Callback', ...
                       sprintf('%s(''context_csize'', %d);', mfilename, volhandle), ...
                       'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item6 = uimenu(item0, 'Label', 'Erosion/Dilation threshold', 'Callback', ...
                       sprintf('%s(''context_erothresh'', %d);', mfilename, volhandle), ...
                       'Visible', 'off', 'Tag', sprintf('ROI_1_%d', volhandle));
        item7 = uimenu(item0, 'Label', 'Edit tool', ...
                       'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle), ...
                       'Separator','on');
        item7_1a = uimenu(item7, 'Label', 'Box tool', 'Callback', ...
                          sprintf('%s(''context_edit'', %d,''box'');', mfilename, volhandle), ...
                          'Tag',sprintf('ROI_EDIT_%d', volhandle), 'Checked','on');
        item7_1b = uimenu(item7, 'Label', 'Polygon tool', 'Callback', ...
                          sprintf('%s(''context_edit'', %d,''poly'');', mfilename, volhandle), ...
                          'Tag',sprintf('ROI_EDIT_%d', volhandle));
        item7_1c = uimenu(item7, 'Label', 'Connected cluster', 'Callback', ...
                          sprintf('%s(''context_edit'', %d,''connect'');', mfilename, volhandle), ...
                          'Tag',sprintf('ROI_EDIT_%d', volhandle));
        item8 = uimenu(item0, 'Label', 'Threshold', 'Callback', ...
                       sprintf('%s(''context_thresh'', %d);', mfilename, volhandle), ...
                       'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item10 = uimenu(item0, 'Label', 'Cleanup clusters', 'Callback', ...
                        sprintf('%s(''cleanup'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item11 = uimenu(item0, 'Label', 'Erode/Dilate', 'Callback', ...
                        sprintf('%s(''erodilate'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item12 = uimenu(item0, 'Label', 'Invert', 'Callback', ...
                        sprintf('%s(''invert'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item13 = uimenu(item0, 'Label', 'Clear', 'Callback', ...
                        sprintf('%s(''clear'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item14 = uimenu(item0, 'Label', 'Add ROI from file(s)', 'Callback', ...
                        sprintf('%s(''addfile'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item15 = uimenu(item0, 'Label', 'Save', 'Callback', ...
                        sprintf('%s(''save'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle),...
                        'Separator','on');
        item16 = uimenu(item0, 'Label', 'Save As', 'Callback', ...
                        sprintf('%s(''saveas'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item17 = uimenu(item0, 'Label', 'Quit', 'Callback', ...
                        sprintf('%s(''context_quit'', %d);', mfilename, volhandle), ...
                        'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        item18 = uimenu(item0, 'Label', 'Help', 'Callback', ...
                        sprintf('spm_help(''%s'');', mfilename));
        % add some stuff outside ROI tool menu
        iorient = findobj(get(st.vols{volhandle}.ax{1}.cm,'Children'), 'flat', 'Label', 'Orientation');
        item19 =  uimenu(iorient, 'Label', 'ROI Space', 'Callback', ...
                         sprintf('%s(''context_space'', %d);', mfilename, volhandle), ...
                         'Visible','off', 'Tag',sprintf('ROI_1_%d', volhandle));
        return;

    case 'context_init'
        Finter = spm_figure('FindWin', 'Interactive');
        spm_input('!DeleteInputObj',Finter);
        loadasroi = spm_input('Initial ROI','!+1','m',{'ROI image',...
                            'ROI space definition', 'SPM result'},[1 0 2],1);
        xyz = [];
        switch loadasroi
            case 1,
                imfname = spm_select(1, 'image', 'Select ROI image');
            case 0,
                imfname = spm_select(1, 'image', 'Select image defining ROI space');
            case 2,
                [SPM, xSPM] = spm_getSPM;
                xyz = xSPM.XYZ;
                imfname = SPM.Vbeta(1).fname;
        end;
        spm_input('!DeleteInputObj',Finter);
        feval('spm_ov_roi','init',volhandle,imfname,loadasroi,xyz);
        return;

    case 'context_selection'
        st.vols{volhandle}.roi.mode = varargin{3};
        obj = findobj(0, 'Tag', ['ROI_SELECTION_', num2str(volhandle)]);
        set(obj, 'Checked', 'off');
        set(gcbo, 'Checked', 'on');
        return;

    case 'context_edit'
        st.vols{volhandle}.roi.tool = varargin{3};
        obj = findobj(0, 'Tag', ['ROI_EDIT_', num2str(volhandle)]);
        set(obj, 'Checked', 'off');
        set(gcbo, 'Checked', 'on');
        return;

    case 'context_box'
        Finter = spm_figure('FindWin', 'Interactive');
        spm_input('!DeleteInputObj',Finter);
        box = spm_input('Selection size {vx vy vz}','!+1','e', ...
                        num2str(st.vols{volhandle}.roi.box), [1 3]);
        spm_input('!DeleteInputObj',Finter);
        st.vols{volhandle}.roi.box = box;
        return;

    case 'context_polyslices'
        Finter = spm_figure('FindWin', 'Interactive');
        spm_input('!DeleteInputObj',Finter);
        polyslices = spm_input('Polygon: slices around current slice','!+1','e', ...
                               num2str(st.vols{volhandle}.roi.polyslices), [1 1]);
        spm_input('!DeleteInputObj',Finter);
        st.vols{volhandle}.roi.polyslices = polyslices;
        return;

    case 'context_csize'
        Finter = spm_figure('FindWin', 'Interactive');
        spm_input('!DeleteInputObj',Finter);
        csize = spm_input('Minimum cluster size (#vx)','!+1','e', ...
                          num2str(st.vols{volhandle}.roi.csize), [1 1]);
        spm_input('!DeleteInputObj',Finter);
        st.vols{volhandle}.roi.csize = csize;
        return;

    case 'context_erothresh'
        Finter = spm_figure('FindWin', 'Interactive');
        spm_input('!DeleteInputObj',Finter);
        erothresh = spm_input('Erosion/Dilation threshold','!+1','e', ...
                              num2str(st.vols{volhandle}.roi.erothresh), [1 1]);
        spm_input('!DeleteInputObj',Finter);
        st.vols{volhandle}.roi.erothresh = erothresh;
        return;

    case 'context_space'
        spm_orthviews('space', volhandle, ...
            st.vols{volhandle}.roi.Vroi.mat, ...
            st.vols{volhandle}.roi.Vroi.dim(1:3));
        iorient = get(findobj(get(st.vols{volhandle}.ax{1}.cm,'Children'), ...
            'flat', 'Label', 'Orientation'),'children');
        if iscell(iorient)
            iorient = cell2mat(iorient);
        end
        set(iorient, 'Checked', 'Off');
        ioroi = findobj(iorient, 'Label','ROI Space');
        set(ioroi, 'Checked', 'On');
        return;

    case 'context_thresh'
        Finter = spm_figure('FindWin', 'Interactive');
        spm_input('!DeleteInputObj',Finter);
        thresh = spm_input('Threshold  {min max}','!+1','e', ...
                           num2str(st.vols{volhandle}.roi.thresh), [1 2]);
        spm_input('!DeleteInputObj',Finter);
        st.vols{volhandle}.roi.thresh = thresh;
        feval('spm_ov_roi', 'thresh', volhandle);
        return;

    case 'context_quit'
        obj = findobj(0, 'Tag',  sprintf('ROI_1_%d', volhandle));
        set(obj, 'Visible', 'off');
        obj = findobj(0, 'Tag',  sprintf('ROI_0_%d', volhandle));
        set(obj, 'Visible', 'on');
        spm_orthviews('rmblobs', volhandle);
        for k=1:3
            set(st.vols{volhandle}.ax{k}.ax,'ButtonDownFcn', st.vols{volhandle}.roi.cb{k});
        end;
        st.vols{volhandle} = rmfield(st.vols{volhandle}, 'roi');
        spm_orthviews('redraw');
        return;

    case 'bdfcn'
        if strcmpi(get(gcf,'SelectionType'), 'extend')
            spm_orthviews('roi','edit',volhandle);
        else
            for k=1:3
                if isequal(st.vols{volhandle}.ax{k}.ax, gca)
                    break;
                end
            end
            feval(st.vols{volhandle}.roi.cb{k});
        end
    otherwise
        fprintf('spm_orthviews(''roi'', ...): Unknown action %s', cmd);
        return;
end;

if update_roi
    if ~isempty(tochange) % change state according to mode
        if strcmp(st.vols{volhandle}.roi.mode,'set')
            toset = tochange;
        else
            toclear = tochange;
        end;
    end;
    % clear first, then set (needed for connect operation)
    if ~isempty(toclear)
        itoclear = sub2ind(st.vols{volhandle}.roi.Vroi.dim(1:3), ...
                           toclear(1,:), toclear(2,:), toclear(3,:));
        st.vols{volhandle}.roi.roi(itoclear) = false;
        if ~isempty(st.vols{volhandle}.roi.xyz)
            st.vols{volhandle}.roi.xyz = setdiff(st.vols{volhandle}.roi.xyz',toclear','rows')';
        else
            st.vols{volhandle}.roi.xyz = [];
        end;
    end;

    if ~isempty(toset)
        % why do we need this round()?? I don't know, but Matlab thinks
        % it is necessary
        itoset = round(sub2ind(st.vols{volhandle}.roi.Vroi.dim(1:3), ...
                               toset(1,:), toset(2,:), toset(3,:)));
        st.vols{volhandle}.roi.roi(itoset) = true;
        if ~isempty(st.vols{volhandle}.roi.xyz)
            st.vols{volhandle}.roi.xyz = union(st.vols{volhandle}.roi.xyz',toset','rows')';
        else
            st.vols{volhandle}.roi.xyz = toset;
        end;
    end;

    if isfield(st.vols{volhandle}, 'blobs')
        nblobs=length(st.vols{volhandle}.blobs);
        if nblobs>1
            blobstmp(1:st.vols{volhandle}.roi.hroi-1) = st.vols{volhandle}.blobs(1:st.vols{volhandle}.roi.hroi-1);
            blobstmp(st.vols{volhandle}.roi.hroi:nblobs-1) = st.vols{volhandle}.blobs(st.vols{volhandle}.roi.hroi+1:nblobs);
            st.vols{volhandle}.blobs=blobstmp;
        else
            if isempty(st.vols{volhandle}.roi.hframe) % save frame
                st.vols{volhandle}=rmfield(st.vols{volhandle},'blobs');
            end;
        end;
    end;
    if isfield(st.vols{volhandle}, 'blobs')
        st.vols{volhandle}.roi.hroi = numel(st.vols{volhandle}.blobs)+1;
    else
        st.vols{volhandle}.roi.hroi = 1;
    end;
    if isempty(st.vols{volhandle}.roi.xyz)  % initialised with empty roi
        spm_orthviews('addcolouredblobs', volhandle, ...
                      [1; 1; 1], 0, st.vols{volhandle}.roi.Vroi.mat,[1 3 1]);
    else
        spm_orthviews('addcolouredblobs', volhandle, ...
                      st.vols{volhandle}.roi.xyz, ones(size(st.vols{volhandle}.roi.xyz,2),1), ...
                      st.vols{volhandle}.roi.Vroi.mat,[1 3 1]); % use color that is more intense than standard rgb range
    end;
    st.vols{volhandle}.blobs{st.vols{volhandle}.roi.hroi}.max=2;
    spm_orthviews('redraw');
end;

spm('pointer','arrow');

function varargout = stack(cmd, varargin)

switch cmd
    case 'init'
        varargout{1}.st = cell(varargin{1},1);
        varargout{1}.top= 0;

    case 'isempty'
        varargout{1} = (varargin{1}.top==0);

    case 'push'
        stck = varargin{1};
        if (stck.top < size(stck.st,1))
            stck.top = stck.top + 1;
            stck.st{stck.top} = varargin{2};
            varargout{1}=stck;
        else
            error('Stack overflow\n');
        end;

    case 'pop'
        if stack('isempty',varargin{1})
            error('Stack underflow\n');
        else
            varargout{2} = varargin{1}.st{varargin{1}.top};
            varargin{1}.top = varargin{1}.top - 1;
            varargout{1} = varargin{1};
        end;
end;

author: spm


text:
function [logStruct] = getLogInfo(logFilePath,varargin)
	% Opens a log file (e.g. Inscopix recording file) and outputs a structure containing the field information for the structure.
	% Biafra Ahanonu
	% started: 2014.01.30
	% based on R code by Biafra Ahanonu started: 2013.09.18
	% inputs
		%
	% outputs
		%

	[logStruct] = ciapkg.io.getLogInfo(logFilePath,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% =========================
%  Robust Bayesian allocation in the stock market, as described in
%  Meucci, A., (2005) "Robust Bayesian Allocation"
%
%  Most recent version of article and code available at http://symmys.com/node/102
%  =========================

clear; clc; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inputs
p_m=.1; % robustness parameter location
p_s=.1; % robustness parameter scatter
load SectorsSnP500

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute weekly returns
Ps=P(1:5:end,:);
R=Ps(2:end,:)./Ps(1:end-1,:)-1;
Dates_P=DP(1:5:end);
Dates_R=Dates_P(2:end);
[Ttot,N]=size(R);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimation
W=52; % rolling estimation period
NumPortf=10;
Ret_hat=[];
Ret_rB=[];
Dates=[];
for t=W+1:Ttot-1
    Ttot-t+2
    Rets=R(t-W:t,:);

    % sample estimate
    m_hat=mean(Rets)';
    S_hat=cov(Rets);
    [de_hat,ds_hat,w_hat] = EfficientFrontier(NumPortf, S_hat, m_hat);

    % Bayesian prior
    S0=diag(diag(S_hat));
    m0=.5*S0*ones(N,1)/N;
    T=size(Rets,1);
    T0=2*T;
    nu0=2*T;

    % Bayesian posterior parameters
    T1=T+T0;
    m1=1/T1*(m_hat*T+m0*T0);
    nu1=T+nu0;
    S1=1/nu1*( S_hat*T + S0*nu0 + (m_hat-m0)*(m_hat-m0)'/(1/T+1/T0)  );
    [d,d,w1] = EfficientFrontier(NumPortf, S1, m1);

    % robustness parameters
    q_m2=chi2inv(p_m,N);
    g_m=sqrt(q_m2/T1*nu1/(nu1-2));
    q_s2=chi2inv(p_s,N*(N+1)/2);
    PickVol=round(.8*NumPortf);
    v=(ds_hat(PickVol))^2;
    g_s=v/(  nu1/(nu1+N+1)+sqrt( 2*nu1*nu1*q_s2/((nu1+N+1)^3)));

    Target=[];

    wu=w_hat(PickVol,:)';
    Ret_hat=[Ret_hat R(t+1,:)*wu];

    for k=1:NumPortf-1
        NewTarget=-(10^10);
        if wu'*S1*wu <= g_s
            NewTarget = m1'*wu-g_m*sqrt(wu'*S1*wu);
        end
        Target=[Target NewTarget];
    end

    [Best,k]=max(Target);
    wu=w1(k,:)';
    Ret_rB=[Ret_rB R(t+1,:)*wu];

    Dates=[Dates Dates_R(t+1)];
end

NAV_hat=cumprod(1+Ret_hat);
NAV_rB=cumprod(1+Ret_rB);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plots

figure
subplot(2,1,1)
h=plot(Dates,Ret_hat);
xlim([Dates(1) Dates(end)])
YLim=get(gca,'ylim');
Datetick('x','mmmyy','keeplimits','keepticks');
grid on
subplot(2,1,2)
h=plot(Dates,Ret_rB);
set(gca,'ylim',YLim,'xlim',[Dates(1) Dates(end)]);
Datetick('x','mmmyy','keeplimits','keepticks');
grid on

figure
subplot(2,1,1)
h=plot(Dates,NAV_hat);
xlim([Dates(1) Dates(end)])
YLim=get(gca,'ylim');
Datetick('x','mmmyy','keeplimits','keepticks');
grid on
subplot(2,1,2)
h=plot(Dates,NAV_rB);
set(gca,'ylim',YLim,'xlim',[Dates(1) Dates(end)]);
Datetick('x','mmmyy','keeplimits','keepticks');
grid on

author: Sable

repo: mcbench-benchmarks


text:
targetWorkCount = 100;
barWidth= int32( targetWorkCount/3 );
p =  TimedProgressBar( targetWorkCount, barWidth, ...
    'Computing, wait for ', ', completed ', 'Concluded in ' );
parfor i=1:targetWorkCount      % could be just a for cycle
    pause(rand);                % Replace with real code
    p.progress;                 % Also percent = p.progress;
end
p.stop;                         % Also percent = p.stop;

repo: CategoryShapes


text:
%random_int Returns random integer is specified range.
%
%  y = random_int(from, to)

function y = random_int(from, to)

y = floor(from + (1 + to - from)*rand);

author: strawlab

repo: MultiCamSelfCal


text:
function title = p00_title ( problem, option )

%*****************************************************************************80
%
%% P00_TITLE sets the title for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Output, string TITLE, the title of the problem.
%    TITLE will never be longer than 80 characters.
%
  if (  problem == 1 )
    title = p01_title ( option );
  elseif ( problem == 2 )
    title = p02_title ( option );
  elseif ( problem == 3 )
    title = p03_title ( option );
  elseif ( problem == 4 )
    title = p04_title ( option );
  elseif ( problem == 5 )
    title = p05_title ( option );
  elseif ( problem == 6 )
    title = p06_title ( option );
  elseif ( problem == 7 )
    title = p07_title ( option );
  elseif ( problem == 8 )
    title = p08_title ( option );
  elseif ( problem == 9 )
    title = p09_title ( option );
  elseif ( problem == 10 )
    title = p10_title ( option );
  elseif ( problem == 11 )
    title = p11_title ( option );
  elseif ( problem == 12 )
    title = p12_title ( option );
  elseif ( problem == 13 )
    title = p13_title ( option );
  elseif ( problem == 14 )
    title = p14_title ( option );
  elseif ( problem == 15 )
    title = p15_title ( option );
  elseif ( problem == 16 )
    title = p16_title ( option );
  elseif ( problem == 17 )
    title = p17_title ( option );
  elseif ( problem == 18 )
    title = p18_title ( option );
  elseif ( problem == 19 )
    title = p19_title ( option );
  elseif ( problem == 20 )
    title = p20_title ( option );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_TITLE - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized problem number = %d\n', problem );
    error ( 'P00_TITLE - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [tm,im] = mapvert(tr,pi)
%MAPVERT find the tree-to-vertex mappings.
%   [TM,IM] = MAPVERT(TR,PI) returns the tree-to-vertex and
%   vertex-to-tree mappings for a given aabb-tree TR and a
%   collection of query vertices PI.
%
%   The tree-to-item mapping TM is a structure representing
%   the intersection of the items PI with the tree TR. TM.II
%   is an M-by-1 array of tree indices and TM.LL is an
%   M-by-1 cell array of item lists. Specifically, items in
%   the list TM.LL{JJ} intersect with the node TM.II(JJ).
%
%   The item-to-tree mapping IM is a structure representing
%   the inverse mapping. IM.II is an N-by-1 array of item
%   indices and IM.LL is an N-by-1 cell array of node lists.
%   Specifically, nodes in the list IM.LL{JJ} intersect with
%   the item IM.II(JJ).
%
%   See also QUERYSET, MAPRECT, MAKETREE

%   Darren Engwirda : 2014 --
%   Email           : de2363@columbia.edu
%   Last updated    : 06/04/2017

%----------------------- call SCANTREE to do the actual work
    if (nargout == +1)
       [tm   ] = scantree(tr,pi,@partvert);
    else
       [tm,im] = scantree(tr,pi,@partvert);
    end

end

function [j1,j2] = partvert(pi,b1,b2)
%PARTVERT partition points between boxes B1,B2 for SCANTREE.

    j1 = true(size(pi,1),1);
    j2 = true(size(pi,1),1);

    nd = size(b1,2) / +2;

    for ax = +1 : nd
%--------------- remains TRUE if inside bounds along axis AX
    j1 = j1 & pi(:,ax)>=b1(ax+nd*0) ...
            & pi(:,ax)<=b1(ax+nd*1) ;

%--------------- remains TRUE if inside bounds along axis AX
    j2 = j2 & pi(:,ax)>=b2(ax+nd*0) ...
            & pi(:,ax)<=b2(ax+nd*1) ;
    end

end

author: CHLNDDEV


text:
% RMART - Remove eye artifacts from EEG data using regression with
%           multiple time lags. Each channel is first made mean-zero.
%           After JL Kenemans et al., Psychophysiology 28:114-21, 1991.
%
%   Usage: >> rmart('datafile','outfile',nchans,chanlist,eogchan,[threshold])
% Example: >> rmart('noisy.floats','clean.floats',31,[2:31],7)
%
%   Input:   datafile - input float data file, multiplexed by channel
%            outfile  - name of output float data file
%            nchans   - number of channels in datafile
%            chanlist - indices of EEG channel(s) to process (1,...,nchans)
%            eogchan  - regressing channel indices(s) (1,...,nchans)
%            threshold- abs threshold value to trigger regression {def|0 -> 80}
%
% Output: Writes [length(chanlist),size(data,2)] floats to 'outfile'
%
% Note: Regression epoch length and number of lags are set in the script.
%       Some machines may require a new byte_order value in the script.
%       note that RUNICA -> ICAPROJ should give better results! See
%       Jung et al., Psychophysiology 111:1745-58, 2000.
%
% Author: Tzyy-Ping Jung, SCCN/INC/UCSD, La Jolla, 1997

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

% 2-22-97  Tzyy-Ping Jung  CNL/Salk Institute, La Jolla, CA
% 2-24-97  Formatted for ICA package release -Scott Makeig
% 12-10-97 Changed name from rmartifact to rmart for toolbox inclusion -sm
% 12-11-97 Adapted to read/write a float matrix -sm & sw
% 09-14-00 Added comments and help -sm
% 01-25-02 reformated help & license -ad

function rmart(datafile,outfile,nchans,chanlist,eogchan,threshold)

if nargin < 5
   help rmart
   return
end

%
% The following parameters may be fine-tuned for a data set
%
DEF_THRESHOLD = 80; % trigger regression on blocks exceeding this (default)
epoch         = 80; % remove artifacts in successive blocks of this length
nlags         = 40; % perform multiple regression filtering of this length

byte_order    = 'b';% (machine-dependent) byte order code for fopen();
MAKE_MEAN_ZERO = 1  % 1/0 flag removing mean offset from each channel

fprintf('Performing artifact regression on data in %s.\n',datafile);

if nargin<6
   threshold = 0;
end
if threshold == 0,
   threshold = DEF_THRESHOLD;
end
fprintf('Regression threshold %g.\n',threshold);

%
% Read the input data
%
[fid,msg]=fopen(datafile,'r',byte_order); % open datafile
if fid < 3,
  fprintf('rmart() - could not open data file: %s\n',msg);
  exit 1
end
data=(fread(fid,'float'))';
status=fclose('all');
if rem(length(data),nchans) == 0 % check length
   fprintf('rmart() - data length not divisible by %d chans.\n',nchans);
    return
end

data = reshape(data,nchans,length(data)/nchans);
[chans,frames] = size(data);
fprintf('Data of size [%d,%d] read.\n',chans,frames);
eog = data(eogchan,:);
data = data(chanlist,:);
procchans = length(chanlist);

fprintf('Regression epoch length %d frames.\n',epoch);
fprintf('Using %d regression lags.\n',nlags);
if length(eogchan)> 1
  fprintf('Processing %d of %d channels using %d EOG channels.\n',...
                        procchans,chans,length(eogchan));
else
  fprintf('Processing %d of %d channels using EOG channel %d.\n',...
                        procchans,chans,eogchan);
end

%
% Process the data
%
for i=1:procchans
  chan = chanlist(i);
  idx=[];
  frame=1+epoch/2+nlags/2;
  if MAKE_MEAN_ZERO
    data(chan,:) = data(chan,:) - mean(data(chan,:)); % make mean-zero
  end

  % Search the EOG & EEG records for values above threshold,
  % Selected frame numbers are registered in the variable "idx".
  % The entries in "idx" are at least epoch apart to avoid double
  % compensation (regression) on the same portion of the EEG data.

  while frame <= length(eog)-epoch/2-nlags/2,  % foreach epoch in channel
    stop = min(frame+epoch-1,eogframes);
    tmp= ...
        find( abs(eog(frame:stop)) >= threshold ...
             | abs(data(chan,frame:stop)) >= threshold);
                                           % find beyond-threshold values
    if length(tmp) ~= 0
      mark = tmp(1)+frame-1;
      if  length(idx) ~= 0
        if mark-idx(length(idx))  < epoch,
           idx=[idx idx(length(idx))+epoch]; % To guarantee idx(i) & idx(i-1)
                                             % are at least EPOCH points apart
           frame = idx(length(idx))+epoch/2;
        else
           idx=[idx mark];
           frame = mark + epoch/2;
        end
      else
        idx=[idx mark];
        frame = mark + epoch/2;
      end
    else
      frame=frame+epoch;
    end
  end % while

  % For each registered frame, take "epoch" points
  % surrounding it from the EEG, and "epoch + lag" points
  % from the EOG channel. Then perform multivariate
  % linear regression on EEG channel.

  for j=1:length(idx);
     art=ones(1,epoch);
     eogtmp=eog(idx(j)-epoch/2-nlags/2:idx(j)+epoch/2-1+nlags/2);

     % Collect EOG data from lag/2 points before to lag/2 points
     % after the regression window.
     for J=nlags:-1:1,
       art=[art ; eogtmp(J:J+epoch-1)];
     end
     eegtmp=data(chan,idx(j)-epoch/2:idx(j)+epoch/2-1);

     eegeog=eegtmp*art';        % perform the regression here
     eogeog=art*art';
     b=eegeog/eogeog;
     eegtmp=eegtmp-b*art;
     data(chan,idx(j)-epoch/2:idx(j)+epoch/2-1)=eegtmp;
   end % j
end % i

%
% Write output file
%
[fid,msg]=fopen(outfile,'w',byte_order);
if fid < 3
  fprintf('rmart() - could not open output file: %s\n',msg);
  return
end
count = fwrite(fid,data,'float');
if count == procchans*frames,
  fprintf('Output file "%s" written, size = [%d,%d] \n\n',...
             outfile,procchans,frames);
else
  fprintf('rmart(): Output file "%s" written, SIZE ONLY [%d,%g]\n',...
             outfile,procchans,count/procchans);
end
fclose('all');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% mario_huggins.m - Nick Clark 28/9/2007
%
% USE HEADPHONES
% WAIT A FEW SECONDS FOR PROCESSING
% LISTEN CAREFULLY AND ENJOY!
%
% Just run this script with the huggins.m in the same directory or in your
% path.  This is a cool litle demo to show how huggins binaural pitch can
% be made into chords. Thanks to James Humes on matlabcentral for inspiring
% me to do this to the mario theme song . . man - this is probably the most
% tragic way I have ever spent a Friday night.
%
% This demo only uses the alto and tenor parts as Huggins pitch is a
% pretty weak percept at the best of times.  You'll need a pretty good set
% of headphones to hear this too - it will not work properly over
% loudspeakers.

t = 0.15; %single beat duration
base_pitch = 440; %tuning parameter 440Hz default
introONLY = 0; %set to one if you only want top hear the intro
sr = 16e3;

%%%%% Notes and timing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The transcription here was originally done by Stewart Bozarth and the
% painstaking programming of the score into MATLAB code was taken from
% James Humes' fileexchange submission
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if introONLY
    keysT = [ 56 56 0 56 0 52 56 0 59 0 0 47 0 0 ];
    keysA=  [ 46 46 0 46 0 46 46 0 51 0 0 47 0 0 ];
    durA =  [ t t t t t t t t t t 2*t t t 2*t];
else %whole song
    keysT = [  56 56 0 56 0 52 56 0 59 0 0 47 0 0 52 0 47 0 44 0 0 49 0 51 0 50 49 0 47 0 56 0 59 0 61 0 57 59 0 56 0 52 54 51 0 52 0 47 0 44 0 0 49 0 51 0 50 49 0 47 0 56 0 59 0 61 0 57 59 0 56 0 52 54 51 0 0 59 58 57 55 0 56 0 48 49 52 0 49 52 54 0 59 58 57 55 0 56 0 64 0 64 64 0 0 0 59 58 57 55 0 56 0 48 49 52 0 49 52 54 0 55 0 0 54 0 52 0 0 0 0 59 58 57 55 0 56 0 48 49 52 0 49 52 54 0 59 58 57 55 0 56 0 64 0 64 64 0 0 0 59 58 57 55 0 56 0 48 49 52 0 49 52 54 0 55 0 0 54 0 52 0 0 0 52 52 0 52 0 52 54 0 56 52 0 49 47 0 0 52 52 0 52 0 52 54 56 0 0 52 52 0 52 0 52 54 0 56 52 0 49 47 0 0 56 56 0 56 0 52 56 0 59 0 0 47 0 0 56 52 0 47 0 48 0 49 57 0 57 49 0 0 51 0 61 0 61 0 61 0 59 0 57 0 56 52 0 49 47 0 0 56 52 0 47 0 48 0 49 57 0 57 49 0 0 47 57 0 57 57 0 56 0 54 0 52 0 0 0 52 0 47 0 44 0 49 51 49 48 50 48 47];
    keysA = [  46 46 0 46 0 46 46 0 51 0 0 47 0 0 44 0 40 0 35 0 0 40 0 42 0 41 40 0 40 0 47 0 51 0 52 0 49 51 0 49 0 56 45 42 0 44 0 40 0 35 0 0 40 0 42 0 41 40 0 40 0 47 0 51 0 52 0 49 51 0 49 0 56 45 42 0 0 56 55 54 51 0 52 0 44 45 47 0 40 44 45 0 56 55 54 51 0 52 0 59 0 59 59 0 0 0 56 55 54 51 0 52 0 44 45 47 0 40 44 45 0 48 0 0 45 0 44 0 0 0 0 56 55 54 51 0 52 0 44 45 47 0 40 44 45 0 56 55 54 51 0 52 0 59 0 59 59 0 0 0 56 55 54 51 0 52 0 44 45 47 0 40 44 45 0 48 0 0 45 0 44 0 0 0 48 48 0 48 0 48 50 0 47 44 0 44 40 0 0 48 48 0 48 0 48 50 47 0 0 48 48 0 48 0 48 50 0 47 44 0 44 40 0 0 46 46 0 46 0 46 46 0 51 0 0 47 0 0 52 49 0 44 0 44 0 45 52 0 52 45 0 0 47 0 57 0 57 0 57 0 56 0 54 0 52 49 0 45 44 0 0 52 49 0 44 0 44 0 45 52 0 52 45 0 0 47 54 0 54 54 0 52 0 51 0 47 44 0 44 40 0 0 44 0 40 49 51 49 48 50 48 47];
    durA =  [  t t t t t t t t t t 2*t t t 2*t t 2*t t 2*t t t t t t t t t t t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t t t t t t t t t t t 2*t t 2*t t 2*t t t t t t t t t t t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t t t t t t t t t t t 2*t 2*t t t t t t t t t t t t t t t 2*t t t t t t t t t t t t t 2*t 2*t t t t t t t t t t t t t t t 2*t t t t t 2*t t t 2*t 4*t 2*t t t t t t t t t t t t t t t 2*t t t t t t t t t t t t t 2*t 2*t t t t t t t t t t t t t t t 2*t t t t t 2*t t t 2*t 4*t t t t t t t t t t t t t t t 2*t t t t t t t t t 4*t 4*t t t t t t t t t t t t t t t 2*t t t t t t t t t t t 2*t t t 2*t t t t t 2*t t t t t t t t t 2*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t t t t t t t 2*t t t t t 2*t t t t t t t t t 2*t t t t t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t (2/3)*t t t t t t t 2*t t 2*t t 2*t t t 4*t 4*t t 6*t];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert key values and make song (may take some time - please wait)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
freqA = base_pitch*2.^((keysA-49)/12);
freqT = base_pitch*2.^((keysT-49)/12);

song = [0 0];
for n = 1:length(keysT)
    song = [song; huggins([freqA(n) freqT(n)],sr,durA(n))];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add a 2.5kHz low-pass filter to make things a little less painful
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[b,a] = butter(4,2500*(sr^-1),'low');
song = filter(b,a,song);

soundsc(song, sr)


text:
% Minimum-phase spectrum of a given spectrum using the complex cepstrum
%
% Description
%  Compute the minimum-phase spectrum of a given spectrum X by dropping the
%  negative quefrencies of X.
%
% Input
%  X   : The spectrum to convert (full DFT length)
%
% Output
%  M   : The corresponding minimum-phase spectrum
%  lM  : The log minimum-phase spectrum
%
% References
%  [1] Alan V. Oppenheim and Ronald W. Schafer, "Digital Signal Processing",
%      Prentice-Hall, 2nd edition, 1978.
%      Note that this 2nd edition contains a chapter about complex cepstrum which
%      has been removed in the 3rd edition (and possibly came back in the 4th ed.)
%
% Copyright (c) 2011 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [M lM] = spec2minphasespec(X)

    X = X(:).';

    rcc = ifft(log(abs(X))); % Compute the real cepstrum

    if mod(length(X),2)==0
        % For even DFT length
        lM = fft([rcc(1),2*rcc(2:end/2),rcc(end/2+1)], length(X));% [1]
    else
        % For odd DFT length
        lM = fft([rcc(1),2*rcc(2:(end-1)/2+1)], length(X));       % [1]
    end

    M = exp(lM);

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [rmse_train, rmse_test] = aistats2012_mohsst5_rmse(F, testset, seed)

if nargin < 3
  [Ytrain, Ytest, Itrain, Itest, data, sea] = aistats2012_mohsst5_sets(testset);
else
  [Ytrain, Ytest, Itrain, Itest, data, sea] = aistats2012_mohsst5_sets(testset, ...
                                                    seed);
end

[M,N] = size(F);

% Error measures
[LON,LAT] = meshgrid(data.longitude, data.latitude);
weights = cosd(LAT(:));
weights = weights(sea);
weights = repmat(weights, [1, N]);
rmse_train = rmsew(Ytrain(Itrain)-F(Itrain),weights(Itrain));
rmse_test = rmsew(Ytest(Itest)-F(Itest),weights(Itest));
rmse_zero = rmsew(Ytest(Itest),weights(Itest));

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [F, TimeVector,DisplayUnits] = in_fread(sFile, ChannelMat, iEpoch, SamplesBounds, iChannels, ImportOptions)
% IN_FREAD: Read a block a data in any recordings file previously opened with in_fopen().
%
% USAGE:  [F, TimeVector, DisplayUnits] = in_fread(sFile, ChannelMat, iEpoch, SamplesBounds, iChannels, ImportOptions);
%         [F, TimeVector, DisplayUnits] = in_fread(sFile, ChannelMat, iEpoch, SamplesBounds, iChannels);                 : Do not apply any pre-preprocessings
%         [F, TimeVector, DisplayUnits] = in_fread(sFile, ChannelMat, iEpoch, SamplesBounds);                            : Read all channels
%
% INPUTS:
%     - sFile         : Structure for importing files in Brainstorm. Created by in_fopen()
%     - iEpoch        : Indice of the epoch to read (only one value allowed)
%     - SamplesBounds : [smpStart smpStop], First and last sample to read in epoch #iEpoch
%     - iChannels     : Array of indices of the channels to import
%     - ImportOptions : Structure created by interface window panel_import_data.m  (look in db_template.m for a description of all the fields).
%
% OUTPUTS:
%     - F          : [nChannels x nTimes], block of recordings
%     - TimeVector : [1 x nTime], time values in seconds
%     - DisplayUnits : char, unit of the recording

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

%% ===== PARSE INPUTS =====
if (nargin < 6)
    ImportOptions = [];
end
if (nargin < 5)
    iChannels = [];
end
TimeVector = [];
% Read channel ranges for faster access
isChanRange = ismember(sFile.format, {'CTF', 'CTF-CONTINUOUS', 'KDF', 'EEG-EDF', 'EEG-BDF', 'BST-BIN', 'EEG-CURRY', 'EEG-DELTAMED', 'EEG-COMPUMEDICS-PFS', 'EEG-MICROMED', 'EEG-NEURONE', 'EEG-NK', 'EEG-OEBIN'});
if isChanRange
    if isempty(iChannels)
        ChannelRange = [];
        iChanRemove = [];
    else
        if ~isequal(iChannels, sort(iChannels))
            error('You need to sort the channels indices before calling in_fread().');
        end
        ChannelRange = [iChannels(1), iChannels(end)];
        iChanRemove = setdiff(ChannelRange(1):ChannelRange(2), iChannels) - ChannelRange(1) + 1;
    end
end

%% ===== OPEN FILE =====
% Open file (for some formats, it is open in the low-level function)
if ismember(sFile.format, {'FIF', 'CTF', 'KIT', 'RICOH', 'BST-DATA', 'SPM-DAT', 'EEG-ANT-CNT', 'EEG-AXION', 'EEG-EEGLAB', 'EEG-GTEC', 'EEG-NEURONE', 'EEG-NEURALYNX', 'EEG-NICOLET', 'EEG-BLACKROCK', 'EEG-RIPPLE', 'EYELINK', 'NIRS-BRS', 'EEG-EGI-MFF', 'MNE-PYTHON'})
    sfid = [];
else
    sfid = fopen(sFile.filename, 'r', sFile.byteorder);
end

% Check whether optional field precision is available
if ~isempty(ImportOptions) && isfield(ImportOptions, 'Precision')
    precision = ImportOptions.Precision;
else
    precision = [];
end

%% ===== READ RECORDINGS BLOCK =====
DisplayUnits = [];
switch (sFile.format)
    case 'FIF'
        [F,TimeVector] = in_fread_fif(sFile, iEpoch, SamplesBounds, iChannels);
    case {'CTF', 'CTF-CONTINUOUS'}
        isContinuous = strcmpi(sFile.format, 'CTF-CONTINUOUS');
        F = in_fread_ctf(sFile, iEpoch, SamplesBounds, ChannelRange, isContinuous);
    case '4D'
        F = in_fread_4d(sFile, sfid, iEpoch, SamplesBounds, iChannels);
    case 'KIT'
        F = in_fread_kit(sFile, iEpoch, SamplesBounds, iChannels);
    case 'RICOH'
        F = in_fread_ricoh(sFile, iEpoch, SamplesBounds, iChannels);
    case 'KDF'
        F = in_fread_kdf(sFile, sfid, SamplesBounds, ChannelRange);
    case 'ITAB'
        F = in_fread_itab(sFile, sfid, SamplesBounds, iChannels);
    case 'MEGSCAN-HDF5'
        F = in_fread_megscan(sFile, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-ADICHT'
        F = in_fread_adicht(sFile, iEpoch, iChannels, SamplesBounds);
    case 'EEG-ANT-CNT'
        F = in_fread_ant(sFile, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-ANT-MSR'
        F = in_fread_msr(sFile, sfid, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-AXION'
        F = in_fread_axion(sFile, SamplesBounds, iChannels, precision);
    case 'EEG-BCI2000'
        F = in_fread_bci2000(sFile, SamplesBounds);
    case {'EEG-BLACKROCK', 'EEG-RIPPLE'}
        F = in_fread_blackrock(sFile, SamplesBounds, iChannels, precision);
    case 'EEG-BRAINAMP'
        F = in_fread_brainamp(sFile, sfid, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-CURRY'
        F = in_fread_curry(sFile, sfid, iEpoch, SamplesBounds, ChannelRange);
    case 'EEG-DELTAMED'
        F = in_fread_deltamed(sFile, sfid, SamplesBounds, ChannelRange);
    case 'EEG-COMPUMEDICS-PFS'
        F = in_fread_compumedics_pfs(sFile, sfid, SamplesBounds, ChannelRange);
    case {'EEG-EDF', 'EEG-BDF'}
        F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelRange);
    case 'EEG-EEGLAB'
        F = in_fread_eeglab(sFile, iEpoch, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-EGI-RAW'
        F = in_fread_egi(sFile, sfid, iEpoch, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-GTEC'
        F = in_fread_gtec(sFile, iEpoch, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-MANSCAN'
        F = in_fread_manscan(sFile, sfid, iEpoch, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-EGI-MFF'
        F = in_fread_mff(sFile, iEpoch, SamplesBounds, ImportOptions);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-MICROMED'
        F = in_fread_micromed(sFile, sfid, SamplesBounds, ChannelRange);
    case 'EEG-NEURALYNX'
        F = in_fread_neuralynx(sFile, SamplesBounds, iChannels);
    case 'EEG-NEURONE'
        F = in_fread_neurone(sFile, SamplesBounds, ChannelRange);
    case 'EEG-NEUROSCAN-CNT'
        F = in_fread_cnt(sFile, sfid, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-NEUROSCAN-EEG'
        F = in_fread_eeg(sFile, sfid, iEpoch, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-NEUROSCAN-AVG'
        F = in_fread_avg(sFile, sfid, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-NEUROSCOPE'
        F = in_fread_neuroscope(sFile, sfid, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'EEG-NICOLET'
        F = in_fread_nicolet(sFile, iEpoch, SamplesBounds, iChannels);
    case 'EEG-OEBIN'
        F = in_fread_oebin(sFile, sfid, SamplesBounds, ChannelRange);
    case 'EEG-NK'
        F = in_fread_nk(sFile, sfid, iEpoch, SamplesBounds, ChannelRange);
    case 'EEG-SMR'
        F = in_fread_smr(sFile, sfid, SamplesBounds, iChannels);
    case 'EEG-SMRX'
        F = in_fread_smrx(sFile, SamplesBounds, iChannels);
    case 'EYELINK'
        [F, TimeVector] = in_fread_eyelink(sFile, iEpoch, SamplesBounds, iChannels);
    case 'NIRS-BRS'
        F = in_fread_nirs_brs(sFile, SamplesBounds);
        if ~isempty(iChannels)
            F = F(iChannels,:);
        end
    case 'SPM-DAT'
        F = in_fread_spm(sFile, SamplesBounds, iChannels);
    case 'BST-BIN'
        F = in_fread_bst(sFile, sfid, SamplesBounds, ChannelRange);
    case 'BST-DATA'
        if ~isempty(SamplesBounds)
            fileSamples = round(sFile.prop.times * sFile.prop.sfreq);
            iTimes = (SamplesBounds(1):SamplesBounds(2)) - fileSamples(1) + 1;
        else
            iTimes = 1:size(sFile.header.F,2);
        end
        if isempty(iChannels)
            iChannels = 1:size(sFile.header.F,1);
        end
        F = sFile.header.F(iChannels, iTimes);
        % Load display units
        if ~isempty(sFile.filename)
            DataMat = load(sFile.filename, 'DisplayUnits');
            if isfield(DataMat, 'DisplayUnits') && ~isempty(DataMat.DisplayUnits)
                DisplayUnits = DataMat.DisplayUnits;
            end
        end

    case 'EEG-INTAN'
        F = in_fread_intan(sFile, SamplesBounds, iChannels, precision);
    case 'EEG-PLEXON'
        F = in_fread_plexon(sFile, SamplesBounds, iChannels, precision);
    case 'EEG-TDT'
        F = in_fread_tdt(sFile, SamplesBounds, iChannels);
    case {'NWB', 'NWB-CONTINUOUS'}
        isContinuous = strcmpi(sFile.format, 'NWB-CONTINUOUS');
        F = in_fread_nwb(sFile, iEpoch, SamplesBounds, iChannels, isContinuous);
    case 'MNE-PYTHON'
        [F, TimeVector] = in_fread_mne(sFile, ChannelMat, iEpoch, SamplesBounds, iChannels);
    otherwise
        error('Cannot read data from this file');
end

% Force the recordings to be in double precision
if ~isempty(precision) && strcmp(precision, 'single')
    F = single(F);
else
    F = double(F);
end
% Remove channels that were not supposed to be read
if isChanRange && ~isempty(iChanRemove)
    F(iChanRemove,:) = [];
end

%% ===== CLOSE FILE =====
if ~isempty(sfid) && ~isempty(fopen(sfid))
    fclose(sfid);
end

%% ===== TIME =====
% If TimeVector was not defined by the reading functions
if isempty(TimeVector)
    if ~isempty(SamplesBounds)
        TimeVector = (SamplesBounds(1) : SamplesBounds(2)) ./ sFile.prop.sfreq;
    elseif ~isempty(iEpoch) && ~isempty(ImportOptions) && strcmpi(ImportOptions.ImportMode, 'Epoch') && ~isempty(sFile.epochs)
        epochSamples = round(sFile.epochs(iEpoch).times * sFile.prop.sfreq);
        TimeVector = (epochSamples(1) : epochSamples(2)) / sFile.prop.sfreq;
    else
        fileSamples = round(sFile.prop.times * sFile.prop.sfreq);
        TimeVector = (fileSamples(1) : fileSamples(2)) / sFile.prop.sfreq;
    end
end
% If epoching the recordings (ie. reading by events): Use imported time window
if ~isempty(ImportOptions) && strcmpi(ImportOptions.ImportMode, 'Event')
    % TimeVector = TimeVector - TimeVector(1) + ImportOptions.EventsTimeRange(1);
    evtOffset = round(ImportOptions.EventsTimeRange(1) * sFile.prop.sfreq) / sFile.prop.sfreq;
    TimeVector = TimeVector - TimeVector(1) + evtOffset;
end

%% ===== GRADIENT CORRECTION =====
% 3rd-order gradient correction
if ~isempty(ImportOptions) && ImportOptions.UseCtfComp && ~strcmpi(sFile.format, 'BST-DATA') && ~isempty(ChannelMat) && ~isempty(ChannelMat.MegRefCoef) && ~isempty(sFile.prop.currCtfComp) && ~isequal(sFile.prop.currCtfComp, sFile.prop.destCtfComp)
    AllSensorTypes = unique({ChannelMat.Channel(iChannels).Type});
    if isempty(AllSensorTypes) || any(ismember(AllSensorTypes, {'MEG','MEG REF','MEG GRAD','MEG MAG'}))
        iMeg = good_channel(ChannelMat.Channel,[],'MEG');
        iRef = good_channel(ChannelMat.Channel,[],'MEG REF');
        if ~isempty(iChannels) && (length(iChannels) ~= length(ChannelMat.Channel))
            error('CTF compensators require that you read all the channels at the same time.');
        else
            F(iMeg,:) = F(iMeg,:) - ChannelMat.MegRefCoef * F(iRef,:);
        end
    end
end

%% ===== SSP PROJECTORS =====
if ~isempty(ImportOptions) && ImportOptions.UseSsp && ~strcmpi(sFile.format, 'BST-DATA') && ~isempty(ChannelMat) && ~isempty(ChannelMat.Projector)
    % Build projector matrix
    Projector = process_ssp2('BuildProjector', ChannelMat.Projector, 1);
    % Apply projector
    if ~isempty(Projector)
        % Get bad channels
        iBadChan = find(sFile.channelflag == -1);
        % Remove bad channels from the projector (similar as in process_megreg)
        if ~isempty(iBadChan)
            Projector(iBadChan,:) = 0;
            Projector(:,iBadChan) = 0;
            Projector(iBadChan,iBadChan) = eye(length(iBadChan));
        end
        % Apply projector
        if ~isempty(iChannels)
            % Channels that are modified by projector.
            isProjected = sum(Projector ~= 0, 2) > 1;
            if any(isProjected(iChannels))
                % If there are projectors involved and only subselection of channels:
                % We must have all data needed to apply the projector, otherwise it doesn't make sense
                missingChannels = setdiff(find(any(Projector(iChannels,:), 1)), iChannels);
                if ~isempty(missingChannels)
                    bst_report('Warning', 'process_import_data_raw', [], ['Missing channels in order to apply existing SSP/ICA projectors. To read the corrected values for channel "' ChannelMat.Channel(iChannels(1)).Name '", first apply the existing projectors with the process Artifacts > Apply SSP and CTF compensation']);
                else
                    F = Projector(iChannels, iChannels) * F;
                end
            end
        else
            F = Projector * F;
        end
    end
end

%% ===== REMOVE BASELINE ======
if ~isempty(ImportOptions) && ~isempty(ImportOptions.RemoveBaseline)
    % Get times to compute the baseline
    switch (ImportOptions.RemoveBaseline)
        case 'all'
            iTimesBl = 1:length(TimeVector);
        case 'time'
            iTimesBl = find((TimeVector >= ImportOptions.BaselineRange(1)) & (TimeVector <= ImportOptions.BaselineRange(2)));
        case 'no'
            iTimesBl = [];
    end
    % Remove baseline
    if ~isempty(iTimesBl)
        % Exclude system channels from the baseline correction
        if ~isempty(ChannelMat) && ~isempty(ChannelMat.Channel)
            iChanBl = find(~ismember(lower({ChannelMat.Channel.Type}), {'stim','video','sysclock'}));
        else
            iChanBl = 1:size(F,1);
        end
        % Sensor selection for the baseline correction
        if ~isempty(ImportOptions) && isfield(ImportOptions, 'BaselineSensorType') && ~isempty(ImportOptions.BaselineSensorType)
            % Select channels based on sensor types (or names) for baseline correction
            iChanBlSensorType = channel_find(ChannelMat.Channel, ImportOptions.BaselineSensorType);
            % Selected channels minus excluded
            iChanBl = intersect(iChanBlSensorType, iChanBl);
        end
        % Compute baseline
        blValue = mean(F(iChanBl,iTimesBl), 2);
        % Remove from recordings
        % F(iChanBl,:) = F(iChanBl,:) - repmat(blValue, [1,size(F,2)]);
        F(iChanBl,:) = bst_bsxfun(@minus, F(iChanBl,:), blValue);
    end
end

%% ===== RESAMPLE =====
if ~isempty(ImportOptions) && ImportOptions.Resample && (size(F,2) > 1) && (abs(ImportOptions.ResampleFreq - sFile.prop.sfreq) > 0.05)
    [F, TimeVector] = process_resample('Compute', F, TimeVector, ImportOptions.ResampleFreq);
end

author: brainstorm-tools


text:
function [mu, varsigma] = fgplvmPosteriorMeanVar(model, X);

% FGPLVMPOSTERIORMEANVAR Mean and variances of the posterior at points given by X.
% FORMAT
% DESC returns the posterior mean and variance for a given set of
% points.
% ARG model : the model for which the posterior will be computed.
% ARG x : the input positions for which the posterior will be
% computed.
% RETURN mu : the mean of the posterior distribution.
% RETURN sigma : the variances of the posterior distributions.
%
% SEEALSO : gpPosteriorMeanVar, fgplvmCreate
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% FGPLVM

[mu, varsigma] = gpPosteriorMeanVar(model, X);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [L_l,R_l] = lindemannInh(L,R,fs,c_inh,dim)
%LINDEMANNINH Signal pre-processing for Lindemann's cross-correlation
%
%   [L_L,R_L] = IOSR.AUDITORY.LINDEMANNINH(L,R,FS) pre-processes left L and
%   right R signals for the cross-correlation function based on Lindemann's
%   precedence model [1,2]. A crucial parameter for Lindemann's model is
%   the "operating point", which controls the amount of inhibition. The
%   parameter is actually a function of the input related to its RMS level.
%   After half-wave rectifying and filtering the input signals L and R
%   (sampled at FS Hz) along the first non-singleton dimension, this
%   function applies a gain related to the inhibition parameter C_INH
%   (default is 0.3). The gain is identical for all rows, columns, etc.
%   Lastly, values outside of the interval [0,1] are not permitted in
%   Lindemann's model and hence these values are clipped.
%
%   [L_L,R_L] = IOSR.AUDITORY.LINDEMANNINH(L,R,FS,C_INH) uses the specified
%   inhibition parameter C_INH. The value must be in the interval [0,1].
%
%   [L_L,R_L] = IOSR.AUDITORY.LINDEMANNINH(L,R,FS,C_INH,DIM) pre-processes
%   L and R along the dimension DIM.
%
%   References
%
%   [1] Lindemann, W. (1986), Extension of a binaural cross-correlation
%       model by contralateral inhibition. I. Simulation of lateralization
%       for stationary signals, The Journal of the Acoustical Society of
%       America 80, 6, 1608-1622.
%
%   [2] Lindemann, W. (1986), Extension of a binaural cross-correlation
%       model by contralateral inhibition. II. The law of the first wave
%       front, The Journal of the Acoustical Society of America 80, 6,
%       1623-1630.
%
%   Further reading
%
%   Hummersone, C., Mason, R., Brookes, T. (2013), A comparison of
%       computational precedence models for source separation in
%       reverberant environments, The Journal of the Audio Engineering
%       Society 61, 7/8, 508-520.
%
%   See also IOSR.AUDITORY.XCORRLINDEMANN.

%   Copyright 2016 University of Surrey.

    %% check input

    assert(isequal(size(L),size(R)), 'iosr:lindemannInh:invalidInputs', 'L and R arrays must be the same size')
    assert(isscalar(fs), 'iosr:lindemannInh:invalidFs', 'FS must be a scalar')

    % default dim
    if nargin<5
        dim = find(size(L)>1,1,'first');
    end

    %% process

    % half-wave rectify
    L_l = hwr(L);
    R_l = hwr(R);

    % filter
    cutofffreq=800;
    [b,a] = butter(1,cutofffreq*2/fs);
    L_l = filter(b,a,L_l,[],dim);
    R_l = filter(b,a,R_l,[],dim);

    % inhibition parameter
    if nargin < 4
        c_inh = .3;
    else
        assert(isscalar(c_inh) & isnumeric(c_inh), 'iosr:lindemannInh:invalidCinh', 'c_inh must be a scalar');
        assert(c_inh>=0 || c_inh<=1, 'iosr:lindemannInh:invalidX', 'c_inh must be in the interval (0,1].')
    end

    % gain
    c_gamma_L = calc_gamma(L,dim);
    c_gamma_R = calc_gamma(R,dim);
    c_gamma = max([c_gamma_L(:); c_gamma_R(:)]);

    % apply parameters
    L_l = apply_gamma(L_l,c_inh,c_gamma);
    R_l = apply_gamma(R_l,c_inh,c_gamma);

    % restrict range
    L_l = clip(L_l);
    R_l = clip(R_l);

end

function y = hwr(x)
%HWR half-wave rectify

    y = max(x,0);

end

function gamma = calc_gamma(x,dim)
%CALC_GAMMA calculate the gamma

    gamma = sqrt(2).*iosr.dsp.rms(x,dim);

end

function y = apply_gamma(x,c_inh,c_gamma)
%APPLY_GAMMA apply gamma to achieve inhibition parameter

    y = (c_inh/c_gamma).*x;

end

function y = clip(x)
%CLIP clip input data to [0,1] interval

    y = x;
    y(y<0) = 0;
    y(y>1) = 1;

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function barycentric_interp_1d_test ( )

%*****************************************************************************80
%
%% BARYCENTRIC_INTERP_1D_TEST tests the BARYCENTRIC_INTERP_1D library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../r8lib' )
  addpath ( '../test_interp_1d' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BARYCENTRIC_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the BARYCENTRIC_INTERP_1D library.\n' );
  fprintf ( 1, '  The R8LIB library is needed.\n' );
  fprintf ( 1, '  These test needs the TEST_INTERP_1D library.\n' );

  prob_num = p00_prob_num ( );
  for prob = 1 : prob_num
    for nd = [ 4, 8, 16, 32, 64, 1000 ]
      barycentric_interp_1d_test01 ( prob, nd );
    end
  end

  prob_num = p00_prob_num ( );
  for prob = 1 : prob_num
    for nd = [ 4, 8, 16, 32, 64, 1000 ]
      barycentric_interp_1d_test02 ( prob, nd );
    end
  end

  prob_num = p00_prob_num ( );
  for prob = 1 : prob_num
    for nd = [ 4, 8, 16, 32, 64, 1000 ]
      barycentric_interp_1d_test03 ( prob, nd );
    end
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BARYCENTRIC_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../r8lib' )
  rmpath ( '../test_interp_1d' )

  return
end


text:
function v = krawtchouk ( n, p, x, m )

%*****************************************************************************80
%
%% KRAWTCHOUK evaluates the Krawtchouk polynomials at X.
%
%  Discussion:
%
%    The polynomial has a parameter P, which must be strictly between
%    0 and 1, and a parameter M which must be a nonnegative integer.
%
%    The Krawtchouk polynomial of order N, with parameters P and M,
%    evaluated at X, may be written K(N,P,X,M).
%
%    The first two terms are:
%
%      K(0,P,X,M) = 1
%      K(1,P,X,M) = X - P * M
%
%    and the recursion, for fixed P and M is
%
%                             ( N + 1 ) * K(N+1,P,X,M) =
%        ( X - ( N + P * ( M - 2 * N))) * K(N,  P,X,M)
%       - ( M - N + 1 ) * P * ( 1 - P ) * K(N-1,P,X,M)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Walter Gautschi,
%    Orthogonal Polynomials: Computation and Approximation,
%    Oxford, 2004,
%    ISBN: 0-19-850672-4,
%    LC: QA404.5 G3555.
%
%  Parameters:
%
%    Input, integer N, the highest order polynomial to evaluate.
%    0 <= N.
%
%    Input, real P, the parameter.  0 < P < 1.
%
%    Input, real X, the evaluation parameter.
%
%    Input, integer M, the parameter.  0 <= M.
%
%    Output, real V(1:N+1), the values of the Krawtchouk polynomials
%    of orders 0 through N at X.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KRAWTCHOUK - Fatal error!\n' );
    fprintf ( 1, '  0 <= N is required.\n' );
    error ( 'KRAWTCHOUK - Fatal error!' );
    return
  end

  if ( p <= 0.0 || 1.0 <= p )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KRAWTCHOUK - Fatal error!\n' );
    fprintf ( 1, '  0 < P < 1 is required.\n' );
    error ( 'KRAWTCHOUK - Fatal error!' );
    return
  end

  if ( m < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KRAWTCHOUK - Fatal error!\n' );
    fprintf ( 1, '  0 <= M is required.\n' );
    error ( 'KRAWTCHOUK - Fatal error!' );
    return
  end

  OFFSET = 1;

  v = zeros(n+1,1);
  v(0+OFFSET) = 1.0;

  if ( 1 <= n )
    v(1+OFFSET) = x - p * m;
  end

  for i = 1 : n - 1
    v(i+1+OFFSET) = ( ( x - ( i + p * ( m - 2 * i ) ) ) * v(i+OFFSET) ...
                        - ( m - i + 1 ) * p * ( 1 - p ) * v(i-1+OFFSET)  ) ...
                    / ( i + 1 );
  end

  return
end


text:
% Test how accurately we can orthogonalize ill-conditioned bases of tangent
% vectors. We test three metrics:
%  1) Is Q actually orthonormal (the figure displays log(|Q'Q|); since Q'Q
%     should be identity, we hope to see 0 on diagonal and -inf everywhere
%     else.)
%  2) Is A = QR? This is in the title of the plots.
%  3) How much time does it take to compute?
%
% Nicolas Boumal, Oct. 5, 2017

clear all;
clc;

M = spherefactory(100);

% M = stiefelcomplexfactory(200, 56);
% M = productmanifold(struct('S', spherefactory(5), 'R', rotationsfactory(3, 10)));

x = M.rand();

% Create a poorly conditioned basis
A = cell(M.dim()-5, 1);
A{1} = M.randvec(x);
for k = 2 : numel(A)
    A{k} = M.lincomb(x, 1, A{k-1}, 1e-6, M.randvec(x));
end

t1 = tic();
[Q1, R1] = orthogonalize(M, x, A);
t1 = toc(t1);

t2 = tic();
[Q2, R2] = orthogonalizetwice(M, x, A);
t2 = toc(t2);

t3 = tic();
[Q3, R3] = orthogonalize_legacy(M, x, A);
t3 = toc(t3);

% To check if Q is really orthonormal
G1 = grammatrix(M, x, Q1);
G2 = grammatrix(M, x, Q2);
G3 = grammatrix(M, x, Q3);

% To check of QR = A
distsq = 0;
for k = 1 : numel(A)
    distsq = distsq + M.norm(x, M.lincomb(x, 1, lincomb(M, x, Q1(1:k), R1(1:k, k)), -1, A{k}))^2;
end
dist1 = sqrt(distsq);

distsq = 0;
for k = 1 : numel(A)
    distsq = distsq + M.norm(x, M.lincomb(x, 1, lincomb(M, x, Q2(1:k), R2(1:k, k)), -1, A{k}))^2;
end
dist2 = sqrt(distsq);

distsq = 0;
for k = 1 : numel(A)
    distsq = distsq + M.norm(x, M.lincomb(x, 1, lincomb(M, x, Q3(1:k), R3(1:k, k)), -1, A{k}))^2;
end
dist3 = sqrt(distsq);

subplot(1, 3, 1);
imagesc(log10(abs(G1))); axis equal, axis tight; colorbar; set(gca, 'CLim', [-20, 0]);
title(sprintf('MGS: CPU time: %g\n||A - QR|| = %g', t1, dist1));
subplot(1, 3, 2);
imagesc(log10(abs(G2))); axis equal, axis tight; colorbar; set(gca, 'CLim', [-20, 0]);
title(sprintf('MGS-twice: CPU time: %g\n||A - QR|| = %g', t2, dist2));
subplot(1, 3, 3);
imagesc(log10(abs(G3))); axis equal, axis tight; colorbar; set(gca, 'CLim', [-20, 0]);
title(sprintf('Legacy: CPU time: %g\n||A - QR|| = %g', t3, dist3));

author: NicolasBoumal

repo: manopt


text:
function permindx = permute_setupperms(n,nperms)
% Set up permutations x observations matrix of observation indices for
% permutation test
%
% :Usage:
% ::
%
%    % approximate test with nperms obs.
%    permindx = permute_setupperms(n,nperms)
%    permindx = permute_setupperms(n,[])
%
%    % exact test (all permutations) of n observations
%    permindx = permute_setupperms(n)
%
% ..
%    Uses code from SnPM3b by Tom Nichols
%    adapted for this function by Tor Wager
% ..

  if nargin < 2 || isempty(nperms) || nperms == 0
    OUT.meth = 'exact';
    nperms = gamma(n+1);
    use_approximate = 0;
  else
    nperms = min(gamma(n+1),nperms);
    use_approximate = 1;

    if nperms > 50000
      warning('More than 50,000 permutations! May cause time/memory problems');
    end
  end

				% interactive
				% % nperms = gamma(n+1);
% % use_approximate = input(sprintf('%d Perms. Enter return for exact test, or num. of perms in random subset: ',nperms));
% % if isempty(use_approximate)
% %     use_approximate = 0;
% % else
% %     nperms = min(use_approximate,nperms);
% %     use_approximate = 1;
% % end

  if use_approximate
      % taken from snpm3b code, by Tom Nichols
      %-Approximate test :
      % Build up random subset of all (within n) permutations
      %===============================================================
    rand('seed',sum(100*clock))	%-Initialise random number generator
    permindx      = zeros(nperms,n);
    permindx(1,:) = 1+rem(0:(n-1), n);

    for i = 2:nperms+1
      permindx(i,:) = randperm(n);
    end
    b = unique(permindx, 'rows');
    if(size(b,1) ~= size(permindx, 1))
      while(size(b,1) ~= size(permindx, 1))
        permindx = b;
        for i = (size(permindx,1)+1):nperms+1
          permindx(i,:) = randperm(n);
        end
        b = unique(permindx, 'rows');
      end
      permindx = [b(1,:); b(randperm(size(b,1)-1)+1,:)];
    end

		       % null permutation is not included in perm indx
    permindx = permindx(2:end,:);
  else
      % taken from snpm3b code, by Tom Nichols
      %
      %-Full permutation test :
      % Build up exhaustive matrix of permutations
      %===============================================================
      %-Compute permutations for a single exchangability block
      %---------------------------------------------------------------
      %-Initialise Xblkpermindx & remaining numbers
    Xblkpermindx = [];
    lef = [1:n]';
%-Loop through numbers left to add to permutations, accumulating permindx
    for i = n:-1:1
				%-Expand Xblkpermindx & lef
      tmp = round(exp(gammaln(n+1)-gammaln(i+1)));
      Exp = meshgrid(1:tmp,1:i); Exp = Exp(:)';
      if ~isempty(Xblkpermindx), Xblkpermindx = Xblkpermindx(:,Exp); end
      lef = lef(:,Exp);
				%-Work out sampling for lef
      tmp1 = round(exp(gammaln(n+1)-gammaln(i+1)));
      tmp2 = round(exp(gammaln(n+1)-gammaln(i)));
      sam = 1+rem(0:i*tmp1-1,i) + ([1:tmp2]-1)*i;
			      %-Add samplings from lef to Xblkpermindx
      Xblkpermindx   = [Xblkpermindx; lef(sam)];
		      %-Delete sampled items from lef & condition size
      lef(sam) = [];
      tmp = round(exp(gammaln(n+1)-gammaln((i-1)+1)));
      lef = reshape(lef,(i-1),tmp);
		%NB:gamma(n+1)/gamma((i-1)+1) == size(Xblkpermindx,2);
    end
    clear lef Exp sam i
				%-Reorient so permutations are in rows
    permindx = Xblkpermindx';
  end

			     %-Check, condition and randomise permindx
%-----------------------------------------------------------------------
%-Check permindxs sum within Xblks to sum to 1
  if ~all(all(sum(permindx,2) == (n+1)*n/2 ))
    error('Invalid permindx computed!'), end

	%-Convert to full permutations from permutations within blocks
  nperms = size(permindx,1);

%-Randomise order of permindxs (except first) to allow interim analysis
  rand('seed',sum(100*clock))	%-Initialise random number generator
  permindx=[permindx(1,:);permindx(randperm(nperms-1)+1,:)];

			 %-Check first permutation is null permutation
  if ~use_approximate
    if ~all(permindx(1,:)==[1:n])
      error('permindx(1,:)~=[1:n]');
    end
  end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function test_bug2502

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_freqstatistics ft_selectdata

% see bug 2502 for more info

% make some data
data1.someparameter = randn(100,10,15);
data1.dimord = 'subj_chan_freq';
data1.freq   = 1:15;
for k = 1:10
  data1.label{k} = ['chan',num2str(k)];
end
data2 = data1;
data2.someparameter = randn(100,10,15);

cfg = [];
cfg.method = 'montecarlo';
cfg.parameter = 'someparameter';
cfg.statistic = 'ft_statfun_indepsamplesT';
cfg.numrandomization = 0;
cfg.design = [ones(1,100) ones(1,100)*2];
stat = ft_freqstatistics(cfg, data1, data2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function obj = updateErrorbar(obj, plotIndex)

	%-------------------------------------------------------------------------%

	%-INITIALIZATION-%

	%-get data structures-%
	plotData = get(obj.State.Plot(plotIndex).Handle);

	%-get error data-%
	yPositiveDelta = date2NumData(plotData.YPositiveDelta);
	yNegativeDelta = date2NumData(plotData.YNegativeDelta);

	xPositiveDelta = date2NumData(plotData.XPositiveDelta);
	xNegativeDelta = date2NumData(plotData.XNegativeDelta);

	%-------------------------------------------------------------------------%

	%-set trace-%

	updateLineseries(obj, plotIndex);

	%-errorbar visible-%
	obj.data{plotIndex}.error_y.visible = true;
	obj.data{plotIndex}.error_x.visible = true;

	%-errorbar type-%
	obj.data{plotIndex}.error_y.type = 'data';
	obj.data{plotIndex}.error_x.type = 'data';

	%-errorbar symmetry-%
	obj.data{plotIndex}.error_y.symmetric = false;

	%-------------------------------------------------------------------------%

	%-set errorbar data-%

	obj.data{plotIndex}.error_y.array = yPositiveDelta;
	obj.data{plotIndex}.error_x.array = xPositiveDelta;

	obj.data{plotIndex}.error_x.arrayminus = xNegativeDelta;
	obj.data{plotIndex}.error_y.arrayminus = yNegativeDelta;

	%-------------------------------------------------------------------------%

	%-errorbar thickness-%
	obj.data{plotIndex}.error_y.thickness = plotData.LineWidth;
	obj.data{plotIndex}.error_x.thickness = plotData.LineWidth;

	%-errorbar width-%
	obj.data{plotIndex}.error_y.width = obj.PlotlyDefaults.ErrorbarWidth;
	obj.data{plotIndex}.error_x.width = obj.PlotlyDefaults.ErrorbarWidth;

	%-errorbar color-%
	errorColor = getStringColor(255*plotData.Color);
	obj.data{plotIndex}.error_y.color = errorColor;
	obj.data{plotIndex}.error_x.color = errorColor;

	%-------------------------------------------------------------------------%

end

author: plotly


text:
%% Copyright (C) 2014-2016, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym atand (@var{x})
%% Symbolic inverse tan function with output in degrees.
%%
%% Example:
%% @example
%% @group
%% atand (sqrt (sym (3)))
%%   @result{} (sym) 60
%%
%% syms x
%% y = atand (x)
%%   @result{} y = (sym)
%%       180⋅atan(x)
%%       ───────────
%%            π
%% @end group
%% @end example
%%
%% @seealso{@@sym/tand, @@sym/atan}
%% @end defmethod

function y = atand(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('lambda a: deg(atan(a))', x);
end

%!error atand (sym(1), 2)
%!assert (isequaln (atand (sym(nan)), sym(nan)))

%!test
%! f1 = atand (sym(1)/2);
%! f2 = atand (1/2);
%! assert (double (f1), f2, -eps)

%!test
%! D = [1 2; 3 4]/4;
%! A = sym([1 2; 3 4])/4;
%! f1 = atand (A);
%! f2 = atand (D);
%! assert (double (f1), f2, -eps)

repo: octsympy


text:
function resampScanNumV = getAssocResampledScanNum(scanNumV,planC)
% Return index of resampled scan created from input scanNum.
% ------------------------------------------------------------------------
% scanNumV   : Vector of scan nos.
% planC
% ------------------------------------------------------------------------
% AI 5/19/22

indexS = planC{end};

scanUIDc = {planC{indexS.scan}(scanNumV).scanUID};

resampScanNumV = nan(1,length(scanNumV)); % Initialize
for nScan = 1:length(planC{indexS.scan})
    baseScanUID = planC{indexS.scan}(nScan).assocBaseScanUID;
    assocIdxV = ismember(scanUIDc,baseScanUID);
    if any(assocIdxV)
        resampScanNumV(scanNumV==find(assocIdxV)) = nScan;
    end
end
resampScanNumV = resampScanNumV(~isnan(resampScanNumV));

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [f, g] = obj_infometric(L, sX, sY, tX, lambda)

% Input:
% L: D*d transformation matrix
% sX: instance matrix for the source domain data. size: N*D
% sY: label vector for the source domain data. size: N*1
% tX: instance matrix for the target domain data. size: M*D
% lambda: regularization parameter
%
% Output:
% f: objective function value
% g: gradient on L
%
% yuanshi@usc.edu
% 2013/9/5
%

% mutual info in the target domain
% do "discriminative clustering" on the target domain
[mi, grad] = compute_mutual_infoL(L, sX, sY, sX);

% mutual info in both domains
% make it hard to distinguish the two domains
all_data = [sX; tX];
all_label = [ones( size(sX,1), 1); zeros( size(tX,1), 1)];  % domain label

if lambda > 0
    [mi2, grad2] = compute_mutual_infoL(L, all_data, all_label, all_data);

    f = -(mi - lambda*mi2);
    g = -(grad - lambda*grad2);
else
    f = -mi;
    g = -grad;
end

author: viggin

repo: domain-adaptation-toolbox


text:
function title = p16_title ( )

%*****************************************************************************80
%
%% P16_TITLE returns a title for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Beale Function.';

  return
end


text:
function net = our_vl_simplenn_move(net, destination)
%VL_SIMPLENN_MOVE  Move a SimpleNN network between CPU and GPU.
%   NET = VL_SIMPLENN_MOVE(NET, 'gpu') moves the network to the
%   current GPU device. NET = VL_SIMPLENN_MOVE(NET, 'cpu') moves the
%   network to the CPU.
%
%   See also: VL_SIMPLENN().

% Copyright (C) 2014-15 Andrea Vedaldi.
% 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).

switch destination
  case 'gpu', moveop = @(x) gpuArray(x) ;
  case 'cpu', moveop = @(x) gather(x) ;
  otherwise, error('Unknown destination ''%s''.', destination) ;
end
for l=1:numel(net.layers)
  switch net.layers{l}.type
    case {'conv', 'convt', 'bnorm','conv_mask'}
      for f = {'filters', 'biases', 'filtersMomentum', 'biasesMomentum'}
        f = char(f) ;
        if isfield(net.layers{l}, f)
          net.layers{l}.(f) = moveop(net.layers{l}.(f)) ;
        end
      end
      for f = {'weights', 'momentum'}
        f = char(f) ;
        if isfield(net.layers{l}, f)
          for j=1:numel(net.layers{l}.(f))
            net.layers{l}.(f){j} = moveop(net.layers{l}.(f){j}) ;
          end
        end
      end
      if(strcmp(net.layers{l}.type,'conv_mask'))
          net.layers{l}.posTemp{1}=moveop(net.layers{l}.posTemp{1});
          net.layers{l}.posTemp{2}=moveop(net.layers{l}.posTemp{2});
      end
    otherwise
      % nothing to do ?
  end
end

repo: interpretableCNN


text:
function TorqueOverTime(nlp, sys)

    domains = sys.Gamma.Nodes.Domain;
    T  = SymVariable('t',[2, 1]);
    for i=1:numel(domains)
        domain = domains{i};

        u = domain.Inputs.Control.u;
        u2r = tovector(norm(u).^2)./(T(2)-T(1));
        u2r_fun = SymFunction(['torqueOverTime_' sys.Gamma.Nodes.Name{i}],u2r,{u,T});
        rs_phase = getPhaseIndex(nlp,sys.Gamma.Nodes.Name{i});
        addRunningCost(nlp.Phase(rs_phase),u2r_fun,{'u', 'T'});

    end
    nlp.update;

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function lA = longAxis(grains,varargin)
% long axis of a grain
%
% the long axis is the direction of the largest
% <grain2d.principalComponents.html,principal component> of a grain
%
% Syntax
%   lA = grains.longAxis
%
% Input
%  grains - @grain2d
%
% Output
%  lA - @vector3d direction of the longest elongation
%
% See also
% grain2d/principalComponents

omega = principalComponents(grains);

lA = vector3d.byPolar(pi/2,omega,'antipodal');

author: mtex-toolbox

repo: mtex

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


text:
function v = makemodel(p,ng)
%MAKEMODEL Helper function to make model matrix
%    P = max model term size, NG = number of grouping variables
% or
%    P = vector of term codes

% We want a matrix with one row per term.  Each row has a 1 for
% the variables participating in the term.  There will be rows
% summing up to 1, 2, ..., p.

if numel(p)==1
   % Create model matrix from a scalar max order value
   vgen = 1:ng;
   v = eye(ng);                      % linear terms
   for j=2:min(p,ng)
      c = nchoosek(vgen,j);          % generate column #'s with 1's
      nrows = size(c,1);             % generate row #'s
      r = repmat((1:nrows)',1,j);    %    and make it conform with c
      m = zeros(nrows,ng);           % create a matrix to get new rows
      m(r(:)+nrows*(c(:)-1)) = 1;    % fill in 1's
      v = [v; m];                    % append rows
   end
else
   % Create model matrix from terms encoded as bit patterms
   nterms = length(p);
   v = zeros(nterms,ng);
   for j=1:nterms
      tm = p(j);
      while(tm)
         % Get last-numbered effect remaining
         lne = 1 + floor(log2(tm));
         tm = bitset(tm, lne, 0);
         v(j,lne) = 1;
      end
   end
end