gin commit from platypus

Modified files: 5

code/analyses/heterogeneous_populations.py

@@ -10,22 +10,28 @@ spike_buffer = {}
 stimulus_name = "gwn300Hz10s0.3.dat"
 
 
-def population_spikes(df, data_location, population_size, delay=0.0):
-    """_summary_
+def population_spikes(df, data_location, population_size, subtract_delay=True):
+    """Assembles the spike times for a 'population' of neurons.
 
     Parameters
     ----------
-    df : _type_
-        _description_
-    population_size : _type_
-        _description_
-    delay : float, optional
-        _description_, by default 0.0
+    df : pd.DataFrame
+        The DataFrame containing the information about where to find the recorded white noise responses.
+    population_size : int
+        The size of the population, i.e. the number of randomly selected trials.
+    subtract_delay : bool, optional
+        Defines whether the neuronal delay should be subtracted from the spike times, by default True
 
     Returns
     -------
-    _type_
-        _description_
+    list:
+        list of spike times of the randomly selected trials.
+    list:
+        list containing whether or not the response needs to be inverted.
+    list:
+        the respective cell ids
+    list:
+        the respective trial ids.
     """
     cell_indices = np.arange(len(df.dataset_id.unique()))
     np.random.shuffle(cell_indices)
@@ -50,7 +56,8 @@ def population_spikes(df, data_location, population_size, delay=0.0):
         trial_index = np.random.randint(min(cell_trials.index), max(cell_trials.index)+1)
         trial = cell_trials[cell_trials.index == trial_index]
         trial_spikes = all_spikes[(all_spikes >= trial.start_time.values[0]) & (all_spikes < trial.end_time.values[0])] - trial.start_time.values[0]
-        trial_spikes -= trial.delay.values[0]
+        if subtract_delay:
+            trial_spikes -= trial.delay.values[0]  # compensate for the delay between stimulus and response
         spikes.append(trial_spikes[trial_spikes > 0])
         inversion_needed.append(trial.inverted.values[0])
         cell_ids.append(c)
@@ -60,6 +67,36 @@ def population_spikes(df, data_location, population_size, delay=0.0):
 
 def mutual_info_heterogeneous_population(df, data_location, stim_location, population_size=10, delay=0.01, repeats=10,
                                          kernel_sigma=0.00125, saving=False, result_location=".", delay_type="equal"):
+    """_summary_
+
+    Parameters
+    ----------
+    df : _type_
+        _description_
+    data_location : _type_
+        _description_
+    stim_location : _type_
+        _description_
+    population_size : int, optional
+        _description_, by default 10
+    delay : float, optional
+        _description_, by default 0.01
+    repeats : int, optional
+        _description_, by default 10
+    kernel_sigma : float, optional
+        _description_, by default 0.00125
+    saving : bool, optional
+        _description_, by default False
+    result_location : str, optional
+        _description_, by default "."
+    delay_type : str, optional
+        _description_, by default "equal"
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
     coherences = []
     frequency = []
     _, stimulus = load_white_noise_stim(os.path.join(stim_location, stimulus_name))
@@ -69,8 +106,9 @@ def mutual_info_heterogeneous_population(df, data_location, stim_location, popul
         r = {"pop_size": population_size, "delay": delay, "true_delay": 0.0, "kernel_sigma": kernel_sigma, "num_inversions": 0, "population_rate": 0.0, "rate_modulation": 0.0,
              "snr": 0.0, "mi":0.0, "mi_100": 0.0, "mi_200": 0.0, "mi_300": 0.0, "cell_ids":[], "mtag_ids":[], "trial_ids":[], "result_file":"" }
         print("population size: %i, delay: %.5fs, kernel: %.5f, repeat: %i" % (population_size, delay, kernel_sigma, i), end="\r" )
-        pop_spikes, inversion_needed, cell_ids, trial_ids = population_spikes(df, data_location, population_size, delay)
-        f, c, mis, rates, true_delay = mutual_info(pop_spikes, delay, inversion_needed, stimulus, freq_bin_edges=[0, 100, 200, 300], kernel_sigma=kernel_sigma)
+        subtract_delay = delay > 0.0
+        pop_spikes, inversion_needed, cell_ids, trial_ids = population_spikes(df, data_location, population_size, subtract_delay)
+        f, c, mis, rates, true_delay = mutual_info(pop_spikes, delay, inversion_needed, stimulus, freq_bin_edges=[0, 100, 200, 300], kernel_sigma=kernel_sigma, delay_type=delay_type)
 
         coherences.append(c)
         frequency = f
@@ -99,7 +137,7 @@ def mutual_info_heterogeneous_population(df, data_location, stim_location, popul
         r["peak_coh"] = max_coh
         r["upper_cutoff"] = upper_cutoff
         r["lower_cutoff"] = lower_cutoff
-        r["peak_freq"] = peak_f   
+        r["peak_freq"] = peak_f
         r["snr"] = snr if snr is not None else 0.0
         r["variability"] = variability if variability is not None else 0.0
         r["mi"] = mis[0]
@@ -119,6 +157,34 @@ def run_heterogeneous_population_analysis(whitenoise_trial_file, data_location,
                                           delays=[0.0, 0.0005, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.0125, 0.015],
                                           kernels=[0.000125, 0.00025, 0.0005, 0.001, 0.002, 0.003, 0.005, 0.01, 0.015, 0.02, 0.025],
                                           repeats=50, saving=False, num_cores=1):
+    """Runs the population analysis for heterogeneous populations.
+
+    Parameters
+    ----------
+    whitenoise_trial_file : str
+        The dataframe containing the information of the white noise recordings.
+    data_location : str
+        The path to the raw data.
+    stim_location : str
+        The path to the folder containing the stimuli.
+    population_sizes : list of ints, optional
+        List of population sizes, by default [2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30]
+    delays : list, optional
+        list of artificial delays added to the spike times. -1 indicates that the original delay should not be changed, by default [0.0, 0.0005, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.0125, 0.015]
+    kernels : list, optional
+        list of Gaussian kernel standard deviations used for firing rate estimation, by default [0.000125, 0.00025, 0.0005, 0.001, 0.002, 0.003, 0.005, 0.01, 0.015, 0.02, 0.025]
+    repeats : int, optional
+        Maximum number of repetitions for each combination of population sizes, delays, kernels. For each repetition a new set of responses will be used. By default 50
+    saving : bool, optional
+        Flag that indicates whether or not coherence spectra, firing rates etc should be saved, by default False
+    num_cores : int, optional
+        number of parallel jobs spawned to run the analyses, by default 1
+
+    Returns
+    -------
+    pd.DataFrame
+        The DataFrame containing the analysis results.
+    """
     result_dicts = []
     whitenoise_trials = pd.read_csv(whitenoise_trial_file, sep=";", index_col=0)
     whitenoise_trials = whitenoise_trials[(whitenoise_trials.stimfile == stimulus_name) & (whitenoise_trials.duration == 10)]  

code/plots/figure4.py

@@ -5,6 +5,7 @@ import numpy as np
 import pandas as pd
 import matplotlib as mplt
 import matplotlib.pyplot as plt
+plt.style.use("./code/plots/pnas_onecolumn.mplstyle")
 
 from .figure_style import subfig_labelsize, subfig_labelweight, despine
 
@@ -116,7 +117,7 @@ def compare_homogeneous_with_similar_heterogeneous(homogeneous, heterogeneous, a
             hom_populations = homogeneous[(homogeneous.pop_size == ps)]
             het_populations = heterogeneous[(heterogeneous.pop_size == ps) & (heterogeneous.delay == 0) & (heterogeneous.kernel_sigma == 0.001)]  
             het_mis = np.zeros(len(het_populations))
-            
+
             for count, (_, row) in enumerate(het_populations.iterrows()):
                 feat_value = row[f]    
                 het_mis[count] = row["mi_100"]
@@ -133,7 +134,7 @@ def compare_homogeneous_with_similar_heterogeneous(homogeneous, heterogeneous, a
 
 
 def layout_figure():
-    fig = plt.figure(figsize=(4.3, 4.2))
+    fig = plt.figure(figsize=(3.42, 4.2))
     shape = (20, 7)
     ax_total = plt.subplot2grid(shape, (0, 0), rowspan=8, colspan=3)
     ax_100 = plt.subplot2grid(shape, (0, 4), rowspan=6, colspan=3 )
@@ -193,14 +194,16 @@ def plot_delay_effect(args):
         despine(ax, ["top", "right"], False)  
     axes[0].set_ylabel("mutual information [bit/s]")
     axes[0].set_xlabel("population size")
-    # ax_total.xaxis.set_label_coords(0.5, -0.15)
+    axes[0].yaxis.set_label_coords(-0.275, 0.5)
   
     axes[1].set_xticklabels([])
     axes[2].set_xticklabels([])
     axes[2].set_ylabel("mutual information [bit/s]")
+    axes[2].yaxis.set_label_coords(-0.225, 0.5)
+
     axes[3].set_xlabel("population size", ha="center") 
     axes[3].xaxis.set_label_coords(0.5, -0.35)
-    
+
     delays[0] += 0.00025  # just for the log scale
     for col, title, cmap, marker in zip(columns[1:], titles[1:], cmaps[1:], markers[1:]):
         avg, yerr = get_mutual_info_delay(df, 16, column=col, error="std", kernel=args.kernel)
@@ -209,7 +212,7 @@ def plot_delay_effect(args):
                           markeredgecolor="white", zorder=2)
         axes[-1].plot(delays*1000, avg/avg[0], c=cmap(0.75), lw=1.0, ls="--", zorder=1)
        
-    axes[-1].legend(ncol=1, handletextpad=0.1, loc="lower left", columnspacing=.5, fontsize=6, frameon=True,
+    axes[-1].legend(ncol=1, handletextpad=0.1, loc="lower left", columnspacing=.5, frameon=True,
                       fancybox=False, shadow=False, markerscale=0.9, borderaxespad=0.05)
     axes[-1].set_xlabel(r"$\sigma_{delay}$ [ms]")
     axes[-1].set_ylabel(r"mutual information [rel.]")
@@ -223,7 +226,8 @@ def plot_delay_effect(args):
     for i in range(1, len(yticklabels), 2): 
         yticklabels[i] = "" 
     yticklabels[yticklabels.index("7.0")] = ""
-    axes[-1].set_xticklabels(yticklabels, rotation=45, fontsize=7)
+    axes[-1].set_xticklabels(yticklabels, rotation=45)
+    axes[-1].yaxis.set_label_coords(-0.275, 0.5)
     despine(axes[-1], ["top", "right"], False)
 
     if args.nosave:

code/plots/figure5.py

@@ -7,6 +7,7 @@ import matplotlib.image as mpimg
 
 from .figure_style import subfig_labelsize, subfig_labelweight, despine
 
+plt.style.use("./code/plots/pnas_onecolumn.mplstyle")
 
 def plot_info_errorbar(axis, pop_size, info, label="", ls="-", color="tab:blue"):
     info_mean = np.mean(info, axis=1)
@@ -32,27 +33,30 @@ def plot_info_per_band(ax1, ax2, ax3, population_size, info_no_delay, info1_dela
     plot_info_errorbar(ax3, population_size, info3_delay, label=labels[3], ls=ls, color=colors[3])
 
 
-def pimp_lif_sim_axes(axis, set_ylabel=True):
+def pimp_lif_sim_axes(axis, set_ylabel=True, set_xlabel=True):
     despine(axis, ["top", "right"], False)
 
     axis.set_xlim([0, 210])
     axis.set_xticks(np.arange(0, 201, 100))
     axis.set_xticks(np.arange(0, 201, 50), minor=True)
     axis.set_xticklabels(np.arange(0, 201, 100))
-    axis.set_xlabel("population size", labelpad=1.5)
+    if set_xlabel:
+        axis.set_xlabel("population size", labelpad=1.5)
 
     axis.set_ylim([0, 750])
     axis.set_yticks(np.arange(0, 801, 200))
     axis.set_yticks(np.arange(0, 801, 100), minor=True)
     if set_ylabel:
         axis.set_ylabel("mutual information [bits/s]")
-
-    axis.legend(frameon=True, fontsize=5, loc=2)
+    else:
+        axis.set_yticklabels([])
+    if set_xlabel:
+        axis.legend(bbox_to_anchor=(0.3, -1.075, 2.1, 0.25), ncol=2, frameon=True)
 
 
-def add_additional_xaxes(fig, axis, population_size, density, conduction_velocity, delay_dx=1):
+def add_additional_xaxes(fig, axis, population_size, density, conduction_velocity, delay_dx=1, add_label=True):
     pos = axis.get_position().bounds
-    ax2 = fig.add_axes((pos[0], 0.2, pos[2], 0.0))
+    ax2 = fig.add_axes((pos[0], 0.325, pos[2], 0.0))
     ax2.yaxis.set_visible(False) # hide the yaxis
     max_rf = max(population_size) / density * 1000
     rf_ticklabels = np.arange(0, np.ceil(max_rf)+1, np.ceil(max_rf/2), dtype=int)
@@ -62,21 +66,23 @@ def add_additional_xaxes(fig, axis, population_size, density, conduction_velocit
     ax2.set_xticks(rf_ticks)
     ax2.set_xticklabels(rf_ticklabels)
     ax2.set_xticks(rf_minorticks, minor=True)
-    ax2.set_xlabel("receptive field size [mm]", labelpad=1.52)
+    if add_label:
+        ax2.set_xlabel("receptive field size [mm]", labelpad=1.52)
     ax2.set_xlim([0, 210])
 
     pos = axis.get_position().bounds
-    ax3 = fig.add_axes((pos[0], 0.11, pos[2], 0.0))
+    ax3 = fig.add_axes((pos[0], 0.2, pos[2], 0.0))
     ax3.yaxis.set_visible(False) # hide the yaxis
     max_delay = np.ceil(max(population_size) / density / conduction_velocity * 1000)
-    delay_xticklabels = np.arange(0.0, max_delay+delay_dx, delay_dx)
+    delay_xticklabels = np.arange(0, max_delay+delay_dx, delay_dx, dtype=int)
     delay_xticks  = delay_xticklabels /1000 * conduction_velocity * density
     minor_delay_xticklabels = np.arange(0.0, max_delay+1, delay_dx/5)
     minor_delay_xticks = minor_delay_xticklabels / 1000 * conduction_velocity * density
     ax3.set_xticks(delay_xticks)
     ax3.set_xticklabels(delay_xticklabels)
     ax3.set_xticks(minor_delay_xticks, minor=True)
-    ax3.set_xlabel("maximum delay [ms]", labelpad=1.5)
+    if add_label:
+        ax3.set_xlabel("maximum delay [ms]", labelpad=1.5)
     ax3.set_xlim([0, 210])
 
 
@@ -101,28 +107,28 @@ def add_model_sketches(fig):
 
 
 def layout_figure(titles, colors):
-    fig = plt.figure(figsize=(7, 5))
+    fig = plt.figure(figsize=(3.42, 3.5))#, constrained_layout=True)
     axes = []
 
     fig_grid = (8, 8)
     colors = ["tab:blue", "tab:red", "tab:orange", "tab:green"]
     # low velocity
-    axes.append(plt.subplot2grid(fig_grid, (3, 0), 4, 2))
+    axes.append(plt.subplot2grid(fig_grid, (0, 0), 4, 2))
     # squid velocity
-    axes.append(plt.subplot2grid(fig_grid, (3, 3), 4, 2))
+    axes.append(plt.subplot2grid(fig_grid, (0, 3), 4, 2))
     # efish velocity
-    axes.append(plt.subplot2grid(fig_grid, (3, 6), 4, 2))
+    axes.append(plt.subplot2grid(fig_grid, (0, 6), 4, 2))
 
-    axes[0].text(0.5, 1.05, titles[0], transform=axes[0].transAxes, ha="center",fontsize=6, color=colors[1])
-    axes[0].text(-0.5, 1.95, "A", transform=axes[0].transAxes, fontsize=subfig_labelsize, weight=subfig_labelweight, color="k")
-    axes[1].text(0.5, 1.05, titles[1], transform=axes[1].transAxes, ha="center",fontsize=6, color=colors[2])
-    axes[1].text(-0.5, 1.95, "B", transform=axes[1].transAxes, fontsize=subfig_labelsize, weight=subfig_labelweight)
-    axes[2].text(0.5, 1.05, titles[2], transform=axes[2].transAxes, ha="center", fontsize=6, color=colors[3])
-    axes[2].text(-0.5, 1.95, "C", transform=axes[2].transAxes,fontsize=subfig_labelsize, weight=subfig_labelweight)
+    axes[0].text(0.5, 0.95, titles[0], transform=axes[0].transAxes, ha="center",fontsize=8, color=colors[1])
+    axes[0].text(-0.5, 1.1, "A", transform=axes[0].transAxes, fontsize=subfig_labelsize, weight=subfig_labelweight, color="k")
+    axes[1].text(0.5, 0.95, titles[1], transform=axes[1].transAxes, ha="center",fontsize=8, color=colors[2])
+    axes[1].text(-0.35, 1.1, "B", transform=axes[1].transAxes, fontsize=subfig_labelsize, weight=subfig_labelweight)
+    axes[2].text(0.5, 0.95, titles[2], transform=axes[2].transAxes, ha="center", fontsize=8, color=colors[3])
+    axes[2].text(-0.35, 1.1, "C", transform=axes[2].transAxes,fontsize=subfig_labelsize, weight=subfig_labelweight)
 
-    fig.subplots_adjust(left=0.10, bottom=0.20, top=0.90, right=0.95, hspace=0.02)
+    fig.subplots_adjust(left=0.15, bottom=0.0, top=0.9, right=0.975, wspace=-0.2)
 
-    add_model_sketches(fig)
+    #add_model_sketches(fig)
     return fig, axes
 
 
@@ -136,10 +142,10 @@ def lif_simulations(args):
     titles = [r"%.1f$m/s$" % 7.0,
               r"%.1f$m/s$" % 25.,
               r"%.1f$m/s$ " % 50.]
-    labels = ["w/o delay; 0 - 100 Hz",
-              "w delay; 0 - 100 Hz", 
-              "w delay; 100 - 200 Hz",
-              "w delay; 200 - 300 Hz"]
+    labels = ["without delay; 0 - 100 Hz",
+              "with delay; 0 - 100 Hz", 
+              "with delay; 100 - 200 Hz",
+              "with delay; 200 - 300 Hz"]
     line_styles = ["--", "-.", ":"]
     colors = ["tab:blue", "tab:red", "tab:orange", "tab:green"]
 
@@ -165,11 +171,11 @@ def lif_simulations(args):
             plot_info_per_band(axes[0], axes[1], axes[2], pop_size, None, info_slow, info_squid, info_efish, colors=colors, ls=line_styles[i], labels=lbls)
 
     for i, a in enumerate(axes):
-        pimp_lif_sim_axes(a, i == 0)
+        pimp_lif_sim_axes(a, i == 0, i==1)
 
-    add_additional_xaxes(fig, axes[0], pop_size, density, vel_low, 5)
-    add_additional_xaxes(fig, axes[1], pop_size, density, vel_med, 1.0)
-    add_additional_xaxes(fig, axes[2], pop_size, density, vel_high, .5)
+    add_additional_xaxes(fig, axes[0], pop_size, density, vel_low, 5, add_label=False)
+    add_additional_xaxes(fig, axes[1], pop_size, density, vel_med, 2.0)
+    add_additional_xaxes(fig, axes[2], pop_size, density, vel_high, 1.0, add_label=False)
 
     if args.nosave:
         plt.show()

code/util.py

@@ -170,7 +170,6 @@ def get_temporal_shift(sta_time, sta):
 
 def load_stim(filename):
     """
-
     Loads a data file saved by relacs. Returns a tuple of dictionaries
     containing the data and the header information
 
@@ -278,33 +277,71 @@ def firing_rate(spikes, duration, sigma=0.005, dt=1./20000.):
 
 
 def coherence_rate(coherence, frequency, start_freq, stop_freq):
+    """Estimates the coherence rate as a lower bound estimation of the mutual information between stimulus and response.
+
+    Parameters
+    ----------
+    coherence : _type_
+        _description_
+    frequency : _type_
+        _description_
+    start_freq : _type_
+        _description_
+    stop_freq : _type_
+        _description_
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
     deltaf = np.mean(np.diff(frequency))
     spectrum = -np.log2((np.ones(coherence.shape) - coherence)) * deltaf
-    
+
     info = np.sum(spectrum[(frequency >= start_freq) & (frequency < stop_freq)])
     return info
 
 
+def mutual_info(spike_responses, artificial_delay, inversion_needed, stimulus, freq_bin_edges, kernel_sigma=0.00125, delay_type="equal", stepsize=1./20000, trial_duration=10.):
+    """Estimates a the lower bound mutual information carried by the responses about the stimulus. The MI is estimated using the stimulus response coherence. See Borst & Theunissen, 2001.
 
-def mutual_info(spike_responses, delay, inversion_needed, stimulus, freq_bin_edges, kernel_sigma=0.00125, delay_type="equal", stepsize=1./20000, trial_duration=10.):
-    """[summary]
-
-    Args:
-        spike_responses (list of list of spike times): [description]
-        delay (float): [description]
-        inversion_needed ([type]): [description]
-        stimulus ([type]): [description]
-        freq_bin_edges ([type]): [description]
-        kernel_sigma (float, optional): [description]. Defaults to 0.00125.
+    Parameters
+    ----------
+    spike_responses : list of np.ndarrays of float
+        The spike times
+    artificial_delay : float
+        The artificial delay that is added to the spike times. If delay is less or equal 0.0 no delay will be added
+    inversion_needed : bool
+        Whether or not the firing rate is inverted, may happen due to the receptor location, stimulus geometry
+    stimulus : np.ndarray
+        The stimulus waveform.
+    freq_bin_edges : list of float
+        The edges of the frequency bands to be analysed.
+    kernel_sigma : float, optional
+        Standard deviation of the Gaussian kernel used for firing rate estimation, by default 0.00125
+    delay_type : str, optional
+        Type of destribution of the delays, by default "equal", alternatively "gaussian"
+    stepsize : float, optional
+        Temporal stepsize of the stimulus trace, by default 1./20000
+    trial_duration : float, optional
+        The duration of a recorded trial, by default 10. The stimulus waveform may be longer than the responses are actually recorded.
 
-    Returns:
-        np.array: array containing the frequency axis of the coherence function
-        np.array: the average coherence function
-        np.array: mutual information values, shape is len of freq_bin_edges + 1
-        np.array: avg_rate; the average firing rate
-        np.array: rate_error; std of firing rates as function of time
-        float: true_delay; avg of the real delays
-        np.array: variability; the time-averaged rate 
+    Returns
+    -------
+    np.array: 
+        array containing the frequency axis of the coherence function
+    np.array:
+        the average coherence function
+    np.array:
+        mutual information values, shape is len of freq_bin_edges + 1
+    np.array:
+        avg_rate; the average firing rate
+    np.array:
+        rate_error; std of firing rates as function of time
+    float: 
+        true_delay; avg of the real delays
+    np.array:
+        variability; the time-averaged rate 
     """
     population_rates = None
     rng = np.random.default_rng()
@@ -313,10 +350,11 @@ def mutual_info(spike_responses, delay, inversion_needed, stimulus, freq_bin_edg
         if isinstance(sr, list):
             sr = np.array(sr)
         dt = 0.0
-        if "equal" in delay_type:
-            dt = rng.random() * delay * 2
-        elif "gaussian" in delay_type:
-            dt = rng.standard_random() * delay
+        if artificial_delay > 0.0:
+            if "equal" in delay_type:
+                dt = rng.random() * artificial_delay * 2
+            elif "gaussian" in delay_type:
+                dt = rng.standard_random() * artificial_delay
         delays[i] = dt
         sr += dt
         rate = firing_rate(sr, trial_duration, kernel_sigma, dt=stepsize)
@@ -337,7 +375,6 @@ def mutual_info(spike_responses, delay, inversion_needed, stimulus, freq_bin_edg
     return f, c, mis, population_rates, np.mean(delays)
 
 
-
 class LIF(object):
 
     def __init__(self, stepsize=0.0001, offset=1.6, tau_m=0.010, tau_a=0.02, da=0.0, D=3.5):
@@ -416,24 +453,3 @@ class LIF(object):
 
     def __repr__(self):
         return self.__str__()
-        
-
-if __name__ == "__main__":
-   # stim = "gwn300Hz10s0.3.dat"
-   # time, stimulus = load_white_noise_stim(stim, 10, 40000)
-   # np.savez("noise_stimulus.npz", time=time, stim=stim)
-   
-   import matplotlib.pyplot as plt
-   from response_features import despine
-   fig = plt.figure(figsize=(1.5, 1.5))
-   ax = fig.add_subplot(111)
-   k = gaussKernel(0.125, 0.0001)
-   ax.plot(k)
-   ax.set_yticklabels([])
-   ax.set_xlim([0, 10000])
-   ax.set_xticks([0, 2500, 5000, 7500, 10000])
-   ax.set_xticklabels([r"-2$\pi$", r"-$\pi$" , "0", r"$\pi$", r"2$\pi$"])
-   despine(ax, ["top", "right"], False)
-   fig.subplots_adjust(bottom=0.2)
-   fig.savefig("gaussian.pdf") 
-   plt.close()

plot_figures.py

@@ -8,7 +8,7 @@ from code.plots.figure5 import command_line_parser as fig5_parser
 from code.plots.figure6 import command_line_parser as fig6_parser
 from code.plots.supp_figure4 import command_line_parser as supfig4_parser
 from code.plots.supp_figure5 import command_line_parser as supfig5_parser
-
+from code.plots.figure1new import command_line_parser as fig1n_parser
 
 def create_parser():
     parser = argparse.ArgumentParser(description="Tool for plotting figures of the Hladnik & Grewe population coding project.")
@@ -17,6 +17,7 @@ def create_parser():
                                        description="")
 
     fig1_parser(subparsers)
+    fig1n_parser(subparsers)
     fig2_parser(subparsers)
     fig3_parser(subparsers)
     fig4_parser(subparsers)