From 2d3832c46f69aac28df70619c7fcab361603bd6d Mon Sep 17 00:00:00 2001 From: BalazsErdos Date: Tue, 30 Dec 2025 21:09:11 +0100 Subject: [PATCH 1/9] added example notebook for replicability analysis --- docs/references.bib | 11 +- examples/plot_replicability.py | 262 +++++++++++++++++++++++++++++++++ 2 files changed, 272 insertions(+), 1 deletion(-) create mode 100644 examples/plot_replicability.py diff --git a/docs/references.bib b/docs/references.bib index 8e44ff4..e87a6d2 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -303,4 +303,13 @@ @INPROCEEDINGS{chatzis2023timeaware volume={}, number={}, pages={1-6}, - doi={10.1109/MLSP55844.2023.10285943}} \ No newline at end of file + doi={10.1109/MLSP55844.2023.10285943}} + +@article {erdos2025extracting, + author = {Erd{\H o}s, Bal{\'a}zs and Chatzis, Christos and Thorsen, Jonathan and Stokholm, Jakob and Smilde, Age K. and Rasmussen, Morten A. and Acar, Evrim}, + title = {Extracting host-specific developmental signatures from longitudinal microbiome data}, + elocation-id = {2025.11.22.689760}, + year = {2025}, + doi = {10.1101/2025.11.22.689760}, + journal = {bioRxiv} +} \ No newline at end of file diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py new file mode 100644 index 0000000..af7f23e --- /dev/null +++ b/examples/plot_replicability.py @@ -0,0 +1,262 @@ +""" +.. _replicability: + +Determine the number of components through replicability analysis +---------------- + +This example shows how to select the number of components for PARAFAC2 models by checking if patterns are replicable :cite:p:`erdos2025extracting`. The process involves fitting the model to different subsets of your data to see if the results stay consistent (i.e. replicable across data subsets). To maximize explanatory power, typically, we select the highest number of components that remains replicable across the data subsets. +""" + +############################################################################### +# Imports and utilities +# ^^^^^^^^^^^^^^^^^^^^^ + +import matplotlib.pyplot as plt +import numpy as np +import tensorly as tl +from tensorly.metrics import congruence_coefficient +from matcouply.decomposition import parafac2_aoadmm +from matcouply.coupled_matrices import CoupledMatrixFactorization +import sklearn +from sklearn.model_selection import RepeatedKFold + +import tlviz + +rng = np.random.default_rng(1) + +############################################################################### +# To fit PARAFAC2 models, we need to solve a non-convex optimization problem, possibly with local minima. It is +# therefore useful to fit several models with the same number of components using many different random +# initialisations. + + +def fit_many_parafac2(X, num_components, num_inits=5): + + best_err = np.inf + decomposition = None + for i in range(num_inits): + trial_decomposition, trial_errs = parafac2_aoadmm( + matrices=X, + rank=num_components, + return_errors=True, + non_negative=[True, True, True], + n_iter_max=500, + absolute_tol=1e-4, + feasibility_tol=1e-4, + inner_tol=1e-4, + inner_n_iter_max=5, + feasibility_penalty_scale=5, + tol=1e-5, + random_state=i, + verbose=0, + ) + + if best_err > trial_errs.rec_errors[-1]: + best_err = trial_errs.rec_errors[-1] + decomposition = trial_decomposition # note, with real data, convergence should be checked + + (est_weights, (est_C, est_B, est_A)) = decomposition + est_B = np.asarray(est_B) + + # normalize factors + As = np.empty(est_A.shape) + Bs = np.empty(est_B.shape) + Cs = np.empty(est_C.shape) + + K = Cs.shape[0] + for r in range(num_components): + norm_Ar = tl.norm(est_A[:, r]) + norm_Cr = tl.norm(est_C[:, r]) + As[:, r] = est_A[:, r] / norm_Ar + Cs[:, r] = est_C[:, r] / norm_Cr + + for k in range(K): + norm_Brk = tl.norm(est_B[k][:, r]) + Bs[k,:,r] = est_B[k,:,r] / norm_Brk + + # calculate scaled B + cB = np.empty(Bs.shape) + for r in range(num_components): + for k in range(Cs.shape[0]): + cB[k,:,r] = Cs[k,r] * Bs[k,:,r] + + + return (est_weights, (As, cB)) + + +############################################################################### +# Creat simulated data +# ^^^^^^^^^^^^^^^^^^^^^^^ +# +# Simulate noisy data with 2 components. + + +def truncated_normal(size): + x = rng.standard_normal(size=size) + x[x < 0] = 0 + return tl.tensor(x) + +I, J, K = 35, 20, 25 +rank = 2 + +A = rng.uniform(size=(I, rank)) + 0.1 # Add 0.1 to ensure that there is signal for all components for all slices +A = tl.tensor(A) + +B_blueprint = truncated_normal(size=(J, rank)) +B_is = [np.roll(B_blueprint, i, axis=0) for i in range(I)] +B_is = [tl.tensor(B_i) for B_i in B_is] + +C = rng.uniform(size=(K, rank)) +C = tl.tensor(C) + +dataset = CoupledMatrixFactorization((None, (A, B_is, C))) + +dataset = dataset.to_tensor() +eta = 0.3 # noise level +noise = np.random.normal(0, 1, dataset.shape) +dataset = dataset + tl.norm(dataset) * eta * noise / tl.norm(noise) +dataset = dataset / tl.norm(dataset) + +############################################################################### +# The replicability analysis boils down to the following steps: +# +# 1. Split the data in a (user-chosen) mode into :math:`N` folds (user-chosen). +# 2. Create :math:`N` subsets by subtracting each fold from the complete dataset. +# 3. Fit multiple initializations to each subset and choose the *best* run +# according to lowest loss (total of :math:`N` *best* runs). +# 4. Compare, in terms of FMS, the best runs across the different subsets +# to evaluate the replicability of the uncovered patterns (:math:`\binom{N}{2}` comparisons). +# 5. Repeat the above process :math:`M` times (user-chosen), to find a total of +# :math:`M \binom{N}{2}` comparisons. + + +############################################################################### +# Split the data and fit PARAFAC2 on each data subset +# ^^^^^^^^^^^^^^^^^^ + +splits = 5 # N +repeats = 5 # M + +models = {} +split_indices = {} # Keeps track of which indices are used in each subset + +for rank in [1, 2, 3, 4]: + + print(f"{rank} components") + + rskf = RepeatedKFold(n_splits=splits, n_repeats=repeats, random_state=1) + + models[rank] = [[] for _ in range(repeats)] + split_indices[rank] = [[] for _ in range(repeats)] + + for split_no, (train_index, _) in enumerate(rskf.split(dataset)): + repeat_no = split_no // splits + + # Append indices to the current repeat + split_indices[rank][repeat_no].append(train_index) + + train = dataset[train_index] + train = train / tl.norm(train) + + current_model = fit_many_parafac2(train, rank) + + # Append model to the current repeat + models[rank][repeat_no].append(current_model) + + +############################################################################### +# Often, the mode we will be splitting within refers to different samples, +# depending on the use-case, it might be reasonable to retain the +# distributions of some properties in each subset. For this goal, +# `RepeatedStratifiedKFold `_ +# can be used. +# +# If pre-processing is used, it is important to apply it to +# each subset in isolation to avoid leaking information from the omitted part of the data. +# For example, in this case we normalize each subset to unit norm independently. +# Note, that ``for split_no, (train_index, _) in enumerate(rskf.split(dataset)):`` may be run in parallel +# for efficiency. + +############################################################################### +# Compute and assess replicability I. +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Since we are subsetting the data on ``mode=0``, and the ``mode=1`` factors of PARAFAC2 are specific +# to the corresponding level in ``mode=0``, only the shared factor matrix can be compared using FMS: + +replicability_stability = {} +for rank in models.keys(): + replicability_stability[rank] = [] + for repeat, current_models in enumerate(models[rank]): + for i, m_i in enumerate(current_models): + for j, m_j in enumerate(current_models): + if i >= j: # include every pair only once and omit i == j + continue + fms = congruence_coefficient(m_i[1][0], m_j[1][0])[0] + replicability_stability[rank].append(fms) + +ranks = sorted(replicability_stability.keys()) +data = [np.ravel(replicability_stability[r]) for r in ranks] + +fig, ax = plt.subplots() +ax.boxplot(data,whis=(0.95,0.05), positions=ranks) +ax.set_xlabel("Number of components") +ax.set_ylabel("FMS_A") +plt.show() + +############################################################################### +# Here, we observe that over-estimating the number of components +# results in a loss of replicable of the patterns, indicated by low FMS. + +############################################################################### +# Compute and assess replicability II. +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# There is an alternative way to estimate the replicability of the uncovered patterns, +# including the factors corresponding to ``mode=0``, and ``mode=1`` :cite:p:`erdos2025extracting`. +# By using only the indices present in both subsets (e.g. the factors +# corresponding to the subjects' data included in both factorizations) + +replicability_alt = {} +for rank in models.keys(): + replicability_alt[rank] = [] + for repeat in range(repeats): + for i, cp_i in enumerate(models[rank][repeat]): + for j, cp_j in enumerate(models[rank][repeat]): + if i >= j: # include every pair only once and omit i == j + continue + weights_i, (A_i, cB_i) = cp_i + weights_j, (A_j, cB_j) = cp_j + + indices_subset_i = sorted(split_indices[rank][repeat][i]) + indices_subset_j = sorted(split_indices[rank][repeat][j]) + + common_indices = sorted(list(set(indices_subset_i).intersection(set(indices_subset_j)))) + indices2use_i = [] + indices2use_j = [] + + for common_idx in common_indices: + indices2use_i.append(indices_subset_i.index(common_idx)) + indices2use_j.append(indices_subset_j.index(common_idx)) + + cB_i = cB_i[indices2use_i, :, :] + cB_j = cB_j[indices2use_j, :, :] + fms = tlviz.factor_tools.factor_match_score( + (weights_i, (A_i, cB_i.reshape(-1, cB_i.shape[2]))), (weights_j, (A_j, cB_j.reshape(-1, cB_j.shape[2]))), consider_weights=False + ) + replicability_alt[rank].append(fms) + +ranks = sorted(replicability_alt.keys()) +data = [np.ravel(replicability_alt[r]) for r in ranks] + +fig, ax = plt.subplots() +ax.boxplot(data, positions=ranks) +ax.set_xlabel("Number of components") +ax.set_ylabel("FMS_CB") +plt.show() + +############################################################################### +# ``common_indices`` contains the indices (e.g. subjects/samples) present in both subsets, +# but since the position of each index can change (e.g. sample no 3 is not guaranteeed at +# the third position in all subsets as the first and second samples might be omitted) we need to +# utilize the indices in the original tensor input. +# +# Similar results can be observed with this approach in terms of the replicability of the patterns. \ No newline at end of file From b3ad5c294918fb9df4ee0ff3ea3df2f6c3522569 Mon Sep 17 00:00:00 2001 From: BalazsErdos Date: Wed, 31 Dec 2025 16:13:25 +0100 Subject: [PATCH 2/9] added explanation in comment --- examples/plot_replicability.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py index af7f23e..aee18b8 100644 --- a/examples/plot_replicability.py +++ b/examples/plot_replicability.py @@ -74,7 +74,8 @@ def fit_many_parafac2(X, num_components, num_inits=5): norm_Brk = tl.norm(est_B[k][:, r]) Bs[k,:,r] = est_B[k,:,r] / norm_Brk - # calculate scaled B + # calculate scaled B; + # since the loadings in B are specific to levels of C, we absorb C into the corresponding B for each component cB = np.empty(Bs.shape) for r in range(num_components): for k in range(Cs.shape[0]): From d7172a528a3b0c56d926e31cc73fd39a83fa68d9 Mon Sep 17 00:00:00 2001 From: berdos <22996725+blzserdos@users.noreply.github.com> Date: Sat, 10 Jan 2026 18:21:24 +0100 Subject: [PATCH 3/9] Update examples/plot_replicability.py fix notation Co-authored-by: Marie Roald --- examples/plot_replicability.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py index aee18b8..6fce895 100644 --- a/examples/plot_replicability.py +++ b/examples/plot_replicability.py @@ -55,7 +55,7 @@ def fit_many_parafac2(X, num_components, num_inits=5): best_err = trial_errs.rec_errors[-1] decomposition = trial_decomposition # note, with real data, convergence should be checked - (est_weights, (est_C, est_B, est_A)) = decomposition + (est_weights, (est_A, est_B, est_C)) = decomposition est_B = np.asarray(est_B) # normalize factors From f50cf44f61b042e732dc0c094a2b5023bd9ddf1d Mon Sep 17 00:00:00 2001 From: berdos <22996725+blzserdos@users.noreply.github.com> Date: Sat, 10 Jan 2026 18:22:05 +0100 Subject: [PATCH 4/9] Update examples/plot_replicability.py fix notation Co-authored-by: Marie Roald --- examples/plot_replicability.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py index 6fce895..47f8e9f 100644 --- a/examples/plot_replicability.py +++ b/examples/plot_replicability.py @@ -220,12 +220,12 @@ def truncated_normal(size): for rank in models.keys(): replicability_alt[rank] = [] for repeat in range(repeats): - for i, cp_i in enumerate(models[rank][repeat]): - for j, cp_j in enumerate(models[rank][repeat]): + for i, model_i in enumerate(models[rank][repeat]): + for j, model_j in enumerate(models[rank][repeat]): if i >= j: # include every pair only once and omit i == j continue - weights_i, (A_i, cB_i) = cp_i - weights_j, (A_j, cB_j) = cp_j + weights_i, (aB_i, C_i) = model_i + weights_j, (aB_j, C_j) = model_j indices_subset_i = sorted(split_indices[rank][repeat][i]) indices_subset_j = sorted(split_indices[rank][repeat][j]) From 491c9df37655c1e6f42b347034e87f94f1928e5b Mon Sep 17 00:00:00 2001 From: berdos <22996725+blzserdos@users.noreply.github.com> Date: Sat, 10 Jan 2026 18:23:08 +0100 Subject: [PATCH 5/9] Update examples/plot_replicability.py improve readability Co-authored-by: Marie Roald --- examples/plot_replicability.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py index 47f8e9f..dd7adc6 100644 --- a/examples/plot_replicability.py +++ b/examples/plot_replicability.py @@ -188,8 +188,8 @@ def truncated_normal(size): for rank in models.keys(): replicability_stability[rank] = [] for repeat, current_models in enumerate(models[rank]): - for i, m_i in enumerate(current_models): - for j, m_j in enumerate(current_models): + for i, model_i in enumerate(current_models): + for j, model_j in enumerate(current_models): if i >= j: # include every pair only once and omit i == j continue fms = congruence_coefficient(m_i[1][0], m_j[1][0])[0] From f246b9af2f55802fe579f5d62bf06be94470e784 Mon Sep 17 00:00:00 2001 From: BalazsErdos Date: Sat, 10 Jan 2026 19:50:56 +0100 Subject: [PATCH 6/9] address PR comments --- .../plot_replicability.codeobj.json | 1348 +++++++++++++++++ docs/auto_examples/plot_replicability.ipynb | 168 ++ docs/auto_examples/plot_replicability.py | 265 ++++ docs/auto_examples/plot_replicability.py.md5 | 1 + docs/auto_examples/plot_replicability.rst | 414 +++++ docs/auto_examples/plot_replicability.zip | Bin 0 -> 23609 bytes 6 files changed, 2196 insertions(+) create mode 100644 docs/auto_examples/plot_replicability.codeobj.json create mode 100644 docs/auto_examples/plot_replicability.ipynb create mode 100644 docs/auto_examples/plot_replicability.py create mode 100644 docs/auto_examples/plot_replicability.py.md5 create mode 100644 docs/auto_examples/plot_replicability.rst create mode 100644 docs/auto_examples/plot_replicability.zip diff --git a/docs/auto_examples/plot_replicability.codeobj.json b/docs/auto_examples/plot_replicability.codeobj.json new file mode 100644 index 0000000..7079e4a --- /dev/null +++ b/docs/auto_examples/plot_replicability.codeobj.json @@ -0,0 +1,1348 @@ +{ + "A": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "B_blueprint": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "B_is": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "C": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "C_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "C_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "CoupledMatrixFactorization": [ + { + "is_class": true, + "is_explicit": false, + "module": "matcouply.coupled_matrices", + "module_short": "matcouply.coupled_matrices", + "name": "CoupledMatrixFactorization" + }, + { + "is_class": true, + "is_explicit": false, + "module": "matcouply", + "module_short": "matcouply", + "name": "CoupledMatrixFactorization" + }, + { + "is_class": true, + "is_explicit": false, + "module": "tensorly._factorized_tensor", + "module_short": "tensorly._factorized_tensor", + "name": "FactorizedTensor" + }, + { + "is_class": true, + "is_explicit": false, + "module": "tensorly", + "module_short": "tensorly", + "name": "FactorizedTensor" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections.abc", + "module_short": "collections.abc", + "name": "Mapping" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections", + "module_short": "collections", + "name": "Mapping" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections.abc", + "module_short": "collections.abc", + "name": "Collection" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections", + "module_short": "collections", + "name": "Collection" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections.abc", + "module_short": "collections.abc", + "name": "Sized" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections", + "module_short": "collections", + "name": "Sized" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections.abc", + "module_short": "collections.abc", + "name": "Iterable" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections", + "module_short": "collections", + "name": "Iterable" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections.abc", + "module_short": "collections.abc", + "name": "Container" + }, + { + "is_class": true, + "is_explicit": false, + "module": "collections", + "module_short": "collections", + "name": "Container" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matcouply.coupled_matrices", + "module_short": "matcouply.coupled_matrices", + "name": "CoupledMatrixFactorization" + } + ], + "I": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "J": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "K": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "RepeatedKFold": [ + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "RepeatedKFold" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection._split", + "name": "_UnsupportedGroupCVMixin" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "_UnsupportedGroupCVMixin" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "_UnsupportedGroupCVMixin" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection._split", + "name": "_RepeatedSplits" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "_RepeatedSplits" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "_RepeatedSplits" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.utils._metadata_requests", + "module_short": "sklearn.utils._metadata_requests", + "name": "_MetadataRequester" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn.utils", + "module_short": "sklearn.utils", + "name": "_MetadataRequester" + }, + { + "is_class": true, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "_MetadataRequester" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold" + } + ], + "_": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "aB_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "aB_i.reshape": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "aB_i.shape": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "tuple" + } + ], + "aB_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "aB_j.reshape": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "aB_j.shape": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "tuple" + } + ], + "ax": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes" + } + ], + "ax.boxplot": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.boxplot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.boxplot" + } + ], + "ax.set_xlabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.set_xlabel" + } + ], + "ax.set_ylabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.set_ylabel" + } + ], + "common_idx": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "int64" + } + ], + "common_indices": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "congruence_coefficient": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "tensorly.metrics", + "module_short": "tensorly.metrics", + "name": "congruence_coefficient" + } + ], + "current_model": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "tuple" + } + ], + "current_models": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "data": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "dataset": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "dataset.shape": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "tuple" + } + ], + "dataset.to_tensor": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "eta": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "float" + } + ], + "fig": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "Figure" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Figure" + } + ], + "fit_many_parafac2": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + } + ], + "fms": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "float64" + } + ], + "i": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "indices2use_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "indices2use_i.append": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "indices2use_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "indices2use_j.append": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "indices_subset_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "indices_subset_i.index": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "indices_subset_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "indices_subset_j.index": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "j": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "model_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "tuple" + } + ], + "model_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "tuple" + } + ], + "models": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "dict" + } + ], + "models.keys": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "noise": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "np.asarray": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "asarray" + } + ], + "np.empty": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "empty" + } + ], + "np.inf": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "float" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "inf" + } + ], + "np.random.default_rng": [ + { + "is_class": false, + "is_explicit": false, + "module": "_cython_3_2_1", + "module_short": "_cython_3_2_1", + "name": "cython_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "default_rng" + } + ], + "np.random.normal": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "module.normal" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "normal" + } + ], + "np.ravel": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ravel" + } + ], + "np.roll": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "roll" + } + ], + "parafac2_aoadmm": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matcouply.decomposition", + "module_short": "matcouply.decomposition", + "name": "parafac2_aoadmm" + } + ], + "plt.show": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "show" + } + ], + "plt.subplots": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "subplots" + } + ], + "rank": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "ranks": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "repeat": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "repeat_no": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "repeats": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "replicability_alt": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "dict" + } + ], + "replicability_alt.keys": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "replicability_stability": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "dict" + } + ], + "replicability_stability.keys": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + } + ], + "rng": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random._generator", + "module_short": "numpy.random", + "name": "Generator" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "Generator" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "Generator" + } + ], + "rng.standard_normal": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random._generator", + "module_short": "numpy.random", + "name": "Generator.standard_normal" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "Generator.standard_normal" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "Generator.standard_normal" + } + ], + "rng.uniform": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random._generator", + "module_short": "numpy.random", + "name": "Generator.uniform" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "Generator.uniform" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "Generator.uniform" + } + ], + "rskf": [ + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "RepeatedKFold" + } + ], + "rskf.split": [ + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "RepeatedKFold.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "RepeatedKFold.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection._split", + "name": "_UnsupportedGroupCVMixin.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "_UnsupportedGroupCVMixin.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "_UnsupportedGroupCVMixin.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection._split", + "module_short": "sklearn.model_selection._split", + "name": "_RepeatedSplits.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.model_selection", + "module_short": "sklearn.model_selection", + "name": "_RepeatedSplits.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "_RepeatedSplits.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.utils._metadata_requests", + "module_short": "sklearn.utils._metadata_requests", + "name": "_MetadataRequester.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.utils", + "module_short": "sklearn.utils", + "name": "_MetadataRequester.split" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn", + "module_short": "sklearn", + "name": "_MetadataRequester.split" + } + ], + "split_indices": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "dict" + } + ], + "split_no": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "splits": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "tl.norm": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "tensorly", + "module_short": "tensorly", + "name": "norm" + } + ], + "tl.tensor": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "tensorly", + "module_short": "tensorly", + "name": "tensor" + } + ], + "tlviz.factor_tools.factor_match_score": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "tlviz.factor_tools", + "module_short": "tlviz.factor_tools", + "name": "factor_match_score" + } + ], + "train": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "train_index": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "truncated_normal": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + } + ], + "weights_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "NoneType" + } + ], + "weights_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "NoneType" + } + ] +} \ No newline at end of file diff --git a/docs/auto_examples/plot_replicability.ipynb b/docs/auto_examples/plot_replicability.ipynb new file mode 100644 index 0000000..70bbe4d --- /dev/null +++ b/docs/auto_examples/plot_replicability.ipynb @@ -0,0 +1,168 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n\n# Determine the number of components through replicability analysis\n\nThis example shows how to select the number of components for PARAFAC2 models by checking if patterns are replicable :cite:p:`erdos2025extracting`. The process involves fitting the model to different subsets of your data to see if the results stay consistent (i.e. replicable across data subsets). To maximize explanatory power, typically, we select the highest number of components that remains replicable across the data subsets.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports and utilities\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\nimport numpy as np\nimport tensorly as tl\nfrom tensorly.metrics import congruence_coefficient\nfrom matcouply.decomposition import parafac2_aoadmm\nfrom matcouply.coupled_matrices import CoupledMatrixFactorization\nimport sklearn\nfrom sklearn.model_selection import RepeatedKFold\n\nimport tlviz\n\nrng = np.random.default_rng(1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To fit PARAFAC2 models, we need to solve a non-convex optimization problem, possibly with local minima. It is\ntherefore useful to fit several models with the same number of components using many different random\ninitialisations.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def fit_many_parafac2(X, num_components, num_inits=5):\n \n best_err = np.inf\n decomposition = None\n for i in range(num_inits):\n trial_decomposition, trial_errs = parafac2_aoadmm(\n matrices=X,\n rank=num_components,\n return_errors=True,\n non_negative=[True, True, True],\n n_iter_max=500,\n absolute_tol=1e-4,\n feasibility_tol=1e-4,\n inner_tol=1e-4,\n inner_n_iter_max=5,\n feasibility_penalty_scale=5,\n tol=1e-5,\n random_state=i,\n verbose=0,\n )\n \n if best_err > trial_errs.rec_errors[-1]:\n best_err = trial_errs.rec_errors[-1]\n decomposition = trial_decomposition # note, with real data, convergence should be checked\n\n (est_weights, (est_A, est_B, est_C)) = decomposition\n est_B = np.asarray(est_B)\n\n # normalize factors\n As = np.empty(est_A.shape)\n Bs = np.empty(est_B.shape)\n Cs = np.empty(est_C.shape)\n\n K = As.shape[0]\n for r in range(num_components):\n norm_Ar = tl.norm(est_A[:, r])\n norm_Cr = tl.norm(est_C[:, r])\n As[:, r] = est_A[:, r] / norm_Ar\n Cs[:, r] = est_C[:, r] / norm_Cr\n\n for k in range(K):\n norm_Brk = tl.norm(est_B[k][:, r])\n Bs[k,:,r] = est_B[k,:,r] / norm_Brk\n\n # calculate scaled B; \n # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component\n aB = np.empty(Bs.shape)\n for r in range(num_components):\n for k in range(As.shape[0]):\n aB[k,:,r] = As[k,r] * Bs[k,:,r] \n\n\n return (est_weights, (Cs, aB))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Creat simulated data\n\nSimulate noisy data with 2 components.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def truncated_normal(size):\n x = rng.standard_normal(size=size)\n x[x < 0] = 0\n return tl.tensor(x)\n\nI, J, K = 25, 20, 35\nrank = 2\n\nA = rng.uniform(size=(I, rank)) + 0.1 # Add 0.1 to ensure that there is signal for all components for all slices\nA = tl.tensor(A)\n\nB_blueprint = truncated_normal(size=(J, rank))\nB_is = [np.roll(B_blueprint, i, axis=0) for i in range(I)]\nB_is = [tl.tensor(B_i) for B_i in B_is]\n\nC = rng.uniform(size=(K, rank))\nC = tl.tensor(C)\n\ndataset = CoupledMatrixFactorization((None, (A, B_is, C)))\n\ndataset = dataset.to_tensor()\neta = 0.3 # noise level\nnoise = np.random.normal(0, 1, dataset.shape)\ndataset = dataset + tl.norm(dataset) * eta * noise / tl.norm(noise)\ndataset = dataset / tl.norm(dataset)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The replicability analysis boils down to the following steps:\n\n1. Split the data in a (user-chosen) mode into $N$ folds (user-chosen).\n2. Create $N$ subsets by subtracting each fold from the complete dataset.\n3. Fit multiple initializations to each subset and choose the *best* run\n according to lowest loss (total of $N$ *best* runs).\n4. Compare, in terms of FMS, the best runs across the different subsets\n to evaluate the replicability of the uncovered patterns ($\\binom{N}{2}$ comparisons).\n5. Repeat the above process $M$ times (user-chosen), to find a total of\n $M \\binom{N}{2}$ comparisons.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Split the data and fit PARAFAC2 on each data subset\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "splits = 5 # N\nrepeats = 5 # M\n\nmodels = {}\nsplit_indices = {} # Keeps track of which indices are used in each subset\n\nfor rank in [1, 2, 3, 4]:\n\n print(f\"{rank} components\")\n\n rskf = RepeatedKFold(n_splits=splits, n_repeats=repeats, random_state=1)\n\n models[rank] = [[] for _ in range(repeats)]\n split_indices[rank] = [[] for _ in range(repeats)]\n\n for split_no, (train_index, _) in enumerate(rskf.split(dataset)):\n repeat_no = split_no // splits\n\n # Append indices to the current repeat\n split_indices[rank][repeat_no].append(train_index)\n \n train = dataset[train_index]\n train = train / tl.norm(train)\n\n current_model = fit_many_parafac2(train, rank)\n \n # Append model to the current repeat\n models[rank][repeat_no].append(current_model)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Often, the mode we will be splitting within refers to different samples,\ndepending on the use-case, it might be reasonable to retain the\ndistributions of some properties in each subset. For this goal,\n[RepeatedStratifiedKFold](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RepeatedStratifiedKFold.html#sklearn.model_selection.RepeatedStratifiedKFold)\ncan be used.\n\nIf pre-processing is used, it is important to apply it to\neach subset in isolation to avoid leaking information from the omitted part of the data.\nFor example, in this case we normalize each subset to unit norm independently.\nNote, that ``for split_no, (train_index, _) in enumerate(rskf.split(dataset)):`` may be run in parallel\nfor efficiency.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compute and assess replicability I.\nSince we are subsetting the data on ``mode=0``, and the ``mode=1`` factors of PARAFAC2 are specific \nto the corresponding level in ``mode=0``, only the shared factor matrix can be compared using FMS: \n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "replicability_stability = {}\nfor rank in models.keys():\n replicability_stability[rank] = []\n for repeat, current_models in enumerate(models[rank]):\n for i, model_i in enumerate(current_models):\n for j, model_j in enumerate(current_models):\n if i >= j: # include every pair only once and omit i == j\n continue\n weights_i, (C_i, _) = model_i\n weights_j, (C_j, _) = model_j\n fms = congruence_coefficient(C_i, C_j)[0]\n replicability_stability[rank].append(fms)\n\nranks = sorted(replicability_stability.keys())\ndata = [np.ravel(replicability_stability[r]) for r in ranks]\n\nfig, ax = plt.subplots()\nax.boxplot(data,whis=(0.95,0.05), positions=ranks)\nax.set_xlabel(\"Number of components\")\nax.set_ylabel(\"FMS_C\")\nplt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, we observe that over-estimating the number of components\nresults in a loss of replicable of the patterns, indicated by low FMS.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compute and assess replicability II.\nThere is an alternative way to estimate the replicability of the uncovered patterns, \nincluding the factors corresponding to ``mode=0``, and ``mode=1`` :cite:p:`erdos2025extracting`. \nBy using only the indices present in both subsets (e.g. the factors \ncorresponding to the subjects' data included in both factorizations)\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "replicability_alt = {}\nfor rank in models.keys():\n replicability_alt[rank] = []\n for repeat in range(repeats):\n for i, model_i in enumerate(models[rank][repeat]):\n for j, model_j in enumerate(models[rank][repeat]):\n if i >= j: # include every pair only once and omit i == j\n continue\n weights_i, (C_i, aB_i) = model_i\n weights_j, (C_j, aB_j) = model_j\n\n indices_subset_i = sorted(split_indices[rank][repeat][i])\n indices_subset_j = sorted(split_indices[rank][repeat][j])\n \n common_indices = sorted(list(set(indices_subset_i).intersection(set(indices_subset_j))))\n indices2use_i = []\n indices2use_j = []\n\n for common_idx in common_indices:\n indices2use_i.append(indices_subset_i.index(common_idx))\n indices2use_j.append(indices_subset_j.index(common_idx))\n\n aB_i = aB_i[indices2use_i, :, :]\n aB_j = aB_j[indices2use_j, :, :]\n fms = tlviz.factor_tools.factor_match_score(\n (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False\n )\n replicability_alt[rank].append(fms)\n\nranks = sorted(replicability_alt.keys())\ndata = [np.ravel(replicability_alt[r]) for r in ranks]\n\nfig, ax = plt.subplots()\nax.boxplot(data, positions=ranks)\nax.set_xlabel(\"Number of components\")\nax.set_ylabel(\"FMS_aB\")\nplt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "``common_indices`` contains the indices (e.g. subjects/samples) present in both subsets,\nbut since the position of each index can change (e.g. sample no 3 is not guaranteeed at\nthe third position in all subsets as the first and second samples might be omitted) we need to\nutilize the indices in the original tensor input.\n\nSimilar results can be observed with this approach in terms of the replicability of the patterns.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/docs/auto_examples/plot_replicability.py b/docs/auto_examples/plot_replicability.py new file mode 100644 index 0000000..4fff486 --- /dev/null +++ b/docs/auto_examples/plot_replicability.py @@ -0,0 +1,265 @@ +""" +.. _replicability: + +Determine the number of components through replicability analysis +---------------- + +This example shows how to select the number of components for PARAFAC2 models by checking if patterns are replicable :cite:p:`erdos2025extracting`. The process involves fitting the model to different subsets of your data to see if the results stay consistent (i.e. replicable across data subsets). To maximize explanatory power, typically, we select the highest number of components that remains replicable across the data subsets. +""" + +############################################################################### +# Imports and utilities +# ^^^^^^^^^^^^^^^^^^^^^ + +import matplotlib.pyplot as plt +import numpy as np +import tensorly as tl +from tensorly.metrics import congruence_coefficient +from matcouply.decomposition import parafac2_aoadmm +from matcouply.coupled_matrices import CoupledMatrixFactorization +import sklearn +from sklearn.model_selection import RepeatedKFold + +import tlviz + +rng = np.random.default_rng(1) + +############################################################################### +# To fit PARAFAC2 models, we need to solve a non-convex optimization problem, possibly with local minima. It is +# therefore useful to fit several models with the same number of components using many different random +# initialisations. + + +def fit_many_parafac2(X, num_components, num_inits=5): + + best_err = np.inf + decomposition = None + for i in range(num_inits): + trial_decomposition, trial_errs = parafac2_aoadmm( + matrices=X, + rank=num_components, + return_errors=True, + non_negative=[True, True, True], + n_iter_max=500, + absolute_tol=1e-4, + feasibility_tol=1e-4, + inner_tol=1e-4, + inner_n_iter_max=5, + feasibility_penalty_scale=5, + tol=1e-5, + random_state=i, + verbose=0, + ) + + if best_err > trial_errs.rec_errors[-1]: + best_err = trial_errs.rec_errors[-1] + decomposition = trial_decomposition # note, with real data, convergence should be checked + + (est_weights, (est_A, est_B, est_C)) = decomposition + est_B = np.asarray(est_B) + + # normalize factors + As = np.empty(est_A.shape) + Bs = np.empty(est_B.shape) + Cs = np.empty(est_C.shape) + + K = As.shape[0] + for r in range(num_components): + norm_Ar = tl.norm(est_A[:, r]) + norm_Cr = tl.norm(est_C[:, r]) + As[:, r] = est_A[:, r] / norm_Ar + Cs[:, r] = est_C[:, r] / norm_Cr + + for k in range(K): + norm_Brk = tl.norm(est_B[k][:, r]) + Bs[k,:,r] = est_B[k,:,r] / norm_Brk + + # calculate scaled B; + # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component + aB = np.empty(Bs.shape) + for r in range(num_components): + for k in range(As.shape[0]): + aB[k,:,r] = As[k,r] * Bs[k,:,r] + + + return (est_weights, (Cs, aB)) + + +############################################################################### +# Creat simulated data +# ^^^^^^^^^^^^^^^^^^^^^^^ +# +# Simulate noisy data with 2 components. + + +def truncated_normal(size): + x = rng.standard_normal(size=size) + x[x < 0] = 0 + return tl.tensor(x) + +I, J, K = 25, 20, 35 +rank = 2 + +A = rng.uniform(size=(I, rank)) + 0.1 # Add 0.1 to ensure that there is signal for all components for all slices +A = tl.tensor(A) + +B_blueprint = truncated_normal(size=(J, rank)) +B_is = [np.roll(B_blueprint, i, axis=0) for i in range(I)] +B_is = [tl.tensor(B_i) for B_i in B_is] + +C = rng.uniform(size=(K, rank)) +C = tl.tensor(C) + +dataset = CoupledMatrixFactorization((None, (A, B_is, C))) + +dataset = dataset.to_tensor() +eta = 0.3 # noise level +noise = np.random.normal(0, 1, dataset.shape) +dataset = dataset + tl.norm(dataset) * eta * noise / tl.norm(noise) +dataset = dataset / tl.norm(dataset) + +############################################################################### +# The replicability analysis boils down to the following steps: +# +# 1. Split the data in a (user-chosen) mode into :math:`N` folds (user-chosen). +# 2. Create :math:`N` subsets by subtracting each fold from the complete dataset. +# 3. Fit multiple initializations to each subset and choose the *best* run +# according to lowest loss (total of :math:`N` *best* runs). +# 4. Compare, in terms of FMS, the best runs across the different subsets +# to evaluate the replicability of the uncovered patterns (:math:`\binom{N}{2}` comparisons). +# 5. Repeat the above process :math:`M` times (user-chosen), to find a total of +# :math:`M \binom{N}{2}` comparisons. + + +############################################################################### +# Split the data and fit PARAFAC2 on each data subset +# ^^^^^^^^^^^^^^^^^^ + +splits = 5 # N +repeats = 5 # M + +models = {} +split_indices = {} # Keeps track of which indices are used in each subset + +for rank in [1, 2, 3, 4]: + + print(f"{rank} components") + + rskf = RepeatedKFold(n_splits=splits, n_repeats=repeats, random_state=1) + + models[rank] = [[] for _ in range(repeats)] + split_indices[rank] = [[] for _ in range(repeats)] + + for split_no, (train_index, _) in enumerate(rskf.split(dataset)): + repeat_no = split_no // splits + + # Append indices to the current repeat + split_indices[rank][repeat_no].append(train_index) + + train = dataset[train_index] + train = train / tl.norm(train) + + current_model = fit_many_parafac2(train, rank) + + # Append model to the current repeat + models[rank][repeat_no].append(current_model) + + +############################################################################### +# Often, the mode we will be splitting within refers to different samples, +# depending on the use-case, it might be reasonable to retain the +# distributions of some properties in each subset. For this goal, +# `RepeatedStratifiedKFold `_ +# can be used. +# +# If pre-processing is used, it is important to apply it to +# each subset in isolation to avoid leaking information from the omitted part of the data. +# For example, in this case we normalize each subset to unit norm independently. +# Note, that ``for split_no, (train_index, _) in enumerate(rskf.split(dataset)):`` may be run in parallel +# for efficiency. + +############################################################################### +# Compute and assess replicability I. +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Since we are subsetting the data on ``mode=0``, and the ``mode=1`` factors of PARAFAC2 are specific +# to the corresponding level in ``mode=0``, only the shared factor matrix can be compared using FMS: + +replicability_stability = {} +for rank in models.keys(): + replicability_stability[rank] = [] + for repeat, current_models in enumerate(models[rank]): + for i, model_i in enumerate(current_models): + for j, model_j in enumerate(current_models): + if i >= j: # include every pair only once and omit i == j + continue + weights_i, (C_i, _) = model_i + weights_j, (C_j, _) = model_j + fms = congruence_coefficient(C_i, C_j)[0] + replicability_stability[rank].append(fms) + +ranks = sorted(replicability_stability.keys()) +data = [np.ravel(replicability_stability[r]) for r in ranks] + +fig, ax = plt.subplots() +ax.boxplot(data,whis=(0.95,0.05), positions=ranks) +ax.set_xlabel("Number of components") +ax.set_ylabel("FMS_C") +plt.show() + +############################################################################### +# Here, we observe that over-estimating the number of components +# results in a loss of replicable of the patterns, indicated by low FMS. + +############################################################################### +# Compute and assess replicability II. +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# There is an alternative way to estimate the replicability of the uncovered patterns, +# including the factors corresponding to ``mode=0``, and ``mode=1`` :cite:p:`erdos2025extracting`. +# By using only the indices present in both subsets (e.g. the factors +# corresponding to the subjects' data included in both factorizations) + +replicability_alt = {} +for rank in models.keys(): + replicability_alt[rank] = [] + for repeat in range(repeats): + for i, model_i in enumerate(models[rank][repeat]): + for j, model_j in enumerate(models[rank][repeat]): + if i >= j: # include every pair only once and omit i == j + continue + weights_i, (C_i, aB_i) = model_i + weights_j, (C_j, aB_j) = model_j + + indices_subset_i = sorted(split_indices[rank][repeat][i]) + indices_subset_j = sorted(split_indices[rank][repeat][j]) + + common_indices = sorted(list(set(indices_subset_i).intersection(set(indices_subset_j)))) + indices2use_i = [] + indices2use_j = [] + + for common_idx in common_indices: + indices2use_i.append(indices_subset_i.index(common_idx)) + indices2use_j.append(indices_subset_j.index(common_idx)) + + aB_i = aB_i[indices2use_i, :, :] + aB_j = aB_j[indices2use_j, :, :] + fms = tlviz.factor_tools.factor_match_score( + (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False + ) + replicability_alt[rank].append(fms) + +ranks = sorted(replicability_alt.keys()) +data = [np.ravel(replicability_alt[r]) for r in ranks] + +fig, ax = plt.subplots() +ax.boxplot(data, positions=ranks) +ax.set_xlabel("Number of components") +ax.set_ylabel("FMS_aB") +plt.show() + +############################################################################### +# ``common_indices`` contains the indices (e.g. subjects/samples) present in both subsets, +# but since the position of each index can change (e.g. sample no 3 is not guaranteeed at +# the third position in all subsets as the first and second samples might be omitted) we need to +# utilize the indices in the original tensor input. +# +# Similar results can be observed with this approach in terms of the replicability of the patterns. \ No newline at end of file diff --git a/docs/auto_examples/plot_replicability.py.md5 b/docs/auto_examples/plot_replicability.py.md5 new file mode 100644 index 0000000..5b41514 --- /dev/null +++ b/docs/auto_examples/plot_replicability.py.md5 @@ -0,0 +1 @@ +f429dbdf7a4cf9ed22965979507994ba \ No newline at end of file diff --git a/docs/auto_examples/plot_replicability.rst b/docs/auto_examples/plot_replicability.rst new file mode 100644 index 0000000..e3dc081 --- /dev/null +++ b/docs/auto_examples/plot_replicability.rst @@ -0,0 +1,414 @@ + +.. DO NOT EDIT. +.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. +.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: +.. "auto_examples/plot_replicability.py" +.. LINE NUMBERS ARE GIVEN BELOW. + +.. only:: html + + .. note:: + :class: sphx-glr-download-link-note + + :ref:`Go to the end ` + to download the full example code. + +.. rst-class:: sphx-glr-example-title + +.. _sphx_glr_auto_examples_plot_replicability.py: + + +.. _replicability: + +Determine the number of components through replicability analysis +---------------- + +This example shows how to select the number of components for PARAFAC2 models by checking if patterns are replicable :cite:p:`erdos2025extracting`. The process involves fitting the model to different subsets of your data to see if the results stay consistent (i.e. replicable across data subsets). To maximize explanatory power, typically, we select the highest number of components that remains replicable across the data subsets. + +.. GENERATED FROM PYTHON SOURCE LINES 11-13 + +Imports and utilities +^^^^^^^^^^^^^^^^^^^^^ + +.. GENERATED FROM PYTHON SOURCE LINES 13-27 + +.. code-block:: Python + + + import matplotlib.pyplot as plt + import numpy as np + import tensorly as tl + from tensorly.metrics import congruence_coefficient + from matcouply.decomposition import parafac2_aoadmm + from matcouply.coupled_matrices import CoupledMatrixFactorization + import sklearn + from sklearn.model_selection import RepeatedKFold + + import tlviz + + rng = np.random.default_rng(1) + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 28-31 + +To fit PARAFAC2 models, we need to solve a non-convex optimization problem, possibly with local minima. It is +therefore useful to fit several models with the same number of components using many different random +initialisations. + +.. GENERATED FROM PYTHON SOURCE LINES 31-88 + +.. code-block:: Python + + + + def fit_many_parafac2(X, num_components, num_inits=5): + + best_err = np.inf + decomposition = None + for i in range(num_inits): + trial_decomposition, trial_errs = parafac2_aoadmm( + matrices=X, + rank=num_components, + return_errors=True, + non_negative=[True, True, True], + n_iter_max=500, + absolute_tol=1e-4, + feasibility_tol=1e-4, + inner_tol=1e-4, + inner_n_iter_max=5, + feasibility_penalty_scale=5, + tol=1e-5, + random_state=i, + verbose=0, + ) + + if best_err > trial_errs.rec_errors[-1]: + best_err = trial_errs.rec_errors[-1] + decomposition = trial_decomposition # note, with real data, convergence should be checked + + (est_weights, (est_A, est_B, est_C)) = decomposition + est_B = np.asarray(est_B) + + # normalize factors + As = np.empty(est_A.shape) + Bs = np.empty(est_B.shape) + Cs = np.empty(est_C.shape) + + K = As.shape[0] + for r in range(num_components): + norm_Ar = tl.norm(est_A[:, r]) + norm_Cr = tl.norm(est_C[:, r]) + As[:, r] = est_A[:, r] / norm_Ar + Cs[:, r] = est_C[:, r] / norm_Cr + + for k in range(K): + norm_Brk = tl.norm(est_B[k][:, r]) + Bs[k,:,r] = est_B[k,:,r] / norm_Brk + + # calculate scaled B; + # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component + aB = np.empty(Bs.shape) + for r in range(num_components): + for k in range(As.shape[0]): + aB[k,:,r] = As[k,r] * Bs[k,:,r] + + + return (est_weights, (Cs, aB)) + + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 89-93 + +Creat simulated data +^^^^^^^^^^^^^^^^^^^^^^^ + +Simulate noisy data with 2 components. + +.. GENERATED FROM PYTHON SOURCE LINES 93-121 + +.. code-block:: Python + + + + def truncated_normal(size): + x = rng.standard_normal(size=size) + x[x < 0] = 0 + return tl.tensor(x) + + I, J, K = 25, 20, 35 + rank = 2 + + A = rng.uniform(size=(I, rank)) + 0.1 # Add 0.1 to ensure that there is signal for all components for all slices + A = tl.tensor(A) + + B_blueprint = truncated_normal(size=(J, rank)) + B_is = [np.roll(B_blueprint, i, axis=0) for i in range(I)] + B_is = [tl.tensor(B_i) for B_i in B_is] + + C = rng.uniform(size=(K, rank)) + C = tl.tensor(C) + + dataset = CoupledMatrixFactorization((None, (A, B_is, C))) + + dataset = dataset.to_tensor() + eta = 0.3 # noise level + noise = np.random.normal(0, 1, dataset.shape) + dataset = dataset + tl.norm(dataset) * eta * noise / tl.norm(noise) + dataset = dataset / tl.norm(dataset) + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 122-132 + +The replicability analysis boils down to the following steps: + +1. Split the data in a (user-chosen) mode into :math:`N` folds (user-chosen). +2. Create :math:`N` subsets by subtracting each fold from the complete dataset. +3. Fit multiple initializations to each subset and choose the *best* run + according to lowest loss (total of :math:`N` *best* runs). +4. Compare, in terms of FMS, the best runs across the different subsets + to evaluate the replicability of the uncovered patterns (:math:`\binom{N}{2}` comparisons). +5. Repeat the above process :math:`M` times (user-chosen), to find a total of + :math:`M \binom{N}{2}` comparisons. + +.. GENERATED FROM PYTHON SOURCE LINES 135-137 + +Split the data and fit PARAFAC2 on each data subset +^^^^^^^^^^^^^^^^^^ + +.. GENERATED FROM PYTHON SOURCE LINES 137-168 + +.. code-block:: Python + + + splits = 5 # N + repeats = 5 # M + + models = {} + split_indices = {} # Keeps track of which indices are used in each subset + + for rank in [1, 2, 3, 4]: + + print(f"{rank} components") + + rskf = RepeatedKFold(n_splits=splits, n_repeats=repeats, random_state=1) + + models[rank] = [[] for _ in range(repeats)] + split_indices[rank] = [[] for _ in range(repeats)] + + for split_no, (train_index, _) in enumerate(rskf.split(dataset)): + repeat_no = split_no // splits + + # Append indices to the current repeat + split_indices[rank][repeat_no].append(train_index) + + train = dataset[train_index] + train = train / tl.norm(train) + + current_model = fit_many_parafac2(train, rank) + + # Append model to the current repeat + models[rank][repeat_no].append(current_model) + + + + + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + 1 components + 2 components + 3 components + 4 components + + + + +.. GENERATED FROM PYTHON SOURCE LINES 169-180 + +Often, the mode we will be splitting within refers to different samples, +depending on the use-case, it might be reasonable to retain the +distributions of some properties in each subset. For this goal, +`RepeatedStratifiedKFold `_ +can be used. + +If pre-processing is used, it is important to apply it to +each subset in isolation to avoid leaking information from the omitted part of the data. +For example, in this case we normalize each subset to unit norm independently. +Note, that ``for split_no, (train_index, _) in enumerate(rskf.split(dataset)):`` may be run in parallel +for efficiency. + +.. GENERATED FROM PYTHON SOURCE LINES 182-186 + +Compute and assess replicability I. +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Since we are subsetting the data on ``mode=0``, and the ``mode=1`` factors of PARAFAC2 are specific +to the corresponding level in ``mode=0``, only the shared factor matrix can be compared using FMS: + +.. GENERATED FROM PYTHON SOURCE LINES 186-209 + +.. code-block:: Python + + + replicability_stability = {} + for rank in models.keys(): + replicability_stability[rank] = [] + for repeat, current_models in enumerate(models[rank]): + for i, model_i in enumerate(current_models): + for j, model_j in enumerate(current_models): + if i >= j: # include every pair only once and omit i == j + continue + weights_i, (C_i, _) = model_i + weights_j, (C_j, _) = model_j + fms = congruence_coefficient(C_i, C_j)[0] + replicability_stability[rank].append(fms) + + ranks = sorted(replicability_stability.keys()) + data = [np.ravel(replicability_stability[r]) for r in ranks] + + fig, ax = plt.subplots() + ax.boxplot(data,whis=(0.95,0.05), positions=ranks) + ax.set_xlabel("Number of components") + ax.set_ylabel("FMS_C") + plt.show() + + + + +.. image-sg:: /auto_examples/images/sphx_glr_plot_replicability_001.png + :alt: plot replicability + :srcset: /auto_examples/images/sphx_glr_plot_replicability_001.png + :class: sphx-glr-single-img + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 210-212 + +Here, we observe that over-estimating the number of components +results in a loss of replicable of the patterns, indicated by low FMS. + +.. GENERATED FROM PYTHON SOURCE LINES 214-220 + +Compute and assess replicability II. +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +There is an alternative way to estimate the replicability of the uncovered patterns, +including the factors corresponding to ``mode=0``, and ``mode=1`` :cite:p:`erdos2025extracting`. +By using only the indices present in both subsets (e.g. the factors +corresponding to the subjects' data included in both factorizations) + +.. GENERATED FROM PYTHON SOURCE LINES 220-259 + +.. code-block:: Python + + + replicability_alt = {} + for rank in models.keys(): + replicability_alt[rank] = [] + for repeat in range(repeats): + for i, model_i in enumerate(models[rank][repeat]): + for j, model_j in enumerate(models[rank][repeat]): + if i >= j: # include every pair only once and omit i == j + continue + weights_i, (C_i, aB_i) = model_i + weights_j, (C_j, aB_j) = model_j + + indices_subset_i = sorted(split_indices[rank][repeat][i]) + indices_subset_j = sorted(split_indices[rank][repeat][j]) + + common_indices = sorted(list(set(indices_subset_i).intersection(set(indices_subset_j)))) + indices2use_i = [] + indices2use_j = [] + + for common_idx in common_indices: + indices2use_i.append(indices_subset_i.index(common_idx)) + indices2use_j.append(indices_subset_j.index(common_idx)) + + aB_i = aB_i[indices2use_i, :, :] + aB_j = aB_j[indices2use_j, :, :] + fms = tlviz.factor_tools.factor_match_score( + (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False + ) + replicability_alt[rank].append(fms) + + ranks = sorted(replicability_alt.keys()) + data = [np.ravel(replicability_alt[r]) for r in ranks] + + fig, ax = plt.subplots() + ax.boxplot(data, positions=ranks) + ax.set_xlabel("Number of components") + ax.set_ylabel("FMS_aB") + plt.show() + + + + +.. image-sg:: /auto_examples/images/sphx_glr_plot_replicability_002.png + :alt: plot replicability + :srcset: /auto_examples/images/sphx_glr_plot_replicability_002.png + :class: sphx-glr-single-img + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 260-265 + +``common_indices`` contains the indices (e.g. subjects/samples) present in both subsets, +but since the position of each index can change (e.g. sample no 3 is not guaranteeed at +the third position in all subsets as the first and second samples might be omitted) we need to +utilize the indices in the original tensor input. + +Similar results can be observed with this approach in terms of the replicability of the patterns. + + +.. rst-class:: sphx-glr-timing + + **Total running time of the script:** (8 minutes 4.271 seconds) + + +.. _sphx_glr_download_auto_examples_plot_replicability.py: + +.. only:: html + + .. container:: sphx-glr-footer sphx-glr-footer-example + + .. container:: sphx-glr-download sphx-glr-download-jupyter + + :download:`Download Jupyter notebook: plot_replicability.ipynb ` + + .. container:: sphx-glr-download sphx-glr-download-python + + :download:`Download Python source code: plot_replicability.py ` + + .. container:: sphx-glr-download sphx-glr-download-zip + + :download:`Download zipped: plot_replicability.zip ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ diff --git a/docs/auto_examples/plot_replicability.zip b/docs/auto_examples/plot_replicability.zip new file mode 100644 index 0000000000000000000000000000000000000000..efd37bc7a6057968fee99897b3f4a05aa20af9e6 GIT binary patch literal 23609 zcmeHP-)|d9a^Bn}LBQt=kjDmPAf%5P$@1#ql;Z-D4#qi#9cS%_14u13B%9Jm!x>^` zXp3tY$UWv?Nb-39N&c1~?|Bc9@2l>fALLLXWy!k-gp0Mvp04WZ>guZM>W_T;>X(1{ zt8e)FuYZ5qKltav|NP(fH{bC8zeU?DjI%+aXJIh%hd~%*3vae)Hk%F4Q`PpPM&oxn z)5$c5w8}2EissXyPE(( zvrwz_GX9t<{G+m1r8?B3Z2eG^I8kqpzdwF)e0rd!@mPnc8ZOl6Qje~K=t2b(HS;q_ z9Hq)nbS@IcIT{6-KAIh!>tr0K2m1$y`X)>KQ3kqmPrb+NW=TBKX{v(gLmYn4shR{C zAB#LP8f6>@lZj3+bCu4AsfGfejzv6A)Y#8_GYd^gsIWw*^AKIq%!in91chWIZ3Uk8 z%HsK>B*xIf!t!lHs#s0^n_wFJsG+b~2&HCmvQV@5qfR<1Tg6BdI<|kpc_9Lt9eF440MY2AL?F> zMj+jxkBkN!218ga-zh&;voLedVWP7I+oG9kg2vN038gg)8!aN$AMY%k`UkPnf5b1e)S>_ z$GOt7@I&xpqmjUtyHJdmKK6hsqgty)@Fy@At}2Xuu2 zqg9-=u3~9tTDF$+4T9a)FXOfm*&S1NL=+r)k0N356N!EAN z**u9j$2dv5@3GD58Nw0~x`l;gIvo%Bm`1b_bzv&+C@7MS7hiM$nGd;-S zu=`B!ep}ydqJ0>eM8-RM22q4DZ|7$EL?!X=9(AS>J+s9ig^la(VLgEzcFj{`ssjX` zOm~BN4%idJIMv-6b+uPkFmFI;%Wd_$l0A8e9@)k1?LIp*rtJ(Vr|7Py$VVVi#cW}! zHO8!V;50G~Z$gKVXzXtWt_}k&rglLu&irB?jv=wh5;UBmlc~iW_faEOW@$&sr{j)d zZEslRwpb;PSW0W=ctpMS&7lte6&)4TOr1d!sZ>#L+n2#`*^HQ z{=v9<^4d`dFtATXp%4iJq75@obs~vtI@6;7$p_2=S?Y&mq_P1VjDEmXj9agPG}x zH?rBl{H)P9tyRmbTn0Wbn4Dq~TyC8G(7rx>Z?%}LA+TWEauywhcv7LYT)ZqE*fnf7 zH0%gV>^+z07NqbF-2V-lnZ}bMd$*>9HGb9+`-4Ev`xo-LcfZbEw{2--`*EP{TNYCmhEgBRWR9z6oq1{z$ingY_&u zlD+xNQ}0mTFogt)%C+`Y3ul64cXWxACu+-4z_`#OWVn|{=Wou*ZJbv7dD!a*p79W< zBII4&;UDD(ymgfX<95kLnPaASXBIjT1+3D@@wum7KtA}JEMOJHohQth2Q?M!J34SP ziX4$4Hw2eSKV>%cRAK9Y6<+=b5kP!D=%5Zf#UGRtT3MW7n<8`+f)~VTo1%P+$>DT@ zatZbm^onwi1e+JH-{G_X`8gQJ!b#l}DcHDD1u>FP?hk%AXV7D}>N-LWHo_Oi$j@;6 z%ZoKFD_?&YMDg^~o6nyPKA($Q@sl7$-9$t^^prW@3LAe2E_u<$GI)KiaCXG$y0X}q z2+x5iG_XKpRAK~nWNPimJlj2d#96XW)<7x;Wt0`f0K|HY6+qzC%1EP;($eW759zqy zG@x<>v%KMTqhSxdUG?d+>5Y@(SV|CVBGoI6y#rtDkFIGfA1?#gp6kRy8g|Z@)=@Ib zMniH1#uhy2!3P~6+H};n$owThiHmPdnxE*&KbK-e6HZVo*nOXKY7=vW5cP24mK7QjiY$=v5w`X}O%TkD>!dJD z**^wIl5hl-RhE)6Q$-r4MlD&~eaWFpVLa_%zmH-0CaH%%5Z9UN-4QYmxN1ZXo?OWm z*(-dck1`qk95QH_B)gLX3eHI|M12)blzTgkQSm}eZKe~}Tg6_BNQ1Bh^M+`tx`_Rd zGMs;JON#HXN*N9%cBiV}US=60?%rNH3a*1}SBjKgoLuZ7m4z&OnEgD|=^hT^vK_{I zYl@fNnjyT)Y#Q!tN%Y2o^Fthrk)!$M0wQ8wZ64+cKs?1gf<(I{;}UQoxv?01 zscw4}V!|&~&lRTGm<&&f|HTq2Fe5VDdCP&@0B6Ol%g9qxltyNDgmZ936AXr?WC;h0IM0=5DCBs? zC3@F-k+y93SwmZdLF46xi)?=s)l9T8UflSNJkFUB40S2MksL=fd1?+^Iwz^JD-khb zknqY8UU3c^5|SA8$v}PARaZw0omlX24lm3L09>5;0j@8p^O(gKE-!Z{2)ZE1BM#ZP z7mA56k-4e5yGH)C)rkRS*gECEi2hwS*%elQQ&__x%1b>Ind@Gk79?lQ^@1|IEl zInPOU>&m*|5B7+Cv=0u0auVDV)1zW%74~H9#+(PyTTbFWj-=}ddMF9X3cBjzHP6_S z;DU#BFu-M08;9qvfbOjs1mzv+gJdq^-+M6sh*{ZTb|pTeK|hZr2K3Pm%jly-QQ2< z_ZsGeua8(40tKz8Vytp-@fiW7y&fH6zgQ9(IR?2C#7cAYyh*Hd%(Q&%kjMLQ; zrS%EkUDCt>Qgu;7uQ-3! zb>aNbWkp_Er;~YnLpye#jD89`c3ex!GlRR}xjZvW42S=g<2SyknZ+`*t2NB7>X@xe zhW3Ve@=wo7V?^vI+{hiRYKZn`2(BswNJHugO#UIT5zn|8+;YcZ>OLTjjV^IL4Wvaq zXrSftojg0Dw74RcOU~9V%0<#3AA1L9Sb8Kzjwnx(DEKnbmA%WfL`ag-V2o>cOWb|o zhv{mzwrW`}sG5uZ-0!)S*&s&(&8FFru>YdGi1?QWN`w59N6&-K&)t47`wIsrCengA zV{Ef#Yoj6DoJ6`k;Bin6Puo1J*F{ejBv8yK-JLqEFfo#Iy6i?dIgKt^_H#oRxFL#d z&zWifIG`@(Ky1M&5s2hx;7q|#`$$ljwz{UH;02cp-oTYfE)gWCi8EP%Kf|>vaYU56 zD9P9}e_M|62I+(W1m*@bKv@GAfggeTx(IluW^RHoFiyUvbb+hJAoOvbvtUA7+prnt z*jzhG%?4$QBxb3AOkE+L+$*-eoI`qVbG+z(|Hsq*zx<8=Ki@w4<{SS1clCJDV77>c zpWv68BOQim^GNl?NwJA9X&wNcsPU*d^^ zqkhzncGMT6P|4=~|HDrIoM5N@sA;Xq#4~Fr4{u^SYoMsKnJq*Zp=eZCxW3UNDY^ps zH;*zHbuoAE%XkjcyIhSKU<4}eR;r5YvPS!QpN!(A}1_Ul>N$*}#l?JK@NZ|So ztj{e&f!Fk0#RAuNybBDxhN=X&1sb@%6M=M9IB-PZ%7u2X(4lL5LiViLhQG*XIIIIT<7k$jH+0as7J4aM58a6IoJr1v6Z*q{{ia6EsW(8nWAQ!(S;y!%?CRN=z9y z!wa+L7P9`l&K%5sYVo|zPuvZ4HfveNI#*V=0O8&s67{2(HE?HS>%8{{-kCkVI^x+3 zSqOS8f|gLrkgQfoCD3y}LamfHv-dzf7qo>8;I;~TrgrWCKAR2&pHrEO#e_k;ZL!ZB zvY@u_%IN2QL_~58{J9@(g@3LY2H18?gb^x-;@ts&u47Vw+^qENEq3b{9sX8|cwD^b(K+qJp0)OsD4}^Xe33yk1 zFYxn|ei?oIgjZ9*`&is_Ke`R}{N(kMO##m$CP)Xe>Z+XqxnY(WIz)zq1QRT8%+6Za74_ce7B~fE+M75xT*XBc} zTNbh1UxjS#*ByK=#bh%kIhxZh}`b768Yq8KBz>&Nb7k!7&SB+Q81SG8~7m;+;H)I&GM$ z0Jv9snchWklh65zlXC^G2fnvsxQh{rcex&9!E=aY1;bwSE?Y&-^ef_gtr+}Jz-}ZK zSaaNq3T#M>0qBM;aNS_=Ww7pk)N}iH3Dtep`gx+dVT@%+H#gY55#7i#OK@s9!v_O& z_oEkZF%O65?nhq<&TWE23C&Fff4xlL^lqFtPtb(yJ}!0}*B*}#-QH5SKd&SP;i?WL zjF7;DrJQ0I7w%x^#r8N0fx>CZhisf5mD%<_&op!9r;z|S43(4o~0kL_npq2wk>8Do0S#eVjKgn-3yV^Lt{cTaJ=0F4pebxQm7P~f~ zrwl(~1S250K*eto#0U>AC=kc0`jGh;{HepgZ;qtDlhx)>`YNmHC(4L=-=l!(S4XO# z#Y+ZCDcq~v$H3|9^h#TEa@PXtR|n?+W#)m)25e)#Ml~S#ySZ%xsjr(5tq&UEp9)4E z36@~wS0$kJ8H##T{VK&au=-VM>Mw`L=XoCmUWUjARfv2}lxKeweVJ$@T>Yw4IGIzh z@~fHx<_^&Y3wPg-wuar)avmLZzfxcU5LiIHtLWh7HNF&yAA!SH%n6xyC}LCr@teSl zL4;p1MSJt@BD=&F*q;|LWc*?o-eKY)_tW8o1{WhuV+h5pa2_ywwY5`=nypJT>PjhQY@{vve&mc|eI zxbG5V=(cyWyTm^qc%^-a6Eud#FRR`V@2nH*X}Q;IaWTbLHPb-|*tYbr6DVo2kiYrd zd-mM>jnx+>BNJ&xLxe;iBr$;eu4vvcz>n0%snL{LOB@{`5!w w%P&P&e~}w{DYK~ix&{gP`qLlxpVqRFFK@hk^;dt*UcbP9|BAo=*5jxD1N*b3vj6}9 literal 0 HcmV?d00001 From 1003ae0a7c770e5167ee8b12593c8ce31c4266c8 Mon Sep 17 00:00:00 2001 From: BalazsErdos Date: Sat, 10 Jan 2026 19:53:14 +0100 Subject: [PATCH 7/9] address PR comments --- examples/plot_replicability.py | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py index dd7adc6..4fff486 100644 --- a/examples/plot_replicability.py +++ b/examples/plot_replicability.py @@ -63,7 +63,7 @@ def fit_many_parafac2(X, num_components, num_inits=5): Bs = np.empty(est_B.shape) Cs = np.empty(est_C.shape) - K = Cs.shape[0] + K = As.shape[0] for r in range(num_components): norm_Ar = tl.norm(est_A[:, r]) norm_Cr = tl.norm(est_C[:, r]) @@ -75,14 +75,14 @@ def fit_many_parafac2(X, num_components, num_inits=5): Bs[k,:,r] = est_B[k,:,r] / norm_Brk # calculate scaled B; - # since the loadings in B are specific to levels of C, we absorb C into the corresponding B for each component - cB = np.empty(Bs.shape) + # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component + aB = np.empty(Bs.shape) for r in range(num_components): - for k in range(Cs.shape[0]): - cB[k,:,r] = Cs[k,r] * Bs[k,:,r] + for k in range(As.shape[0]): + aB[k,:,r] = As[k,r] * Bs[k,:,r] - return (est_weights, (As, cB)) + return (est_weights, (Cs, aB)) ############################################################################### @@ -97,7 +97,7 @@ def truncated_normal(size): x[x < 0] = 0 return tl.tensor(x) -I, J, K = 35, 20, 25 +I, J, K = 25, 20, 35 rank = 2 A = rng.uniform(size=(I, rank)) + 0.1 # Add 0.1 to ensure that there is signal for all components for all slices @@ -192,7 +192,9 @@ def truncated_normal(size): for j, model_j in enumerate(current_models): if i >= j: # include every pair only once and omit i == j continue - fms = congruence_coefficient(m_i[1][0], m_j[1][0])[0] + weights_i, (C_i, _) = model_i + weights_j, (C_j, _) = model_j + fms = congruence_coefficient(C_i, C_j)[0] replicability_stability[rank].append(fms) ranks = sorted(replicability_stability.keys()) @@ -201,7 +203,7 @@ def truncated_normal(size): fig, ax = plt.subplots() ax.boxplot(data,whis=(0.95,0.05), positions=ranks) ax.set_xlabel("Number of components") -ax.set_ylabel("FMS_A") +ax.set_ylabel("FMS_C") plt.show() ############################################################################### @@ -224,8 +226,8 @@ def truncated_normal(size): for j, model_j in enumerate(models[rank][repeat]): if i >= j: # include every pair only once and omit i == j continue - weights_i, (aB_i, C_i) = model_i - weights_j, (aB_j, C_j) = model_j + weights_i, (C_i, aB_i) = model_i + weights_j, (C_j, aB_j) = model_j indices_subset_i = sorted(split_indices[rank][repeat][i]) indices_subset_j = sorted(split_indices[rank][repeat][j]) @@ -238,10 +240,10 @@ def truncated_normal(size): indices2use_i.append(indices_subset_i.index(common_idx)) indices2use_j.append(indices_subset_j.index(common_idx)) - cB_i = cB_i[indices2use_i, :, :] - cB_j = cB_j[indices2use_j, :, :] + aB_i = aB_i[indices2use_i, :, :] + aB_j = aB_j[indices2use_j, :, :] fms = tlviz.factor_tools.factor_match_score( - (weights_i, (A_i, cB_i.reshape(-1, cB_i.shape[2]))), (weights_j, (A_j, cB_j.reshape(-1, cB_j.shape[2]))), consider_weights=False + (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False ) replicability_alt[rank].append(fms) @@ -251,7 +253,7 @@ def truncated_normal(size): fig, ax = plt.subplots() ax.boxplot(data, positions=ranks) ax.set_xlabel("Number of components") -ax.set_ylabel("FMS_CB") +ax.set_ylabel("FMS_aB") plt.show() ############################################################################### From 0e0bd6595f24e8804479aeb12e401477e3a5937d Mon Sep 17 00:00:00 2001 From: cchatzis Date: Mon, 26 Jan 2026 17:23:03 +0100 Subject: [PATCH 8/9] addressed PR comments --- .../plot_replicability.codeobj.json | 68 ++++----------- docs/auto_examples/plot_replicability.ipynb | 10 +-- docs/auto_examples/plot_replicability.py | 54 +++++------- docs/auto_examples/plot_replicability.py.md5 | 2 +- docs/auto_examples/plot_replicability.rst | 82 ++++++++---------- docs/auto_examples/plot_replicability.zip | Bin 23609 -> 23277 bytes 6 files changed, 82 insertions(+), 134 deletions(-) diff --git a/docs/auto_examples/plot_replicability.codeobj.json b/docs/auto_examples/plot_replicability.codeobj.json index 7079e4a..2b2edd7 100644 --- a/docs/auto_examples/plot_replicability.codeobj.json +++ b/docs/auto_examples/plot_replicability.codeobj.json @@ -284,9 +284,9 @@ { "is_class": false, "is_explicit": false, - "module": "numpy", - "module_short": "numpy", - "name": "ndarray" + "module": "builtins", + "module_short": "builtins", + "name": "list" } ], "aB_i": [ @@ -298,24 +298,6 @@ "name": "ndarray" } ], - "aB_i.reshape": [ - { - "is_class": false, - "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "builtin_function_or_method" - } - ], - "aB_i.shape": [ - { - "is_class": false, - "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "tuple" - } - ], "aB_j": [ { "is_class": false, @@ -325,24 +307,6 @@ "name": "ndarray" } ], - "aB_j.reshape": [ - { - "is_class": false, - "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "builtin_function_or_method" - } - ], - "aB_j.shape": [ - { - "is_class": false, - "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "tuple" - } - ], "ax": [ { "is_class": false, @@ -822,20 +786,20 @@ "name": "asarray" } ], - "np.empty": [ + "np.concatenate": [ { "is_class": false, "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "builtin_function_or_method" + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" }, { "is_class": false, "is_explicit": false, "module": "numpy", "module_short": "numpy", - "name": "empty" + "name": "concatenate" } ], "np.inf": [ @@ -858,8 +822,8 @@ { "is_class": false, "is_explicit": false, - "module": "_cython_3_2_1", - "module_short": "_cython_3_2_1", + "module": "_cython_3_1_2", + "module_short": "_cython_3_1_2", "name": "cython_function_or_method" }, { @@ -1331,18 +1295,18 @@ { "is_class": false, "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "NoneType" + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" } ], "weights_j": [ { "is_class": false, "is_explicit": false, - "module": "builtins", - "module_short": "builtins", - "name": "NoneType" + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" } ] } \ No newline at end of file diff --git a/docs/auto_examples/plot_replicability.ipynb b/docs/auto_examples/plot_replicability.ipynb index 70bbe4d..f72ac30 100644 --- a/docs/auto_examples/plot_replicability.ipynb +++ b/docs/auto_examples/plot_replicability.ipynb @@ -40,14 +40,14 @@ }, "outputs": [], "source": [ - "def fit_many_parafac2(X, num_components, num_inits=5):\n \n best_err = np.inf\n decomposition = None\n for i in range(num_inits):\n trial_decomposition, trial_errs = parafac2_aoadmm(\n matrices=X,\n rank=num_components,\n return_errors=True,\n non_negative=[True, True, True],\n n_iter_max=500,\n absolute_tol=1e-4,\n feasibility_tol=1e-4,\n inner_tol=1e-4,\n inner_n_iter_max=5,\n feasibility_penalty_scale=5,\n tol=1e-5,\n random_state=i,\n verbose=0,\n )\n \n if best_err > trial_errs.rec_errors[-1]:\n best_err = trial_errs.rec_errors[-1]\n decomposition = trial_decomposition # note, with real data, convergence should be checked\n\n (est_weights, (est_A, est_B, est_C)) = decomposition\n est_B = np.asarray(est_B)\n\n # normalize factors\n As = np.empty(est_A.shape)\n Bs = np.empty(est_B.shape)\n Cs = np.empty(est_C.shape)\n\n K = As.shape[0]\n for r in range(num_components):\n norm_Ar = tl.norm(est_A[:, r])\n norm_Cr = tl.norm(est_C[:, r])\n As[:, r] = est_A[:, r] / norm_Ar\n Cs[:, r] = est_C[:, r] / norm_Cr\n\n for k in range(K):\n norm_Brk = tl.norm(est_B[k][:, r])\n Bs[k,:,r] = est_B[k,:,r] / norm_Brk\n\n # calculate scaled B; \n # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component\n aB = np.empty(Bs.shape)\n for r in range(num_components):\n for k in range(As.shape[0]):\n aB[k,:,r] = As[k,r] * Bs[k,:,r] \n\n\n return (est_weights, (Cs, aB))" + "def fit_many_parafac2(X, num_components, num_inits=5):\n \n best_err = np.inf\n decomposition = None\n for i in range(num_inits):\n trial_decomposition, trial_errs = parafac2_aoadmm(\n matrices=X,\n rank=num_components,\n return_errors=True,\n non_negative=[True, True, True],\n n_iter_max=500,\n absolute_tol=1e-4,\n feasibility_tol=1e-4,\n inner_tol=1e-4,\n inner_n_iter_max=5,\n feasibility_penalty_scale=5,\n tol=1e-5,\n random_state=i,\n verbose=0,\n )\n \n if best_err < trial_errs.rec_errors[-1]:\n continue\n \n best_err = trial_errs.rec_errors[-1]\n decomposition = trial_decomposition # note, with real data, convergence should be checked\n\n (est_weights, (est_A, est_B, est_C)) = decomposition\n est_B = np.asarray(est_B)\n\n # Normalize the decomposition:\n A_norm = tl.norm(est_A, axis=0)\n B_norm = tl.norm(est_B[0], axis=0) # This is the same for all B_i because of the PARAFAC2 constraint\n C_norm = tl.norm(est_C, axis=0)\n est_weights = A_norm * B_norm * C_norm # The PARAFAC2 AO-ADMM returns None as the weights\n\n As = est_A / A_norm\n Bs = [est_B_i / B_norm for est_B_i in est_B]\n Cs = est_C / C_norm\n\n # calculate scaled B; \n # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component\n aB = [a_i * B_i for a_i, B_i in zip(As, Bs)]\n\n return (est_weights, (aB, Cs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Creat simulated data\n\nSimulate noisy data with 2 components.\n\n" + "## Generate simulated data\n\nSimulate noisy data with 2 components.\n\n" ] }, { @@ -108,7 +108,7 @@ }, "outputs": [], "source": [ - "replicability_stability = {}\nfor rank in models.keys():\n replicability_stability[rank] = []\n for repeat, current_models in enumerate(models[rank]):\n for i, model_i in enumerate(current_models):\n for j, model_j in enumerate(current_models):\n if i >= j: # include every pair only once and omit i == j\n continue\n weights_i, (C_i, _) = model_i\n weights_j, (C_j, _) = model_j\n fms = congruence_coefficient(C_i, C_j)[0]\n replicability_stability[rank].append(fms)\n\nranks = sorted(replicability_stability.keys())\ndata = [np.ravel(replicability_stability[r]) for r in ranks]\n\nfig, ax = plt.subplots()\nax.boxplot(data,whis=(0.95,0.05), positions=ranks)\nax.set_xlabel(\"Number of components\")\nax.set_ylabel(\"FMS_C\")\nplt.show()" + "replicability_stability = {}\nfor rank in models.keys():\n replicability_stability[rank] = []\n for repeat, current_models in enumerate(models[rank]):\n for i, model_i in enumerate(current_models):\n for j, model_j in enumerate(current_models):\n if i >= j: # include every pair only once and omit i == j\n continue\n weights_i, (_, C_i) = model_i\n weights_j, (_, C_j) = model_j\n fms = congruence_coefficient(C_i, C_j)[0]\n replicability_stability[rank].append(fms)\n\nranks = sorted(replicability_stability.keys())\ndata = [np.ravel(replicability_stability[r]) for r in ranks]\n\nfig, ax = plt.subplots()\nax.boxplot(data,whis=(0.95,0.05), positions=ranks)\nax.set_xlabel(\"Number of components\")\nax.set_ylabel(\"FMS_C\")\nplt.show()" ] }, { @@ -133,7 +133,7 @@ }, "outputs": [], "source": [ - "replicability_alt = {}\nfor rank in models.keys():\n replicability_alt[rank] = []\n for repeat in range(repeats):\n for i, model_i in enumerate(models[rank][repeat]):\n for j, model_j in enumerate(models[rank][repeat]):\n if i >= j: # include every pair only once and omit i == j\n continue\n weights_i, (C_i, aB_i) = model_i\n weights_j, (C_j, aB_j) = model_j\n\n indices_subset_i = sorted(split_indices[rank][repeat][i])\n indices_subset_j = sorted(split_indices[rank][repeat][j])\n \n common_indices = sorted(list(set(indices_subset_i).intersection(set(indices_subset_j))))\n indices2use_i = []\n indices2use_j = []\n\n for common_idx in common_indices:\n indices2use_i.append(indices_subset_i.index(common_idx))\n indices2use_j.append(indices_subset_j.index(common_idx))\n\n aB_i = aB_i[indices2use_i, :, :]\n aB_j = aB_j[indices2use_j, :, :]\n fms = tlviz.factor_tools.factor_match_score(\n (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False\n )\n replicability_alt[rank].append(fms)\n\nranks = sorted(replicability_alt.keys())\ndata = [np.ravel(replicability_alt[r]) for r in ranks]\n\nfig, ax = plt.subplots()\nax.boxplot(data, positions=ranks)\nax.set_xlabel(\"Number of components\")\nax.set_ylabel(\"FMS_aB\")\nplt.show()" + "replicability_alt = {}\nfor rank in models.keys():\n replicability_alt[rank] = []\n for repeat in range(repeats):\n for i, model_i in enumerate(models[rank][repeat]):\n for j, model_j in enumerate(models[rank][repeat]):\n if i >= j: # include every pair only once and omit i == j\n continue\n weights_i, (aB_i, C_i) = model_i\n weights_j, (aB_j, C_j) = model_j\n\n indices_subset_i = sorted(split_indices[rank][repeat][i])\n indices_subset_j = sorted(split_indices[rank][repeat][j])\n \n common_indices = sorted(list(set(indices_subset_i).intersection(set(indices_subset_j))))\n indices2use_i = []\n indices2use_j = []\n\n for common_idx in common_indices:\n indices2use_i.append(indices_subset_i.index(common_idx))\n indices2use_j.append(indices_subset_j.index(common_idx))\n\n aB_i = np.concatenate([aB_i[idx] for idx in indices2use_i])\n aB_j = np.concatenate([aB_j[idx] for idx in indices2use_j])\n fms = tlviz.factor_tools.factor_match_score(\n (weights_i, (C_i, aB_i)), (weights_j, (C_j, aB_j)), consider_weights=False\n )\n replicability_alt[rank].append(fms)\n\nranks = sorted(replicability_alt.keys())\ndata = [np.ravel(replicability_alt[r]) for r in ranks]\n\nfig, ax = plt.subplots()\nax.boxplot(data, positions=ranks)\nax.set_xlabel(\"Number of components\")\nax.set_ylabel(\"FMS_aB\")\nplt.show()" ] }, { @@ -160,7 +160,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.9" + "version": "3.12.3" } }, "nbformat": 4, diff --git a/docs/auto_examples/plot_replicability.py b/docs/auto_examples/plot_replicability.py index 4fff486..eb1ecf8 100644 --- a/docs/auto_examples/plot_replicability.py +++ b/docs/auto_examples/plot_replicability.py @@ -51,42 +51,34 @@ def fit_many_parafac2(X, num_components, num_inits=5): verbose=0, ) - if best_err > trial_errs.rec_errors[-1]: - best_err = trial_errs.rec_errors[-1] - decomposition = trial_decomposition # note, with real data, convergence should be checked + if best_err < trial_errs.rec_errors[-1]: + continue + + best_err = trial_errs.rec_errors[-1] + decomposition = trial_decomposition # note, with real data, convergence should be checked (est_weights, (est_A, est_B, est_C)) = decomposition est_B = np.asarray(est_B) - # normalize factors - As = np.empty(est_A.shape) - Bs = np.empty(est_B.shape) - Cs = np.empty(est_C.shape) - - K = As.shape[0] - for r in range(num_components): - norm_Ar = tl.norm(est_A[:, r]) - norm_Cr = tl.norm(est_C[:, r]) - As[:, r] = est_A[:, r] / norm_Ar - Cs[:, r] = est_C[:, r] / norm_Cr + # Normalize the decomposition: + A_norm = tl.norm(est_A, axis=0) + B_norm = tl.norm(est_B[0], axis=0) # This is the same for all B_i because of the PARAFAC2 constraint + C_norm = tl.norm(est_C, axis=0) + est_weights = A_norm * B_norm * C_norm # The PARAFAC2 AO-ADMM returns None as the weights - for k in range(K): - norm_Brk = tl.norm(est_B[k][:, r]) - Bs[k,:,r] = est_B[k,:,r] / norm_Brk + As = est_A / A_norm + Bs = [est_B_i / B_norm for est_B_i in est_B] + Cs = est_C / C_norm # calculate scaled B; # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component - aB = np.empty(Bs.shape) - for r in range(num_components): - for k in range(As.shape[0]): - aB[k,:,r] = As[k,r] * Bs[k,:,r] - + aB = [a_i * B_i for a_i, B_i in zip(As, Bs)] - return (est_weights, (Cs, aB)) + return (est_weights, (aB, Cs)) ############################################################################### -# Creat simulated data +# Generate simulated data # ^^^^^^^^^^^^^^^^^^^^^^^ # # Simulate noisy data with 2 components. @@ -192,8 +184,8 @@ def truncated_normal(size): for j, model_j in enumerate(current_models): if i >= j: # include every pair only once and omit i == j continue - weights_i, (C_i, _) = model_i - weights_j, (C_j, _) = model_j + weights_i, (_, C_i) = model_i + weights_j, (_, C_j) = model_j fms = congruence_coefficient(C_i, C_j)[0] replicability_stability[rank].append(fms) @@ -226,8 +218,8 @@ def truncated_normal(size): for j, model_j in enumerate(models[rank][repeat]): if i >= j: # include every pair only once and omit i == j continue - weights_i, (C_i, aB_i) = model_i - weights_j, (C_j, aB_j) = model_j + weights_i, (aB_i, C_i) = model_i + weights_j, (aB_j, C_j) = model_j indices_subset_i = sorted(split_indices[rank][repeat][i]) indices_subset_j = sorted(split_indices[rank][repeat][j]) @@ -240,10 +232,10 @@ def truncated_normal(size): indices2use_i.append(indices_subset_i.index(common_idx)) indices2use_j.append(indices_subset_j.index(common_idx)) - aB_i = aB_i[indices2use_i, :, :] - aB_j = aB_j[indices2use_j, :, :] + aB_i = np.concatenate([aB_i[idx] for idx in indices2use_i]) + aB_j = np.concatenate([aB_j[idx] for idx in indices2use_j]) fms = tlviz.factor_tools.factor_match_score( - (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False + (weights_i, (C_i, aB_i)), (weights_j, (C_j, aB_j)), consider_weights=False ) replicability_alt[rank].append(fms) diff --git a/docs/auto_examples/plot_replicability.py.md5 b/docs/auto_examples/plot_replicability.py.md5 index 5b41514..753bd08 100644 --- a/docs/auto_examples/plot_replicability.py.md5 +++ b/docs/auto_examples/plot_replicability.py.md5 @@ -1 +1 @@ -f429dbdf7a4cf9ed22965979507994ba \ No newline at end of file +57c18f67c14eb468e9b98ee466710970 \ No newline at end of file diff --git a/docs/auto_examples/plot_replicability.rst b/docs/auto_examples/plot_replicability.rst index e3dc081..62f0d61 100644 --- a/docs/auto_examples/plot_replicability.rst +++ b/docs/auto_examples/plot_replicability.rst @@ -61,7 +61,7 @@ To fit PARAFAC2 models, we need to solve a non-convex optimization problem, poss therefore useful to fit several models with the same number of components using many different random initialisations. -.. GENERATED FROM PYTHON SOURCE LINES 31-88 +.. GENERATED FROM PYTHON SOURCE LINES 31-80 .. code-block:: Python @@ -88,38 +88,30 @@ initialisations. verbose=0, ) - if best_err > trial_errs.rec_errors[-1]: - best_err = trial_errs.rec_errors[-1] - decomposition = trial_decomposition # note, with real data, convergence should be checked + if best_err < trial_errs.rec_errors[-1]: + continue + + best_err = trial_errs.rec_errors[-1] + decomposition = trial_decomposition # note, with real data, convergence should be checked (est_weights, (est_A, est_B, est_C)) = decomposition est_B = np.asarray(est_B) - # normalize factors - As = np.empty(est_A.shape) - Bs = np.empty(est_B.shape) - Cs = np.empty(est_C.shape) - - K = As.shape[0] - for r in range(num_components): - norm_Ar = tl.norm(est_A[:, r]) - norm_Cr = tl.norm(est_C[:, r]) - As[:, r] = est_A[:, r] / norm_Ar - Cs[:, r] = est_C[:, r] / norm_Cr + # Normalize the decomposition: + A_norm = tl.norm(est_A, axis=0) + B_norm = tl.norm(est_B[0], axis=0) # This is the same for all B_i because of the PARAFAC2 constraint + C_norm = tl.norm(est_C, axis=0) + est_weights = A_norm * B_norm * C_norm # The PARAFAC2 AO-ADMM returns None as the weights - for k in range(K): - norm_Brk = tl.norm(est_B[k][:, r]) - Bs[k,:,r] = est_B[k,:,r] / norm_Brk + As = est_A / A_norm + Bs = [est_B_i / B_norm for est_B_i in est_B] + Cs = est_C / C_norm # calculate scaled B; # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component - aB = np.empty(Bs.shape) - for r in range(num_components): - for k in range(As.shape[0]): - aB[k,:,r] = As[k,r] * Bs[k,:,r] - + aB = [a_i * B_i for a_i, B_i in zip(As, Bs)] - return (est_weights, (Cs, aB)) + return (est_weights, (aB, Cs)) @@ -129,14 +121,14 @@ initialisations. -.. GENERATED FROM PYTHON SOURCE LINES 89-93 +.. GENERATED FROM PYTHON SOURCE LINES 81-85 -Creat simulated data +Generate simulated data ^^^^^^^^^^^^^^^^^^^^^^^ Simulate noisy data with 2 components. -.. GENERATED FROM PYTHON SOURCE LINES 93-121 +.. GENERATED FROM PYTHON SOURCE LINES 85-113 .. code-block:: Python @@ -175,7 +167,7 @@ Simulate noisy data with 2 components. -.. GENERATED FROM PYTHON SOURCE LINES 122-132 +.. GENERATED FROM PYTHON SOURCE LINES 114-124 The replicability analysis boils down to the following steps: @@ -188,12 +180,12 @@ The replicability analysis boils down to the following steps: 5. Repeat the above process :math:`M` times (user-chosen), to find a total of :math:`M \binom{N}{2}` comparisons. -.. GENERATED FROM PYTHON SOURCE LINES 135-137 +.. GENERATED FROM PYTHON SOURCE LINES 127-129 Split the data and fit PARAFAC2 on each data subset ^^^^^^^^^^^^^^^^^^ -.. GENERATED FROM PYTHON SOURCE LINES 137-168 +.. GENERATED FROM PYTHON SOURCE LINES 129-160 .. code-block:: Python @@ -244,7 +236,7 @@ Split the data and fit PARAFAC2 on each data subset -.. GENERATED FROM PYTHON SOURCE LINES 169-180 +.. GENERATED FROM PYTHON SOURCE LINES 161-172 Often, the mode we will be splitting within refers to different samples, depending on the use-case, it might be reasonable to retain the @@ -258,14 +250,14 @@ For example, in this case we normalize each subset to unit norm independently. Note, that ``for split_no, (train_index, _) in enumerate(rskf.split(dataset)):`` may be run in parallel for efficiency. -.. GENERATED FROM PYTHON SOURCE LINES 182-186 +.. GENERATED FROM PYTHON SOURCE LINES 174-178 Compute and assess replicability I. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Since we are subsetting the data on ``mode=0``, and the ``mode=1`` factors of PARAFAC2 are specific to the corresponding level in ``mode=0``, only the shared factor matrix can be compared using FMS: -.. GENERATED FROM PYTHON SOURCE LINES 186-209 +.. GENERATED FROM PYTHON SOURCE LINES 178-201 .. code-block:: Python @@ -278,8 +270,8 @@ to the corresponding level in ``mode=0``, only the shared factor matrix can be c for j, model_j in enumerate(current_models): if i >= j: # include every pair only once and omit i == j continue - weights_i, (C_i, _) = model_i - weights_j, (C_j, _) = model_j + weights_i, (_, C_i) = model_i + weights_j, (_, C_j) = model_j fms = congruence_coefficient(C_i, C_j)[0] replicability_stability[rank].append(fms) @@ -304,12 +296,12 @@ to the corresponding level in ``mode=0``, only the shared factor matrix can be c -.. GENERATED FROM PYTHON SOURCE LINES 210-212 +.. GENERATED FROM PYTHON SOURCE LINES 202-204 Here, we observe that over-estimating the number of components results in a loss of replicable of the patterns, indicated by low FMS. -.. GENERATED FROM PYTHON SOURCE LINES 214-220 +.. GENERATED FROM PYTHON SOURCE LINES 206-212 Compute and assess replicability II. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -318,7 +310,7 @@ including the factors corresponding to ``mode=0``, and ``mode=1`` :cite:p:`erdos By using only the indices present in both subsets (e.g. the factors corresponding to the subjects' data included in both factorizations) -.. GENERATED FROM PYTHON SOURCE LINES 220-259 +.. GENERATED FROM PYTHON SOURCE LINES 212-251 .. code-block:: Python @@ -331,8 +323,8 @@ corresponding to the subjects' data included in both factorizations) for j, model_j in enumerate(models[rank][repeat]): if i >= j: # include every pair only once and omit i == j continue - weights_i, (C_i, aB_i) = model_i - weights_j, (C_j, aB_j) = model_j + weights_i, (aB_i, C_i) = model_i + weights_j, (aB_j, C_j) = model_j indices_subset_i = sorted(split_indices[rank][repeat][i]) indices_subset_j = sorted(split_indices[rank][repeat][j]) @@ -345,10 +337,10 @@ corresponding to the subjects' data included in both factorizations) indices2use_i.append(indices_subset_i.index(common_idx)) indices2use_j.append(indices_subset_j.index(common_idx)) - aB_i = aB_i[indices2use_i, :, :] - aB_j = aB_j[indices2use_j, :, :] + aB_i = np.concatenate([aB_i[idx] for idx in indices2use_i]) + aB_j = np.concatenate([aB_j[idx] for idx in indices2use_j]) fms = tlviz.factor_tools.factor_match_score( - (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False + (weights_i, (C_i, aB_i)), (weights_j, (C_j, aB_j)), consider_weights=False ) replicability_alt[rank].append(fms) @@ -373,7 +365,7 @@ corresponding to the subjects' data included in both factorizations) -.. GENERATED FROM PYTHON SOURCE LINES 260-265 +.. GENERATED FROM PYTHON SOURCE LINES 252-257 ``common_indices`` contains the indices (e.g. subjects/samples) present in both subsets, but since the position of each index can change (e.g. sample no 3 is not guaranteeed at @@ -385,7 +377,7 @@ Similar results can be observed with this approach in terms of the replicability .. rst-class:: sphx-glr-timing - **Total running time of the script:** (8 minutes 4.271 seconds) + **Total running time of the script:** (12 minutes 12.010 seconds) .. _sphx_glr_download_auto_examples_plot_replicability.py: diff --git a/docs/auto_examples/plot_replicability.zip b/docs/auto_examples/plot_replicability.zip index efd37bc7a6057968fee99897b3f4a05aa20af9e6..5b9e2b5b5ac51a3245975bc19792a80dc0d5b28d 100644 GIT binary patch delta 1588 zcmbu9TTc@~6vwwMt?POt-oVD1F~riXwN&Ie`WfnV2B>2$q zrV?V}JgW~z6CO>~NaCCFY1G7z;LL8>m862=_Odf`F8_0Oe>0n(d_O+;I;O@0Z6_sh zeSDlQoShE#%&L+!`>x}T>>FV(-GQ`nl^&%VlA^_*E-FsgV30jrF_94eU$ zizdYuHc46m%?s*;yMa+)Rg+>AL)6k?0vOAfj;2)CFVo5|qosKZW&X+ze2E1 zjRgee%^EPORglp!xQ9x{5=FqAw;iTn9$tspUQjX#wQ2??mMfQSC6{%zTt89mLxk^_ zEy{anHWTM+#H{nS;-Hzxt_=ZK65hH5$0jF14OvSyLP5?X2pCS`&3uY00!427xqzOg zhHJ|Rn8g} zXLc*bIV1$5D{A?{HNo%Rb&G=Qt!J577JH{p>vRoj7 z*aG2rC;a^OSvJ^E|ASpAcVa}`9+-QM|H^qCfqii9>Foi!=lBQG>&=4yhx7vRWLmK9 z8C$?UAWjdl$+x$F3%eKAiuR4cjqbmsoqZoV7jjrumx!n)}?NBlV+3drifmm z2alfK&5uK+mmVbudhj5{g9nv*6A$7+Jb4j37M$7L>~yn9V+WGmo%zq3|9dk(=J%_% zH?Q0Jt|h%4Js#42K8 zy4jx|&D!x-iyG39q9RRFZ_5#Nqadyn3#Fnet7V-{*g6t}BqDMej#OTYlMmB!IuS!! zmRh!@f+4N!(+Z8WoNo5wxJC)AuXDBa93BLx*b0a&v7j*yHOa0kPt!{R7wW!Y;O%mMb%zZy28k?e2S)rEmAHaW0SUGz*WEap<3!V(N2+NKB-J!0GNV?g ztd2=K(BhC#Au)Y0l`d_iAevk@z|tvh#idi-3e1dZ>KBM-(&|okl1(@R5=~1|VQ8s) z_&?x;kJH<;ZY=4F2bCPYU+G6Tz-RCqo!=r_59cQb1~!C-P~ws^e8WHH6o znEX2b*l#mgNS!g5e381^o@~k_&5Kl#w3bk=FXmjp9R|q1i$gx|MEpFQT;_s~ixFRr ey~cHE{%yhC Date: Mon, 26 Jan 2026 17:31:04 +0100 Subject: [PATCH 9/9] Addressed PR --- docs/sg_execution_times.rst | 52 ++++++++++++++++++++++++++++++++ examples/plot_replicability.py | 54 +++++++++++++++------------------- 2 files changed, 75 insertions(+), 31 deletions(-) create mode 100644 docs/sg_execution_times.rst diff --git a/docs/sg_execution_times.rst b/docs/sg_execution_times.rst new file mode 100644 index 0000000..e1c02e5 --- /dev/null +++ b/docs/sg_execution_times.rst @@ -0,0 +1,52 @@ + +:orphan: + +.. _sphx_glr_sg_execution_times: + + +Computation times +================= +**12:12.010** total execution time for 6 files **from all galleries**: + +.. container:: + + .. raw:: html + + + + + + + + .. list-table:: + :header-rows: 1 + :class: table table-striped sg-datatable + + * - Example + - Time + - Mem (MB) + * - :ref:`sphx_glr_auto_examples_plot_replicability.py` (``../examples/plot_replicability.py``) + - 12:12.010 + - 0.0 + * - :ref:`sphx_glr_auto_examples_plot_bikesharing.py` (``../examples/plot_bikesharing.py``) + - 00:00.000 + - 0.0 + * - :ref:`sphx_glr_auto_examples_plot_custom_penalty.py` (``../examples/plot_custom_penalty.py``) + - 00:00.000 + - 0.0 + * - :ref:`sphx_glr_auto_examples_plot_examining_different_number_of_components.py` (``../examples/plot_examining_different_number_of_components.py``) + - 00:00.000 + - 0.0 + * - :ref:`sphx_glr_auto_examples_plot_semiconductor_etch_analysis.py` (``../examples/plot_semiconductor_etch_analysis.py``) + - 00:00.000 + - 0.0 + * - :ref:`sphx_glr_auto_examples_plot_simulated_nonnegative.py` (``../examples/plot_simulated_nonnegative.py``) + - 00:00.000 + - 0.0 diff --git a/examples/plot_replicability.py b/examples/plot_replicability.py index 4fff486..eb1ecf8 100644 --- a/examples/plot_replicability.py +++ b/examples/plot_replicability.py @@ -51,42 +51,34 @@ def fit_many_parafac2(X, num_components, num_inits=5): verbose=0, ) - if best_err > trial_errs.rec_errors[-1]: - best_err = trial_errs.rec_errors[-1] - decomposition = trial_decomposition # note, with real data, convergence should be checked + if best_err < trial_errs.rec_errors[-1]: + continue + + best_err = trial_errs.rec_errors[-1] + decomposition = trial_decomposition # note, with real data, convergence should be checked (est_weights, (est_A, est_B, est_C)) = decomposition est_B = np.asarray(est_B) - # normalize factors - As = np.empty(est_A.shape) - Bs = np.empty(est_B.shape) - Cs = np.empty(est_C.shape) - - K = As.shape[0] - for r in range(num_components): - norm_Ar = tl.norm(est_A[:, r]) - norm_Cr = tl.norm(est_C[:, r]) - As[:, r] = est_A[:, r] / norm_Ar - Cs[:, r] = est_C[:, r] / norm_Cr + # Normalize the decomposition: + A_norm = tl.norm(est_A, axis=0) + B_norm = tl.norm(est_B[0], axis=0) # This is the same for all B_i because of the PARAFAC2 constraint + C_norm = tl.norm(est_C, axis=0) + est_weights = A_norm * B_norm * C_norm # The PARAFAC2 AO-ADMM returns None as the weights - for k in range(K): - norm_Brk = tl.norm(est_B[k][:, r]) - Bs[k,:,r] = est_B[k,:,r] / norm_Brk + As = est_A / A_norm + Bs = [est_B_i / B_norm for est_B_i in est_B] + Cs = est_C / C_norm # calculate scaled B; # since the loadings in B are specific to levels of A, we absorb A into the corresponding B for each component - aB = np.empty(Bs.shape) - for r in range(num_components): - for k in range(As.shape[0]): - aB[k,:,r] = As[k,r] * Bs[k,:,r] - + aB = [a_i * B_i for a_i, B_i in zip(As, Bs)] - return (est_weights, (Cs, aB)) + return (est_weights, (aB, Cs)) ############################################################################### -# Creat simulated data +# Generate simulated data # ^^^^^^^^^^^^^^^^^^^^^^^ # # Simulate noisy data with 2 components. @@ -192,8 +184,8 @@ def truncated_normal(size): for j, model_j in enumerate(current_models): if i >= j: # include every pair only once and omit i == j continue - weights_i, (C_i, _) = model_i - weights_j, (C_j, _) = model_j + weights_i, (_, C_i) = model_i + weights_j, (_, C_j) = model_j fms = congruence_coefficient(C_i, C_j)[0] replicability_stability[rank].append(fms) @@ -226,8 +218,8 @@ def truncated_normal(size): for j, model_j in enumerate(models[rank][repeat]): if i >= j: # include every pair only once and omit i == j continue - weights_i, (C_i, aB_i) = model_i - weights_j, (C_j, aB_j) = model_j + weights_i, (aB_i, C_i) = model_i + weights_j, (aB_j, C_j) = model_j indices_subset_i = sorted(split_indices[rank][repeat][i]) indices_subset_j = sorted(split_indices[rank][repeat][j]) @@ -240,10 +232,10 @@ def truncated_normal(size): indices2use_i.append(indices_subset_i.index(common_idx)) indices2use_j.append(indices_subset_j.index(common_idx)) - aB_i = aB_i[indices2use_i, :, :] - aB_j = aB_j[indices2use_j, :, :] + aB_i = np.concatenate([aB_i[idx] for idx in indices2use_i]) + aB_j = np.concatenate([aB_j[idx] for idx in indices2use_j]) fms = tlviz.factor_tools.factor_match_score( - (weights_i, (C_i, aB_i.reshape(-1, aB_i.shape[2]))), (weights_j, (C_j, aB_j.reshape(-1, aB_j.shape[2]))), consider_weights=False + (weights_i, (C_i, aB_i)), (weights_j, (C_j, aB_j)), consider_weights=False ) replicability_alt[rank].append(fms)