diff --git a/docs/auto_examples/plot_replicability.codeobj.json b/docs/auto_examples/plot_replicability.codeobj.json new file mode 100644 index 0000000..2b2edd7 --- /dev/null +++ b/docs/auto_examples/plot_replicability.codeobj.json @@ -0,0 +1,1312 @@ +{ + "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": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "aB_i": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "aB_j": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "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.concatenate": [ + { + "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": "concatenate" + } + ], + "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_1_2", + "module_short": "_cython_3_1_2", + "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": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "weights_j": [ + { + "is_class": false, + "is_explicit": false, + "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 new file mode 100644 index 0000000..f72ac30 --- /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 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": [ + "## Generate 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, (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()" + ] + }, + { + "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.12.3" + } + }, + "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..eb1ecf8 --- /dev/null +++ b/docs/auto_examples/plot_replicability.py @@ -0,0 +1,257 @@ +""" +.. _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]: + 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 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 + + 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 = [a_i * B_i for a_i, B_i in zip(As, Bs)] + + return (est_weights, (aB, Cs)) + + +############################################################################### +# Generate 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, (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]) + + 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 = 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)), (weights_j, (C_j, aB_j)), 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..753bd08 --- /dev/null +++ b/docs/auto_examples/plot_replicability.py.md5 @@ -0,0 +1 @@ +57c18f67c14eb468e9b98ee466710970 \ 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..62f0d61 --- /dev/null +++ b/docs/auto_examples/plot_replicability.rst @@ -0,0 +1,406 @@ + +.. 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-80 + +.. 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]: + 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 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 + + 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 = [a_i * B_i for a_i, B_i in zip(As, Bs)] + + return (est_weights, (aB, Cs)) + + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 81-85 + +Generate simulated data +^^^^^^^^^^^^^^^^^^^^^^^ + +Simulate noisy data with 2 components. + +.. GENERATED FROM PYTHON SOURCE LINES 85-113 + +.. 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 114-124 + +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 127-129 + +Split the data and fit PARAFAC2 on each data subset +^^^^^^^^^^^^^^^^^^ + +.. GENERATED FROM PYTHON SOURCE LINES 129-160 + +.. 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 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 +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 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 178-201 + +.. 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 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 206-212 + +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 212-251 + +.. 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, (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]) + + 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 = 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)), (weights_j, (C_j, aB_j)), 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 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 +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:** (12 minutes 12.010 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 0000000..5b9e2b5 Binary files /dev/null and b/docs/auto_examples/plot_replicability.zip differ diff --git a/docs/references.bib b/docs/references.bib index d63b1d3..fc78fa6 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -305,6 +305,14 @@ @INPROCEEDINGS{chatzis2023timeaware pages={1-6}, 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} +} @article{chatzis2025dcmf, title={dCMF: Learning interpretable evolving patterns from temporal multiway data}, author={Chatzis, Christos and Schenker, Carla and Cohen, J{\'e}r{\'e}my E and Acar, Evrim}, 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 new file mode 100644 index 0000000..eb1ecf8 --- /dev/null +++ b/examples/plot_replicability.py @@ -0,0 +1,257 @@ +""" +.. _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]: + 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 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 + + 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 = [a_i * B_i for a_i, B_i in zip(As, Bs)] + + return (est_weights, (aB, Cs)) + + +############################################################################### +# Generate 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, (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]) + + 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 = 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)), (weights_j, (C_j, aB_j)), 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