[noisesplit] almost done

modelsusceptovern.py

@@ -93,10 +93,10 @@ def plot_overn(ax, s, files, nmax=1e6, title=False):
     ns = ns[indx]
     stats = stats[indx]
     ax.set_visible(True)
-    ax.plot(ns, stats[:, 7], '0.5', lw=1, zorder=50, label='99.8\\%')
+    ax.plot(ns, stats[:, 7], zorder=50, label='99.8\\%', **s.lsMax)
     ax.fill_between(ns, stats[:, 2], stats[:, 6], fc='0.85', zorder=40, label='5--95\\%')
     ax.fill_between(ns, stats[:, 3], stats[:, 5], fc='0.5', zorder=45, label='25-75\\%')
-    ax.plot(ns, stats[:, 4], zorder=50, label='median', **s.lsSpine)
+    ax.plot(ns, stats[:, 4], zorder=50, label='median', **s.lsMedian)
     #ax.plot(ns, stats[:, 8], '0.0')
     if title:
         if 'noise_frac' in data:

noisesplit.py

@@ -1,6 +1,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from pathlib import Path
+from spectral import whitenoise
 from plotstyle import plot_style
 
 
@@ -33,9 +34,11 @@ def plot_chi2(ax, s, freqs, chi2, nsegs):
     for fac, delta in zip([1, 2, 3, 4, 6, 8, 10],
                           [0.5, 1, 1, 2, 3, 4, 5]):
         if fac*ten >= vmax:
-            vmax = fac*ten
+            vmax = prev_fac*ten
-            ten *= delta
+            ten *= prev_delta
             break
+        prev_fac = fac
+        prev_delta = delta
     pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax,
                        rasterized=True)
     ax.set_xlim(0, 300)
@@ -54,11 +57,60 @@ def plot_chi2(ax, s, freqs, chi2, nsegs):
     cb.outline.set_linewidth(0)
     cax.set_ylabel(r'$|\chi_2|$ [Hz]')
     cax.set_yticks_delta(ten)
+    return cax
 
 
-def plot_chi2_contrast(ax1, ax2, s, cell_name, contrast, nsmall, nlarge):
+def plot_overn(ax, s, files, nmax=1e6):
-    data_files = sims_path.glob(f'chi2-noisen-{cell_name}-{1000*contrast:03.0f}-*.npz')
+    ns = []
-    files, nums = sort_files(cell_name, data_files, 2)
+    stats = []
+    for fname in files:
+        data = np.load(fname)
+        n = data['n']
+        if nmax is not None and n > nmax:
+            continue
+        alpha = data['alpha']
+        freqs = data['freqs']
+        pss = data['pss']
+        dt_fix = 1 # 0.0005
+        chi2 = np.abs(data['prss'])/dt_fix*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
+        ns.append(n)
+        i0 = np.argmin(freqs < 0)
+        i1 = np.argmax(freqs > 300)
+        if i1 == 0:
+            i1 = len(freqs)
+        chi2 = chi2[i0:i1, i0:i1]
+        stats.append(np.quantile(chi2, [0, 0.001, 0.05, 0.25, 0.5,
+                                        0.75, 0.95, 0.998, 1.0]))
+    ns = np.array(ns)
+    stats = np.array(stats)
+    indx = np.argsort(ns)
+    ns = ns[indx]
+    stats = stats[indx]
+    ax.set_visible(True)
+    ax.plot(ns, stats[:, 7], zorder=50, label='99.8\\%', **s.lsMax)
+    ax.fill_between(ns, stats[:, 2], stats[:, 6], fc='0.85', zorder=40, label='5--95\\%')
+    ax.fill_between(ns, stats[:, 3], stats[:, 5], fc='0.5', zorder=45, label='25-75\\%')
+    ax.plot(ns, stats[:, 4], zorder=50, label='median', **s.lsMedian)
+    #ax.plot(ns, stats[:, 8], '0.0')
+    ax.set_xlim(1e2, nmax)
+    ax.set_xscale('log')
+    ax.set_yscale('log')
+    ax.set_yticks_log(numticks=5)
+    if nmax > 1e6:
+        ax.set_ylim(3e-1, 5e3)
+        ax.set_minor_yticks_log(numticks=5)
+        ax.set_xticks_log(numticks=4)
+        ax.set_minor_xticks_log(numticks=8)
+    else:
+        ax.set_ylim(1e0, 1.3e3)
+        #ax.set_minor_yticks_log(numticks=5)
+        ax.set_xticks_log(numticks=5)
+        #ax.set_minor_xticks_log(numticks=6)
+    ax.set_xlabel('segments')
+    ax.set_ylabel('$|\\chi_2|$ [Hz]')
+
+
+def plot_chi2_contrast(ax1, ax2, s, files, nums, nsmall, nlarge):
     for ax, n in zip([ax1, ax2], [nsmall, nlarge]):
         i = nums.index(n)
         data = np.load(files[i])
@@ -67,12 +119,11 @@ def plot_chi2_contrast(ax1, ax2, s, cell_name, contrast, nsmall, nlarge):
         freqs = data['freqs']
         pss = data['pss']
         chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
-        plot_chi2(ax, s, freqs, chi2, n)
+        cax = plot_chi2(ax, s, freqs, chi2, n)
+        cax.set_ylabel('')
 
 
-def plot_chi2_split(ax1, ax2, s, cell_name, nsmall, nlarge):
+def plot_chi2_split(ax1, ax2, s, files, nums, nsmall, nlarge):
-    data_files = sims_path.glob(f'chi2-split-{cell_name}-*.npz')
-    files, nums = sort_files(cell_name, data_files, 1)
     for ax, n in zip([ax1, ax2], [nsmall, nlarge]):
         i = nums.index(n)
         data = np.load(files[i])
@@ -82,11 +133,17 @@ def plot_chi2_split(ax1, ax2, s, cell_name, nsmall, nlarge):
         freqs = data['freqs']
         pss = data['pss']
         chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
-        plot_chi2(ax, s, freqs, chi2, n)
+        cax = plot_chi2(ax, s, freqs, chi2, n)
+        cax.set_ylabel('')
     return alpha, noise_frac
 
         
 def plot_chi2_data(ax, s, cell_name, run):
+    data_file = data_path / f'{cell_name}-baseline.npz'
+    data = np.load(data_file)
+    eodf = float(data['eodf'])
+    ratebase = float(data['ratebase/Hz'])
+    cvbase = float(data['cvbase'])
     data_file = data_path / f'{cell_name}-spectral-s{run:02d}.npz'
     data = np.load(data_file)
     n = data['n']
@@ -94,57 +151,138 @@ def plot_chi2_data(ax, s, cell_name, run):
     freqs = data['freqs']
     pss = data['pss']
     chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
-    print(f'Measured cell {data_file.name} at {100*alpha:.1f}% contrast')
+    print(f'Measured cell {data_file.name} at {100*alpha:.1f}% contrast: r={ratebase:3.0f}Hz, CV={cvbase:4.2f}')
     plot_chi2(ax, s, freqs, chi2, n)
-    return alpha
+    return alpha, ratebase, eodf
 
     
-def plot_noise_split(ax, contrast, noise_contrast, noise_frac):
+def plot_ram(ax, contrast, eodf, wtime, wnoise):
-    axr, axs, axn = ax.subplots(3, 1, hspace=0.1)
     tmax = 50
+    am = 1 + contrast*wnoise
+    eod = np.sin(2*np.pi*eodf*wtime)*am
+
+    ax.show_spines('l')
+    ax.plot(1e3*wtime, eod, clip_on=False, **s.lsEOD)
+    ax.plot(1e3*wtime, +am, clip_on=False, **s.lsAM)
+    ax.plot(1e3*wtime, -am, clip_on=False, **s.lsAM)
+    ax.set_xlim(0, tmax)
+    ax.set_ylim(-1.3, 1.3)
+    ax.set_yticks_delta(1)
+    ax.set_ylabel('EOD')
+    ax.text(1, 1, f'RAM ($c={100*contrast:.0f}$\\,\\%)', ha='right',
+            transform=ax.transAxes, color=s.lsAM['color'])
+
     
+def plot_noise_split(ax, contrast, noise_contrast, noise_frac,
+                     wtime, wnoise):
+    axr, axs, axn = ax.subplots(3, 1, hspace=0.2)
+    cmax = 26
+    cdelta = 20
+    tmax = 50
+
     axr.show_spines('l')
+    axr.axhline(0, **s.lsGrid)
+    axr.plot(1e3*wtime, 100*contrast*wnoise, clip_on=False, **s.lsAM)
     axr.set_xlim(0, tmax)
-    axr.set_ylim(-8, 8)
+    axr.set_ylim(-cmax, cmax)
-    axr.set_yticks_delta(6)
+    axr.set_yticks_delta(cdelta)
     axr.set_ylabel('\\%')
+    axr.text(1, 1, f'RAM ($c={100*contrast:.0f}$\\,\\%)', ha='right',
+             transform=axr.transAxes, color=s.lsAM['color'])
     
     axs.show_spines('l')
+    axs.axhline(0, **s.lsGrid)
+    axs.plot(1e3*wtime, 100*noise_contrast*wnoise, clip_on=False, **s.lsAMsplit)
     axs.set_xlim(0, tmax)
-    axs.set_ylim(-8, 8)
+    axs.set_ylim(-cmax, cmax)
-    axs.set_yticks_delta(6)
+    axs.set_yticks_delta(cdelta)
     axs.set_ylabel('\\%')
-    
+    if noise_contrast > 0:
-    axn.set_ylim(-6, 6)
+        axs.text(1, 1, f'$s_{{\\xi}}(t)$ ($c={100*noise_contrast:.0f}$\\,\\%)',
+                 ha='right', transform=axs.transAxes,
+                 color=s.lsAMsplit['color'])
+
+    ntime = np.linspace(0, 1e-3*tmax, 800)
+    rng = np.random.default_rng(45432)
+    nnoise = rng.normal(size=len(ntime))
+    axn.show_spines('l')
+    axn.axhline(0, **s.lsGrid)
+    axn.plot(1e3*ntime, noise_frac*nnoise, clip_on=False, **s.lsNoise)
+    axn.set_ylim(-2, 2)
     axn.set_xlim(0, tmax)
-    axn.set_yticks_delta(6)
+    axn.set_yticks_delta(5)
     axn.set_yticks_blank()
-    axn.set_xticks_delta(25)
+    #axn.set_xticks_delta(25)
-    axn.set_xlabel('Time', 'ms')
+    #axn.set_xlabel('Time', 'ms')
+    y = 0.8 if noise_frac < 1 else 1.2
+    axn.text(1, y, f'Intrinsic noise (${100*noise_frac:.0f}$\\,\\%)',
+             ha='right', transform=axn.transAxes, color=s.lsNoise['color'])
+    if noise_frac < 1:
+        axn.xscalebar(1, 0, 10, 'ms', ha='right')
+        
+    return axr
 
     
 if __name__ == '__main__':
-    cell_name = '2012-07-03-ak-invivo-1'
+    #cell_name = ['2012-07-03-ak-invivo-1', 0]
+    cell_name = ['2017-07-18-ai-invivo-1', 1]  # Take this! at 3% model, 5% data
     nsmall = 100
     nlarge = 1000000
     contrast = 0.03
+
+    wdt = 0.0001
+    wnoise = whitenoise(0, 300, wdt, 0.05, rng=np.random.default_rng(51234))
+    wtime = np.arange(len(wnoise))*wdt
                  
     s = plot_style()
-    fig, axs = plt.subplots(2, 4, cmsize=(s.plot_width, 0.4*s.plot_width),
+    fig, axs = plt.subplots(3, 4, cmsize=(s.plot_width, 0.7*s.plot_width),
-                            width_ratios=[1, 0, 1, 1, 1])
+                            width_ratios=[1, 0, 1, 1, 0.15, 1])
-    fig.subplots_adjust(leftm=7, rightm=8, topm=2, bottomm=3.5,
+    fig.subplots_adjust(leftm=8, rightm=1.5, topm=3, bottomm=4,
-                        wspace=0.4, hspace=0.6)
+                        wspace=0.25, hspace=0.8)
-    axs[1, 0].set_visible(False)
+    axs[0, 2].set_visible(False)
-    data_contrast = plot_chi2_data(axs[0, 0], s, cell_name[:13], 0)
+    axs[0, 3].set_visible(False)
-    plot_noise_split(axs[0, 1], data_contrast, 0, 1)
+    
-    plot_chi2_contrast(axs[0, 2], axs[0, 3], s, cell_name, contrast, nsmall, nlarge)
+    # data:
-    noise_contrast, noise_frac = plot_chi2_split(axs[1, 2], axs[1, 3], s,
+    axs[0, 1].text(-2.42, 1.2, 'P-unit data', fontsize='large',
-                                                 cell_name, nsmall, nlarge)
+                   transform=axs[0, 1].transAxes, color=s.punit_color1)
-    plot_noise_split(axs[1, 1], contrast, noise_contrast, noise_frac)
+    data_contrast, ratebase, eodf = plot_chi2_data(axs[0, 1], s, cell_name[0],
+                                                   cell_name[1])
+    plot_ram(axs[0, 0], data_contrast, eodf, wtime, wnoise)
+    axs[0, 1].text(-1.5, 1.2, f'$r={ratebase:.0f}$\\,Hz',
+                   transform=axs[0, 1].transAxes, fontsize='large')
+    
+    # model:
+    data_files = sims_path.glob(f'chi2-noisen-{cell_name[0]}-{1000*contrast:03.0f}-*.npz')
+    files, nums = sort_files(cell_name[0], data_files, 2)
+    axs[1, 1].text(-2.42, 1.2, 'P-unit model', fontsize='large',
+                   transform=axs[1, 1].transAxes, color=s.model_color1)
+    plot_chi2_contrast(axs[1, 1], axs[1, 2], s, files, nums, nsmall, nlarge)
+    axr1 = plot_noise_split(axs[1, 0], contrast, 0, 1, wtime, wnoise)
+    plot_overn(axs[1, 3], s, files, nmax=1e6)
+    axs[1, 3].legend(loc='lower center', bbox_to_anchor=(0.5, 1.1),
+                     markerfirst=False, title='$|\\chi_2|$ percentiles')
+
+    # model noise split:
+    data_files = sims_path.glob(f'chi2-split-{cell_name[0]}-*.npz')
+    files, nums = sort_files(cell_name[0], data_files, 1)
+    axs[2, 1].text(-2.42, 1.2, 'P-unit model', fontsize='large',
+                   transform=axs[2, 1].transAxes, color=s.model_color1)
+    axs[2, 1].text(-1.5, 1.2, f'(noise split)', fontsize='large',
+                   transform=axs[2, 1].transAxes)
+    noise_contrast, noise_frac = plot_chi2_split(axs[2, 1], axs[2, 2], s,
+                                                 files, nums, nsmall, nlarge)
+    axr2 = plot_noise_split(axs[2, 0], 0, noise_contrast, noise_frac,
+                            wtime, wnoise)
+    plot_overn(axs[2, 3], s, files, nmax=1e6)
+    
+    fig.common_xticks(axs[:, 1])
     fig.common_xticks(axs[:, 2])
     fig.common_xticks(axs[:, 3])
-    fig.common_yticks(axs[0, 2:])
+    fig.common_yticks(axs[1, 1:3])
-    fig.common_yticks(axs[1, 2:])
+    fig.common_yticks(axs[2, 1:3])
-    #fig.tag(axs, xoffs=-4.5, yoffs=1.8)
+    fig.tag([axs[0, :2],
+             [axr1] + axs[1, 1:].tolist(),
+             [axr2] + axs[2, 1:].tolist()],
+            xoffs=-4.5, yoffs=2)
     fig.savefig()
     print()

plotstyle.py

@@ -55,7 +55,11 @@ def plot_style():
     pt.make_line_styles(ns, 'ls', 'Line', '', palette['black'], '-',
                         lwthin)
     pt.make_line_styles(ns, 'ls', 'EOD', '', palette['gray'], '-', lwthin)
-    pt.make_line_styles(ns, 'ls', 'AM', '', palette['red'], '-', lwthick)
+    pt.make_line_styles(ns, 'ls', 'AM', '', palette['red'], '-', lwmid)
+    pt.make_line_styles(ns, 'ls', 'AMsplit', '', palette['orange'], '-', lwmid)
+    pt.make_line_styles(ns, 'ls', 'Noise', '', palette['gray'], '-', lwmid)
+    pt.make_line_styles(ns, 'ls', 'Median', '', palette['black'], '-', lwthick)
+    pt.make_line_styles(ns, 'ls', 'Max', '', palette['black'], '-', lwmid)
 
     ns.lsStim = dict(color='gray', lw=ns.lwmid)
     ns.lsRaster = dict(color='black', lw=ns.lwthin)

punitexamplecell.py

@@ -14,9 +14,9 @@ run2 = 1
 example_cells = [
                  ['2021-06-18-ae-invivo-1', 3],    #  98Hz, 1%, ok
                  ['2012-03-30-ah', 2],             # 177Hz, 2.5%, 2.0, nice
-                 ##['2012-07-03-ak', 0],             # 120Hz, 2.5%, 1.8, broader
+                 ##['2012-07-03-ak', 0],             # 120Hz, 2.5%, 1.8, broader, the one model cell, nice triangle up to 1%!
-                 ##['2012-12-20-ac', 0],             # 213Hz, 2.5%, 2.1, ok
+                 ##['2012-12-20-ac', 0],             # 213Hz, 2.5%, 2.1, ok, model cell, weak triangle up to 1%!
-                 #['2017-07-18-ai-invivo-1', 1],    #  78Hz, 5%, 2.3, weak
+                 #['2017-07-18-ai-invivo-1', 1],    #  78Hz, 5%, 2.3, weak, nice model cell with clear triangle up to 10%!
                  ##['2019-06-28-ae', 0],             # 477Hz, 10%, 2.6, weak
                  ##['2020-10-27-aa-invivo-1', 4],    # 259Hz, 0.5%, 2.0, ok
                  ##['2020-10-27-ae-invivo-1', 4],    # 375Hz, 0.5%, 4.3, nice, additional low freq line