Scripts producing figures for kCSD paper (#104)

* corrected algorithm for distribution of basis sources in 2D

* updated documentation: validation part

* Added script that produces plot for unitary potential propagation

* code refactoring

* Renamed script that produces reliability map plot

* Renamed script that produces example of kCSD reconstruction with reliability map in 2D

* Code and documentation refactoring

* Code refactoring, changed scripts names

* Common scale on the plot

* Common scale on the plot, error formulation as in the paper

* Fixed typo

* Removed redundant variables.

* Fixed round problem - generating electrodes in 3D

* Fixed potential unit

* Electrodes markers on the X-axis

Author: m-kowalska

figures/kCSD_properties/figure_Tbasis.py

@@ -88,8 +88,6 @@ def make_subplot(ax, true_csd, est_csd, estm_x, title=None, ele_pos=None,
         l3 = ax.plot(estm_x, est_csd_LC, '.', label='kCSD L-Curve', lw=2.)
     else:
         l2 = ax.plot(estm_x, est_csd, '--', label='kCSD', lw=2.)
-    s1 = ax.scatter(ele_pos, np.zeros(len(ele_pos)), 17, 'k', label='Electrodes')
-#    ax.legend(fontsize=10)
     ax.set_xlim([0, 1])
     if xlabel:
         ax.set_xlabel('Depth ($mm$)')
@@ -99,8 +97,23 @@ def make_subplot(ax, true_csd, est_csd, estm_x, title=None, ele_pos=None,
         ax.set_title(title)
     if np.max(est_csd) < 1.2:
         ax.set_ylim(-0.2, 1.2)
-    elif np.max(est_csd) > 500:
-        ax.set_yticks([-5000, 0, 5000], [-5000, 0, 5000])
+        s1 = ax.scatter(ele_pos, np.zeros(len(ele_pos))-0.2, 17, 'k', label='Electrodes')
+        s1.set_clip_on(False)
+    elif np.max(est_csd) < 1.7:
+        ax.set_ylim(-10000, 10000)
+        s3 = ax.scatter(ele_pos, np.zeros(len(ele_pos))-10000, 17, 'k', label='Electrodes')
+        s3.set_clip_on(False)
+        ax.set_yticks([-7000, 0, 7000])
+    if np.max(est_csd) > 500:
+        ax.set_ylim(-7000, 7000)
+        s3 = ax.scatter(ele_pos, np.zeros(len(ele_pos))-7000, 17, 'k', label='Electrodes')
+        s3.set_clip_on(False)
+        ax.set_yticks([-5000, 0, 5000])
+    elif np.max(est_csd) > 50:
+        ax.set_ylim(-100, 100)
+        s2 = ax.scatter(ele_pos, np.zeros(len(ele_pos))-100, 17, 'k', label='Electrodes')
+        s2.set_clip_on(False)
+        ax.set_yticks([-70, 0, 70])
     ax.set_xticks([0, 0.5, 1])
     set_axis(ax, letter=letter)
     # ax.legend(frameon=False, loc='upper center', ncol=3)
@@ -493,10 +506,10 @@ def generate_figure_CVLC(R, MU, n_src, true_csd_xlims, total_ele, save_path,
     MU = 0.25
     method = 'cross-validation'  # L-curve
 #    method = 'L-curve'
-    Rs = np.arange(0.1, 0.4, 0.05)
-#    Rs = np.array([0.2])
+#    Rs = np.arange(0.1, 0.4, 0.05)
+    Rs = np.array([0.2])
     lambdas = np.zeros(1)
-#    generate_figure(R, MU, N_SRC, TRUE_CSD_XLIMS, TOTAL_ELE, SAVE_PATH,
-#                    method=method, Rs=Rs, lambdas=lambdas, noise=0)
-    generate_figure_CVLC(R, MU, N_SRC, TRUE_CSD_XLIMS, TOTAL_ELE, SAVE_PATH,
-                         Rs=Rs, lambdas=None, noise=10)
+    generate_figure(R, MU, N_SRC, TRUE_CSD_XLIMS, TOTAL_ELE, SAVE_PATH,
+                    method=method, Rs=Rs, lambdas=lambdas, noise=0)
+#    generate_figure_CVLC(R, MU, N_SRC, TRUE_CSD_XLIMS, TOTAL_ELE, SAVE_PATH,
+#                         Rs=Rs, lambdas=None, noise=10)

…res/kCSD_properties/figure_vectors_v2.py …D_properties/figure_eigensources_M_1D.pyfigures/kCSD_properties/figure_vectors_v2.py renamed to figures/kCSD_properties/figure_eigensources_M_1D.py figures/kCSD_properties/figure_vectors_v2.py renamed to figures/kCSD_properties/figure_eigensources_M_1D.py

@@ -4,16 +4,14 @@
 import os
 from os.path import expanduser
 import numpy as np
+from numpy.linalg import LinAlgError
 import matplotlib.pyplot as plt
 from figure_properties import *
 import matplotlib.gridspec as gridspec
-from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
-from mpl_toolkits.axes_grid1.inset_locator import mark_inset
-from matplotlib.ticker import FuncFormatter
 import datetime
 import time
 
-from kcsd import SpectralStructure, KCSD1D
+from kcsd import KCSD1D
 import targeted_basis as tb
 
 __abs_file__ = os.path.abspath(__file__)
@@ -23,38 +21,24 @@ def _html(r, g, b):
     return "#{:02X}{:02X}{:02X}".format(r, g, b)
 
 
-def stability_M(csd_profile, n_src, ele_lims, true_csd_xlims,
-                total_ele, ele_pos, pots,
-                method='cross-validation', Rs=None, lambdas=None):
+def stability_M(n_src, total_ele, ele_pos, pots, R_init=0.23):
     """
     Investigates stability of reconstruction for different number of basis
     sources
 
     Parameters
     ----------
-    csd_profile: function
-        Function to produce csd profile.
     n_src: int
         Number of basis sources.
-    ele_lims: list
-        Boundaries for electrodes placement.
-    true_csd_xlims: list
-        Boundaries for ground truth space.
     total_ele: int
         Number of electrodes.
     ele_pos: numpy array
         Electrodes positions.
     pots: numpy array
         Values of potentials at ele_pos.
-    method: string
-        Determines the method of regularization.
-        Default: cross-validation.
-    Rs: numpy 1D array
-        Basis source parameter for crossvalidation.
-        Default: None.
-    lambdas: numpy 1D array
-        Regularization parameter for crossvalidation.
-        Default: None.
+    R_init: float
+        Initial value of R parameter - width of basis source
+        Default: 0.23.
 
     Returns
     -------
@@ -70,14 +54,19 @@ def stability_M(csd_profile, n_src, ele_lims, true_csd_xlims,
     for i, value in enumerate(n_src):
         pots = pots.reshape((len(ele_pos), 1))
         obj = KCSD1D(ele_pos, pots, src_type='gauss', sigma=0.3, h=0.25,
-                     gdx=0.01, n_src_init=n_src[i], ext_x=0, xmin=0, xmax=1)
-        if method == 'cross-validation':
-            obj.cross_validate(Rs=Rs, lambdas=lambdas)
-        elif method == 'L-curve':
-            obj.L_curve(Rs=Rs, lambdas=lambdas)
-        ss = SpectralStructure(obj)
-        eigenvectors[i], eigenvalues[i] = ss.evd()
-
+                     gdx=0.01, n_src_init=n_src[i], ext_x=0, xmin=0, xmax=1,
+                     R_init=R_init)
+        try:
+            eigenvalue, eigenvector = np.linalg.eigh(obj.k_pot +
+                                                     obj.lambd *
+                                                     np.identity
+                                                     (obj.k_pot.shape[0]))
+        except LinAlgError:
+            raise LinAlgError('EVD is failing - try moving the electrodes'
+                              'slightly')
+        idx = eigenvalue.argsort()[::-1]
+        eigenvalues[i] = eigenvalue[idx]
+        eigenvectors[i] = eigenvector[:, idx]
         obj_all.append(obj)
     return obj_all, eigenvalues, eigenvectors
 
@@ -112,8 +101,7 @@ def set_axis(ax, x, y, letter=None):
 
 
 def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
-                    save_path, method='cross-validation', Rs=None,
-                    lambdas=None, noise=0):
+                    save_path, noise=0, R_init=0.23):
     """
     Generates figure for spectral structure decomposition.
 
@@ -135,18 +123,12 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
         Electrodes limits.
     save_path: string
         Directory.
-    method: string
-        Determines the method of regularization.
-        Default: cross-validation.
-    Rs: numpy 1D array
-        Basis source parameter for crossvalidation.
-        Default: None.
-    lambdas: numpy 1D array
-        Regularization parameter for crossvalidation.
-        Default: None.
     noise: float
         Determines the level of noise in the data.
         Default: 0.
+    R_init: float
+        Initial value of R parameter - width of basis source
+        Default: 0.23.
 
     Returns
     -------
@@ -159,19 +141,16 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
                                                             noise=noise)
 
     n_src_M = [2, 4, 8, 16, 32, 64, 128, 256, 512]
-    OBJ_M, eigenval_M, eigenvec_M = stability_M(csd_profile, n_src_M,
-                                                ele_lims, true_csd_xlims,
+    OBJ_M, eigenval_M, eigenvec_M = stability_M(n_src_M,
                                                 total_ele, ele_pos, pots,
-                                                method=method, Rs=Rs,
-                                                lambdas=lambdas)
+                                                R_init=R_init)
 
     plt_cord = [(2, 0), (2, 2), (2, 4),
                 (3, 0), (3, 2), (3, 4),
                 (4, 0), (4, 2), (4, 4),
                 (5, 0), (5, 2), (5, 4)]
 
-
-    letters = ['C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'O']
+    letters = ['C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N']
 
     BLACK = _html(0, 0, 0)
     ORANGE = _html(230, 159, 0)
@@ -199,10 +178,8 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
                 linestyle=linestyles[indx], color=colors[indx],
                 marker=markers[indx], label='M='+str(n_src_M[i]),
                 markersize=10)
-#    ax.set_title(' ', fontsize=12)
     ht, lh = ax.get_legend_handles_labels()
     set_axis(ax, -0.05, 1.05, letter='A')
-#    ax.legend(loc='lower left')
     ax.set_xlabel('Number of components')
     ax.set_ylabel('Eigenvalues')
     ax.set_yscale('log')
@@ -213,7 +190,6 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
     ax = fig.add_subplot(gs[0, 3:])
     ax.plot(n_src_M, eigenval_M[:, 0], marker='s', color='k', markersize=5,
             linestyle=' ')
-    #ax.set_title(' ', fontsize=12)
     set_axis(ax, -0.05, 1.05, letter='B')
     ax.set_xlabel('Number of basis sources')
     ax.set_xscale('log')
@@ -229,13 +205,9 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
                     eigenvec_M[j, :, i]),
                     linestyle=linestyles[idx], color=colors[idx],
                     label='M='+str(n_src_M[j]), lw=2)
-            #ax.set_title(r"$\tilde{K}*v_{{%(i)d}}$" % {'i': i+1})
-            ax.text(0.5, 1., r"$\tilde{K}*v_{{%(i)d}}$" % {'i': i+1},
-                    horizontalalignment='center', transform=ax.transAxes, fontsize=20)
-#            ax.locator_params(axis='y', nbins=3)
-
-#            ax.set_xlabel('Depth (mm)', fontsize=12)
-#            ax.set_ylabel('CSD (mA/mm)', fontsize=12)
+            ax.text(0.5, 1., r"$\tilde{K}\cdot{v_{{%(i)d}}}$" % {'i': i+1},
+                    horizontalalignment='center', transform=ax.transAxes,
+                    fontsize=20)
             set_axis(ax, -0.10, 1.1, letter=letters[i])
             if i < 9:
                 ax.get_xaxis().set_visible(False)
@@ -245,23 +217,12 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
             if i % 3 == 0:
                 ax.set_ylabel('CSD ($mA/mm$)')
                 ax.yaxis.set_label_coords(-0.18, 0.5)
-#            ax.yaxis.get_major_formatter().set_powerlimits((0, 1))
-#            ax.tick_params(direction='out', pad=10)
-#            ax.yaxis.get_major_formatter(FormatStrFormatter('%.2f'))
             ax.ticklabel_format(style='sci', axis='y', scilimits=((0.0, 0.0)))
             ax.spines['right'].set_visible(False)
             ax.spines['top'].set_visible(False)
-#     ht, lh = ax.get_legend_handles_labels()
-
-#     ax = fig.add_subplot(gs[3, :])
-#     ax.legend(ht, lh,  fancybox=False, shadow=False, ncol=len(src_idx),
-#               loc='upper center', frameon=False, bbox_to_anchor=(0.5, 0.0))
-#     ax.axis('off')
-
-#    plt.tight_layout()
     fig.legend(ht, lh, loc='lower center', ncol=5, frameon=False)
-    fig.savefig(os.path.join(save_path, 'vectors_' + method +
-                             '_noise_' + str(noise) + '.png'), dpi=300)
+    fig.savefig(os.path.join(save_path, 'vectors_' + '_noise_' +
+                             str(noise) + 'R0_2' + '.png'), dpi=300)
 
     plt.show()
 
@@ -281,6 +242,6 @@ def generate_figure(csd_profile, R, MU, true_csd_xlims, total_ele, ele_lims,
     CSD_PROFILE = tb.csd_profile
     R = 0.2
     MU = 0.25
+    R_init = 0.2
     generate_figure(CSD_PROFILE, R, MU, TRUE_CSD_XLIMS, TOTAL_ELE, ELE_LIMS,
-                    SAVE_PATH, method='cross-validation',
-                    Rs=np.arange(0.1, 0.5, 0.05), noise=None)
+                    SAVE_PATH, noise=None, R_init=R_init)