kESI electrode potential fixing, fixes with Kuba, equal conductivity vs kCSD, linear solver benchmark

[mdovgialo]

Need to experiment and figure out what happening with MFEM correctly solved electrode potential (far boundary). Why the results of kESI are uniform...

TO use:

• far boundary
• far boundary + regrounding at ground electrode-
• far boundary + grounding electrode?
• regroundiung kCSD?

[mdovgialo]

For some parts it seems to work. When we use purely Numerical leadfield from simulated point sources at electrodes with far boundary.

For example 4 sphere model, point dipole in it, modelled using MFEM with far boundary and electrode in the brain as ground:

Biased kCSD solution:

Biased kESI solution:

But it seems like regularisation parameter selection is very unstable. More than that, analytical solution or unbiased solution makes kESI fall apart completely. Althout with manual lambda selection it can work kinda and have some sorts of sphere barely visible:

Unbiased kESI solution with manal ok lambda selection:

Unbiased kCSD solution to compare:

My suspicion is instability due to large numbers in leadfields?. Also np.linalg.solv takes ages for some reason for more than 200 channels.

[mdovgialo]

experimenting with replacing solvers

Trying to calculate good kesi or kcsd solution at sensible regularisation value, calculating leave one out crossvalidation speed.

Solving the kernel * x = potential (i think) equations for the x

We are changing the solver in this line

kESI/src/kesi/_engine.py

Line 103 in 17fffaf

def _leave_one_out_estimate(self, KERNEL, X, i, IDX):

kCSD

solver	speed	result	notes
numpy.linalg.solve	5 it/s	~1e-5	was used until now
scipy.linalg.solve assume_a gen	5 it/s	~1e-5	equal to numpy!
scipy.linalg.solve assume_a pos	500 it/s	~1e-5	equal to numpy!
scipy cholesky decomp	500 it/s	~1e-5	equal to numpy!
numpy.linalg.pinv * rhs	30 it/s	~1e-5	equal to numpy!
scipy.linalg.lstsq * rhs	30 it/s	~1e-5	equal to numpy!

kCSD linear equation system solving is stable across solvers

kESI

solver	speed	result	notes
numpy.linalg.solve	5 it/s	~1e-5	was used until now
scipy.linalg.solve assume_a gen	5 it/s	~1e-5	NOT EQUAL TO NUMPY, ILL CONDITIONED MATRIX warning
scipy.linalg.solve assume_a pos	400 it/s	~1e-5	NOT EQUAL TO NUMPY, ILL CONDITIONED MATRIX warning
scipy cholesky decomp	crash	N/A	Singular minor erorr, cannot decompose
numpy.linalg.pinv * rhs	30 it/s	~1e-5	NOT EQUAL TO NUMPY
scipy.linalg.lstsq * rhs	30 it/s	~1e-5	NOT EQUAL TO NUMPY

results review

We can look at the results of leave one out error, per electrode left out. In kCSD it looks like this:

In kESI:

255 electrodes, 100 samples of signal

inverse and least squares remove some part of the result (which makes it actually better, as it means lower error), and the whole image is quite noisy, unlike kCSD. Signal used in the simulation had no noise.

[mdovgialo]

kCSD solution using scipy.linalg.solve assume pos vs numpy.linalg.solve

on set of 198 electrodes:

numpy took 2 minutes, scipy, 30 seconds.

Set of 255 electrodes:

scipy takes 1 minute 26 seconds, numpy - 2 hours

Solutions crossection:

Seems like it's exactly the same, and we've got a tremendous speedup.

[danek8317]

Re: #38 (comment)

The cross-validation errors seem comparable for kCSD and kESI, going up to 0.075, which is good. The equation you mention is to compute the inverse $K$, which is needed to compute CSD. I would be interested in seeing $K$ and $\tilde{K}$ matrices for kCSD and kESI compared.

How to visualize them:

• as images, square for $K$, rectangle for $\tilde{K}$ for a modest subsection of reconstruction points,
• same as above but showing logarithm of absolute value of elements
• sorted by absolute value of matrix elements, probably on a log scale

The size of these matrices should be the same for kCSD and kESI. The images would help us understand differences

I would do the tests for small matrices first to get the gist without waiting for results. Even 40 channels gives 1600 elements of K..., so perhaps start with 20 electrodes? Just to identify the differences between the matrices. Do you have 1/N factors in the definition of $K$ and $\tilde{K}$?

[mdovgialo]

Visualising KERNEL from kernel constructor:

kCSD

kESI, far boundary

kESI, far boundary, regrounded at a point between FCz and Cz

it seems like some electrodes are quite big outliers, which could make the solution numerically unstable.

And it seems some Ox_L1 electrodes are outside of the 4 sheres. Maybe, when I was scaling their locations down in excel I managed to miscalculate and only scale their X, Z axis, without Y.... Will try again without them.

Todo: fix the layered electrodes...

[mdovgialo]

Comparing solution of electrodes:

Regrounded artificially in the grounding neck electrode:

Biggest difference is that forcing BC to be 0 in neck and infinity is overall smoother potential with smaller differences at material boundaries.

Biased

Sources only in the brain material

Ground truth - dipole in the brain, off center
Potentials at electrodes calculated using FEM with grounding at infinity (100 meters) and the "neck".

kernel matrices

kCSD

kESI

Electrode potentials calculated in 4 spheres with grounding at infinity

kESI regrounded

Electrode potentials calculated in 4 spheres with grounding at infinity, then potential at a point between FCz and Cz was subtracted for every leadfield

kESI electrode

Electrode potentials calculated in 4 spheres with grounding at infinity and a disk simulating grounding through the neck

CSD

Cross marks the dipole position

kCSD

kESI

kESI regrounded

kESI electrode

Seems to have the least artifacts? 6 orders of magnitude higher CSD than kCSD...

UnBiased

Sources in the whole domain

Ground truth - dipole in the brain, off center
Potentials at electrodes calculated using FEM with grounding at infinity (100 meters) and the "neck".

kernel matrices

kCSD

kESI

kESI regrounded

kESI electrode

CSD

kCSD

kESI

kESI regrounded

kESI electrode

Actually recreates the spherical domain even in ubniased mode

Conclusion?

Seems like solution with 0 potential at infinity (100 meters) and in the neck area makes the potential not explode to giant values, and makes the kernel matrix "behave" and look like kCSD kernel matrix. Just regrounding the potential also seems to do the trick but only in biased mode, which never samples the negative potential. Seems like having negative potential from positive monopolar makes the kernel solution misbehave...

[mdovgialo]

Analytical

Potentials at electrodes calculated using FEM with grounding at infinity (100 meters) and the "neck". Only 2 layers of electrodes instead of 3, because the inner one is below the dipole, and is unsampleable (for now)

Biased

kCSD

kESI with grounding electrode

Unbiased

kCSD

,

kESI with grounding electrode

[mdovgialo]

Trying to precalculate kCSD electrode potentials

I've taken a 0.5 and 1 mm grid meshes, of 4 spheres "MRI", treat them the same way, as we would treat real MRIs, but instead of MFEM solver use 1/4 pi sigma R "solver". Thus creating electrode correction vtk file, which we can use with kESI, to check, if we'll get the same results as kESI - kCSD.

Comparing with kESI regrounded and kESI with electrode

regrounded at z = 0.0795

Just the potential from 1/4 pi R and kESI with far boundary and grounding electrode:

Regrounded at z = 0.0795

kESI solutions using kCSD precalculated leadfields

kernels

Kernels look remarkably similar, the only difference seems to be scaling - kESI kCSD implementation has the matrix scaled by 1e6?

Let's zoom in into the most inner electrode layer

and on one "Line"

If we try to minimize np.sqrt(np.sum((kcsd['KERNEL']/10**i - prec['KERNEL'])**2))) - difference between kcsd kernel and precalculated kCSD kernel we get:

lowest error with kCSD kernel rescaled by 1e-8, sum error equal to 0.5661129733749277,

imshow(kcsd['KERNEL']*1e-8 - prec['KERNEL'])
seems to mostly invert the numbers, so it's not exactly 1e-8...

Somehow the best scale is:

np.sqrt(np.sum((kcsd['KERNEL']*x - prec['KERNEL'])**2))

x = 1.2222122212221224e-08
with error = 0.0008348384586058734

error looks like:

imshow(kcsd['KERNEL']*best_x - prec['KERNEL'])

kcsd['KERNEL']*best_x - prec['KERNEL'] distribution:

It's seems, like the kernel matrices match pretty closely except of the scale. 99.7% of the matrix elements are close enough to each other at 0.1% tolerance. At 1% tolerance, all matrix elements are "equal":

In [120]: np.sum(np.isclose(kcsd['KERNEL']*best_x, prec['KERNEL'], rtol=1e-3))/np.size(kcsd['KERNEL'])
Out[120]: 0.9966166858900423

In [121]: np.sum(np.isclose(kcsd['KERNEL']*best_x, prec['KERNEL'], rtol=1e-2))/np.size(kcsd['KERNEL'])
Out[121]: 1.0

kCSD Regrounding

if we reground precalculated kCSD at a point between FCz and Cz:

it becomes kinda similar to regrounded far boundary FEM solution, but not quite...

adding constant to precalculated kCSD

Reminder: potential looks like this:

Then we add a constant to the potential

• 0
  

• 10 V
  

• 20 V
  

• 40 V

Average of the matrix goes up, standard deviation goes down. Kinda makes, sense as we enforce similarity between electrodes?

CSD

potential calulated from point dipole using FEM:

Normal kESI kCSD

kESI kCSD but leadfield precalculated on 1 mm mesh grid

kESI kCSD but leadfield precalculated on 1 mm mesh grid and 10 V constant added

kESI kCSD but leadfield precalculated on 1 mm mesh grid and 20 V constant added

kESI kCSD but leadfield precalculated on 1 mm mesh grid and 40 V constant added

kESI kCSD but leadfield precalculated on 1 mm mesh grid and regrounded between FCz and Cz

animated cycling through them:

Interestingly enough, the added constant CSDs look the same, but they aren't, they differ quite significantly (up to percents)

CSD value at the marker, at maximum of the signal depending on the constant added to the leadfields:

[mdovgialo]

TODO:

Why there is 1.2222122212221224e-08 scaling difference between kCSD and kCSD-precalculated???

kCSD - precalculated CSD solution is on a similar scale as kESI deep electrode + boundary solution.

This makes sense, as the shape and order of magnitude of leadfields are similar:

Although kernels differ:

[mdovgialo]

It seems that, by following kESI-kCSD tutorial, the code was modelling both basis sources and electrodes as R_init sized spheres. Is that correct even? In kESI using MFEM or FeNiCs electrodes were simulated as point sources.

Trying to replace electrode model in:

pbf_kcsd = pbf.Analytical(convolver_interface,
                                  potential=electrode_model_src.potential)

results in this intertesting difference, diagonals becomes slimmer, but the values stay at similar order of magnitude.

Worth to note, that for some reason, pbf.Analytical object has this:

    @_sum_of_not_none
    def _potential_basis_functions(self, electrode):
        return (self._potential_divided_by_relative_conductivity_if_available(
                                                                     electrode),
                super()._potential_basis_functions(electrode))

    def _potential_divided_by_relative_conductivity_if_available(self,
                                                                 electrode):
        return self._divide_by_relative_conductivity_if_available(
                                       self.potential(electrode.x - self.SRC_X,
                                                      electrode.y - self.SRC_Y,
                                                      electrode.z - self.SRC_Z),
                                       electrode)

    def _divide_by_relative_conductivity_if_available(self,
                                                      potential,
                                                      electrode):
        try:
            factor = 1.0 / electrode.conductivity

        except AttributeError:
            return potential

        return factor * potential

This might double apply the conductivity factor, that does not change the shape of electrode potential though, only scales everything...

Question is, why is precalculated kCSD kernel the same as non precalculated but scaled, but the point electrode model is slimmer? Maybe the source and electrode models are reversed???

Question 2: Maybe it's not slimmer, but the diagonal values are just higher (due to point source having higher max potential) and on new color scale it look slimmer? Need to look at crossections of this matrix

Q3: Simulating basis sources in kESI in infinite space???

Identical but scaled? Weird...

Distances to layers don't that much??

kCSD precalcualted seems exactly the same as not precalculated with point electrode model or shperical electrode model but scaled?

Overlaying all to compare visually in post...

Seems like kcsd, kcsd precalculated, kcsd with point source electrode K matrices all have overlapping shapes, but different scaling and maybe offset. While kESI K matrices from far boundary and ground electrode in the neck are a bit different.

• electrode in the neck is just also a boundary condition same as far boundary - 0 potential.

What is interesting, that kESI and kCSD rely on second order spatial differentiation, which should kill any offset, but it seems that offset does matter for the end result. Why?

[mdovgialo]

After consultation with Jakub, some bugs with quadrature weights selection were fixed, and correct estimation of conductance was estimated.

As a results kCSD and kCSD-precalculated have almost exactly the same results:

Kernels:

CSD:

kESI solution has kernels a bit scaled up (makes sense, as "leadfield" is offset compared to theoretical):

But CSD solution is on the same scale, but slightly different, more adhereing to the spherical geometry:

Overall solution value at the maximum point:

[mdovgialo]

Calculating kESI solution when all conductances are equal 0.33:

Leadfield for Cz_L2

Green - kCSD - precalculated using 1/4 pi sigma r
Blue - boundary far and electrode, realistic 4 spheres + air
orange - 4 spheres + air all equal conductivity, far boundary

Solution looks exactly the same as kCSD and kCSD precalc (has additional noise/oscillations):

adding deep grounding electrode to equal conductances:

Creates, basically the same leadfield

regrounded at z=0.075

CSD looks like kCSD:

Kernels overall look like:

kCSD, kCSD precalculated, kESI on realistic spheres with all conductivities set to 0.33 with or without ground electrode, look almost exactly the same.

While kESI on realistic spheres with realistic conductivities and ground electrode has higher K matrix values and different shape, but results in similar valued CSD, but of significantly different shape

TODO: might need testing increasing the amount of MFEM solver iterations to get rid of artifacts, increasing resolution of the mesh might also help