pyKCSD - implementation step by step¶

Implementation of 1-dimensional kCSD method in Python. This is a brutal port with minimal changes to the existing version.

Data¶

Let's work with two gaussian sources - generate CSD and calculate potentials.

Potentials measured by an electrode at the height H from the tissue:

$$\phi(x_e) = \frac{1}{4\pi\sigma}\int_{-\infty}^{\infty} \frac{I(x) dx}{\sqrt{(x-x_e)^2 + H^2}}$$

Sampling¶

Lets now insert 10 uniformly distributed electrodes and plot their readings.

If we were now to calculate CSD using finite difference approximation, we would get:

$$c \approx -\frac{\phi_{i-1} + \phi_{i+1} - 2\phi_{i}}{(x_i - x_{i-1})^2}$$

Now time for kCSD¶

Travel through the existing Matlab code with partial Python adaptation.

Firstly we are setting some basic parameters.

Calculating contributions of sources at different distances.

We are integrating such a function:

$$f(x)= \frac{1}{2 \sigma}(\sqrt{(arg - \text{curr_pos})^2+R^2}-\left|{arg-\text{curr_pos}}\right|) \cdot \text{gauss_rescale(src, curr_pos, h)}$$

where Gauss rescale $g(x,\mu,3\sigma)$:

$$\sigma_{n2} = \left(\frac{3\sigma}{3}\right)^2$$

$$g(x) = \frac{1}{\sqrt{2\pi\sigma_{n2}}} \cdot \exp\left( -\frac{(x-\mu)^2}{2 \sigma_{n2}} \cdot (\left| x-\mu \right|<3\sigma)\right)$$

It is important to state here that $\sigma$ in the $3 \sigma$ doesn't mean tissue conductivity. It is just a scaling constant.

$$y = \int_{src-4h}^{src+4h} \frac{1}{2 \sigma}(\sqrt{(x - \text{curr_pos})^2+R^2}-\left|{arg-\text{curr_pos}}\right|) \cdot \text{gauss_rescale(src, curr_pos, h)} \text{ dcurr_pos}$$

Using the functions above we are creating a table with the potential generated by a single source base elements as a function of distance.

We can adjust the number of points used to calculate the integrals in the dimensions where we assume the current distribution. Doing this we can balance between the accuracy and performance.

Note: this is definitely a sparse and potentialy large matrix.

Possibly some memory can be saved here, especially for bigger simulations.

Unfortunately this matrix is dense.

Now give me my pots!

Now its time for CSD

This is the result for $\lambda$ = 0. Not very impressive yet, but the direction is good.

Regularization¶

Now let's calculate the regularization parameter $\lambda$.

The $\lambda$ parameter measures the importance of solution smoothness against the importance of the training data. Higher $\lambda$ means that the smoothness is more important and the training data -- not so much.

def calc_err(lambd, Pot, K, Ind_test, Ind_train, norm_order, subset, subset_size): nr_of_frames = size(Pot,2); if subset == 'all': Time_ind = np.array(range(0,nr_of_frames-1)) if subset == 'part': Time_ind = randperm(nr_of_frames,subset_size); if subset == 'one':
diffs = max(Pot)-min(Pot); Time_ind = find(diffs==max(diffs),1);

K_train = K(Ind_train, Ind_train);
Pot_train = Pot(Ind_train,Time_ind);
Pot_test = Pot(Ind_test,Time_ind);
I = identity(len(Ind_train))
beta = linalg.inv(K_train + lambd*I)*Pot_train;
K_cross = K(Ind_test, Ind_train);
Pot_est = K_cross * beta;
err = norm(Pot_test - Pot_est, norm_order);
return err

def get_choose_lambda_parameters(args): propertyArgIn = args;

while len(propertyArgIn) >= 2:
    prop = propertyArgIn[0];
    val = propertyArgIn[1];
    propertyArgIn = propertyArgIn[2:];

    if prop == 'n_folds':
        n_folds = val;
    elif prop == 'n_iter':
        n_iter = val;
    elif prop == 'sampling':
        sampling = val;
    else:
        print 'no method defined for input: ', prop
        print 'Available inputs: n_folds, n_iter, sampling, choose_R'

if n_folds is None:
    #[n_el, ] = size(pots);
    #n_folds = n_el;
    n_folds = size(pots)
if n_iter is None:
    n_iter = 1;
if sampling is None:
    sampling = 1

return [n_folds, n_iter, sampling]

def cross_validation(lambd, pot, K, n_folds, Ind_perm, norm_order, subset, subset_size):
""" K - fold cross validation preparing subsets of electrodes """

nEl = len(pot);
width = int(floor(nEl/n_folds))
Ind = [0 for i in xrange(n_folds)]
All_ind = range(0,nEl-1)

for i in xrange(0, n_folds-1):
    Ind_set = np.array(range(0,int(width))) + width*(i-1) #check this out in Matlab
    print 'Ind_perm:', Ind_perm
    print 'Ind_set', Ind_set
    #Ind_set_perm = [Ind_perm[z] for z in Ind_set]
    Ind_set_perm = Ind_perm.copy()
    Ind[i] = Ind_set_perm

last = max(Ind_set)

if last < nEl:
    Ind[n_folds-1] = [Ind[n_folds-1], Ind_perm[last:nEl-1]]

errors = np.zeros(n_folds);

for i in xrange(0,n_folds-1):
    Ind_test = Ind[i];
    Ind_train = All_ind;
    Ind_train[Ind[i]] =[];
    errors[i] = calc_err(lambd, Pot, K, Ind_test, Ind_train, norm_order, subset, subset_size);

err = mean(errors);
return err

def calcCvError(varargin): """A cross-validation estimator of the error of the estimation. Used in methods that choose parameters through cross-validation. """ if varargin==None: n_folds = len(elec_pos); else: n_folds = varargin[0];

Ind_perm = np.random.permutation(range(len(elec_pos)))
err = cross_validation(lambd, pots, KPot, n_folds, Ind_perm, norm_order, subset, subset_size)

return err

def lambda_sampling_1(n_folds, n_iter): n = len(lambdas); errors = zeros(n); errors_iter = zeros(n_iter);

for i in xrange(n-1):
    lambd = lambdas[i];
    for j in xrange(0,n_iter-1):
        errors_iter[j] = calcCvError([n_folds])
    errors[i] = mean(errors_iter)
value = lambdas[np.where(errors == min(errors))[0][0]]
return value

def chooseLambda(varargin): """ Chooses the regularisation lambda parameter for ridge regression. The user can enter options by providing 'property_name', property_value pairs:

    'n_folds'     number of folds to perform Cross validation (CV)
    'n_iter'      number of iterations for the CV procedure
    'sampling'    ways of looking for the optimal lambda:
                    1 - simple sampling
                    2 - using fminbnd function
"""
[n_folds, n_iter, sampling] = get_choose_lambda_parameters(varargin)
lambd = lambda_sampling_1(n_folds, n_iter)
return lambd

import sklearn from sklearn.cross_validation import KFold sklearn.linear_model import RidgeCV

def chooseLambda(): validator = RidgeCV(alphas=lambdas, fit_intercept=True, normalize=False, scoring=None, score_func=None, loss_func=None, cv=None, gcv_mode=None, store_cv_values=True) lambd = 0.0 return lambd

norm_order = 2 subset = 'part' subset_size = 500 lambd = chooseLambda(['n_folds', None, 'n_iter', 2, 'sampling', None]) print 'lambda = ', lambd csd = estimate_csd() plot(x,true_csd) plot(linspace(min(elec_pos), max(elec_pos), len(output)),csd) lambd = 0.0 csd = estimate_csd() plot(linspace(min(elec_pos), max(elec_pos), len(output)),csd) legend((r'kCSD estimated CSD, $\lambda$=4.0', r'true CSD',r'kCSD estimated CSD, $\lambda$=0.0'), shadow = True, loc = (0.01, 0.75))

import sklearn

from sklearn.cross_validation import KFold X = np.array([[0., 0.], [1., 1.], [-1., -1.], [2., 2.]]) Y = np.array([0, 1, 0, 1])

kf = KFold(len(Y), n_folds=2, indices=False) print(kf)

sklearn.cross_validation.KFold(n=4, n_folds=2)

for train, test in kf:

... print("%s %s" % (train, test))

TODO¶

• Split into statistical part (independent of dimensionality) and basis function management part (dependent on dimensionality)
• Currently used cross-validation: leave-one-out $\rightarrow$ room for improvement, other method may be used (maybe SciPy has something interesting to offer)