test meeg chan neighborhood

function test_suite = test_meeg_chan_neighborhood()
    % tests for cosmo_meeg_chan_neighborhood
    %
    % #   For CoSMoMVPA's copyright information and license terms,   #
    % #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions = localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_neighbors()
    % note: this tests assumes that meeg_chan_neighbors works properly
    if cosmo_skip_test_if_no_external('fieldtrip')
        return
    end

    % only test every test_step-th channel
    test_step = 10;

    ds = cosmo_synthetic_dataset('type', 'meeg', 'size', 'big');
    msk = ds.fa.chan < 10 | mod(ds.fa.chan, 5) == 0 | ds.fa.chan > 250;
    imsk = find(msk);
    rp = randperm(numel(imsk));
    ds = cosmo_slice(ds, imsk(rp), 2);
    ds = cosmo_dim_prune(ds);
    n = numel(ds.a.fdim.values{1});
    ds.a.fdim.values{1} = ds.a.fdim.values{1}(randperm(n));

    nbrs = cosmo_meeg_chan_neighbors(ds, 'chantype', 'meg_planar', 'radius', 0);
    nh = cosmo_meeg_chan_neighborhood(ds, nbrs);

    assertEqual(nh.a.fdim.values, {{nbrs.label}});
    assertEqual(nh.a.fdim.labels, {'chan'});

    n = numel(nh.neighbors);
    for k = 1:test_step:n
        idx = nh.neighbors{k};
        chan_label = ds.a.fdim.values{1}(ds.fa.chan(idx));
        assert(all(cosmo_match(chan_label, nh.a.fdim.values{1}{k})));
    end

    % test correspondence with neighbors
    args = {'chantype', 'all', 'delaunay', true};
    nbrs = cosmo_meeg_chan_neighbors(ds, args{:});
    nh = cosmo_meeg_chan_neighborhood(ds, args{:});
    nh2 = cosmo_meeg_chan_neighborhood(ds, nbrs);
    assertEqual(nh, nh2);

    ds_label = ds.a.fdim.values{1};

    n = numel(nh.neighbors);
    for k = 1:test_step:n
        idx = nh.neighbors{k};

        chan_label = ds.a.fdim.values{1}(ds.fa.chan(idx));

        % ensure both arguments for intersect are column vectors,
        % because Octave behaves differently than Matlab
        assertEqual(intersect(ds_label(:), nbrs(k).neighblabel(:)), ...
                    unique(chan_label'));

    end

    % try with dataset labels
    args = {'chantype', 'all', 'count', 5, 'label', 'dataset'};
    nbrs = cosmo_meeg_chan_neighbors(ds, args{:});
    nh = cosmo_meeg_chan_neighborhood(ds, nbrs);
    n = numel(nh.neighbors);
    assertEqual(n, numel(ds.a.fdim.values{1}));
    assertEqual(1:n, nh.fa.chan);
    ds_label = ds.a.fdim.values{1};
    nh_label = nh.a.fdim.values{1};

    for k = 1:test_step:n
        i = find(cosmo_match(nh_label, ds_label{k}));
        assertEqual(nbrs(i).label, ds_label{k});

        overlap = cosmo_overlap({ds_label(ds.fa.chan(nh.neighbors{i}))}, ...
                                {nbrs(i).neighblabel});
        assertEqual(overlap, 1);
    end

    % test for number of channels
    nbrs = cosmo_meeg_chan_neighbors(ds, 'count', 5, ...
                                     'chantype', 'meg_combined_from_planar');
    nh = cosmo_meeg_chan_neighborhood(ds, nbrs);
    assertEqual(numel(nh.neighbors), 102);
    h = cellfun(@numel, nh.neighbors) / 7;
    assert(all(h >= 5 & h <= 10));
    mh = mean(h);
    assert(mh > 6 && mh < 9);

    nbrs = cosmo_meeg_chan_neighbors(ds, 'count', 5, 'chantype', 'meg_planar');
    nh = cosmo_meeg_chan_neighborhood(ds, nbrs);
    assertEqual(numel(nh.neighbors), 204);
    h = cellfun(@numel, nh.neighbors) / 7;
    assert(all(h == 5));

    % test tiny dataset
    ds = cosmo_synthetic_dataset('type', 'meeg');
    opt = struct();
    opt.delaunay = true;
    opt.label = 'dataset';
    opt.chantype = 'meg_axial';

    nh = cosmo_meeg_chan_neighborhood(ds, opt);
    assertEqual(nh.neighbors, {[1 4]});

    opt.chantype = 'meg_planar';
    nh = cosmo_meeg_chan_neighborhood(ds, opt);
    assertEqual(nh.neighbors, {[2 5 3 6]; [2 5 3 6]});