Replace deprecated scipy integrate.simps with integrate.simpson.

figures/kCSD_properties/L_curve_simulation.py

@@ -1,23 +1,25 @@
-'''
-This script is used to test Current Source Density Estimates, 
+"""
+This script is used to test Current Source Density Estimates,
 using the kCSD method Jan et.al (2012)
 
 This script is in alpha phase.
 
 This was written by :
-Chaitanya Chintaluri, 
+Chaitanya Chintaluri,
 Laboratory of Neuroinformatics,
 Nencki Institute of Exprimental Biology, Warsaw.
-'''
+"""
+
 import os
 import numpy as np
-from scipy.integrate import simps 
+from scipy.integrate import simpson
 import matplotlib.pyplot as plt
 from kcsd import KCSD1D
 from figure_LC import make_plots
 from figure_LCandCVperformance import make_plot_perf
 from figure_LCandCV import plot_surface
 
+
 def do_kcsd(i, ele_pos, pots, **params):
     """
     Function that calls the KCSD1D module
@@ -26,127 +28,197 @@ def do_kcsd(i, ele_pos, pots, **params):
     pots = pots.reshape(num_ele, 1)
     ele_pos = ele_pos.reshape(num_ele, 1)
     k = KCSD1D(ele_pos, pots, **params)
-    noreg_csd = k.values('CSD')
+    noreg_csd = k.values("CSD")
     k.cross_validate(Rs=Rs, lambdas=lambdas)
     errsy = k.errs
-    LandR[1,0,i] = k.lambd
-    LandR[1,1,i] = k.R
-    est_csd_cv = k.values('CSD')
+    LandR[1, 0, i] = k.lambd
+    LandR[1, 1, i] = k.R
+    est_csd_cv = k.values("CSD")
     k.L_curve(lambdas=lambdas, Rs=Rs)
-    LandR[0,0,i] = k.lambd
-    LandR[0,1,i] = k.R
-    est_csd = k.values('CSD')
-    est_pot = k.values('POT')
+    LandR[0, 0, i] = k.lambd
+    LandR[0, 1, i] = k.R
+    est_csd = k.values("CSD")
+    est_pot = k.values("POT")
     return k, est_csd, est_pot, noreg_csd, errsy, est_csd_cv
 
-def generate_csd_1D(src_width, name, srcs=np.array([0]),
-                    start_x=0., end_x=1.,
-                    start_y=0., end_y=1.,
-                    res_x=50, res_y=50):
+
+def generate_csd_1D(
+    src_width,
+    name,
+    srcs=np.array([0]),
+    start_x=0.0,
+    end_x=1.0,
+    start_y=0.0,
+    end_y=1.0,
+    res_x=50,
+    res_y=50,
+):
     """
     Gives CSD profile at the requested spatial location, at 'res' resolution
     """
     csd_y = np.linspace(start_x, end_x, res_x)
     csd_x = np.linspace(start_x, end_x, res_x)
-    src_peak = .7*np.exp(-((csd_x-0.25)**2)/(2.*src_width))
-    src_thr1 = .35*np.exp(-((csd_x-0.45)**2)/(2.*src_width))
-    src_thr2 = .35*np.exp(-((csd_x-0.65)**2)/(2.*src_width))
-    f = src_peak-src_thr1-src_thr2
+    src_peak = 0.7 * np.exp(-((csd_x - 0.25) ** 2) / (2.0 * src_width))
+    src_thr1 = 0.35 * np.exp(-((csd_x - 0.45) ** 2) / (2.0 * src_width))
+    src_thr2 = 0.35 * np.exp(-((csd_x - 0.65) ** 2) / (2.0 * src_width))
+    f = src_peak - src_thr1 - src_thr2
     csd_y = np.linspace(start_x, end_x, res_x)
     return csd_x, csd_y, f
 
+
 def integrate_1D(x0, csd_x, csd, h):
-    m = np.sqrt((csd_x-x0)**2 + h**2) - abs(csd_x-x0)
+    m = np.sqrt((csd_x - x0) ** 2 + h**2) - abs(csd_x - x0)
     y = csd * m
-    I = simps(y, csd_x)
+    I = simpson(y, csd_x)
     return I
 
+
 def calculate_potential_1D(csd, measure_locations, csd_space_x, h):
     sigma = 1
     print(measure_locations.shape, csd_space_x.shape, csd.shape)
     pots = np.zeros(len(measure_locations))
     for ii in range(len(measure_locations)):
         pots[ii] = integrate_1D(measure_locations[ii], csd_space_x, csd, h)
-    pots *= 1/(2.*sigma) #eq.: 26 from Potworowski et al
+    pots *= 1 / (2.0 * sigma)  # eq.: 26 from Potworowski et al
     return pots
 
+
 def transformAndNormalizeLongGrid(D):
-    #D is a dictionary
+    # D is a dictionary
     locations = np.zeros((len(D), 2))
     for i in range(len(D)):
-        locations[i, 0] = D[i]['x']
-        locations[i, 1] = D[i]['y']
+        locations[i, 0] = D[i]["x"]
+        locations[i, 1] = D[i]["y"]
     return locations
 
+
 def generate_electrodes(inpos, lpos, b):
     loc = {}
-    cordx = -(inpos[0] - lpos[0])/b
-    cordy = -(inpos[1] - lpos[1])/b
+    cordx = -(inpos[0] - lpos[0]) / b
+    cordy = -(inpos[1] - lpos[1]) / b
+
     def _addLocToDict(D, chnum, x, y):
         d = {}
-        d['x'] = x
-        d['y'] = y
+        d["x"] = x
+        d["y"] = y
         D[chnum] = d
-    #first add locations just as they suppose to be
+
+    # first add locations just as they suppose to be
     for i in range(b):
-        _addLocToDict(loc, i, inpos[0] + cordx*i, inpos[1] + cordy*i)
+        _addLocToDict(loc, i, inpos[0] + cordx * i, inpos[1] + cordy * i)
     loc_mtrx = transformAndNormalizeLongGrid(loc)
     loc_mtrx = loc_mtrx.T
     return loc_mtrx[0], loc_mtrx[1]
 
+
 def electrode_config(ele_res, true_csd, csd_x, csd_y, inpos, lpos):
     """
     What is the configuration of electrode positions, between what and what positions
     """
     ele_x, ele_y = generate_electrodes(inpos, lpos, ele_res)
     pots = calculate_potential_1D(true_csd, ele_y.T, csd_y, h=1)
-    ele_pos = np.vstack((ele_x, ele_y)).T     #Electrode configs
+    ele_pos = np.vstack((ele_x, ele_y)).T  # Electrode configs
     return ele_pos, pots
 
+
 def main_loop(src_width, total_ele, inpos, lpos, nm, noise=0, srcs=1):
     """
     Loop that decides the random number seed for the CSD profile,
     electrode configurations and etc.
     """
-    #TrueCSD
-    t_csd_x, t_csd_y, true_csd = generate_csd_1D(src_width, nm, srcs=srcs,
-                                                 start_x=0, end_x=1.,
-                                                 start_y=0, end_y=1,
-                                                 res_x=101, res_y=101)
-    if type(noise) ==  float: n_spec = [noise]
-    else: n_spec = noise
+    # TrueCSD
+    t_csd_x, t_csd_y, true_csd = generate_csd_1D(
+        src_width,
+        nm,
+        srcs=srcs,
+        start_x=0,
+        end_x=1.0,
+        start_y=0,
+        end_y=1,
+        res_x=101,
+        res_y=101,
+    )
+    if type(noise) == float:
+        n_spec = [noise]
+    else:
+        n_spec = noise
     for i, noise in enumerate(n_spec):
-        plt.close('all')
+        plt.close("all")
         noise = np.round(noise, 5)
-        print('number of reconstruction: ', i, 'noise level: ', noise)
-        #Electrodes
-        ele_pos, pots = electrode_config(total_ele, true_csd, t_csd_x, t_csd_y, inpos, lpos)
+        print("number of reconstruction: ", i, "noise level: ", noise)
+        # Electrodes
+        ele_pos, pots = electrode_config(
+            total_ele, true_csd, t_csd_x, t_csd_y, inpos, lpos
+        )
         ele_y = ele_pos[:, 1]
         gdX = 0.01
-        x_lims = [0, 1] #CSD estimation place
+        x_lims = [0, 1]  # CSD estimation place
         np.random.seed(srcs)
-        pots_n = pots + (np.random.rand(total_ele)*np.max(abs(pots))-np.max(abs(pots))/2)*noise
-        k, est_csd, est_pot, noreg_csd, errsy, est_csd_cv = do_kcsd(i, ele_y, pots_n, h=1., gdx=gdX, sigma = 1,
-                                                                    xmin=x_lims[0], xmax=x_lims[1], n_src_init=1e3)
-        save_as = nm + '_noise' + str(np.round(noise*100, 1))
+        pots_n = (
+            pots
+            + (np.random.rand(total_ele) * np.max(abs(pots)) - np.max(abs(pots)) / 2)
+            * noise
+        )
+        k, est_csd, est_pot, noreg_csd, errsy, est_csd_cv = do_kcsd(
+            i,
+            ele_y,
+            pots_n,
+            h=1.0,
+            gdx=gdX,
+            sigma=1,
+            xmin=x_lims[0],
+            xmax=x_lims[1],
+            n_src_init=1e3,
+        )
+        save_as = nm + "_noise" + str(np.round(noise * 100, 1))
         m_norm = k.m_norm
         m_resi = k.m_resi
         curve_surf = k.curve_surf
         title = [str(k.lambd), str(k.R), str(noise), nm]
         if plot:
-            make_plots(title, m_norm, m_resi, true_csd, curve_surf, ele_y, pots_n,
-                       pots, k.estm_x, est_pot, est_csd, noreg_csd, save_as)
-            plot_surface(curve_surf, errsy, save_as+'surf')
-        vals_to_save = {'m_norm':m_norm, 'm_resi':m_resi, 'true_csd':true_csd, 
-                        'curve_surf':curve_surf, 'ele_y':ele_y, 'pots_n':pots_n,
-                        'pots':pots, 'estm_x':k.estm_x, 'est_pot':est_pot, 
-                        'est_csd':est_csd, 'noreg_csd':noreg_csd, 'errsy':errsy}
-        np.savez('data_fig4_and_fig13_'+save_as, **vals_to_save)
-        print(true_csd.shape, est_csd[:,0].shape)
-        RMS_wek[0, i] = np.linalg.norm(true_csd/np.linalg.norm(true_csd) - est_csd[:,0]/np.linalg.norm(est_csd[:,0]))
-        RMS_wek[1, i] = np.linalg.norm(true_csd/np.linalg.norm(true_csd) - est_csd_cv[:,0]/np.linalg.norm(est_csd_cv[:,0]))
-
-if __name__=='__main__':
+            make_plots(
+                title,
+                m_norm,
+                m_resi,
+                true_csd,
+                curve_surf,
+                ele_y,
+                pots_n,
+                pots,
+                k.estm_x,
+                est_pot,
+                est_csd,
+                noreg_csd,
+                save_as,
+            )
+            plot_surface(curve_surf, errsy, save_as + "surf")
+        vals_to_save = {
+            "m_norm": m_norm,
+            "m_resi": m_resi,
+            "true_csd": true_csd,
+            "curve_surf": curve_surf,
+            "ele_y": ele_y,
+            "pots_n": pots_n,
+            "pots": pots,
+            "estm_x": k.estm_x,
+            "est_pot": est_pot,
+            "est_csd": est_csd,
+            "noreg_csd": noreg_csd,
+            "errsy": errsy,
+        }
+        np.savez("data_fig4_and_fig13_" + save_as, **vals_to_save)
+        print(true_csd.shape, est_csd[:, 0].shape)
+        RMS_wek[0, i] = np.linalg.norm(
+            true_csd / np.linalg.norm(true_csd)
+            - est_csd[:, 0] / np.linalg.norm(est_csd[:, 0])
+        )
+        RMS_wek[1, i] = np.linalg.norm(
+            true_csd / np.linalg.norm(true_csd)
+            - est_csd_cv[:, 0] / np.linalg.norm(est_csd_cv[:, 0])
+        )
+
+
+if __name__ == "__main__":
     saveDir = "./LCurve/"
     try:
         os.chdir(saveDir)
@@ -158,26 +230,25 @@ def main_loop(src_width, total_ele, inpos, lpos, nm, noise=0, srcs=1):
     src_width = 0.001
     noises = 9
     seeds = 9
-    inpos = [0.5, 0.1]#od dolu
+    inpos = [0.5, 0.1]  # od dolu
     lpos = [0.5, 0.9]
     noise_lvl = np.linspace(0, 0.5, noises)
-    Rs = np.linspace(0.025, 8*0.025, 8)
+    Rs = np.linspace(0.025, 8 * 0.025, 8)
     lambdas = np.logspace(-7, -3, 50)
     sim_results = np.zeros((2, 4, seeds, noises))
     sim_results[:, 3] = noise_lvl
     for src in range(seeds):
-        print('noise seed:', src)
-        RMS_wek = np.zeros((2,len(noise_lvl)))
-        LandR = np.zeros((2,2, len(noise_lvl)))
-        mypath = 'LC' + str(src) + '/'
-        if not os.path.isdir(mypath): os.makedirs(mypath)
+        print("noise seed:", src)
+        RMS_wek = np.zeros((2, len(noise_lvl)))
+        LandR = np.zeros((2, 2, len(noise_lvl)))
+        mypath = "LC" + str(src) + "/"
+        if not os.path.isdir(mypath):
+            os.makedirs(mypath)
         os.chdir(mypath)
-        main_loop(src_width, total_ele, inpos,
-                  lpos, 'LC', noise=noise_lvl, srcs=src)
-        sim_results[:,:2, src] = LandR
+        main_loop(src_width, total_ele, inpos, lpos, "LC", noise=noise_lvl, srcs=src)
+        sim_results[:, :2, src] = LandR
         sim_results[:, 2, src] = RMS_wek
-        os.chdir('..')
-    np.save('sim_results', sim_results)
-    sim_results = np.load('sim_results.npy')
+        os.chdir("..")
+    np.save("sim_results", sim_results)
+    sim_results = np.load("sim_results.npy")
     make_plot_perf(sim_results, noise_lvl)
-