glhmm.utils

Some public useful functions - Gaussian Linear Hidden Markov Model @author: Diego Vidaurre 2023

glhmm.utils.get_FO(Gamma, indices, summation=False)[source]

Calculates the fractional occupancy of each state.

Parameters:

Gammaarray-like, shape (n_samples, n_states)

The state probability time series.

indicesarray-like, shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

summationbool, optional, default=False

If True, the sum of each row is not normalized, otherwise it is.

Returns:

FOarray-like, shape (n_sessions, n_states)

The fractional occupancy of each state per session.

glhmm.utils.get_FO_entropy(Gamma, indices)[source]

Calculates the entropy of each session, if we understand fractional occupancies as probabilities.

Parameters:

Gammaarray-like of shape (n_samples, n_states)

The Gamma represents the state probability timeseries.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

Returns:

entropyarray-like of shape (n_sessions,)

The entropy of each session.

glhmm.utils.get_gamma_similarity(gamma1, gamma2)[source]

Computes a measure of similarity between two sets of state time courses.

These can have different numbers of states, but they must have the same number of time points.

Parameters:

gamma1numpy.ndarray

First set of state time courses with shape (T, K).

gamma2numpy.ndarray

Second set of state time courses with shape (T, K2), where K2 may be different from K.

Returns:

Sfloat

Similarity, measured as the sum of joint probabilities under the optimal state alignment.

assignumpy.ndarray

Optimal state alignment for gamma2 (uses Munkres’ algorithm).

gamma2numpy.ndarray

The second set of state time courses reordered to match gamma1.

glhmm.utils.get_life_times(vpath, indices, threshold=0)[source]

Calculates the average, median and maximum life times for each state.

Parameters:

vpatharray-like of shape (n_samples,)

The viterbi path represents the most likely state sequence.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

thresholdint, optional, default=0

A threshold value used to exclude visits with a duration below this value.

Returns:

meanLFarray-like of shape (n_sessions, n_states)

The average visit duration for each state in each trial/session.

medianLFarray-like of shape (n_sessions, n_states)

The median visit duration for each state in each trial/session.

maxLFarray-like of shape (n_sessions, n_states)

The maximum visit duration for each state in each trial/session.

Notes:

A visit to a state is defined as a contiguous sequence of time points in which the state is active. The duration of a visit is the number of time points in the sequence. This function uses the get_visits function to compute the visits and exclude those below the threshold.

glhmm.utils.get_maxFO(Gamma, indices)[source]

Calculates the maximum fractional occupancy per session.

The first argument can also be a viterbi path (vpath).

Parameters:

Gammaarray-like of shape (n_samples, n_states); or a vpath, array of shape (n_samples,)

The Gamma represents the state probability timeseries and the vpath represents the most likely state sequence.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

Returns:

maxFO: array-like of shape (n_sessions,)

The maximum fractional occupancy across states for each trial/session

Notes:

The maxFO is useful to assess the amount of state mixing. For more information, see [^1].

References:

[^1]: Ahrends, R., et al. (2022). Data and model considerations for estimating time-varying functional connectivity in fMRI. NeuroImage 252, 119026.

https://pubmed.ncbi.nlm.nih.gov/35217207/)

glhmm.utils.get_state_evoked_response(Gamma, indices)[source]

Calculates the state evoked response

The first argument can also be a viterbi path (vpath).

Parameters:

Gammaarray-like of shape (n_samples, n_states), or a vpath array of shape (n_samples,)

The Gamma represents the state probability timeseries and the vpath represents the most likely state sequence.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

Returns:

serarray-like of shape (n_samples, n_states)

The state evoked response matrix.

Raises:

Exception

If the input data violates any of the following conditions: - There is only one trial/session - Not all trials/sessions have the same length.

glhmm.utils.get_state_evoked_response_entropy(Gamma, indices)[source]

Calculates the entropy of each time point, if we understand state evoked responses as probabilities.

Parameters:

Gamma: array-like of shape (n_samples, n_states)

The Gamma represents the state probability timeseries.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

Returns:

entropy: array-like of shape (n_samples,)

The entropy of each time point.

glhmm.utils.get_state_onsets(vpath, indices, threshold=0)[source]

Calculates the state onsets, i.e., the time points when each state activates.

Parameters:

vpatharray-like of shape (n_samples, n_states)

The viterbi path represents the most likely state sequence.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

thresholdint, optional, default=0

A threshold value used to exclude visits with a duration below this value.

Returns:

onsetslist of lists of ints

A list of the time points when each state activates for each trial/session.

Notes:

A visit to a state is defined as a contiguous sequence of time points in which the state is active. This function uses the get_visits function to compute the visits and exclude those below the threshold.

glhmm.utils.get_switching_rate(Gamma, indices)[source]

Calculates the switching rate.

The first argument can also be a viterbi path (vpath).

Parameters:

Gammaarray-like of shape (n_samples, n_states), or a vpath array of shape (n_samples,)

The Gamma represents the state probability timeseries and the vpath represents the most likely state sequence.

indicesarray-like of shape (n_sessions, 2)

The start and end indices of each trial/session in the input data.

Returns:

SRarray-like of shape (n_sessions, n_states)

The switching rate matrix.

glhmm.utils.get_visits(vpath, k, threshold=0)[source]

Computes a list of visits for state k, given a viterbi path (vpath).

Parameters:

vpatharray-like of shape (n_samples,)

The viterbi path represents the most likely state sequence.

kint

The state for which to compute the visits.

thresholdint, optional, default=0

A threshold value used to exclude visits with a duration below this value.

Returns:

lengthslist of floats

A list of visit durations for state k, where each duration is greater than the threshold.

onsetslist of ints

A list of onset time points for each visit.

Notes:

A visit to state k is defined as a contiguous sequence of time points in which state k is active.

glhmm.utils.load_stability_results(save_dir)[source]

Load HMM stability training results from disk.

Handles two cases automatically: - summary_results.pkl present -> loads directly (fast path). - Only individual hmm_K*_rep*.pkl files -> rebuilds the summary from them

by recomputing all N*(N-1)/2 pairwise Gamma similarities across repetitions.

Parameters:

save_dir (str or Path):

Directory where run_stability_training() saved its outputs.

Returns:

results (dict):

Dictionary with K values as keys, each containing: - ‘FE’: list of free energy arrays, one per repetition. - ‘similarity_scores’: list of Gamma similarity floats (all pairwise comparisons).

state_range (list of int):

Sorted list of K values found in the directory.

glhmm.utils.osf_download_data(osf_url, data_dir='data', folder=None)[source]

Download files from an OSF project to a local directory.

Queries the OSF storage API for the given project, optionally navigates into a named sub-folder, and downloads every file that does not yet exist locally. Files already present are silently skipped, so the function is safe to re-run.

Parameters:

osf_url (str):

OSF project URL (e.g. 'https://osf.io/8qcyj/') or bare project identifier (e.g. '8qcyj'). The project ID is extracted automatically so any standard OSF URL format works.

data_dir (str or Path, optional), default=``’data’``:

Local directory to download files into. Created automatically if it does not exist.

folder (str or None, optional), default=None:

Name of a sub-folder inside the project’s OSF Storage to download from. None downloads all files from the storage root level.

Returns:

None

Examples:

Download all files from the root of a project:

utils.osf_download_data('https://osf.io/8qcyj/')

Download files from a specific sub-folder:

utils.osf_download_data('https://osf.io/8qcyj/', folder='Simulation_data_numpy')
glhmm.utils.run_stability_training(Y, indices, state_range, n_repeats, save_dir, log_preproc=None, covtype='full', model_mean='no', options=None)[source]

Train HMMs across a range of K values to assess solution stability.

For each K and random repetition: initialises an HMM; trains with full-batch EM until convergence; saves the model to disk; then computes all N*(N-1)/2 pairwise Gamma similarities across repetitions to measure how reproducible the state solution is across random initialisations.

Parameters:

Y (numpy.ndarray):

Preprocessed data array of shape (n_total_timepoints, n_features), with all subjects concatenated along the time axis.

indices (numpy.ndarray):

Start and end indices for each subject, shape (n_subjects, 2).

state_range (iterable of int):

K values to test, e.g. range(5, 13).

n_repeats (int):

Number of independent random initialisations per K value.

save_dir (str or Path):

Directory to write per-model pickle files and the summary_results.pkl summary.

log_preproc (preprocessing log or None, optional), default=None:

Log returned by preproc.preprocess_data(). Passed as preproclogY to the HMM so that state parameters can be back-transformed to the original space.

covtype (str, optional), default=’full’:

Covariance type passed to glhmm(). Options: 'full' (state-specific FC matrices), 'diag' (diagonal, faster), 'sharedfull' (one shared FC matrix), 'shareddiag'.

model_mean (str, optional), default=’no’:

Whether to model per-state activation means. Use 'no' for standardised data; 'state' if activation levels carry information.

options (dict or None, optional), default=None:

Training options passed to hmm.train(). Defaults to {'cyc': 500, 'min_cyc': 25, 'tol': 1e-5, 'verbose': False}.

Returns:

results (dict):

Dictionary with K values as keys, each containing: - ‘FE’: list of free energy arrays, one per repetition. - ‘similarity_scores’: list of Gamma similarity floats (all pairwise comparisons).

glhmm.utils.run_stability_training_stochastic(files, state_range, n_repeats, save_dir, log_preproc=None, covtype='full', model_mean='no', options=None)[source]

Train HMMs stochastically across a range of K values to assess solution stability.

For each K and random repetition: initialises an HMM; trains with stochastic mini-batch EM; calls hmm.decode() to obtain the Gamma time series (stochastic training returns empty Gamma by design); saves the model to disk; then computes all N*(N-1)/2 pairwise Gamma similarities across repetitions.

Use this function when your dataset is too large to hold in RAM. Data must be split into one .npy or .npz file per subject on disk (see io.save_subjects_file()). For in-memory data, use run_stability_training().

Parameters:

files (list of str or Path):

Paths to per-subject preprocessed data files (one file per subject).

state_range (iterable of int):

K values to test, e.g. range(5, 13).

n_repeats (int):

Number of independent random initialisations per K value.

save_dir (str or Path):

Directory to write per-model pickle files and the summary_results.pkl summary.

log_preproc (preprocessing log or None, optional), default=None:

Log returned by preproc.preprocess_data(). Passed as preproclogY to the HMM so that state parameters can be back-transformed to the original space.

covtype (str, optional), default=’full’:

Covariance type passed to glhmm(). Options: 'full' (state-specific FC matrices), 'diag' (diagonal, faster), 'sharedfull' (one shared FC matrix), 'shareddiag'.

model_mean (str, optional), default=’no’:

Whether to model per-state activation means. Use 'no' for standardised data; 'state' if activation levels carry information.

options (dict or None, optional), default=None:

Training options passed to hmm.train(). stochastic is always set to True. Defaults to {'Nbatch': 20, 'initNbatch': 20, 'initcyc': 50, 'cyc': 500, 'min_cyc': 100, 'forget_rate': 0.5, 'base_weights': 0.75, 'cyc_to_go_under_th': 10, 'deactivate_states': False, 'verbose': False}.

Returns:

results (dict):

Dictionary with K values as keys, each containing: - ‘FE’: list of free energy arrays, one per repetition. - ‘similarity_scores’: list of Gamma similarity floats (all pairwise comparisons).