Fork me on GitHub

pyhrf.vbjde.vem_tools module

TOOLS and FUNCTIONS for VEM JDE Used in different versions of VEM

pyhrf.vbjde.vem_tools.A_Entropy(Sigma_A, M, J)
pyhrf.vbjde.vem_tools.Compute_FreeEnergy(y_tilde, m_A, Sigma_A, mu_Ma, sigma_Ma, m_H, Sigma_H, AuxH, R, R_inv, sigmaH, sigmaG, m_C, Sigma_C, mu_Mc, sigma_Mc, m_G, Sigma_G, AuxG, q_Z, neighboursIndexes, Beta, Gamma, gamma, gamma_h, gamma_g, sigma_eps, XX, W, J, D, M, N, K, hyp, Gamma_X, Gamma_WX, plot=False, bold=False, S=1)
pyhrf.vbjde.vem_tools.H_Entropy(Sigma_H, D)
pyhrf.vbjde.vem_tools.PolyMat(Nscans, paramLFD, tr)

Build polynomial basis

pyhrf.vbjde.vem_tools.Q_Entropy(q_Z, M, J)
pyhrf.vbjde.vem_tools.Q_expectation_Ptilde(q_Z, neighboursIndexes, Beta, gamma, K, M)
pyhrf.vbjde.vem_tools.RF_Entropy(Sigma_RF, D)
pyhrf.vbjde.vem_tools.RF_expectation_Ptilde(m_X, Sigma_X, sigmaX, R, R_inv, D)
pyhrf.vbjde.vem_tools.RL_Entropy(Sigma_RL, M, J)
pyhrf.vbjde.vem_tools.RL_expectation_Ptilde(m_X, Sigma_X, mu_Mx, sigma_Mx, q_Z)
pyhrf.vbjde.vem_tools.Z_Entropy(q_Z, M, J)
pyhrf.vbjde.vem_tools.beta_gradient(beta, labels_proba, labels_neigh, neighbours_indexes, gamma, gradient_method='m1')

Computes the gradient of the beta function.

The maximization of f(\beta^{m}) needs the computation of its derivative with respect to \beta^{m}.

Method 1

\frac{\partial f(\beta^{m})}{\partial \beta^{m}} = -\frac{1}{2} \sum\limits_{j}
\sum\limits_{k \in N(j)} \sum\limits_{i \in \{0,1\}} \left\{ p_{mf_{j}}(i)p_{mf_{k}}(i) -
\widetilde{p}_{q^{m}_{j}}(i)\widetilde{p}_{q^{m}_{k}}(i) \right\} - \lambda_{\beta}

Method 2

\frac{\partial f(\beta^{m})}{\partial \beta^{m}} = - \sum\limits_{j}
\sum\limits_{k \in N(j)} \sum\limits_{i \in \{0,1\}} \widetilde{p}_{q^{m}_{k}} (i) \left\{ p_{mf_{j}} (i) -
\frac{1}{2}\widetilde{p}_{q^{m}_{j}} (i) \right\} - \lambda_{\beta}

where

p_{mf_{j}} (i) = \frac{\exp \left( \beta \sum\limits_{k \in N(j)} \widetilde{p}_{q^{m}_{k}} (i) \right)}
{\sum\limits_{i \in \{0,1\}} \exp \left( \beta \sum\limits_{k \in N(j)} \widetilde{p}_{q^{m}_{k}} (i) \right)}

Parameters:
  • beta (float) –
  • labels_proba (ndarray) –
  • labels_neigh (ndarray) –
  • neighbours_indexes (ndarray) –
  • gamma (float) –
  • gradient_method (str) – for testing purposes
Returns:

gradient – the gradient estimated in beta

Return type:

float

pyhrf.vbjde.vem_tools.beta_maximization(beta, labels_proba, neighbours_indexes, gamma)

Computes the Beta Maximization step of the JDE VEM algorithm.

The maximization over each \beta^{m} corresponds to the M-step obtained for a standard Hiddden MRF model:

\hat{\beta}^{m} = \underset{\beta^{m}}{\mathrm{arg\, max}} f(\beta^{m})

Parameters:
  • beta (ndarray) – initial value of beta
  • labels_proba (ndarray) –
  • neighbours_indexes (ndarray) –
  • gamma (float) –
Returns:

  • beta (float) – the new value of beta
  • success (bool) – True if the maximization has succeeded

Notes

See beta_gradient() function.

pyhrf.vbjde.vem_tools.buildFiniteDiffMatrix(order, size, regularization=None)

Build the finite difference matrix used for the hrf regularization prior.

Parameters:
  • order (int) – difference order (see numpy.diff function)
  • size (int) – size of the matrix
  • regularization (array like, optional) – one dimensional vector of factors used for regularizing the hrf
Returns:

diffMat – the finite difference matrix

Return type:

ndarray, shape (size, size)

pyhrf.vbjde.vem_tools.computeFit(hrf_mean, nrls_mean, X, nb_voxels, nb_scans)

Compute the estimated induced signal by each stimulus.

Parameters:
  • hrf_mean (ndarray) –
  • nrls_mean (ndarray) –
  • X (OrderedDict) –
  • nb_voxels (int) –
  • nb_scans (int) –
Returns:

Return type:

ndarray

pyhrf.vbjde.vem_tools.computeFit_asl(H, m_A, G, m_C, W, XX)

Compute Fit

pyhrf.vbjde.vem_tools.compute_contrasts(condition_names, contrasts, m_A, m_C, Sigma_A, Sigma_C, M, J)
pyhrf.vbjde.vem_tools.compute_mat_X_2(nbscans, tr, lhrf, dt, onsets, durations=None)
pyhrf.vbjde.vem_tools.constraint_norm1_b(Ftilde, Sigma_F, positivity=False, perfusion=None)

Constrain with optimization strategy

pyhrf.vbjde.vem_tools.contrasts_mean_var_classes(contrasts, condition_names, nrls_mean, nrls_covar, nrls_class_mean, nrls_class_var, nb_contrasts, nb_classes, nb_voxels)

Computes the contrasts nrls from the conditions nrls and the mean and variance of the gaussian classes of the contrasts (in the cases of all inactive conditions and all active conditions).

Parameters:
  • contrasts (OrderedDict) – TODO.
  • condition_names (list) – TODO.
  • nrls_mean (ndarray, shape (nb_voxels, nb_conditions)) – TODO.
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels)) – TODO.
  • nrls_class_mean (ndarray, shape (nb_conditions, nb_classes)) – TODO.
  • nrls_class_var (ndarray, shape (nb_conditions, nb_classes)) – TODO.
  • nb_contrasts (int) –
  • nb_classes (int) –
  • nb_voxels (int) –
Returns:

  • contrasts_mean (ndarray, shape (nb_voxels, nb_contrasts))
  • contrasts_var (ndarray, shape (nb_voxels, nb_contrasts))
  • contrasts_class_mean (ndarray, shape (nb_contrasts, nb_classes))
  • contrasts_class_var (ndarray, shape (nb_contrasts, nb_classes))

pyhrf.vbjde.vem_tools.cosine_drifts_basis(nb_scans, param_lfd, tr)

Build cosine drifts basis.

Parameters:
  • nb_scans (int) –
  • param_lfd (int) – TODO
  • tr (float) –
Returns:

drifts_basis – K is determined by the scipy.linalg.orth function and corresponds to the effective rank of the matrix it is applied to (see function’s docstring)

Return type:

ndarray, shape (nb_scans, K)

pyhrf.vbjde.vem_tools.covariance_matrix(order, D, dt)
pyhrf.vbjde.vem_tools.create_conditions(onsets, durations, nb_conditions, nb_scans, hrf_len, tr, dt)

Generate the occurrences matrix.

Parameters:
  • onsets (dict) – dictionary of onsets for each condition.
  • durations (dict) – dictionary of durations for each condition.
  • nb_conditions (int) – number of experimental conditions.
  • nb_scans (int) – number of scans.
  • hrf_len (int) – number of points of the hrf
  • tr (float) – time of repetition
  • dt (float) – hrf temporal precision
Returns:

  • X (dict) – dictionary of the occurence matrix
  • occurence_matrix (ndarray)
  • condition_names (list)

pyhrf.vbjde.vem_tools.create_neighbours(graph)

Transforms the graph list in ndarray. This is for performances purposes. Sets the empty neighbours to -1.

Parameters:graph (list of ndarray) – each graph[i] represents the list of neighbours of the ith voxel
Returns:neighbours_indexes
Return type:ndarray
pyhrf.vbjde.vem_tools.drifts_coeffs_fit(signal, drift_basis)

# TODO

Parameters:
  • signal (ndarray, shape (nb_scans, nb_voxels)) –
  • drift_basis (ndarray, shape (nb_scans, int)) –
Returns:

drift_coeffs

Return type:

ndarray, shape

pyhrf.vbjde.vem_tools.expectation_A_asl(H, G, m_C, W, XX, Gamma, Gamma_X, q_Z, mu_Ma, sigma_Ma, J, y_tilde, Sigma_H, sigma_eps_m)

Expectation-A step:

p_A = argmax_h(E_pc,pq,ph,pg[log p(a|y, h, c, g, q; theta)]) \propto exp(E_pc,ph,pg[log p(y|h, a, c, g; theta)] + E_pq[log p(a|q; mu_Ma, sigma_Ma)])

Returns:
Return type:m_A, Sigma_A of probability distribution p_A of the current iteration
pyhrf.vbjde.vem_tools.expectation_A_ms(m_A, Sigma_A, H, G, m_C, W, XX, Gamma, Gamma_X, q_Z, mu_Ma, sigma_Ma, J, y_tilde, Sigma_H, sigma_eps_m, N, M, D, S)

Expectation-A step:

p_A = argmax_h(E_pc,pq,ph,pg[log p(a|y, h, c, g, q; theta)]) \propto exp(E_pc,ph,pg[log p(y|h, a, c, g; theta)] + E_pq[log p(a|q; mu_Ma, sigma_Ma)])

Returns:
Return type:m_A, Sigma_A of probability distribution p_A of the current iteration
pyhrf.vbjde.vem_tools.expectation_C_asl(G, H, m_A, W, XX, Gamma, Gamma_X, q_Z, mu_Mc, sigma_Mc, J, y_tilde, Sigma_G, sigma_eps_m)

Expectation-C step:

p_C = argmax_h(E_pa,pq,ph,pg[log p(a|y, h, a, g, q; theta)]) \propto exp(E_pa,ph,pg[log p(y|h, a, c, g; theta)] + E_pq[log p(c|q; mu_Mc, sigma_Mc)])

Returns:
Return type:m_C, Sigma_C of probability distribution p_C of the current iteration
pyhrf.vbjde.vem_tools.expectation_C_ms(m_C, Sigma_C, G, H, m_A, W, XX, Gamma, Gamma_X, q_Z, mu_Mc, sigma_Mc, J, y_tilde, Sigma_G, sigma_eps_m, N, M, D, S)

Expectation-C step:

p_C = argmax_h(E_pa,pq,ph,pg[log p(a|y, h, a, g, q; theta)]) \propto exp(E_pa,ph,pg[log p(y|h, a, c, g; theta)] + E_pq[log p(c|q; mu_Mc, sigma_Mc)])

Returns:
Return type:m_C, Sigma_C of probability distribution p_C of the current iteration
pyhrf.vbjde.vem_tools.expectation_G_asl(Sigma_C, m_C, m_A, H, XX, W, WX, Gamma, Gamma_WX, XW_Gamma_WX, J, y_tilde, cov_noise, R_inv, sigmaG, prior_mean_term, prior_cov_term)

Expectation-G step:

p_G = argmax_g(E_pa,pc,ph[log p(g|y, a, c, h; theta)]) \propto exp(E_pa,pc,ph[log p(y|h, a, c, g; theta) + log p(g; sigmaG)])

Returns:
Return type:m_G, Sigma_G of probability distribution p_G of the current iteration
pyhrf.vbjde.vem_tools.expectation_G_ms(Sigma_C, m_C, m_A, H, XX, W, WX, Gamma, Gamma_WX, XW_Gamma_WX, J, y_tilde, cov_noise, R_inv, sigmaG, prior_mean_term, prior_cov_term, N, M, D, S)

Expectation-G step:

p_G = argmax_g(E_pa,pc,ph[log p(g|y, a, c, h; theta)]) \propto exp(E_pa,pc,ph[log p(y|h, a, c, g; theta) + log p(g; sigmaG)])

Returns:
Return type:m_G, Sigma_G of probability distribution p_G of the current iteration
pyhrf.vbjde.vem_tools.expectation_H_asl(Sigma_A, m_A, m_C, G, XX, W, Gamma, Gamma_X, X_Gamma_X, J, y_tilde, cov_noise, R_inv, sigmaH, prior_mean_term, prior_cov_term)

Expectation-H step:

p_H = argmax_h(E_pa,pc,pg[log p(h|y, a, c, g; theta)]) \propto exp(E_pa,pc,pg[log p(y|h, a, c, g; theta) + log p(h; sigmaH)])

Returns:
Return type:m_H, Sigma_H of probability distribution p_H of the current iteration
pyhrf.vbjde.vem_tools.expectation_H_ms(Sigma_A, m_A, m_C, G, XX, W, Gamma, Gamma_X, X_Gamma_X, J, y_tilde, cov_noise, R_inv, sigmaH, prior_mean_term, prior_cov_term, N, M, D, S)

Expectation-H step:

p_H = argmax_h(E_pa,pc,pg[log p(h|y, a, c, g; theta)]) \propto exp(E_pa,pc,pg[log p(y|h, a, c, g; theta) + log p(h; sigmaH)])

Returns:
Return type:m_H, Sigma_H of probability distribution p_H of the current iteration
pyhrf.vbjde.vem_tools.expectation_H_ms_concat(Sigma_A, m_A, m_C, G, XX, W, Gamma, Gamma_X, X_Gamma_X, J, y_tilde, cov_noise, R_inv, sigmaH, prior_mean_term, prior_cov_term, S)

Expectation-H step:

p_H = argmax_h(E_pa,pc,pg[log p(h|y, a, c, g; theta)]) \propto exp(E_pa,pc,pg[log p(y|h, a, c, g; theta) + log p(h; sigmaH)])

Returns:
Return type:m_H, Sigma_H of probability distribution p_H of the current iteration
pyhrf.vbjde.vem_tools.expectation_Ptilde_Likelihood(y_tilde, m_A, Sigma_A, H, Sigma_H, m_C, Sigma_C, G, Sigma_G, XX, W, sigma_eps, Gamma, J, D, M, N, Gamma_X, Gamma_WX)
pyhrf.vbjde.vem_tools.expectation_Q_asl(Sigma_A, m_A, Sigma_C, m_C, sigma_Ma, mu_Ma, sigma_Mc, mu_Mc, Beta, p_q_t, p_Q, neighbours_indexes, graph, M, J, K)
pyhrf.vbjde.vem_tools.expectation_Q_async_asl(Sigma_A, m_A, Sigma_C, m_C, sigma_Ma, mu_Ma, sigma_Mc, mu_Mc, Beta, p_q_t, p_Q, neighbours_indexes, graph, M, J, K)
pyhrf.vbjde.vem_tools.expectation_Q_ms(Sigma_A, m_A, Sigma_C, m_C, sigma_Ma, mu_Ma, sigma_Mc, mu_Mc, Beta, p_q_t, p_Q, neighbours_indexes, graph, M, J, K, S)
pyhrf.vbjde.vem_tools.expectation_ptilde_hrf(hrf_mean, hrf_covar, sigma_h, hrf_regu_prior, hrf_regu_prior_inv, hrf_len)

Expectation with respect to p_tilde hrf.

\mathrm{E}_{\widetilde{p}_{h}}\left[ \log p(h | \sigma_{h}) \right] = -\frac{D+1}{2}\log 2\pi -
\frac{D-1}{2}\log \sigma_{h} - \frac{\log \left| \mathbf{R} \right|}{2} - \frac{m^{t}_{h}\mathbf{R}^{-1}m_{h}
+ \mathrm{tr} \left( \Sigma_{h} \mathbf{R}^{-1} \right)}{2 \sigma_{h}}

pyhrf.vbjde.vem_tools.expectation_ptilde_labels(labels_proba, neighbours_indexes, beta, nb_conditions, nb_classes)

Expectation with respect to p_tilde q (or z).

\mathrm{E}_{\widetilde{p}_{q}} \left[ \log p (q | \beta ) \right] = \sum\limits_{m} & \left\{
  - \sum\limits_{j} \left\{ \log \left( \sum\limits^{1}_{i=0} \exp \left(
  \beta^{m} \sum\limits_{k \in N(j)} \widetilde{p}_{q^{m}_{k}} (i) \right)
  \right) \right\} \right. \\
  & \left. - \beta^{m} \sum\limits_{j}\sum\limits_{k \in N(j)}\sum\limits^{1}_{i=0}
  \left[ p^{MF}_{j}(i) \left( \frac{p^{MF}_{k}(i)}{2} -
  \widetilde{p}_{q^{m}_{k}} (i) \right) -
  \frac{1}{2}\widetilde{p}_{q^{m}_{j}}(i)\widetilde{p}_{q^{m}_{k}}(i) \right]
  \right\}

pyhrf.vbjde.vem_tools.expectation_ptilde_likelihood(data_drift, nrls_mean, nrls_covar, hrf_mean, hrf_covar, occurence_matrix, noise_var, noise_struct, nb_voxels, nb_scans)

Expectation with respect to likelihood.

\mathrm{E}_{\widetilde{p}_{a}\widetilde{p}_{h}\widetilde{p}_{q}}\left[\log p(y | a,h,q; \theta) \right] =
-\frac{NJ}{2} \log 2\pi + \frac{J}{2}\log\left| \Lambda_{j} \right| - N\sum\limits_{j \in J}\log
v_{b_{j}} + \frac{1}{2v_{b_{j}}}\sum\limits_{j \in J}V_{j}

where

V_{j} = \widetilde{m}^{t}_{a_{j}}\mathbf{X}^{t}_{h}\Lambda_{j}\mathbf{X}_{h}\widetilde{m}_{a_{j}}
+ \mathrm{tr}\left( \Sigma_{a_{j}}\mathbf{X}^{t}_{h}\Lambda_{j}\mathbf{X}_h \right) -
2\widetilde{m}^{t}_{a_{j}} \mathbf{X}^{t}_{h}\Lambda_{j}\left( y_{j} - \mathbf{P}\ell_{j} \right)

Parameters:
  • data_drift (ndarray, shape (nb_scans, nb_voxels)) – This is the BOLD data minus the drifts (y_tilde in the paper)
  • nrls_mean (ndarray, shape (nb_voxels, nb_conditions)) –
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels)) –
  • hrf_mean (ndarray, shape (hrf_len,)) –
  • hrf_covar (ndarray, shape (hrf_len, hrf_len)) –
  • occurence_matrix (ndarray, shape (nb_conditions, nb_scans, hrf_len)) –
  • noise_var (ndarray, shape (nb_voxels,)) –
  • noise_struct (ndarray, shape (nb_scans, nb_scans)) –
  • nb_voxels (int) –
  • nb_scans (int) –
Returns:

ptilde_likelihood

Return type:

float

pyhrf.vbjde.vem_tools.expectation_ptilde_nrls(labels_proba, nrls_class_mean, nrls_class_var, nrls_mean, nrls_covar)

Expectation with respect to p_tilde a.

\mathrm{E}_{\widetilde{p}_{a}\widetilde{p}_{q}}[\log p (a | q, \theta_{a})] = \sum\limits_{m}\sum\limits_{j}
& \left\{ \left[1 - \widetilde{p}_{q^{m}_{j}}(1) \right] \left[\log\frac{1}{\sqrt{2\pi\sigma^{2m}_{0}}} -
\frac{\left(m_{a^{m}_{j}} - \mu^{m}_{0} \right)^{2} + \Sigma_{a^{m,m}_{j}}}{2\sigma^{2m}_{0}} \right] +
\right. \\ &  \left. \widetilde{p}_{q^{m}_{j}}(1)  \left[\log\frac{1}{\sqrt{2\pi\sigma^{2m}_{1}}} -
\frac{\left(m_{a^{m}_{j}} - \mu^{m}_{1} \right)^{2} + \Sigma_{a^{m,m}_{j}}}{2\sigma^{2m}_{1}} \right] \right\}

pyhrf.vbjde.vem_tools.fit_hrf_two_gammas(hrf_mean, dt, duration)

Fits the estimated HRF to the standard two gammas model.

Parameters:
  • dt (float) –
  • duration (float) –
  • hrf_mean (ndarray, shape (hrf_len,)) –
Returns:

  • delay_of_response (float)
  • delay_of_undershoot (float)
  • dispersion_of_response (float)
  • dispersion_of_undershoot (float)
  • ratio_resp_under (float)
  • delay (float)

pyhrf.vbjde.vem_tools.free_energy_computation(nrls_mean, nrls_covar, hrf_mean, hrf_covar, hrf_len, labels_proba, data_drift, occurence_matrix, noise_var, noise_struct, nb_conditions, nb_voxels, nb_scans, nb_classes, nrls_class_mean, nrls_class_var, neighbours_indexes, beta, sigma_h, hrf_regu_prior, hrf_regu_prior_inv, gamma, hrf_hyperprior)

Compute the free energy functional.

\mathcal{F}(q, \theta) = \mathrm{E}_{q}\left[ \log p(y, A, H, Z ; \theta) \right] +  \mathcal{G}(q)

where E_{q}[\cdot] denotes the expectation with respect to q and \mathcal{G}(q) is the entropy of q.

Returns:free_energy
Return type:float
pyhrf.vbjde.vem_tools.fun(Beta, p_Q, Qtilde_sumneighbour, neighboursIndexes, gamma)

function to minimize

pyhrf.vbjde.vem_tools.grad_fun(Beta, p_Q, Qtilde_sumneighbour, neighboursIndexes, gamma)

function to minimize

pyhrf.vbjde.vem_tools.hrf_entropy(hrf_covar, hrf_len)

Compute the entropy of the hemodynamic response function. The entropy of a multivariate normal distribution is

\ln\left( \sqrt{(2\pi e)^{n} \left|\Sigma \right|} \right)

where n is the dimensionality of the vector space and \left|\Sigma \right| is the determinant of the covariance matrix.

Parameters:
  • hrf_covar (ndarray, shape (hrf_len, hrf_len)) – Covariance matrix of the HRF
  • hrf_len (int) – size of the HRF
Returns:

entropy

Return type:

float

pyhrf.vbjde.vem_tools.hrf_expectation(nrls_covar, nrls_mean, occurence_matrix, noise_struct, hrf_regu_prior_inv, sigmaH, nb_voxels, y_tilde, noise_var, prior_mean_term=0.0, prior_cov_term=0.0)

Computes the VE-H step of the JDE-VEM algorithm.

m^{(r)}_{H}  &= \Sigma^{(r)}_{H} \left( \sum\limits_{i \in \mathcal{P}} \widetilde{S}_{i}^{t}
\widetilde{y}^{(r)}_{i} \right) \\

\left(\Sigma^{(r)}_{H}\right)^{-1} &= \frac{\mathbf{R}^{-1}}{\sigma^{2(r)}_{h}} + \sum\limits_{i\in
\mathcal{P}} \left( \sum\limits_{m,m'} \sigma^{(r-1)}_{A_{mi}A_{m'i}} \mathbf{X}^{t}_{m} \Gamma^{(r)}_{i}
\mathbf{X}_{m'} + \widetilde{S}_{i}^{t} \Gamma^{(r)}_{i} \widetilde{S}_{i}\right)

where

\widetilde{S}_{i} &= \sum\limits^{M}_{m=1} m^{(r-1)}_{A_{mi}}\mathbf{X}_{m} \\

\widetilde{y}^{(r)}_{i} &= \Gamma^{(r)}_{i} \left( y_{i} - \mathbf{P}\ell^{(r)}_{i} \right)

Here, m^{(r-1)}_{A_{mi}} and \sigma^{(r-1)}_{A_{mi}A_{m'i}} denote the m^{th} and (m,m')^{th} entries of the mean vector and covariance matrix of the current q^{(r-1)}_{A_{i}}, respectively.

Parameters:
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels)) –
  • nrls_mean (ndarray, shape (nb_voxels, nb_conditions)) –
  • occurence_matrix (ndarray, shape (nb_conditions, nb_scans, hrf_len)) –
  • noise_struct (ndarray, shape (nb_scans, nb_scans)) –
  • hrf_regu_prior_inv (ndarray, shape (hrf_len, hrf_len)) – inverse of the hrf regularization prior matrix R
  • sigmaH (float) –
  • nb_voxels (int) –
  • y_tilde (ndarray, shape (nb_scans, nb_voxels)) –
  • noise_var (ndarray, shape (nb_voxels,)) –
  • prior_mean_term (float, optional) –
  • prior_cov_term (float, optional) –
Returns:

  • hrf_mean (ndarray, shape (hrf_len,))
  • hrf_covar (ndarray, shape (hrf_len, hrf_len))

pyhrf.vbjde.vem_tools.labels_entropy(labels_proba)

Compute the labels entropy.

-\sum\limits^{1}_{x=0} p_{Q}(x) \log p_{Q}(x)

Parameters:labels_proba (ndarray, shape (nb_conditions, nb_classes, nb_voxels)) – Probability of each voxel to be in one class
Returns:entropy
Return type:float
pyhrf.vbjde.vem_tools.labels_expectation(nrls_covar, nrls_mean, nrls_class_var, nrls_class_mean, beta, labels_proba, neighbours_indexes, nb_conditions, nb_classes, nb_voxels=None, parallel=True, nans_init=False)

Computes the E-Z (or E-Q) step of the JDE-VEM algorithm.

Using the mean-field approximation, \widetilde{p}^{(r)}_{Q^{m}}(q^{m}) is approximated by a factorized density \widetilde{p}^{(r)}_{Q^{m}}(q^{m})=
\prod\limits_{j \in \mathcal{P}_{\gamma}}\widetilde{p}^{(r)}_{Q^{m}_{j}}(q^{m}_{j}) such that if q^{m}_{j} = i, then \widetilde{p}^{(r)}_{Q^{m}_{j}}(i) \propto \mathcal{N} \left(m^{(r)}_{A^{m}_{j}};
\mu^{(r-1)}_{im}, v^{(r-1)}_{im} \right) f\left(q^{m}_{j} = i \,|\, \widetilde{q}^{m}_{\sim j};
\beta^{(r-1)}_{m}, v^{(r-1)}_{m} \right) where \widetilde{q}^{m} is a particular configuration of q^{m} updated at each iteration according to a specific scheme and

f\left( q^{m}_{j} \,|\, \widetilde{q}^{m}_{\sim j}; \beta^{(r-1)}_{m}, v^{(r-1)}_{m} \right) \propto
\exp\left(\alpha^{m(r)}_{j}(q^{m}_{j}) + \beta^{(r-1)}_{m}\sum\limits_{k \sim j}
I\left(\widetilde{q}^{m}_{k} = q^{m}_{j}\right) \right)

where

\left\{ \alpha^{m(r)}_{j} = \left( -v^{(r)}_{A_{j}^{m}A^{m''}_{j}} \left[ \frac{1}{v^{(r-1)}_{0m}},
\frac{1}{v^{(r-1)}_{1m}} \right] \right)^{t}, j \in \mathcal{P}_{\gamma}\right\}

and v^{(r)}_{A_{j}^{m}A^{m''}_{j}} denotes the (m,m') entries of the covariance matrix \Sigma^{(r)}_{A_{j}}

Notes

The mean-field fixed point equation is defined in:

Celeux, G., Forbes, F., & Peyrard, N. (2003). EM procedures using mean field-like approximations for Markov model-based image segmentation. Pattern Recognition, 36(1), 131–144. https://doi.org/10.1016/S0031-3203(02)00027-4
Parameters:
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels)) –
  • nrls_mean (ndarray, shape (nb_voxels, nb_conditions)) –
  • nrls_class_var (ndarray, shape (nb_conditions, nb_classes)) –
  • nrls_class_mean (ndarray, shape (nb_conditions, nb_classes)) –
  • beta (ndarray, shape) –
  • labels_proba (ndarray, shape (nb_conditions, nb_classes, nb_voxels)) –
  • neighbours_indexes (ndarray, shape (nb_voxels, max(len(a) for a in graph))) – This is the version of graph array where arrays from graph smaller than the maximum ones are filled with -1
  • nb_conditions (int) –
  • nb_classes (int) –
  • nb_voxels (int) –
  • parallel (bool) –
  • nans_init (bool) –
Returns:

labels_proba

Return type:

ndarray, shape (nb_conditions, nb_classes, nb_voxels)

pyhrf.vbjde.vem_tools.maximization_LA_asl(Y, m_A, m_C, XX, WP, W, WP_Gamma_WP, H, G, Gamma)
pyhrf.vbjde.vem_tools.maximization_Mu_asl(H, G, matrix_covH, matrix_covG, sigmaH, sigmaG, sigmaMu, Omega, R_inv)
pyhrf.vbjde.vem_tools.maximization_beta_m2_asl(beta, p_Q, Qtilde_sumneighbour, Qtilde, neighboursIndexes, maxNeighbours, gamma, MaxItGrad, gradientStep)
pyhrf.vbjde.vem_tools.maximization_beta_m2_scipy_asl(Beta, p_Q, Qtilde_sumneighbour, Qtilde, neighboursIndexes, maxNeighbours, gamma, MaxItGrad, gradientStep)

Maximize beta

pyhrf.vbjde.vem_tools.maximization_beta_m4_asl(beta, p_Q, Qtilde_sumneighbour, Qtilde, neighboursIndexes, maxNeighbours, gamma, MaxItGrad, gradientStep)
pyhrf.vbjde.vem_tools.maximization_class_proba(labels_proba, nrls_mean, nrls_covar)

Computes the M-(mu, sigma) step of the JDE-VEM algorithm.

\bar{q}^{(r)}_{mk}   & = \sum \limits_{i \in \mathcal{P}} q^{(r)}_{Z_{mi}} (k) \\
\mu^{(r+1)}_{mk}     & = \frac{\sum \limits_{i \in \mathcal{P}} q^{(r)}_{Z_{mi}} (k) m^{(r)}_{A_{mi}}}{\bar{q}^{(r)}_{mk}} \\
\sigma^{2(r+1)}_{mk} &= \frac{\sum \limits_{i \in \mathcal{P}} q^{(r)}_{Z_{mi}} (k) \left(\left(m^{(r)}_{A_{mi}} - \mu^{(r+1)}_{mk}\right)^{2} + \sigma^{(r)}_{A_{im}A_{im}} \right)}{\bar{q}^{(r)}_{mk}}

pyhrf.vbjde.vem_tools.maximization_drift_coeffs(data, nrls_mean, occurence_matrix, hrf_mean, noise_struct, drift_basis)

Computes the M-(l, Gamma) step of the JDE-VEM algorithm. In the AR(1) case:

\ell^{(r)}_{j} = \left( \bm{P}^{\intercal} \bm{\Lambda}^{(r)}_{j} \bm{P} \right)^{-1} \bm{P}^{\intercal} \bm{\Lambda}^{(r)}_{j} \left( y_{j} - \bm{\widetilde{S}}_{j} m^{(r)}_{H} \right)

pyhrf.vbjde.vem_tools.maximization_mu_sigma_asl(q_Z, m_X, Sigma_X)
pyhrf.vbjde.vem_tools.maximization_mu_sigma_ms(q_Z, m_X, Sigma_X, M, J, S, K)
pyhrf.vbjde.vem_tools.maximization_noise_var(occurence_matrix, hrf_mean, hrf_covar, nrls_mean, nrls_covar, noise_struct, data_drift, nb_scans)

Computes the M-sigma_epsilon step of the JDE-VEM algorithm.

\sigma^{2(r)}_{j} = \frac{1}{N} \left(\mathrm{E}_{\widetilde{p}^{(r)}_{A_{j}}} \left[a^{t}_{j}
\widetilde{\Lambda}^{(r)}_{j} a_{j} \right] - 2 \left(m^{(r)}_{A_{j}}\right)^{t} \widetilde{G}^{(r)}_{j}
y^{(r)}_{j} + \left(y^{(r)}_{j}\right)^{t} \Lambda^{(r)}_{j}y^{(r)}_{j} \right)

where matrix \widetilde{\Lambda}^{(r)}_{j} = \mathrm{E}_{\widetilde{p}^{(r)}_{H}}
\left[ G^{t}\Lambda^{(r)}_{j}G \right] is a M \times M whose element (m, m') is given by

\widetilde{g}^{t}_{m}\Lambda^{(r)}_{j}\widetilde{g}_{m'} + \mathrm{tr}\left(\Lambda^{(r)}_{j} X_{m}
\Sigma^{(r)}_{H} X^{t}_{m'} \right)

pyhrf.vbjde.vem_tools.maximization_sigmaH(D, Sigma_H, R, m_H)

Computes the M-sigma_h step of the JDE-VEM algorithm.

\sigma^{2(r+1)}_{h} = \frac{\mathrm{tr} ((\Sigma^{(r)}_{H} + m^{(r)}_{H} (m^{(r)}_{H})^{\intercal})\mathbf{R}^{-1})}{D-1}

pyhrf.vbjde.vem_tools.maximization_sigmaH_prior(D, Sigma_H, R, m_H, gamma_h)

Computes the M-sigma_h step of the JDE-VEM algorithm with a prior.

\sigma_{h}^{(r)} = \frac{(D-1) + \sqrt{8 \lambda_{\sigma_{h}} C + (D-1)^{2}}}{4\lambda_{\sigma_{h}}}

where

C = \mathrm{tr}\left( \left( \Sigma^{(r)}_{H} + m^{(r)}_{H} (m^{(r)}_{H})^{\intercal} \right) \mathbf{R}^{-1} \right)

pyhrf.vbjde.vem_tools.maximization_sigma_asl(D, Sigma_H, R_inv, m_H, use_hyp, gamma_h)
pyhrf.vbjde.vem_tools.maximization_sigma_noise_asl(XX, m_A, Sigma_A, H, m_C, Sigma_C, G, Sigma_H, Sigma_G, W, y_tilde, Gamma, Gamma_X, Gamma_WX, N)

Maximization sigma_noise

pyhrf.vbjde.vem_tools.maximum(iterable)

Return the maximum and the indice of the maximum of an iterable.

Parameters:iterable (iterable or numpy array) –
Returns:
  • iter_max (the maximum)
  • iter_max_indice (the indice of the maximum)
pyhrf.vbjde.vem_tools.mult(v1, v2)

Multiply two vectors.

The first vector is made vertical and the second one horizontal. The result will be a matrix of size len(v1), len(v2).

Parameters:
  • v1 (ndarray) – unidimensional
  • v2 (ndarray) – unidimensional
Returns:

x

Return type:

ndarray, shape (len(v1), len(v2))

pyhrf.vbjde.vem_tools.norm1_constraint(function, variance)

Returns the function constrained with optimization strategy.

Parameters:
  • function (array_like) – function to optimize under norm1 constraint
  • variance (array_like) – variance of the function, must be the same size
Returns:

optimized_function

Return type:

numpy array

Raises:

ValueError – If len(variance) != len(function)

pyhrf.vbjde.vem_tools.normpdf(x, mu, sigma)
pyhrf.vbjde.vem_tools.nrls_entropy(nrls_covar, nb_conditions)

Compute the entropy of neural response levels. The entropy of a multivariate normal distribution is

\ln\left( \sqrt{(2\pi e)^{n} \left|\Sigma \right|} \right)

where n is the dimensionality of the vector space and \left|\Sigma \right| is the determinant of the covariance matrix.

Parameters:
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels)) – Covariance of the NRLs
  • nb_conditions (int) –
Returns:

entropy

Return type:

float

pyhrf.vbjde.vem_tools.nrls_expectation(hrf_mean, nrls_mean, occurence_matrix, noise_struct, labels_proba, nrls_class_mean, nrls_class_var, nb_conditions, y_tilde, nrls_covar, hrf_covar, noise_var)

Computes the VE-A step of the JDE-VEM algorithm.

\Sigma^{(r)}_{A_{j}} &= \left( \sum\limits^{I}_{i=1} \Delta_{ij} + \widetilde{H}_{j} \right)^{-1} \\
m^{(r)}_{A_{j}} &= \Sigma^{(r)}_{A{j}} \left( \sum\limits^{I}_{i=1} \Delta_{ij} \mu_{i}^{(r-1)} + \widetilde{G}^{t}\Gamma_{j}^{(r-1)}\left(y_{j} - P\ell^{(r-1)} \right)\right)

where:

The mth column of \widetilde{G} is denote by \widetilde{g}_{m}  = X_{m}m^{(r)}_{H} \in \mathbb{R}^{N}

\Delta_{ij} &= \mathrm{diag}_{M}\left[\frac{q^{(r-1)}_{Z_{mj}}(i)}{\sigma^{2(r)}_{mi}}\right] \\

\widetilde{H}_{j} &= \widetilde{g}^{t}_{m} \Gamma^{(r-1)}_{j} \widetilde{g}_{m'} + \mathrm{tr}\left( \Gamma^{(r-1)}_{j}X_{m} \Sigma^{(r)}_{H}X^{t}_{m'} \right)

Parameters:
  • hrf_mean (ndarray, shape (hrf_len,)) –
  • nrls_mean (ndarray, shape (nb_voxels, nb_conditions)) –
  • occurence_matrix (ndarray, shape (nb_conditions, nb_scans, hrf_len)) –
  • noise_struct (ndarray, shape (nb_scans, nb_scans)) –
  • labels_proba (ndarray, shape (nb_conditions, nb_classes, nb_voxels)) –
  • nrls_class_mean (ndarray, shape (nb_conditions, nb_classes)) –
  • nrls_class_var (ndarray, shape (nb_conditions, nb_classes)) –
  • nb_conditions (int) –
  • y_tilde (ndarray, shape (nb_scans, nb_voxels)) – BOLD data minus drifts
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels)) –
  • hrf_covar (ndarray, shape (hrf_len, hrf_len)) –
  • noise_var (ndarray, shape (nb_voxels,)) –
Returns:

  • nrls_mean (ndarray, shape (nb_voxels, nb_conditions))
  • nrls_covar (ndarray, shape (nb_conditions, nb_conditions, nb_voxels))

pyhrf.vbjde.vem_tools.plot_convergence(ni, M, cA, cC, cH, cG, cAH, cCG, SUM_q_Z, mua1, muc1, FE)
pyhrf.vbjde.vem_tools.plot_response_functions_it(ni, NitMin, M, H, G, Mu=None, prior=None)
pyhrf.vbjde.vem_tools.polyFit(signal, tr, order, p)
pyhrf.vbjde.vem_tools.poly_drifts_basis(nb_scans, param_lfd, tr)

Build polynomial drifts basis.

Parameters:
  • nb_scans (int) –
  • param_lfd (int) – TODO
  • tr (float) –
Returns:

drifts_basis – K is determined by the scipy.linalg.orth function and corresponds to the effective rank of the matrix it is applied to (see function’s docstring)

Return type:

ndarray, shape (nb_scans, K)

pyhrf.vbjde.vem_tools.ppm_contrasts(contrasts_mean, contrasts_var, contrasts_class_mean, contrasts_class_var, threshold_a='std_inact', threshold_g=0.95)

Computes the ppm for the given contrast using either the standard deviation of the “all inactive conditions” class gaussian (default) or the intersection of the [all inactive conditions] and [all active conditions] classes gaussians as threshold for the PPM_a and 0.95 (default) for the PPM_g. Be carefull, this computation considers the mean of the inactive class as zero.

Parameters:
  • contrasts_mean (ndarray, shape (nb_voxels, nb_contrasts)) –
  • contrasts_var (ndarray, shape (nb_voxels, nb_contrasts)) –
  • contrasts_class_mean (ndarray, shape (nb_contrasts, nb_classes)) –
  • contrasts_class_var (ndarray, shape (nb_contrasts, nb_classes)) –
  • threshold_a (str, optional) – if “std_inact” (default) uses the standard deviation of the [all inactive conditions] gaussian class as PPM_a threshold, if “intersect” uses the intersection of the [all inactive/all active conditions] gaussian classes
  • threshold_g (float, optional) – the threshold of the PPM_g
Returns:

  • ppm_a_contrasts (ndarray, shape (nb_voxels, nb_contrasts))
  • ppm_g_contrasts (ndarray, shape (nb_voxels, nb_contrasts))

pyhrf.vbjde.vem_tools.ppms_computation(elements_mean, elements_var, class_mean, class_var, threshold_a='std_inact', threshold_g=0.9)

Considering the elements_mean and elements_var from a gaussian distribution, commutes the posterior probability maps considering for the alpha threshold, either the standard deviation of the [all inactive conditions] gaussian class or the intersection of the [all (in)active conditions] gaussian classes; and for the gamma threshold 0.9 (default).

The posterior probability maps (PPM) per experimental condition is computed as p(a^{m}_{j} > \delta | y_{j}) > \alpha. Note that we have to thresholds to set. We set \delta to get a posterior probability distribution, and \alpha is the threshold that we set to see a certain level of significance. As default, we chose a threshold \delta for each experimental condition m as the intersection of the two Gaussian densities of the Gaussian Mixture Model (GMM) that represent active and non-active voxel.

\frac{(\delta - \mu^{m}_{1})^{2}}{v^{m}_{1}} - \frac{(\delta - \mu^{m}_{0})^{2}}{v^{m}_{0}} =
\log\left( \frac{v^{m}_{1}}{v^{m}_{0}} \right)

\mu^{m}_{i} and v^{m}_{i} being the parameters of the GMM in a^{m}_{j} corresponding to active (i=0) and non-active (i=1) voxels for experimental condition m.

\delta = \frac{\mu_{1}\sigma_{0}^{2} \pm \sigma_{0}\sigma_{1}\sqrt{\mu^{2}_{1} + 2 \left(\sigma^{2}_{0} -
\sigma^{2}_{1} \right) \log \left( \frac{\sigma_{0}}{\sigma_{1}} \right)}}{\sigma^{2}_{0} - \sigma^{2}_{1}}

Be careful, this computation considers the mean of the inactive class as zero.

Notes

nb_elements refers either to the number of contrasts (for the PPMs contrasts computation) or for the number of conditions (for the PPMs nrls computation).

Parameters:
  • elements_mean (ndarray, shape (nb_voxels, nb_elements)) –
  • elements_var (ndarray, shape (nb_voxels, nb_elements)) –
  • class_mean (ndarray, shape (nb_elements, nb_classes)) –
  • class_var (ndarray, shape (nb_elements, nb_classes)) –
  • threshold_a (str, optional) – if “std_inact” (default) uses the standard deviation of the [all inactive conditions] gaussian class as PPM_a threshold. If “intersect” uses the intersection of the [all inactive/all active conditions] gaussian classes.
  • threshold_g (float, optional) – the threshold of the PPM_g
Returns:

  • ppm_a (ndarray, shape (nb_voxels, nb_elements))
  • ppm_g (ndarray, shape (nb_voxels, nb_elements))

pyhrf.vbjde.vem_tools.roc_curve(dvals, labels, rocN=None, normalize=True)

Compute ROC curve coordinates and area

  • dvals - a list with the decision values of the classifier
  • labels - list with class labels, in {0, 1}

returns (FP coordinates, TP coordinates, AUC )

pyhrf.vbjde.vem_tools.sum_over_neighbours(neighbours_indexes, array_to_sum)

Sums the array_to_sum over the neighbours in the graph.