diff --git a/CanlabCore/Unit_tests/canlab_run_all_tests.m b/CanlabCore/Unit_tests/canlab_run_all_tests.m index 2257e139..5e013a72 100644 --- a/CanlabCore/Unit_tests/canlab_run_all_tests.m +++ b/CanlabCore/Unit_tests/canlab_run_all_tests.m @@ -63,7 +63,9 @@ case 'include' % run everything end - suite = [suite, TestSuite.fromFile(fpath)]; %#ok + % canlab_safe_suite_from_file warn-skips files that are not valid test + % files instead of letting a NonTestFile error abort the whole run. + suite = [suite, canlab_safe_suite_from_file(fpath)]; %#ok end tag = char(p.Results.Tag); diff --git a/CanlabCore/Unit_tests/canlab_test_runner_robustness.m b/CanlabCore/Unit_tests/canlab_test_runner_robustness.m new file mode 100644 index 00000000..149254dd --- /dev/null +++ b/CanlabCore/Unit_tests/canlab_test_runner_robustness.m @@ -0,0 +1,76 @@ +function tests = canlab_test_runner_robustness +%CANLAB_TEST_RUNNER_ROBUSTNESS Discovery must not abort on a stray non-test file. +% +% canlab_run_all_tests globs every canlab_test_*.m under Unit_tests and +% concatenates the suites. A file that matches the name pattern but is not a +% valid matlab.unittest test (an old plain-assert script, or a function whose +% internal name does not match the filename) makes TestSuite.fromFile throw +% MATLAB:unittest:TestSuite:NonTestFile. Without guarding, that single error +% aborts discovery of the entire suite. canlab_safe_suite_from_file warn-skips +% such files instead; these tests pin that behavior. + +tests = functiontests(localfunctions); +end + + +function test_nontest_file_is_warn_skipped(tc) %#ok<*DEFNU> +% A canlab_test_*.m whose function name != filename and that never calls +% functiontests is not a valid test file: skip it, warn, return empty. +folder = tc.applyFixture( ... + matlab.unittest.fixtures.TemporaryFolderFixture).Folder; +fpath = local_write(folder, 'canlab_test_bogus_nontest.m', { ... + 'function some_other_name()' + 'assert(true);' + 'end'}); + +% Use lastwarn rather than verifyWarning so the check does not depend on the +% ambient warning display state (a caller may have warnings off; lastwarn +% still records the id even when the warning event is not displayed). +lastwarn('', ''); +[suite, ok] = canlab_safe_suite_from_file(fpath); +[~, warn_id] = lastwarn; + +tc.verifyFalse(ok, 'expected ok=false for a non-test file'); +tc.verifyEmpty(suite, 'expected an empty suite for a non-test file'); +tc.verifyEqual(warn_id, 'canlab_run_all_tests:skippedNonTestFile', ... + 'expected a skippedNonTestFile warning'); +end + + +function test_valid_test_file_is_loaded(tc) +% A well-formed functiontests file loads normally (ok=true, non-empty suite). +folder = tc.applyFixture( ... + matlab.unittest.fixtures.TemporaryFolderFixture).Folder; +fpath = local_write(folder, 'canlab_test_valid_dummy.m', { ... + 'function tests = canlab_test_valid_dummy' + 'tests = functiontests(localfunctions);' + 'end' + '' + 'function test_trivial(tc)' + 'tc.verifyTrue(true);' + 'end'}); + +[suite, ok] = canlab_safe_suite_from_file(fpath); +tc.verifyTrue(ok, 'expected ok=true for a valid test file'); +tc.verifyNotEmpty(suite); +tc.verifyEqual(numel(suite), 1); +end + + +function fpath = local_write(folder, name, lines) +% Write a cellstr of lines to folder/name and return the full path. Uses +% fopen/fprintf (not writelines) so the test runs on older MATLAB too. +fpath = fullfile(folder, name); +fid = fopen(fpath, 'w'); +tc_assert_open(fid, fpath); +cleanup = onCleanup(@() fclose(fid)); +fprintf(fid, '%s\n', lines{:}); +end + + +function tc_assert_open(fid, fpath) +if fid < 0 + error('canlab_test_runner_robustness:cannotWrite', ... + 'Could not open %s for writing.', fpath); +end +end diff --git a/CanlabCore/Unit_tests/helpers/canlab_get_dpsp_hot_warm.m b/CanlabCore/Unit_tests/helpers/canlab_get_dpsp_hot_warm.m new file mode 100644 index 00000000..e94fdd88 --- /dev/null +++ b/CanlabCore/Unit_tests/helpers/canlab_get_dpsp_hot_warm.m @@ -0,0 +1,27 @@ +function [hot, warm, ok] = canlab_get_dpsp_hot_warm() +%CANLAB_GET_DPSP_HOT_WARM Load the DPSP single-subject Hot/Warm sample maps. +% +% [hot, warm, ok] = canlab_get_dpsp_hot_warm() +% +% Returns the single-subject Hot and Warm condition maps from +% Sample_datasets/DPSP_pain_rejection_participant_maps as fmri_data objects. +% Used by the predictive_model / xval_* tests so the load semantics live in +% one place. ok is false (and hot/warm are []) if the sample files are not +% on the path, so callers can assume/skip gracefully. + +sample_dir = fullfile(fileparts(fileparts(which('fmri_data'))), ... + 'Sample_datasets', 'DPSP_pain_rejection_participant_maps'); +hot_file = fullfile(sample_dir, 'DPSP_single_subject_images_hot.mat'); +warm_file = fullfile(sample_dir, 'DPSP_single_subject_images_warm.mat'); + +ok = exist(hot_file, 'file') == 2 && exist(warm_file, 'file') == 2; +hot = []; +warm = []; + +if ok + H = load(hot_file); + W = load(warm_file); + hot = H.single_subject_images_hot; + warm = W.single_subject_images_warm; +end +end diff --git a/CanlabCore/Unit_tests/helpers/canlab_safe_suite_from_file.m b/CanlabCore/Unit_tests/helpers/canlab_safe_suite_from_file.m new file mode 100644 index 00000000..b50b122f --- /dev/null +++ b/CanlabCore/Unit_tests/helpers/canlab_safe_suite_from_file.m @@ -0,0 +1,42 @@ +function [suite, ok] = canlab_safe_suite_from_file(fpath) +%CANLAB_SAFE_SUITE_FROM_FILE Build a test suite from one file without aborting on bad files. +% +% [suite, ok] = canlab_safe_suite_from_file(fpath) +% +% Wraps matlab.unittest.TestSuite.fromFile so that a file which is not a +% valid test (e.g. a stray canlab_test_*.m that is an old plain-assert +% script or whose internal function name does not match the filename) +% does not throw and abort discovery of the whole suite. Such a file is +% warn-skipped and an empty Test array is returned instead. +% +% This matters because canlab_run_all_tests globs every canlab_test_*.m +% under Unit_tests and concatenates the results; a single NonTestFile +% error from fromFile would otherwise take down the entire run. +% +% :Inputs: +% **fpath:** absolute path to a candidate .m test file. +% +% :Outputs: +% **suite:** a matlab.unittest.Test array (empty if the file is not a +% valid test file). +% **ok:** logical, true if fromFile succeeded, false if the file was +% skipped. +% +% :See also: canlab_run_all_tests, matlab.unittest.TestSuite + +ok = true; +suite = matlab.unittest.Test.empty; + +try + suite = matlab.unittest.TestSuite.fromFile(fpath); +catch ME + ok = false; + if strcmp(ME.identifier, 'MATLAB:unittest:TestSuite:NonTestFile') + warning('canlab_run_all_tests:skippedNonTestFile', ... + 'Skipping %s: not a valid matlab.unittest test file.', fpath); + else + warning('canlab_run_all_tests:fromFileError', ... + 'Skipping %s: %s (%s)', fpath, ME.message, ME.identifier); + end +end +end diff --git a/CanlabCore/Unit_tests/image_vector/canlab_test_check_roi_extraction.m b/CanlabCore/Unit_tests/image_vector/canlab_test_check_roi_extraction.m new file mode 100644 index 00000000..3d609903 --- /dev/null +++ b/CanlabCore/Unit_tests/image_vector/canlab_test_check_roi_extraction.m @@ -0,0 +1,58 @@ +function tests = canlab_test_check_roi_extraction +%CANLAB_TEST_CHECK_ROI_EXTRACTION Multi-region ROI extraction + write/reload round-trip. +% +% Complements canlab_test_extract_roi (single-mask extraction) by exercising +% the two things that test does not: +% 1. Multi-region extraction with 'unique_mask_values' (one average per +% integer label in a parcellation), using atlas_labels_combined.img. +% 2. Invariance of the extracted region averages across a write-to-disk and +% reload cycle - a regression guard on the NIfTI I/O path. +% +% Converted from the old standalone script Unit_tests/old_to_integrate/ +% check_roi_extraction.m, which printed PASS/FAIL but never asserted. + +tests = functiontests(localfunctions); +end + + +function test_unique_mask_values_roundtrip(tc) %#ok<*DEFNU> +mask_file = which('atlas_labels_combined.img'); +tc.assumeNotEmpty(mask_file, 'atlas_labels_combined.img not on path'); + +mask_image = fmri_data(mask_file, 'noverbose'); + +% Synthetic per-region timeseries: region i carries the signal ts*sqrt(i), +% so every region has a distinct, known mean trajectory across images. +wh_region = mask_image.dat; +regions = unique(wh_region); +nimgs = 20; +ts = linspace(-10, 10, nimgs); + +dat = mask_image; +dat.dat = zeros(size(mask_image.dat, 1), nimgs); +for i = 1:numel(regions) + idx = wh_region == regions(i); + dat.dat(idx, :) = repmat(ts .* sqrt(i), sum(idx), 1); +end + +cl1 = extract_roi_averages(dat, mask_image, 'unique_mask_values'); +all_reg1 = cat(2, cl1(:).dat); + +% One row per image, one column per non-zero region label. +tc.verifyEqual(size(all_reg1, 1), nimgs); +tc.verifyGreaterThan(size(all_reg1, 2), 1); + +% Round-trip through disk in a scratch folder, then re-extract. +tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +dat.fullpath = fullfile(pwd, 'test_roi_image.nii'); +write(dat, 'overwrite'); +reloaded = fmri_data(dat.fullpath, 'noverbose'); + +cl2 = extract_roi_averages(reloaded, mask_image, 'unique_mask_values'); +all_reg2 = cat(2, cl2(:).dat); + +tc.verifyEqual(size(all_reg2), size(all_reg1)); +% Small absolute slack absorbs single-precision NIfTI write rounding (the +% original script used a 1e-3 relative threshold for the same reason). +tc.verifyEqual(all_reg2, all_reg1, 'AbsTol', 1e-2, 'RelTol', 1e-3); +end diff --git a/CanlabCore/Unit_tests/image_vector/canlab_test_resampling_pattern_expression.m b/CanlabCore/Unit_tests/image_vector/canlab_test_resampling_pattern_expression.m new file mode 100644 index 00000000..79ebfebb --- /dev/null +++ b/CanlabCore/Unit_tests/image_vector/canlab_test_resampling_pattern_expression.m @@ -0,0 +1,48 @@ +function tests = canlab_test_resampling_pattern_expression +%CANLAB_TEST_RESAMPLING_PATTERN_EXPRESSION Pattern expression is stable across resampling. +% +% apply_mask(..., 'pattern_expression') resamples the weight map into the +% data space when the two differ. This test checks that the resulting +% pattern-expression scores barely change regardless of which interpolation +% method does that resampling (default vs nearest vs spline) - i.e. the +% dot-product readout is not an artifact of the interpolation choice. +% +% Uses the SIIPS signature against the emotionreg sample. Skipped if the +% SIIPS image set is not available on the path (it ships with +% Neuroimaging_Pattern_Masks). +% +% Converted from the old standalone visual script +% Unit_tests/old_to_integrate/resampling_pattern_expression_unit_test1.m, +% which only plotted the variants and asserted nothing. + +tests = functiontests(localfunctions); +end + + +function test_pattern_expression_stable_across_interp(tc) %#ok<*DEFNU> +obj = canlab_get_sample_fmri_data(); + +try + siips = load_image_set('siips', 'noverbose'); +catch ME + tc.assumeFail(['SIIPS signature not available on this runner: ' ME.message]); +end +tc.assumeNotEmpty(siips.dat, 'SIIPS loaded but empty'); + +n = size(obj.dat, 2); + +% Default path lets apply_mask resample internally; the other two pre-resample +% the weight map with an explicit interpolation method. +pe_default = apply_mask(obj, siips, 'pattern_expression'); +pe_nearest = apply_mask(obj, resample_space(siips, obj, 'nearest'), 'pattern_expression'); +pe_spline = apply_mask(obj, resample_space(siips, obj, 'spline'), 'pattern_expression'); + +for v = {pe_default, pe_nearest, pe_spline} + tc.verifyEqual(numel(v{1}), n); + tc.verifyTrue(all(isfinite(v{1})), 'pattern expression returned non-finite values'); +end + +% Empirically these correlate at >= 0.9999; 0.99 is a safe regression floor. +tc.verifyGreaterThan(corr(pe_default, pe_nearest), 0.99); +tc.verifyGreaterThan(corr(pe_default, pe_spline), 0.99); +end diff --git a/CanlabCore/Unit_tests/old_to_integrate/check_roi_extraction.m b/CanlabCore/Unit_tests/old_to_integrate/check_roi_extraction.m deleted file mode 100644 index b25734b2..00000000 --- a/CanlabCore/Unit_tests/old_to_integrate/check_roi_extraction.m +++ /dev/null @@ -1,64 +0,0 @@ -function check_roi_extraction() - -% test ROI extraction -mask_image = fmri_data(which('atlas_labels_combined.img')); -dat = mask_image; -wh_region = dat.dat; - -% create timeseries data -- add data "after" 1st image -nimgs = 50; - -n = unique(wh_region); -ts = linspace(-10, 10, nimgs); -v = size(dat.dat, 1); -dat.dat = zeros(v, nimgs); - -for i = 1:length(n) - my_indx = wh_region == n(i); - - dat.dat(my_indx, :) = repmat(ts .* sqrt(i), sum(my_indx), 1); -end - -fprintf('\nDATA GENERATED\n\n'); - -cl1 = extract_roi_averages(dat, mask_image, 'unique_mask_values'); - -% extracted average values for all the regions in mask_image. -all_reg1 = cat(2, cl1(:).dat); - -dat.fullpath = fullfile(pwd, 'test_image.img'); -write(dat); -fprintf('\nSAVED DATA TO %s\n\n', dat.fullpath); - -% reload -fprintf('\nLOADING DATA FROM %s\n\n', dat.fullpath); -dat = fmri_data(dat.fullpath, which('gray_matter_mask.img')); - -% region obj -> extract_data -% test extract_roi_averages -cl2 = extract_roi_averages(dat, mask_image, 'unique_mask_values'); - -% extracted average values for all the regions in mask_image. -all_reg2 = cat(2, cl2(:).dat); - -% uniq_reg = sort(unique(mask_image.dat)); - -% compare roi averages of initially data generated and data loaded from file -if is_equal(all_reg1, all_reg2, 0.001) - fprintf('\nExtracted ROI averages are equal with an error threshold of 0.001\nPASS\n'); -else - fprintf('\nExtracted ROI averages are not equal with an error threshold of 0.001\nFAIL\n'); -end - - -% ------------------------------------------------------ -% INLINE FUNCTION -% ------------------------------------------------------ - - function res = is_equal(A, B, error_threshold) - err = abs((A-B)./B) <= error_threshold; - - res = all(err(:)); - end % is_equal - -end % main function \ No newline at end of file diff --git a/CanlabCore/Unit_tests/old_to_integrate/resampling_pattern_expression_unit_test1.m b/CanlabCore/Unit_tests/old_to_integrate/resampling_pattern_expression_unit_test1.m deleted file mode 100644 index b7f654e0..00000000 --- a/CanlabCore/Unit_tests/old_to_integrate/resampling_pattern_expression_unit_test1.m +++ /dev/null @@ -1,40 +0,0 @@ -function resampling_pattern_expression_unit_test1 - -obj = load_image_set('emotionreg'); % test dataset -siips = load_image_set('siips'); % sampled to higher-res mask on loading -siips_resamp = resample_space(siips, obj, 'spline'); -siips_orig = fmri_data(siips.image_names); % original pattern - -val = []; -val(:, 1) = apply_mask(obj, siips_orig, 'pattern_expression'); -val(:, 2) = apply_mask(obj, resample_space(siips_orig, obj), 'pattern_expression'); -val(:, 3) = apply_mask(obj, resample_space(siips_orig, obj, 'nearest'), 'pattern_expression'); -val(:, 4) = apply_mask(obj, resample_space(siips_orig, obj, 'spline'), 'pattern_expression'); - -val(:, 5) = apply_mask(obj, siips, 'pattern_expression'); -val(:, 6) = apply_mask(obj, resample_space(siips, obj), 'pattern_expression'); -val(:, 7) = apply_mask(obj, resample_space(siips, obj, 'nearest'), 'pattern_expression'); -val(:, 8) = apply_mask(obj, resample_space(siips, obj, 'spline'), 'pattern_expression'); - -val(:, 9) = apply_mask(obj, siips, 'pattern_expression'); -val(:, 10) = apply_mask(obj, resample_space(siips, obj), 'pattern_expression'); -val(:, 11) = apply_mask(obj, resample_space(siips, obj, 'nearest'), 'pattern_expression'); -val(:, 12) = apply_mask(obj, resample_space(siips, obj, 'spline'), 'pattern_expression'); - -val(:, 13) = apply_mask(obj, siips_resamp, 'pattern_expression'); -sdat = apply_siips(obj); -val(:, 14) = sdat{1}; - -m = mean(val); -s = std(val); -c = corr(val); - -create_figure('siips comparison', 1, 3); -errorbar(m, s); title('means and standard errors') -subplot(1, 3, 2); plot(m); title('means (zooming in on differences') -subplot(1, 3, 3); imagesc(c); colorbar -title('Correlations among variants') - -set(gca, 'YDir', 'reverse'); axis tight - -end diff --git a/CanlabCore/Unit_tests/predictive_model/canlab_test_predictive_model_api.m b/CanlabCore/Unit_tests/predictive_model/canlab_test_predictive_model_api.m new file mode 100644 index 00000000..515bf5ac --- /dev/null +++ b/CanlabCore/Unit_tests/predictive_model/canlab_test_predictive_model_api.m @@ -0,0 +1,311 @@ +function tests = canlab_test_predictive_model_api +%CANLAB_TEST_PREDICTIVE_MODEL_API End-to-end test of the @predictive_model public API. +% +% Converted from Unit_tests/predictive_model_unit_test.m. Exercises the +% sklearn-style @predictive_model object surface on the DPSP Hot vs Warm +% dataset: construction, fit/predict/score, crossval (cv_splitter), bootstrap, +% weight_map_object, permutation_test, calibrate/predict_proba, +% select_features, grid_search, stability_selection, @pipeline, +% fmri_data.predict 'newapi' routing, and the regression/report/summary path. +% Visualisation methods (plot/confusionchart/montage) are smoke-tested and +% skipped on a headless runner. Complements the xval_* wrapper tests, which +% cover the legacy entry points rather than the object's own methods. +% +% The shared fit -> crossval -> bootstrap chain is computed once in setupOnce +% and cached; method-specific tests build their own small models. Small +% nboot/nperm/grid values keep this to a few minutes; it is still much slower +% than the rest of the suite. Skipped if DPSP_hotwarm is not on the path. + +tests = functiontests(localfunctions); +end + + +% ===================================================================== +% Fixtures +% ===================================================================== + +function setupOnce(tc) %#ok<*DEFNU> +tc.assumeNotEmpty(which('load_image_set'), 'load_image_set not on path'); +try + hw_obj = load_image_set('DPSP_hotwarm', 'noverbose'); +catch ME + tc.assumeFail(['DPSP_hotwarm sample not available: ' ME.message]); + return +end + +% Restrict to a sparse gray-matter ROI so the many cross-validations below +% are fast: mask to gray matter, then deterministically thin to every 8th +% voxel (~20k features instead of ~195k). hw_obj and X stay in lockstep, so +% weight_map_object / montage still back-project consistently. Falls back to +% the full volume if the gray-matter mask is not on the path. +gm_file = which('gray_matter_mask.img'); +if ~isempty(gm_file) + hw_obj = remove_empty(apply_mask(hw_obj, gm_file)); + thin = get_wh_image(hw_obj, 1); + thin.dat = zeros(size(hw_obj.dat, 1), 1); + thin.dat(1:8:end) = 1; + hw_obj = remove_empty(apply_mask(hw_obj, thin)); +end + +X = double(hw_obj.dat'); +Y = hw_obj.Y; +id = grp2idx(hw_obj.metadata_table.subj_id); + +tc.TestData.hw_obj = hw_obj; +tc.TestData.X = X; +tc.TestData.Y = Y; +tc.TestData.id = id; + +% Shared expensive chain, computed once: in-sample fit, cross-validation, +% and bootstrap on the cross-validated model. +tc.TestData.pm_insample = fit( ... + predictive_model('algorithm', 'svm', 'task', 'classification', 'random_state', 0), ... + X, Y, 'id', id); + +pm_cv = predictive_model('algorithm', 'svm', 'task', 'classification', ... + 'modeloptions', {'KernelFunction', 'linear'}, 'random_state', 0); +pm_cv = crossval(pm_cv, X, Y, ... + 'cv', cv_splitter.stratified_group_kfold(5), 'groups', id, ... + 'scoring', 'balanced_accuracy'); +tc.TestData.pm_cv = pm_cv; + +% Minimum-viable nboot: just enough to compute boot std / z / p and exercise +% the empirical p-floor. This is a smoke test, not real inference. +tc.TestData.pm_boot = bootstrap(pm_cv, X, Y, 'nboot', 5, 'groups', id, 'verbose', false); +end + + +function setup(tc) +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); +end + + +function teardown(tc) +close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end +end + + +% ===================================================================== +% Construction / in-sample fit / predict / score +% ===================================================================== + +function test_construction_hyperparameters(tc) +pm = predictive_model('algorithm', 'svm', 'task', 'classification', 'random_state', 0); +tc.verifyEqual(pm.algorithm, 'svm'); +tc.verifyEqual(pm.task, 'classification'); +tc.verifyFalse(pm.is_fitted, 'should start unfitted'); +tc.verifyTrue(pm.is_classifier); +end + + +function test_fit_insample(tc) +pm = tc.TestData.pm_insample; +tc.verifyTrue(pm.is_fitted); +tc.verifyEqual(pm.fit_type, 'insample'); +tc.verifyLessThanOrEqual(numel(pm.weights.w), size(tc.TestData.X, 2)); +end + + +function test_predict(tc) +[yhat, scores] = predict(tc.TestData.pm_insample, tc.TestData.X); +tc.verifyNumElements(yhat, numel(tc.TestData.Y)); +tc.verifyEqual(size(scores, 1), numel(tc.TestData.Y)); +end + + +function test_score_balanced_accuracy(tc) +pm = tc.TestData.pm_insample; +pm.scorer = cv_scorer.balanced_accuracy(); +v = score(pm, tc.TestData.X, tc.TestData.Y); +tc.verifyGreaterThanOrEqual(v, 0); +tc.verifyLessThanOrEqual(v, 1); +end + + +% ===================================================================== +% crossval / bootstrap / weight map +% ===================================================================== + +function test_crossval_stratified_group_kfold(tc) +pm = tc.TestData.pm_cv; +tc.verifyEqual(pm.fit_type, 'crossval'); +tc.verifyEqual(pm.cv_partition.nfolds, 5); +tc.verifyNumElements(pm.fitted_values.yfit, numel(tc.TestData.Y)); +tc.verifyFalse(isnan(pm.error_metrics.balanced_accuracy.value), 'metric populated'); +end + + +function test_bootstrap_weight_stats(tc) +nboot = 5; % must match setupOnce +pm = tc.TestData.pm_boot; +tc.verifyEqual(size(pm.weights.boot_w, 2), nboot, 'boot_w should have nboot columns'); +tc.verifyNumElements(pm.weights.z, numel(pm.weights.w)); +tc.verifyNumElements(pm.weights.p, numel(pm.weights.w)); +tc.verifyTrue(all(pm.weights.p >= 0 & pm.weights.p <= 1), 'p in [0,1]'); +tc.verifyTrue(all(abs(pm.weights.z) < 1e6), 'z floored, not 1e15'); +p_floor = 2 / (nboot + 1); +tc.verifyGreaterThanOrEqual(min(pm.weights.p), p_floor - eps, 'empirical p floor'); +end + + +function test_weight_map_object(tc) +pm = tc.TestData.pm_boot; +hw = tc.TestData.hw_obj; +[pm_w, si_full] = weight_map_object(pm, hw); +[~, si_thr] = weight_map_object(pm, hw, 'use', 'thresh_fdr'); +tc.verifyClass(si_full, 'statistic_image'); +tc.verifyEqual(size(si_full.dat, 1), size(hw.dat, 1)); +tc.verifyEqual(size(si_thr.dat, 1), size(hw.dat, 1)); +tc.verifyNotEmpty(pm_w.weights.weight_obj, 'weight_obj cached on pm'); +tc.verifyClass(pm_w.weights.weight_obj, 'statistic_image'); +end + + +% ===================================================================== +% permutation / calibration / feature selection +% ===================================================================== + +function test_permutation_test_within_subjects(tc) +pm = predictive_model('algorithm', 'svm', 'task', 'classification', ... + 'modeloptions', {'KernelFunction', 'linear'}, 'random_state', 0); +pm.cv = cv_splitter.stratified_group_kfold(3); +nperm = 3; % minimum-viable: enough to form a null and a defined p-value +pm = permutation_test(pm, tc.TestData.X, tc.TestData.Y, ... + 'nperm', nperm, 'groups', tc.TestData.id, 'verbose', false); +tc.verifyNumElements(pm.permutation_results.null_scores, nperm); +tc.verifyGreaterThan(pm.permutation_results.p_value, 0); +tc.verifyLessThanOrEqual(pm.permutation_results.p_value, 1); +% DPSP is paired (Hot+Warm per subject) -> auto should pick within_subjects. +tc.verifyEqual(pm.permutation_results.permutation, 'within_subjects'); +end + + +function test_calibrate_predict_proba(tc) +pm_cal = calibrate(tc.TestData.pm_cv, tc.TestData.X, tc.TestData.Y); +tc.verifyTrue(isfield(pm_cal.fitted_values, 'calibrator'), 'calibrator present'); +P = predict_proba(pm_cal, tc.TestData.X); +tc.verifyTrue(all(P >= 0 & P <= 1), 'predict_proba in [0,1]'); +end + + +function test_select_features(tc) +pm = predictive_model('algorithm', 'svm', 'task', 'classification'); +pm = select_features(pm, tc.TestData.X, tc.TestData.Y, 'k', 500, 'verbose', false); +tc.verifyEqual(pm.diagnostics.feature_selection.n_selected, 500); +end + + +% ===================================================================== +% grid search / stability selection / pipeline / newapi +% ===================================================================== + +function test_grid_search(tc) +pm = predictive_model('algorithm', 'svm', 'task', 'classification', ... + 'cv', cv_splitter.stratified_group_kfold(3)); +pm = grid_search(pm, tc.TestData.X, tc.TestData.Y, ... + struct('BoxConstraint', [0.1 1]), 'groups', tc.TestData.id, 'verbose', false); +tc.verifyTrue(isfield(pm.diagnostics.grid_search, 'best_score')); +end + + +function test_stability_selection(tc) +pm = stability_selection( ... + predictive_model('algorithm', 'linear_svm', 'task', 'classification'), ... + tc.TestData.X, tc.TestData.Y, 'nboot', 5, 'k', 1000, 'threshold', 0.6, ... + 'groups', tc.TestData.id, 'verbose', false); +ss = pm.diagnostics.stability_selection; +tc.verifyNumElements(ss.selection_freq, size(tc.TestData.X, 2)); +tc.verifyTrue(all(ss.selection_freq >= 0 & ss.selection_freq <= 1), 'freq in [0,1]'); +end + + +function test_pipeline_pca_svm(tc) +est = predictive_model('algorithm', 'svm', 'task', 'classification'); +pipe = pipeline({ {'pca', 'k', 20} }, est); +pipe = crossval(pipe, tc.TestData.X, tc.TestData.Y, ... + 'groups', tc.TestData.id, 'cv', cv_splitter.stratified_group_kfold(5)); +tc.verifyEqual(pipe.fit_type, 'crossval'); +tc.verifyNumElements(pipe.weights.w, size(tc.TestData.X, 2), ... + 'pipeline should back-project weights to full feature space'); +si_pipe = weight_map_object(pipe, tc.TestData.hw_obj); +tc.verifyClass(si_pipe, 'statistic_image'); +end + + +function test_fmri_data_predict_newapi(tc) +[cverr, st, ~, pm_api] = predict(tc.TestData.hw_obj, 'algorithm_name', 'cv_svm', ... + 'nfolds', 5, 'newapi', 'verbose', 0); +tc.verifyClass(pm_api, 'predictive_model'); +tc.verifyEqual(pm_api.fit_type, 'crossval'); +tc.verifyTrue(isfield(st, 'dist_from_hyperplane_xval')); +tc.verifyTrue(isfinite(cverr)); +end + + +% ===================================================================== +% regression metrics / report_accuracy / summary +% ===================================================================== + +function test_regression_metrics_and_reports(tc) +rng(7); +Y = tc.TestData.Y; +Yr = double(Y) + 0.5 * randn(size(Y)); +pm_reg = predictive_model('algorithm', 'pcr', 'task', 'regression'); +pm_reg = crossval(pm_reg, tc.TestData.X, Yr, ... + 'cv', cv_splitter.group_kfold(5), 'groups', tc.TestData.id); + +tc.verifyTrue(isfield(pm_reg.error_metrics, 'predicted_r2')); +tc.verifyTrue(isfield(pm_reg.error_metrics, 'out_of_sample_r2')); + +% predicted_r2 must equal 1 - PRESS/SST computed by hand. +yf = pm_reg.fitted_values.yfit(:); +yo = Yr(:); +pr2_manual = 1 - sum((yo - yf).^2) / sum((yo - mean(yo)).^2); +tc.verifyEqual(pm_reg.error_metrics.predicted_r2.value, pr2_manual, 'AbsTol', 1e-9); + +acc_r = report_accuracy(pm_reg, 'noverbose'); +tc.verifyEqual(acc_r.task, 'regression'); +tc.verifyFalse(isnan(acc_r.predicted_r2)); + +s_reg = summary(pm_reg, 'noverbose'); +tc.verifyEqual(s_reg.provenance.fit_type, 'crossval'); + +% Classification report_accuracy: ROC-derived fields populated. +acc_c = report_accuracy(tc.TestData.pm_cv, 'noverbose'); +tc.verifyEqual(acc_c.task, 'classification'); +tc.verifyFalse(isnan(acc_c.auc)); +tc.verifyFalse(isnan(acc_c.npv)); +end + + +% ===================================================================== +% Visualisation (smoke; skipped on a headless runner) +% ===================================================================== + +function test_rocplot_metrics(tc) +% rocplot('noplot') returns metrics without a display. +ROC = rocplot(tc.TestData.pm_boot, 'noplot'); +tc.verifyTrue(isfield(ROC, 'AUC')); +tc.verifyGreaterThanOrEqual(ROC.AUC, 0); +tc.verifyLessThanOrEqual(ROC.AUC, 1); +end + + +function test_display_methods_smoke(tc) +pm = tc.TestData.pm_boot; +hw = tc.TestData.hw_obj; +try + h = plot(pm); close(h); % classification dispatcher (violin + ROC) + confusionchart(pm); close(gcf); + montage(pm, hw); close all force; % delegates to weight-map montage +catch ME + if strcmp(canlab_classify_environment_error(ME), 'genuine') + rethrow(ME); + end + tc.assumeFail(['display method needs a graphics environment: ' ME.message]); +end +end diff --git a/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_SVM.m b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_SVM.m new file mode 100644 index 00000000..93d2753d --- /dev/null +++ b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_SVM.m @@ -0,0 +1,105 @@ +function tests = canlab_test_xval_SVM +%CANLAB_TEST_XVAL_SVM Cross-validated SVM classification via xval_SVM. +% +% Converted from Unit_tests/xval_SVM_unit_test.m. Drives xval_SVM on the DPSP +% Hot vs Warm single-subject maps (a paired, within-person design) and checks +% the returned @predictive_model object: class/state, canonical-path fields, +% shapes/ranges, fit_type + omitted markers, validate_object, and clone. +% +% The model is fit once in setupOnce and cached; each test reads the cached +% object. Skipped if the DPSP sample data is not on the path. + +tests = functiontests(localfunctions); +end + + +function setupOnce(tc) %#ok<*DEFNU> +[hot, warm, ok] = canlab_get_dpsp_hot_warm(); +tc.assumeTrue(ok, 'DPSP sample data not on path'); + +n_hot = size(hot.dat, 2); +n_warm = size(warm.dat, 2); + +rng(0); +p = size(hot.dat, 1); +keep_vox = randsample(p, min(5000, p)); +X = double([hot.dat(keep_vox, :) warm.dat(keep_vox, :)])'; +Y = [ones(n_hot, 1); -ones(n_warm, 1)]; +id = [(1:n_hot)'; (1:n_warm)']; % each subject contributes both conditions + +tc.TestData.X = X; +tc.TestData.Y = Y; +tc.TestData.id = id; +tc.TestData.pm = xval_SVM(X, Y, id, 'nooptimize', 'norepeats', ... + 'nobootstrap', 'noverbose', 'noplot'); +end + + +function test_class_and_state(tc) +pm = tc.TestData.pm; +tc.verifyClass(pm, 'predictive_model'); +tc.verifyTrue(pm.is_fitted, 'is_fitted should be true'); +tc.verifyTrue(pm.is_classifier, 'binary Y should give is_classifier true'); +end + + +function test_canonical_path_fields_populated(tc) +pm = tc.TestData.pm; +tc.verifyNotEmpty(pm.Y); +tc.verifyNotEmpty(pm.id); +tc.verifyNotEmpty(pm.fitted_values.yfit); +tc.verifyNotEmpty(pm.fitted_values.dist_from_hyperplane_xval); +tc.verifyNotEmpty(pm.weights.w); +tc.verifyNotEmpty(pm.error_metrics.crossval_accuracy.value); +tc.verifyNotEmpty(pm.error_metrics.d_singleinterval.value); +tc.verifyNotEmpty(pm.ml_model); +tc.verifyNotEmpty(pm.cv_partition.trIdx); +tc.verifyNotEmpty(pm.cv_partition.teIdx); +tc.verifyNotEmpty(pm.cv_partition.nfolds); +end + + +function test_shapes_and_ranges(tc) +pm = tc.TestData.pm; +n = numel(tc.TestData.Y); +nfeat = size(tc.TestData.X, 2); + +tc.verifyNumElements(pm.fitted_values.yfit, n); +tc.verifyNumElements(pm.fitted_values.dist_from_hyperplane_xval, n); +tc.verifyEqual(size(pm.weights.w, 1), nfeat, ... + 'weights.w should have one row per feature'); + +cv_acc = pm.error_metrics.crossval_accuracy.value; +tc.verifyGreaterThanOrEqual(cv_acc, 0); +tc.verifyLessThanOrEqual(cv_acc, 100); + +% DPSP is paired, so within-person scoring should be populated. +tc.verifyFalse(isnan(pm.error_metrics.crossval_accuracy_within.value)); +tc.verifyFalse(isnan(pm.error_metrics.d_within.value)); +tc.verifyTrue(pm.diagnostics.mult_obs_within_person, ... + 'should detect multiple observations per id'); +end + + +function test_fit_type_and_omitted_markers(tc) +pm = tc.TestData.pm; +tc.verifyEqual(pm.fit_type, 'crossval'); +tc.verifyClass(pm.omitted_cases, 'logical'); +tc.verifyClass(pm.omitted_features, 'logical'); +tc.verifyNumElements(pm.omitted_cases, numel(tc.TestData.Y)); +tc.verifyNumElements(pm.omitted_features, size(tc.TestData.X, 2)); +end + + +function test_validate_object_accepts(tc) +% validate_object throws on an invalid object; reaching the end is a pass. +tc.TestData.pm.validate_object('noverbose'); +end + + +function test_clone_clears_fitted_state(tc) +pm2 = clone(tc.TestData.pm); +tc.verifyFalse(pm2.is_fitted, 'clone should clear fitted state'); +tc.verifyEqual(pm2.modeloptions, tc.TestData.pm.modeloptions, ... + 'clone should preserve modeloptions'); +end diff --git a/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_SVR.m b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_SVR.m new file mode 100644 index 00000000..c5dafebb --- /dev/null +++ b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_SVR.m @@ -0,0 +1,72 @@ +function tests = canlab_test_xval_SVR +%CANLAB_TEST_XVAL_SVR Cross-validated SVR regression via xval_SVR. +% +% Converted from Unit_tests/xval_SVR_unit_test.m. Predicts a synthetic +% continuous outcome from the DPSP Hot-Warm contrast maps via cross-validated +% linear SVR and checks the returned @predictive_model object. Model fit once +% in setupOnce; skipped if the DPSP sample data is not on the path. + +tests = functiontests(localfunctions); +end + + +function setupOnce(tc) %#ok<*DEFNU> +[hot, warm, ok] = canlab_get_dpsp_hot_warm(); +tc.assumeTrue(ok, 'DPSP sample data not on path'); + +hot_vs_warm = image_math(hot, warm, 'minus'); + +rng(0); +[p, n] = size(hot_vs_warm.dat); +keep_vox = randsample(p, min(3000, p)); +X = double(hot_vs_warm.dat(keep_vox, :))'; + +% Synthetic outcome predictable from a sparse subset of features. +b_true = zeros(size(X, 2), 1); +b_true(1:30) = randn(30, 1); +Y = X * b_true + 0.5 * std(X * b_true) * randn(n, 1); +id = (1:n)'; + +tc.TestData.X = X; +tc.TestData.Y = Y; +tc.TestData.pm = xval_SVR(X, Y, id, 'nooptimize', 'norepeats', ... + 'nobootstrap', 'noverbose', 'noplot'); +end + + +function test_class_and_state(tc) +pm = tc.TestData.pm; +tc.verifyClass(pm, 'predictive_model'); +tc.verifyTrue(pm.is_fitted, 'is_fitted should be true'); +tc.verifyTrue(pm.is_regressor, 'continuous Y should give is_regressor true'); +end + + +function test_canonical_path_fields_populated(tc) +pm = tc.TestData.pm; +tc.verifyNotEmpty(pm.Y); +tc.verifyNotEmpty(pm.fitted_values.yfit); +tc.verifyNotEmpty(pm.weights.w); +tc.verifyNotEmpty(pm.error_metrics.prediction_outcome_r.value); +tc.verifyNotEmpty(pm.ml_model); +end + + +function test_weight_shape(tc) +pm = tc.TestData.pm; +tc.verifyEqual(size(pm.weights.w, 1), size(tc.TestData.X, 2), ... + 'weights.w should have one row per feature'); +end + + +function test_fit_type_and_omitted_markers(tc) +pm = tc.TestData.pm; +tc.verifyEqual(pm.fit_type, 'crossval'); +tc.verifyClass(pm.omitted_cases, 'logical'); +tc.verifyClass(pm.omitted_features, 'logical'); +end + + +function test_validate_object_accepts(tc) +tc.TestData.pm.validate_object('noverbose'); +end diff --git a/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_discriminant_classifier.m b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_discriminant_classifier.m new file mode 100644 index 00000000..2fec72d4 --- /dev/null +++ b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_discriminant_classifier.m @@ -0,0 +1,66 @@ +function tests = canlab_test_xval_discriminant_classifier +%CANLAB_TEST_XVAL_DISCRIMINANT_CLASSIFIER Cross-validated LDA classification. +% +% Converted from Unit_tests/xval_discriminant_classifier_unit_test.m. Drives +% cross-validated LDA (fitcdiscr) on the DPSP Hot vs Warm single-subject maps +% (100 random voxels, so n > p) and checks the returned @predictive_model +% object. Model fit once in setupOnce; skipped if DPSP data is not on the path. + +tests = functiontests(localfunctions); +end + + +function setupOnce(tc) %#ok<*DEFNU> +[hot, warm, ok] = canlab_get_dpsp_hot_warm(); +tc.assumeTrue(ok, 'DPSP sample data not on path'); + +rng(0); +p = size(hot.dat, 1); +keep_vox = randsample(p, 100); % LDA needs n > p +X = double([hot.dat(keep_vox, :) warm.dat(keep_vox, :)])'; +labels = int32([ones(size(hot.dat, 2), 1); 2 * ones(size(warm.dat, 2), 1)]); + +tc.TestData.pm = xval_discriminant_classifier(X, labels, 'nFolds', 5, ... + 'verbose', false, 'doplot', false); +end + + +function test_class_and_state(tc) +pm = tc.TestData.pm; +tc.verifyClass(pm, 'predictive_model'); +tc.verifyTrue(pm.is_fitted, 'is_fitted should be true'); +end + + +function test_canonical_path_fields_populated(tc) +pm = tc.TestData.pm; +tc.verifyNotEmpty(pm.Y); +tc.verifyNotEmpty(pm.fitted_values.yfit); +tc.verifyNotEmpty(pm.fitted_values.predictions); +tc.verifyNotEmpty(pm.fitted_values.Y_per_fold); +tc.verifyNotEmpty(pm.error_metrics.accuracy.value); +tc.verifyNotEmpty(pm.error_metrics.overallAccuracy.value); +tc.verifyNotEmpty(pm.cv_partition.trIdx); +tc.verifyNotEmpty(pm.cv_partition.teIdx); +tc.verifyNotEmpty(pm.cv_partition.nfolds); +end + + +function test_fold_models_match_nfolds(tc) +pm = tc.TestData.pm; +tc.verifyNumElements(pm.fold_models, pm.cv_partition.nfolds, ... + 'fold_models length should equal nfolds'); +end + + +function test_fit_type_and_omitted_markers(tc) +pm = tc.TestData.pm; +tc.verifyEqual(pm.fit_type, 'crossval'); +tc.verifyClass(pm.omitted_cases, 'logical'); +tc.verifyClass(pm.omitted_features, 'logical'); +end + + +function test_validate_object_accepts(tc) +tc.TestData.pm.validate_object('noverbose'); +end diff --git a/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_regression_multisubject.m b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_regression_multisubject.m new file mode 100644 index 00000000..d6bebbe3 --- /dev/null +++ b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_regression_multisubject.m @@ -0,0 +1,74 @@ +function tests = canlab_test_xval_regression_multisubject +%CANLAB_TEST_XVAL_REGRESSION_MULTISUBJECT OLS+PCA multisubject regression. +% +% Converted from Unit_tests/xval_regression_multisubject_unit_test.m. Runs +% xval_regression_multisubject ('ols' + 'pca') on the DPSP Hot-Warm contrast +% maps against a synthetic continuous outcome and checks the returned +% @predictive_model object. Model fit once in setupOnce; skipped if the DPSP +% sample data is not on the path. + +tests = functiontests(localfunctions); +end + + +function setupOnce(tc) %#ok<*DEFNU> +[hot, warm, ok] = canlab_get_dpsp_hot_warm(); +tc.assumeTrue(ok, 'DPSP sample data not on path'); + +hot_vs_warm = image_math(hot, warm, 'minus'); + +rng(0); +[p, n] = size(hot_vs_warm.dat); +keep_vox = randsample(p, min(5000, p)); +X = double(hot_vs_warm.dat(keep_vox, :))'; % n x v + +b_true = zeros(size(X, 2), 1); +b_true(1:50) = randn(50, 1); +Y = X * b_true + 0.5 * std(X * b_true) * randn(n, 1); + +tc.TestData.X = X; +tc.TestData.Y = Y; +tc.TestData.n = n; +tc.TestData.pm = xval_regression_multisubject('ols', {Y}, {X}, ... + 'pca', 'ndims', 10, 'holdout_method', 'balanced4', 'noverbose'); +end + + +function test_class_and_state(tc) +pm = tc.TestData.pm; +tc.verifyClass(pm, 'predictive_model'); +tc.verifyTrue(pm.is_fitted, 'is_fitted should be true'); +end + + +function test_canonical_path_fields_populated(tc) +pm = tc.TestData.pm; +tc.verifyNotEmpty(pm.fitted_values.subjfit); +tc.verifyNotEmpty(pm.weights.subjbetas); +tc.verifyNotEmpty(pm.weights.mean_vox_weights); +tc.verifyNotEmpty(pm.error_metrics.pred_err.value); +tc.verifyNotEmpty(pm.error_metrics.pred_err_null.value); +tc.verifyNotEmpty(pm.error_metrics.var_reduction.value); +tc.verifyNotEmpty(pm.error_metrics.r_each_subject.value); +tc.verifyNotEmpty(pm.Y); +tc.verifyNotEmpty(pm.inputParameters); +end + + +function test_shapes_and_variance_explained(tc) +pm = tc.TestData.pm; +tc.verifyNumElements(pm.fitted_values.subjfit, 1, ... + 'subjfit cell length should equal number of datasets'); +tc.verifyNumElements(pm.fitted_values.subjfit{1}, tc.TestData.n); +tc.verifyEqual(size(pm.weights.mean_vox_weights, 1), size(tc.TestData.X, 2), ... + 'mean_vox_weights should have one row per voxel'); + +r2 = pm.error_metrics.r_squared.value; +tc.verifyNotEmpty(r2); +tc.verifyGreaterThan(r2(1), 0, 'OLS+PCA should explain non-trivial variance'); +end + + +function test_validate_object_accepts(tc) +tc.TestData.pm.validate_object('noverbose'); +end diff --git a/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_regression_multisubject_featureselect.m b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_regression_multisubject_featureselect.m new file mode 100644 index 00000000..08311c82 --- /dev/null +++ b/CanlabCore/Unit_tests/predictive_model/canlab_test_xval_regression_multisubject_featureselect.m @@ -0,0 +1,75 @@ +function tests = canlab_test_xval_regression_multisubject_featureselect +%CANLAB_TEST_XVAL_REGRESSION_MULTISUBJECT_FEATURESELECT Feature-select regression. +% +% Converted from Unit_tests/xval_regression_multisubject_featureselect_unit_test.m. +% Runs xval_regression_multisubject_featureselect ('ols') on the DPSP Hot-Warm +% contrast maps against a synthetic continuous outcome. Beyond the usual +% predictive_model checks, it pins the per-fold weight padding: weights.w_perfold +% must span the FULL feature space, not the selected-feature count (the bug this +% test originally guarded). Skipped if the DPSP sample data is not on the path. + +tests = functiontests(localfunctions); +end + + +function setupOnce(tc) %#ok<*DEFNU> +[hot, warm, ok] = canlab_get_dpsp_hot_warm(); +tc.assumeTrue(ok, 'DPSP sample data not on path'); + +hvw = image_math(hot, warm, 'minus'); + +rng(0); +[p, n] = size(hvw.dat); +keep_vox = randsample(p, min(2000, p)); +X = double(hvw.dat(keep_vox, :))'; + +b_true = zeros(size(X, 2), 1); +b_true(1:30) = randn(30, 1); +Y = X * b_true + 0.5 * std(X * b_true) * randn(n, 1); + +tc.TestData.X = X; +tc.TestData.Y = Y; +tc.TestData.pm = xval_regression_multisubject_featureselect('ols', {Y}, {X}, ... + 'holdout_method', 'balanced4', 'noverbose'); +end + + +function test_class_and_state(tc) +pm = tc.TestData.pm; +tc.verifyClass(pm, 'predictive_model'); +tc.verifyTrue(pm.is_fitted, 'is_fitted should be true'); +end + + +function test_canonical_path_fields_populated(tc) +% mean_vox_weights is only populated when pcsquash is enabled; this no-PCA +% test does not exercise that path, so it is not checked here. +pm = tc.TestData.pm; +tc.verifyNotEmpty(pm.fitted_values.subjfit); +tc.verifyNotEmpty(pm.weights.w_perfold); +tc.verifyNotEmpty(pm.error_metrics.pred_err.value); +tc.verifyNotEmpty(pm.error_metrics.pred_err_null.value); +tc.verifyNotEmpty(pm.error_metrics.var_reduction.value); +tc.verifyNotEmpty(pm.Y); +end + + +function test_w_perfold_spans_full_feature_space(tc) +% Regression guard: w_perfold must have one row per ORIGINAL feature, not per +% selected feature. +pm = tc.TestData.pm; +tc.verifyEqual(size(pm.weights.w_perfold, 1), size(tc.TestData.X, 2), ... + 'weights.w_perfold should span the full feature count'); +end + + +function test_variance_explained(tc) +r2 = tc.TestData.pm.error_metrics.r_squared.value; +tc.verifyNotEmpty(r2); +tc.verifyGreaterThan(r2(1), 0, 'expected r_squared > 0'); +end + + +function test_validate_object_accepts(tc) +tc.TestData.pm.validate_object('noverbose'); +end diff --git a/CanlabCore/Unit_tests/predictive_model_unit_test.m b/CanlabCore/Unit_tests/predictive_model_unit_test.m deleted file mode 100644 index c6891160..00000000 --- a/CanlabCore/Unit_tests/predictive_model_unit_test.m +++ /dev/null @@ -1,195 +0,0 @@ -function predictive_model_unit_test() -% predictive_model_unit_test End-to-end test of the new sklearn-style -% @predictive_model API, using DPSP_hotwarm as the test bed. -% -% Small nboot / nperm values are used so the test runs in ~minutes, -% not hours. The point is to exercise every public method, not to -% produce final-quality inference. -% -% Run: cd to this directory, then `predictive_model_unit_test`. - - fprintf('=== predictive_model_unit_test ===\n'); - - % --- Load DPSP Hot vs Warm --- - hw_obj = load_image_set('DPSP_hotwarm', 'noverbose'); - X = double(hw_obj.dat'); - Y = hw_obj.Y; - id = grp2idx(hw_obj.metadata_table.subj_id); - fprintf(' Loaded DPSP_hotwarm: %d obs, %d features, %d subjects\n', ... - numel(Y), size(X, 2), numel(unique(id))); - - % --- 1. Construction + hyperparameters --- - pm = predictive_model('algorithm', 'svm', ... - 'task', 'classification', ... - 'random_state', 0); - assert(strcmp(pm.algorithm, 'svm'), 'algorithm'); - assert(strcmp(pm.task, 'classification'), 'task'); - assert(~pm.is_fitted, 'starts unfitted'); - assert(pm.is_classifier, 'is_classifier'); - fprintf(' 1. Construction + clone OK\n'); - - % --- 2. fit (in-sample) --- - pm = fit(pm, X, Y, 'id', id); - assert(pm.is_fitted, 'is_fitted after fit'); - assert(strcmp(pm.fit_type, 'insample'), 'fit_type=insample'); - assert(numel(pm.weights.w) <= size(X, 2), 'weights.w shape'); - fprintf(' 2. fit OK (fit_type=%s, |w|=%d)\n', pm.fit_type, numel(pm.weights.w)); - - % --- 3. predict --- - [yhat, scores] = predict(pm, X); - assert(numel(yhat) == numel(Y), 'predict yhat length'); - assert(size(scores, 1) == numel(Y), 'predict scores length'); - fprintf(' 3. predict OK (train acc=%.2f)\n', mean(yhat == Y)); - - % --- 4. score --- - pm.scorer = cv_scorer.balanced_accuracy(); - v = score(pm, X, Y); - assert(v >= 0 && v <= 1, 'score in [0,1]'); - fprintf(' 4. score (balanced_accuracy) = %.3f\n', v); - - % --- 5. crossval with stratified_group_kfold --- - pm = predictive_model('algorithm','svm','task','classification', ... - 'modeloptions', {'KernelFunction','linear'}, ... - 'random_state', 0); - pm = crossval(pm, X, Y, ... - 'cv', cv_splitter.stratified_group_kfold(5), ... - 'groups', id, ... - 'scoring', 'balanced_accuracy'); - assert(strcmp(pm.fit_type, 'crossval'), 'fit_type=crossval'); - assert(pm.cv_partition.nfolds == 5, '5 folds'); - assert(numel(pm.fitted_values.yfit) == numel(Y), 'cv yfit length'); - assert(~isnan(pm.error_metrics.balanced_accuracy.value), 'metric populated'); - fprintf(' 5. crossval (stratified_group_kfold(5)) bal_acc = %.3f\n', ... - pm.error_metrics.balanced_accuracy.value); - - % --- 6. bootstrap (small nboot for speed) --- - pm = bootstrap(pm, X, Y, 'nboot', 25, 'groups', id, 'verbose', false); - assert(size(pm.weights.boot_w, 2) == 25, 'nboot=25 samples'); - assert(numel(pm.weights.z) == numel(pm.weights.w), 'z length matches w'); - assert(numel(pm.weights.p) == numel(pm.weights.w), 'p length matches w'); - assert(all(pm.weights.p >= 0 & pm.weights.p <= 1), 'p in [0,1]'); - assert(all(abs(pm.weights.z) < 1e6), 'z floored not 1e15'); - p_floor = 2 / (25 + 1); - assert(min(pm.weights.p) >= p_floor - eps, 'empirical p floor'); - fprintf(' 6. bootstrap (nboot=25) OK\n'); - fprintf(' z range = [%g .. %g]; p floor = %.4f\n', ... - min(pm.weights.z), max(pm.weights.z), min(pm.weights.p)); - - % --- 7. weight_map_object: pre- and post-bootstrap thresholded, + caching --- - [pm_w, si_full] = weight_map_object(pm, hw_obj); - [~, si_thr] = weight_map_object(pm, hw_obj, 'use', 'thresh_fdr'); - assert(isa(si_full, 'statistic_image'), 'si_full class'); - assert(size(si_full.dat, 1) == size(hw_obj.dat, 1), 'si_full voxel count'); - assert(size(si_thr.dat, 1) == size(hw_obj.dat, 1), 'si_thr voxel count'); - assert(~isempty(pm_w.weights.weight_obj), 'weight_obj cached on pm'); - assert(isa(pm_w.weights.weight_obj, 'statistic_image'), 'cached weight_obj class'); - fprintf(' 7. weight_map_object (full + thresh_fdr + cache) OK\n'); - fprintf(' full-w nonzero: %d, FDR-thresh nonzero: %d\n', ... - sum(si_full.dat ~= 0), sum(si_thr.dat ~= 0)); - - % --- 8. permutation_test (small nperm) --- - pm_for_perm = predictive_model('algorithm','svm','task','classification', ... - 'modeloptions', {'KernelFunction','linear'}, ... - 'random_state', 0); - pm_for_perm.cv = cv_splitter.stratified_group_kfold(3); - pm_for_perm = permutation_test(pm_for_perm, X, Y, ... - 'nperm', 10, 'groups', id, ... - 'verbose', false); - assert(numel(pm_for_perm.permutation_results.null_scores) == 10, 'nperm=10 null scores'); - assert(pm_for_perm.permutation_results.p_value > 0, 'p_value > 0'); - assert(pm_for_perm.permutation_results.p_value <= 1, 'p_value <= 1'); - % DPSP is a paired design (Hot+Warm per subject) — auto should pick within_subjects. - assert(strcmp(pm_for_perm.permutation_results.permutation, 'within_subjects'), ... - 'auto should pick within_subjects on paired DPSP'); - fprintf(' 8. permutation_test (nperm=10, %s) observed=%.3f, p=%.3f\n', ... - pm_for_perm.permutation_results.permutation, ... - pm_for_perm.permutation_results.observed, ... - pm_for_perm.permutation_results.p_value); - - % --- 9. calibrate + predict_proba --- - pm_cal = calibrate(pm, X, Y); - assert(isfield(pm_cal.fitted_values, 'calibrator'), 'calibrator present'); - P = predict_proba(pm_cal, X); - assert(all(P >= 0 & P <= 1), 'predict_proba in [0,1]'); - fprintf(' 9. calibrate + predict_proba: P mean=%.3f\n', mean(P)); - - % --- 10. select_features --- - pm_fs = predictive_model('algorithm','svm','task','classification'); - pm_fs = select_features(pm_fs, X, Y, 'k', 500, 'verbose', false); - assert(pm_fs.diagnostics.feature_selection.n_selected == 500, 'k=500 selected'); - fprintf(' 10. select_features OK (kept 500)\n'); - - % --- 11. visualisation methods (smoke test; no display assertions) --- - ROC = rocplot(pm, 'noplot'); - assert(isfield(ROC, 'AUC') && ROC.AUC >= 0 && ROC.AUC <= 1, 'rocplot AUC in [0,1]'); - h = plot(pm); close(h); % classification dispatcher (violin + ROC) - confusionchart(pm); close(gcf); - fprintf(' 11. rocplot/plot/confusionchart OK (AUC=%.3f)\n', ROC.AUC); - - % --- 12. montage / surface delegates (smoke test) --- - montage(pm, hw_obj); close all force; - fprintf(' 12. montage(pm, source) OK\n'); - - % --- 13. grid_search (tiny grid) --- - pm_gs = predictive_model('algorithm','svm','task','classification', ... - 'cv', cv_splitter.stratified_group_kfold(3)); - pm_gs = grid_search(pm_gs, X, Y, struct('BoxConstraint', [0.1 1]), ... - 'groups', id, 'verbose', false); - assert(isfield(pm_gs.diagnostics.grid_search, 'best_score'), 'grid_search best_score'); - fprintf(' 13. grid_search OK (best %s=%.3f)\n', ... - pm_gs.scorer.name, pm_gs.diagnostics.grid_search.best_score); - - % --- 14. stability_selection (small nboot) --- - pm_ss = stability_selection( ... - predictive_model('algorithm','linear_svm','task','classification'), ... - X, Y, 'nboot', 20, 'k', 1000, 'threshold', 0.6, 'groups', id, 'verbose', false); - ss = pm_ss.diagnostics.stability_selection; - assert(numel(ss.selection_freq) == size(X, 2), 'selection_freq length'); - assert(all(ss.selection_freq >= 0 & ss.selection_freq <= 1), 'selection_freq in [0,1]'); - fprintf(' 14. stability_selection OK (%d stable voxels)\n', ss.n_stable); - - % --- 15. @pipeline: PCA -> svm, leakage-free crossval + weight back-projection --- - est = predictive_model('algorithm','svm','task','classification'); - pipe = pipeline({ {'pca','k',20} }, est); - pipe = crossval(pipe, X, Y, 'groups', id, 'cv', cv_splitter.stratified_group_kfold(5)); - assert(strcmp(pipe.fit_type, 'crossval'), 'pipeline crossval fit_type'); - assert(numel(pipe.weights.w) == size(X, 2), 'pipeline back-projected weight length'); - si_pipe = weight_map_object(pipe, hw_obj); - assert(isa(si_pipe, 'statistic_image'), 'pipeline weight_map_object class'); - fprintf(' 15. pipeline (PCA->svm) OK (cv bal_acc=%.3f)\n', ... - pipe.error_metrics.balanced_accuracy.value); - - % --- 16. fmri_data.predict 'newapi' routing --- - [cverr, st, ~, pm_api] = predict(hw_obj, 'algorithm_name', 'cv_svm', ... - 'nfolds', 5, 'newapi', 'verbose', 0); - assert(isa(pm_api, 'predictive_model'), 'newapi returns predictive_model'); - assert(strcmp(pm_api.fit_type, 'crossval'), 'newapi fit_type=crossval'); - assert(isfield(st, 'dist_from_hyperplane_xval'), 'newapi stats has xval distances'); - fprintf(' 16. fmri_data.predict ''newapi'' OK (cverr=%.3f)\n', cverr); - - % --- 17. predicted_r2 / out_of_sample_r2 + report_accuracy + summary --- - % Regression model on a continuous outcome (signed label scaled to vary). - rng(7); - Yr = double(Y) + 0.5 * randn(size(Y)); - pm_reg = predictive_model('algorithm', 'pcr', 'task', 'regression'); - pm_reg = crossval(pm_reg, X, Yr, 'cv', cv_splitter.group_kfold(5), 'groups', id); - assert(isfield(pm_reg.error_metrics, 'predicted_r2'), 'predicted_r2 present'); - assert(isfield(pm_reg.error_metrics, 'out_of_sample_r2'), 'out_of_sample_r2 present'); - % predicted_r2 must equal 1 - PRESS/SST computed by hand. - yf = pm_reg.fitted_values.yfit(:); yo = Yr(:); - pr2_manual = 1 - sum((yo - yf).^2) / sum((yo - mean(yo)).^2); - assert(abs(pm_reg.error_metrics.predicted_r2.value - pr2_manual) < 1e-9, ... - 'predicted_r2 matches 1 - PRESS/SST'); - acc_r = report_accuracy(pm_reg, 'noverbose'); - assert(strcmp(acc_r.task, 'regression') && ~isnan(acc_r.predicted_r2), 'report_accuracy regression'); - s_reg = summary(pm_reg, 'noverbose'); - assert(strcmp(s_reg.provenance.fit_type, 'crossval'), 'summary provenance'); - % Classification report_accuracy: ROC-derived fields populated. - acc_c = report_accuracy(pm, 'noverbose'); - assert(strcmp(acc_c.task, 'classification') && ~isnan(acc_c.auc) && ~isnan(acc_c.npv), ... - 'report_accuracy classification (auc + npv)'); - fprintf(' 17. predicted_r2 / report_accuracy / summary OK (pred R^2=%.3f)\n', ... - pm_reg.error_metrics.predicted_r2.value); - - fprintf('\npredictive_model_unit_test: PASS\n'); -end diff --git a/CanlabCore/Unit_tests/xval_SVM_unit_test.m b/CanlabCore/Unit_tests/xval_SVM_unit_test.m deleted file mode 100644 index 535e6be4..00000000 --- a/CanlabCore/Unit_tests/xval_SVM_unit_test.m +++ /dev/null @@ -1,104 +0,0 @@ -function xval_SVM_unit_test() -% xval_SVM_unit_test Smoke test of xval_SVM on the DPSP Hot vs. Warm task. -% -% Reproduces the binary-classification setup from -% docs/markdown_tutorials/multivariate_classification_with_SVM/ -% multivariate_decoding_part1_classification_with_SVM.mlx -% and asserts that the wrapper returns a populated @predictive_model -% object whose categorised sub-structs agree with the legacy flat aliases. -% -% Data: Sample_datasets/DPSP_pain_rejection_participant_maps -% Single-subject Hot and Warm condition maps. -% -% Run: cd to this directory, then `xval_SVM_unit_test`. - - fprintf('=== xval_SVM_unit_test ===\n'); - - % --- Load DPSP Hot and Warm --- - canlabcore_dir = fileparts(fileparts(which('fmri_data'))); - sample_dir = fullfile(canlabcore_dir, 'Sample_datasets', ... - 'DPSP_pain_rejection_participant_maps'); - H = load(fullfile(sample_dir, 'DPSP_single_subject_images_hot.mat')); - W = load(fullfile(sample_dir, 'DPSP_single_subject_images_warm.mat')); - hot = H.single_subject_images_hot; - warm = W.single_subject_images_warm; - - n_hot = size(hot.dat, 2); - n_warm = size(warm.dat, 2); - - % Stack into one (subjects+subjects) x voxels matrix - rng(0); - p = size(hot.dat, 1); - keep_vox = randsample(p, min(5000, p)); - X = double([hot.dat(keep_vox, :) warm.dat(keep_vox, :)])'; - Y = [ones(n_hot, 1); -ones(n_warm, 1)]; - id = [(1:n_hot)'; (1:n_warm)']; % each subject contributes both conditions - - % --- Run xval_SVM (fast: no optimize / no repeats / no bootstrap) --- - pmodel_obj = xval_SVM(X, Y, id, 'nooptimize', 'norepeats', 'nobootstrap', ... - 'noverbose', 'noplot'); - - % --- Class invariants --- - assert(isa(pmodel_obj, 'predictive_model'), ... - 'Expected output of class predictive_model, got %s', class(pmodel_obj)); - assert(pmodel_obj.is_fitted, 'Returned object reports is_fitted = false'); - assert(pmodel_obj.is_classifier, 'Y is binary, expected is_classifier true'); - fprintf(' class = %s, is_fitted = %d, is_classifier = %d\n', ... - class(pmodel_obj), pmodel_obj.is_fitted, pmodel_obj.is_classifier); - - % --- Canonical-path access (legacy flat aliases removed) --- - assert(~isempty(pmodel_obj.Y), 'Y empty'); - assert(~isempty(pmodel_obj.id), 'id empty'); - assert(~isempty(pmodel_obj.fitted_values.yfit), 'yfit empty'); - assert(~isempty(pmodel_obj.fitted_values.dist_from_hyperplane_xval), 'dist_from_hyperplane_xval empty'); - % class_probability_xval is xval_SVM-legacy (per-fold fitPosterior). - % The new pipeline stores raw scores in fitted_values.scores; to get - % calibrated probabilities call calibrate(pm, X, Y) + predict_proba. - assert(~isempty(pmodel_obj.weights.w), 'weights.w empty'); - assert(~isempty(pmodel_obj.error_metrics.crossval_accuracy.value), 'crossval_accuracy.value empty'); - assert(~isempty(pmodel_obj.error_metrics.d_singleinterval.value), 'd_singleinterval.value empty'); - assert(~isempty(pmodel_obj.ml_model), 'ml_model empty'); - assert(~isempty(pmodel_obj.cv_partition.trIdx), 'trIdx empty'); - assert(~isempty(pmodel_obj.cv_partition.teIdx), 'teIdx empty'); - assert(~isempty(pmodel_obj.cv_partition.nfolds), 'nfolds empty'); - fprintf(' Canonical-path access OK\n'); - - % --- Shape / sanity --- - n = numel(Y); - assert(numel(pmodel_obj.fitted_values.yfit) == n, 'yfit length mismatch'); - assert(numel(pmodel_obj.fitted_values.dist_from_hyperplane_xval) == n, 'dist_from_hyperplane_xval length mismatch'); - assert(size(pmodel_obj.weights.w, 1) == size(X, 2), 'weights.w length should match number of features'); - cv_acc = pmodel_obj.error_metrics.crossval_accuracy.value; - assert(cv_acc >= 0 && cv_acc <= 100, 'crossval_accuracy out of [0,100]'); - % Within-person scoring fields populated by crossval (DPSP is paired): - assert(~isnan(pmodel_obj.error_metrics.crossval_accuracy_within.value), 'crossval_accuracy_within populated'); - assert(~isnan(pmodel_obj.error_metrics.d_within.value), 'd_within populated'); - assert(pmodel_obj.diagnostics.mult_obs_within_person == true, 'Should detect multiple obs per id'); - fprintf(' Shape/sanity OK: cv_acc = %.1f%%, d_single = %.2f, cv_acc_within = %.1f%%, d_within = %.2f\n', ... - cv_acc, ... - pmodel_obj.error_metrics.d_singleinterval.value, ... - pmodel_obj.error_metrics.crossval_accuracy_within.value, ... - pmodel_obj.error_metrics.d_within.value); - - % --- fit_type + omitted markers (Phase B) --- - assert(strcmp(pmodel_obj.fit_type, 'crossval'), ... - 'fit_type should be ''crossval'', got ''%s''', pmodel_obj.fit_type); - assert(islogical(pmodel_obj.omitted_cases), 'omitted_cases must be logical'); - assert(islogical(pmodel_obj.omitted_features), 'omitted_features must be logical'); - assert(numel(pmodel_obj.omitted_cases) == numel(Y), 'omitted_cases length should match original Y'); - assert(numel(pmodel_obj.omitted_features) == size(X, 2), 'omitted_features length should match original feature count'); - fprintf(' fit_type=%s, omitted_cases=%d, omitted_features=%d\n', ... - pmodel_obj.fit_type, sum(pmodel_obj.omitted_cases), sum(pmodel_obj.omitted_features)); - - % --- validate_object accepts the returned object --- - pmodel_obj.validate_object('noverbose'); - fprintf(' validate_object OK\n'); - - % --- Clone preserves hyperparameters, clears fitted state --- - pmodel_obj2 = clone(pmodel_obj); - assert(~pmodel_obj2.is_fitted, 'clone did not clear fitted state'); - assert(isequal(pmodel_obj2.modeloptions, pmodel_obj.modeloptions), 'clone lost modeloptions'); - fprintf(' clone() OK\n'); - - fprintf('xval_SVM_unit_test: PASS\n'); -end diff --git a/CanlabCore/Unit_tests/xval_SVR_unit_test.m b/CanlabCore/Unit_tests/xval_SVR_unit_test.m deleted file mode 100644 index 7e7c6c81..00000000 --- a/CanlabCore/Unit_tests/xval_SVR_unit_test.m +++ /dev/null @@ -1,62 +0,0 @@ -function xval_SVR_unit_test() -% xval_SVR_unit_test Smoke test of xval_SVR on the DPSP Hot-Warm contrast. -% -% Predicts a synthetic continuous outcome from Hot-Warm contrast maps -% (one per subject) via cross-validated linear SVR. Asserts that the -% wrapper returns a @predictive_model object with consistent -% categorised <-> legacy alias values. -% -% Run: cd to this directory, then `xval_SVR_unit_test`. - - fprintf('=== xval_SVR_unit_test ===\n'); - - canlabcore_dir = fileparts(fileparts(which('fmri_data'))); - sample_dir = fullfile(canlabcore_dir, 'Sample_datasets', ... - 'DPSP_pain_rejection_participant_maps'); - H = load(fullfile(sample_dir, 'DPSP_single_subject_images_hot.mat')); - W = load(fullfile(sample_dir, 'DPSP_single_subject_images_warm.mat')); - hot_vs_warm = image_math(H.single_subject_images_hot, ... - W.single_subject_images_warm, 'minus'); - - rng(0); - [p, n] = size(hot_vs_warm.dat); - keep_vox = randsample(p, min(3000, p)); - X = double(hot_vs_warm.dat(keep_vox, :))'; - % Synthetic predictable continuous outcome. - b_true = zeros(size(X, 2), 1); - b_true(1:30) = randn(30, 1); - Y = X * b_true + 0.5 * std(X * b_true) * randn(n, 1); - id = (1:n)'; - - pmodel_obj = xval_SVR(X, Y, id, 'nooptimize', 'norepeats', ... - 'nobootstrap', 'noverbose', 'noplot'); - - assert(isa(pmodel_obj, 'predictive_model'), ... - 'Expected predictive_model, got %s', class(pmodel_obj)); - assert(pmodel_obj.is_fitted, 'is_fitted false'); - assert(pmodel_obj.is_regressor, 'Y is continuous, expected is_regressor true'); - fprintf(' class = %s, is_fitted = %d, is_regressor = %d\n', ... - class(pmodel_obj), pmodel_obj.is_fitted, pmodel_obj.is_regressor); - - % Canonical-path access (legacy flat aliases removed). - assert(~isempty(pmodel_obj.Y), 'Y empty'); - assert(~isempty(pmodel_obj.fitted_values.yfit), 'yfit empty'); - assert(~isempty(pmodel_obj.weights.w), 'weights.w empty'); - assert(~isempty(pmodel_obj.error_metrics.prediction_outcome_r.value), 'prediction_outcome_r empty'); - assert(~isempty(pmodel_obj.ml_model), 'ml_model empty'); - fprintf(' Canonical-path access OK\n'); - - fprintf(' cv: r = %.3f, d = %.2f\n', ... - pmodel_obj.error_metrics.prediction_outcome_r.value, ... - pmodel_obj.error_metrics.d_singleinterval.value); - - % --- fit_type + omitted markers (Phase B) --- - assert(strcmp(pmodel_obj.fit_type, 'crossval'), 'fit_type should be crossval'); - assert(islogical(pmodel_obj.omitted_cases), 'omitted_cases must be logical'); - assert(islogical(pmodel_obj.omitted_features), 'omitted_features must be logical'); - fprintf(' fit_type=%s, omitted_cases=%d, omitted_features=%d\n', ... - pmodel_obj.fit_type, sum(pmodel_obj.omitted_cases), sum(pmodel_obj.omitted_features)); - - pmodel_obj.validate_object('noverbose'); - fprintf('xval_SVR_unit_test: PASS\n'); -end diff --git a/CanlabCore/Unit_tests/xval_discriminant_classifier_unit_test.m b/CanlabCore/Unit_tests/xval_discriminant_classifier_unit_test.m deleted file mode 100644 index 2a8810e2..00000000 --- a/CanlabCore/Unit_tests/xval_discriminant_classifier_unit_test.m +++ /dev/null @@ -1,62 +0,0 @@ -function xval_discriminant_classifier_unit_test() -% xval_discriminant_classifier_unit_test Smoke test on DPSP Hot vs Warm. -% -% Drives cross-validated LDA via fitcdiscr on the DPSP single-subject -% Hot vs Warm maps. Asserts @predictive_model return and -% categorised <-> legacy alias agreement. - - fprintf('=== xval_discriminant_classifier_unit_test ===\n'); - - canlabcore_dir = fileparts(fileparts(which('fmri_data'))); - sample_dir = fullfile(canlabcore_dir, 'Sample_datasets', ... - 'DPSP_pain_rejection_participant_maps'); - H = load(fullfile(sample_dir, 'DPSP_single_subject_images_hot.mat')); - W = load(fullfile(sample_dir, 'DPSP_single_subject_images_warm.mat')); - hot = H.single_subject_images_hot; - warm = W.single_subject_images_warm; - - % Stratify: 100 random voxels (LDA needs n > p) - rng(0); - p = size(hot.dat, 1); - keep_vox = randsample(p, 100); - X = double([hot.dat(keep_vox, :) warm.dat(keep_vox, :)])'; - labels = int32([ones(size(hot.dat, 2), 1); ... - 2*ones(size(warm.dat, 2), 1)]); - - pmodel_obj = xval_discriminant_classifier(X, labels, 'nFolds', 5, ... - 'verbose', false, 'doplot', false); - - assert(isa(pmodel_obj, 'predictive_model'), ... - 'Expected predictive_model, got %s', class(pmodel_obj)); - assert(pmodel_obj.is_fitted, 'is_fitted false'); - fprintf(' class = %s, is_fitted = %d\n', class(pmodel_obj), pmodel_obj.is_fitted); - - % Categorised <-> legacy alias agreement - % Canonical-path access (legacy flat aliases removed). - assert(~isempty(pmodel_obj.Y), 'Y empty'); - assert(~isempty(pmodel_obj.fitted_values.yfit), 'yfit empty'); - assert(~isempty(pmodel_obj.fitted_values.predictions), 'predictions empty'); - assert(~isempty(pmodel_obj.fitted_values.Y_per_fold), 'Y_per_fold empty'); - assert(~isempty(pmodel_obj.error_metrics.accuracy.value), 'accuracy.value empty'); - assert(~isempty(pmodel_obj.error_metrics.overallAccuracy.value), 'overallAccuracy.value empty'); - assert(~isempty(pmodel_obj.cv_partition.trIdx), 'trIdx empty'); - assert(~isempty(pmodel_obj.cv_partition.teIdx), 'teIdx empty'); - assert(~isempty(pmodel_obj.cv_partition.nfolds), 'nfolds empty'); - assert(numel(pmodel_obj.fold_models) == pmodel_obj.cv_partition.nfolds, ... - 'fold_models length should equal nfolds'); - fprintf(' Canonical-path access OK\n'); - - fprintf(' overall_acc = %.1f%%, mean_fold_acc = %.1f%%\n', ... - pmodel_obj.error_metrics.overallAccuracy.value, ... - mean(pmodel_obj.error_metrics.accuracy.value)); - - % --- fit_type + omitted markers (Phase B) --- - assert(strcmp(pmodel_obj.fit_type, 'crossval'), 'fit_type should be crossval'); - assert(islogical(pmodel_obj.omitted_cases), 'omitted_cases must be logical'); - assert(islogical(pmodel_obj.omitted_features), 'omitted_features must be logical'); - fprintf(' fit_type=%s, omitted_cases=%d, omitted_features=%d\n', ... - pmodel_obj.fit_type, sum(pmodel_obj.omitted_cases), sum(pmodel_obj.omitted_features)); - - pmodel_obj.validate_object('noverbose'); - fprintf('xval_discriminant_classifier_unit_test: PASS\n'); -end diff --git a/CanlabCore/Unit_tests/xval_regression_multisubject_featureselect_unit_test.m b/CanlabCore/Unit_tests/xval_regression_multisubject_featureselect_unit_test.m deleted file mode 100644 index b0f45baf..00000000 --- a/CanlabCore/Unit_tests/xval_regression_multisubject_featureselect_unit_test.m +++ /dev/null @@ -1,59 +0,0 @@ -function xval_regression_multisubject_featureselect_unit_test() -% xval_regression_multisubject_featureselect_unit_test DPSP smoke test. -% -% Verifies that xval_regression_multisubject_featureselect returns a -% populated @predictive_model object after the per-fold vox_weights -% padding fix. Predicts a synthetic continuous outcome from -% Hot - Warm DPSP contrast maps. - - fprintf('=== xval_regression_multisubject_featureselect_unit_test ===\n'); - - canlabcore_dir = fileparts(fileparts(which('fmri_data'))); - sd = fullfile(canlabcore_dir, 'Sample_datasets', 'DPSP_pain_rejection_participant_maps'); - H = load(fullfile(sd, 'DPSP_single_subject_images_hot.mat')); - W = load(fullfile(sd, 'DPSP_single_subject_images_warm.mat')); - hvw = image_math(H.single_subject_images_hot, W.single_subject_images_warm, 'minus'); - - rng(0); - [p, n] = size(hvw.dat); - keep_vox = randsample(p, min(2000, p)); - X = double(hvw.dat(keep_vox, :))'; - b_true = zeros(size(X, 2), 1); - b_true(1:30) = randn(30, 1); - Y = X * b_true + 0.5 * std(X * b_true) * randn(n, 1); - - pm = xval_regression_multisubject_featureselect('ols', {Y}, {X}, ... - 'holdout_method', 'balanced4', 'noverbose'); - - assert(isa(pm, 'predictive_model'), 'class mismatch: %s', class(pm)); - assert(pm.is_fitted, 'is_fitted=false'); - fprintf(' class = %s, is_fitted = %d\n', class(pm), pm.is_fitted); - - % Canonical-path access (legacy flat aliases removed). - % mean_vox_weights is only populated when pcsquash is enabled; this no-PCA test doesn't exercise that path. - assert(~isempty(pm.fitted_values.subjfit), 'subjfit empty'); - assert(~isempty(pm.weights.w_perfold), 'weights.w_perfold empty (was vox_weights)'); - assert(~isempty(pm.error_metrics.pred_err.value), 'pred_err empty'); - assert(~isempty(pm.error_metrics.pred_err_null.value), 'pred_err_null empty'); - assert(~isempty(pm.error_metrics.var_reduction.value), 'var_reduction empty'); - assert(~isempty(pm.Y), 'Y empty'); - fprintf(' Canonical-path access OK\n'); - - % Verify weights.w_perfold matches the full feature space (the bug we fixed - % was that this was the *selected*-features count instead). - assert(size(pm.weights.w_perfold, 1) == size(X, 2), ... - 'weights.w_perfold should span full feature count (%d), got %d', ... - size(X, 2), size(pm.weights.w_perfold, 1)); - fprintf(' w_perfold shape OK: %d x %d (features x folds)\n', ... - size(pm.weights.w_perfold, 1), size(pm.weights.w_perfold, 2)); - - r2 = pm.error_metrics.r_squared.value; - assert(~isempty(r2) && r2(1) > 0, 'expected r_squared > 0; got %.3f', r2(1)); - fprintf(' r_squared = %.3f, pred_err = %.3f, pred_err_null = %.3f\n', ... - r2(1), ... - pm.error_metrics.pred_err.value(1), ... - pm.error_metrics.pred_err_null.value(1)); - - pm.validate_object('noverbose'); - fprintf('xval_regression_multisubject_featureselect_unit_test: PASS\n'); -end diff --git a/CanlabCore/Unit_tests/xval_regression_multisubject_unit_test.m b/CanlabCore/Unit_tests/xval_regression_multisubject_unit_test.m deleted file mode 100644 index 0696e886..00000000 --- a/CanlabCore/Unit_tests/xval_regression_multisubject_unit_test.m +++ /dev/null @@ -1,75 +0,0 @@ -function xval_regression_multisubject_unit_test() -% xval_regression_multisubject_unit_test Smoke test using DPSP sample data. -% -% Verifies that xval_regression_multisubject returns a populated -% @predictive_model object whose categorised sub-structs and legacy -% flat aliases agree, on a real (small) brain dataset. -% -% Data: Sample_datasets/DPSP_pain_rejection_participant_maps -% Single-subject Hot - Warm contrast maps, regressed against a -% synthetic continuous outcome to exercise the regression code path. -% -% Run: cd to this directory, then `xval_regression_multisubject_unit_test`. - - fprintf('=== xval_regression_multisubject_unit_test ===\n'); - - % --- Load DPSP Hot - Warm contrast across subjects --- - canlabcore_dir = fileparts(fileparts(which('fmri_data'))); - sample_dir = fullfile(canlabcore_dir, 'Sample_datasets', ... - 'DPSP_pain_rejection_participant_maps'); - H = load(fullfile(sample_dir, 'DPSP_single_subject_images_hot.mat')); - W = load(fullfile(sample_dir, 'DPSP_single_subject_images_warm.mat')); - hot_obj = H.single_subject_images_hot; - warm_obj = W.single_subject_images_warm; - hot_vs_warm = image_math(hot_obj, warm_obj, 'minus'); - - % Keep the test light: ~5000 random voxels, all subjects. - rng(0); - [p, n] = size(hot_vs_warm.dat); - keep_vox = randsample(p, min(5000, p)); - X = double(hot_vs_warm.dat(keep_vox, :))'; % n x v - - % Synthetic continuous outcome predictable from X with sparse weights. - b_true = zeros(size(X, 2), 1); - b_true(1:50) = randn(50, 1); - Y = X * b_true + 0.5 * std(X * b_true) * randn(n, 1); - - % --- Run xval_regression_multisubject with OLS + PCA on a single dataset --- - pmodel_obj = xval_regression_multisubject('ols', {Y}, {X}, ... - 'pca', 'ndims', 10, 'holdout_method', 'balanced4', 'noverbose'); - - % --- Class invariants --- - assert(isa(pmodel_obj, 'predictive_model'), ... - 'Expected output of class predictive_model, got %s', class(pmodel_obj)); - assert(pmodel_obj.is_fitted, 'Returned object reports is_fitted = false'); - fprintf(' class = %s, is_fitted = %d\n', class(pmodel_obj), pmodel_obj.is_fitted); - - % --- Canonical-path access (legacy flat aliases removed) --- - assert(~isempty(pmodel_obj.fitted_values.subjfit), 'subjfit empty'); - assert(~isempty(pmodel_obj.weights.subjbetas), 'subjbetas empty'); - assert(~isempty(pmodel_obj.weights.mean_vox_weights), 'mean_vox_weights empty'); - assert(~isempty(pmodel_obj.error_metrics.pred_err.value), 'pred_err empty'); - assert(~isempty(pmodel_obj.error_metrics.pred_err_null.value), 'pred_err_null empty'); - assert(~isempty(pmodel_obj.error_metrics.var_reduction.value), 'var_reduction empty'); - assert(~isempty(pmodel_obj.error_metrics.r_each_subject.value), 'r_each_subject empty'); - assert(~isempty(pmodel_obj.Y), 'Y empty'); - assert(~isempty(pmodel_obj.inputParameters), 'inputParameters empty'); - fprintf(' Canonical-path access OK\n'); - - % --- Shape / sanity --- - assert(numel(pmodel_obj.fitted_values.subjfit) == 1, 'subjfit cell length should equal num datasets'); - assert(numel(pmodel_obj.fitted_values.subjfit{1}) == n, 'subjfit{1} length should equal n'); - assert(size(pmodel_obj.weights.mean_vox_weights, 1) == size(X, 2), 'mean_vox_weights should have one row per voxel'); - r2 = pmodel_obj.error_metrics.r_squared.value; - assert(~isempty(r2) && r2(1) > 0, 'OLS+PCA should explain non-trivial variance'); - fprintf(' Shape/sanity OK: r_squared = %.3f, pred_err = %.3f, pred_err_null = %.3f\n', ... - r2(1), ... - pmodel_obj.error_metrics.pred_err.value(1), ... - pmodel_obj.error_metrics.pred_err_null.value(1)); - - % --- validate_object accepts the returned object --- - pmodel_obj.validate_object('noverbose'); - fprintf(' validate_object OK\n'); - - fprintf('xval_regression_multisubject_unit_test: PASS\n'); -end