What is network_control.energies.minimum_energy_fast doing?

[JohannesWiesner]

I've implemented a function called state_to_state_transition in our nict repo and would like to make sure that I don't implement functions that have already been implemented by you. The purpose of this function is to either calculate optimal control energy or minimum control energy for all combinations of input brain states. The function outputs a TxNxM 3D-array where

where T is the number of time points, N is the number of nodes and M denotes each combination of states.

In other words, it's just a 3D array that stores multiple outputs from either minimum_input or optimal_input.

Since we are both interested in combinations (e.g. only AB but not BA) and permutations (e.g. AB but also BA), I've implemented this logic in nict._get_layout. (so the user can choose between 'combinations' and 'permutations'

Now I've discovered minimum_energy_fast but I can't really figure out what this function does? E.g what do you mean with numpy array (N x n_transitions) array? And what is meant by:

Returns:
  E: numpy array (N x n_transitions)
        Regional energy for all state transition pairs.
        Notes,
            np.sum(E, axis=0)
                collapse over regions to yield energy associated with all transitions.
            np.sum(E, axis=0).reshape(n_states, n_states)
                collapse over regions and reshape into a state by state transition matrix.

P.S.: Especially pinging @jastiso here because I promised to share the link to this repo.

[lindenmp]

Hi Johannes,

The minimum_energy_fast function simply estimates energy associated with all possible state transitions simultaneously, eschewing the need to process in serial. It does so by taking advantage of the fact that a few steps in the process need only be estimated once (e.g., computing and inverting the Gramian; see section titled Evaluating Minimum Control Energy at the bottom of this document).

The output E has regions on the rows and state transitions on the columns, so if you sum over regions and then reshape to a square matrix, you end up with a matrix of transition energies. Note, this is only true if you use our utils.expand_states function to define x0_mat and xf_mat (sorry, that's not so clear from the code!). See this notebook for an example.

Thanks, Linden

[JohannesWiesner]

Hi Linden, I never thanked you for the answer which I would like to do now (also thanks for the detailed explanation)! Over the holidays I reworked the nict package. It now contains a better separation between computing control energy for one subject vs. multiple subjects (nict.single_subject, nict.multi_subjects).

I have a few follow-up questions because I think one could speed up computations in the nict package based on what you have said? In my package, I have this state_to_state_transition function. It allows to compute all possible state-to-state-transitions based on an input A and a matrix of states X where the rows denote each state and the columns denote the nodes.

As you can see the function just uses a for-loop to iterate through all combinations of possible source and target states to either compute minimum or optimal control energy.

Couldn't one also use the same logic here as in minimum_energy_fast? And is there also an analog for computing optimal control energy (or is this form of mathematical speed-up only possible when computing minimum control energy)? Lastly minimum_energy_fast is constrained to binary states, but here in Mannheim we are working with with activation maps (beta-images / T-images, etc.). Would such a speed up also be working for dimensional states?

[lindenmp]

Hi Johannes,

No problem, we're glad you're finding the package useful!

Yes, one of the reasons we wanted to implement minimum_energy_fast was to achieve exactly what you're doing in your own state_to_state_transition function but faster! So I would definitely suggest swapping in our minimum_energy_fast to speed things up. Something like the following in place of your for-loop should work:

E = np.zeros((n_states, n_states))
e = minimum_energy_fast(A_norm, T, B, x0_mat, xf_mat)
E[:] = np.sum(e, axis=0).reshape(n_states, n_states)

There are two considerations: 1) You will only get e out, not x or u. 2) You need to define the matrix of state transitions ahead of time (x0_mat, xf_mat). This should work with non-binary states, but I haven't tried it myself. In any case, I think you'll have to make a few changes to your X matrix. We define x0_mat and xf_mat with nodes on rows and transitions on columns. Note, every single transition is stored in its own column, which means the actual states themselves are duplicated many times in each of x0_mat and xf_mat. For example, if the first 3 transitions are from state_0 to state_1, state_0 to state_2, and state_0 to state_3, then the first 3 columns of x0_mat are identical (all representing the node activity levels associated with state_0). So you'll have to duplicate your activation map-derived states.

Sorry, we won't have a fast version of optimal control at the moment.

Hope this helps!

[JohannesWiesner]

Hi Linden, thanks for the explanation, i've written the following test script to ensure that both minimum_input , nict.single_subject.state_to_state_transition and minimum_energy_fast produce the same error. Computing the error for one state-to-state transition using nict gives the same result as manually computing minimum_input using the same source and target state. (which makes sense because nict is simply using minimum_input under the hood). However, minimum_energy_fast does not give the same error (and they are surprisingly large?). Did I do something wrong or is this a bug?:

# -*- coding: utf-8 -*-
"""
Compare minimum_energy_fast from network control with manually computed
error using nict package

@author: Johannes.Wiesner
"""
import sys
sys.path.append('./nict')

from nict.single_subject import state_to_state_transition
from nict.single_subject import _set_transition_order
from network_control.utils import matrix_normalization
from network_control.energies import minimum_energy_fast
from network_control.energies import minimum_input
import numpy as np

# create artifical data
np.random.seed(28)
A = np.random.rand(10,10)
A_norm = matrix_normalization(A,c=1,version='continuous')
X = np.random.rand(4,10) # 4 (dimensional) sample states and 10 nodes
# X = np.random.choice(a=[0,1],size=(4,10),p=[0.5,0.5]) # 4 (boolean) sample states and 10 nodes
B = np.eye(10)
S = np.eye(10)
T = 10 # make T large enough to get big errors
version='continous'
n_states,n_nodes = X.shape

# 1.) first compute all state-to-state-transitions using johannes nict package 
X_out,U_out,E_out = state_to_state_transition(A=A_norm,T=T,B=B,X=X,rho=None,S=None,order='permutations')

# 2.) Now do sanity check by simply recomputing one randomly chosen state-to-state-transition
# and compare that with nict package
x,u,e = minimum_input(A=A_norm,T=T,B=B,x0=X[0,:],xf=X[1,:])
print(np.equal(E_out[0],e))

# 3.) now run minimum_energy_fast from network_control package and check if the
# computational error is the same
# create x0_mat and xf_mat using nict package (we could probably also use 
# expand states here but the outcome is the same)
n_transitions,transition_idxs = _set_transition_order(n_states=n_states,order='permutations')
x0_mat = np.zeros((n_nodes,n_transitions))
xf_mat = np.zeros((n_nodes,n_transitions))

for idx,src,tgt in transition_idxs:
    x0_mat[:,idx] = X[src,:]
    xf_mat[:,idx] = X[tgt,:]

# run minimum energy fast and sum error over nodes to get one value for each
# transition
e_fast = minimum_energy_fast(A=A_norm,T=T,B=B,x0_mat=x0_mat,xf_mat=xf_mat)
e_fast = np.sum(e_fast,axis=0)
print(np.equal(E_out[0],e_fast[0]))

Follow up question regarding:

You will only get e out, not x or u.

Is this something that could be implemented in the future or simply not possible? If it could be implemented I could open a feature request for this? :) Receiving u would definitely be something that I would need for my project

Sorry, we won't have a fast version of optimal control at the moment.

No need to say sorry! I am the one that should come up with improvements (I am just really sure that I can't...) Same as above, I would then open a feature request if it would be theoretically possible to create optimal_energy_fast?

EDIT:

I tried for both boolean and dimensional states: X = np.random.rand(4,10) # 4 (dimensional) sample states and 10 nodes or X = np.random.choice(a=[0,1],size=(4,10),p=[0.5,0.5]) # 4 (boolean) sample states and 10 nodes