<jupyter_start><jupyter_text>Parameter exploration<jupyter_code>parametrization = pp.cartesian_product({
     'mue_ext_mean':[ round(elem, 3) for elem in np.linspace(0, 7.0, 71)],
     'mui_ext_mean':[ round(elem, 3) for elem in np.linspace(0, 7.0, 71)],
    'a':[0.0, 7.5, 15.0, 30.0],
    'b':[0.0, 20.0, 40.0, 80.0],
})

print("{} parameterizations".format(len(parametrization[parametrization.keys()[0]])))<jupyter_output><empty_output><jupyter_text>Run simulations<jupyter_code># ---- initialize pypet environment ----
trajectoryName = 'results' + datetime.datetime.now().strftime("-%Y-%m-%d-%HH-%MM-%SS")
HDF_FILE = os.path.join(paths.HDF_DIR, 'aln-grid-bifurcation-adaptation.hdf')
trajectoryFileName = HDF_FILE

import multiprocessing
ncores = multiprocessing.cpu_count()
print("Number of cores: {}".format(ncores))



env = pp.Environment(trajectory=trajectoryName,filename=trajectoryFileName,
                    file_title='aln bifurcation diagrams with increasing adaptation',
                    large_overview_tables=True,
                    multiproc=True,           
                    ncores=ncores,
                    wrap_mode='QUEUE',
                    log_stdout = False
                    )

# Get the trajectory from the environment
traj = env.v_trajectory
trajectoryName = traj.v_name

pe.add_parameters(traj, params)
traj.f_explore(parametrization)
env.f_run(rm.runModels_stimulus)
env.f_disable_logging()
print("Done.")<jupyter_output><empty_output><jupyter_text>Data processing  Load results from disk<jupyter_code>HDF_FILE = os.path.join(paths.HDF_DIR, 'aln-grid-bifurcation-adaptation.hdf')
trajectoryFileName = HDF_FILE
# ---- load pypet trajectory "trajectoryFileName" ----
print("Analyzing File \"{}\"".format(trajectoryFileName))
print("All Trajectories:")
print(pe.getTrajectoryNameInsideFile(trajectoryFileName))
trajectoryName = pe.getTrajectoryNameInsideFile(trajectoryFileName)[-1]

print("Analyzing trajectory \"{}\".".format(trajectoryName))
trajLoaded = pp.Trajectory(trajectoryName,add_time=False)
trajLoaded.f_load(trajectoryName,filename=trajectoryFileName, force=True)
trajLoaded.v_auto_load = True
print("{} results found".format(len(trajLoaded.f_get_results())))<jupyter_output><empty_output><jupyter_text>Process data<jupyter_code>nResults = len(trajLoaded.f_get_run_names()) # number of results in trajectory 
dt = trajLoaded.f_get_parameters()['parameters.simulation.dt'].f_get() # dt of simulation

# ---- explored parameter list ----
exploredParameters = trajLoaded.f_get_explored_parameters()
niceParKeys = [p.split('.')[-1] for p in exploredParameters.keys()]

# ---- lcreate pandas df with results as keys ----
dfResults = pd.DataFrame(columns=niceParKeys,dtype=object)

# range of parameters
for nicep, p in zip(niceParKeys,exploredParameters.keys()):
    dfResults[nicep] = exploredParameters[p].f_get_range()<jupyter_output><empty_output><jupyter_text>Serial Processing<jupyter_code># ---- make a dictionary with results ----
resultDicts = []
makeDict = True

#measures = ['domfr_exc', 'domfr_power_exc', 'max_exc', 'min_exc', 'updowndiff_exc', 'domfr_inh','domfr_power_inh', 'max_inh', 'min_inh', 'updowndiff_inh'] 
measures = ['domfr_exc', 'max_exc', 'updowndiff_exc', 'domfr_power_exc', 'max_inh']

for rInd in tqdm.tqdm(range(nResults), total=nResults):
    
    for measure in measures:
        dfResults.loc[rInd, measure] = func.analyse_run(measure, trajLoaded.results[rInd].f_to_dict(), dt)    
    if makeDict:
            resultDicts.append(trajLoaded.results[rInd].f_to_dict())
            
print("done.")<jupyter_output><empty_output><jupyter_text>Parallel Processing<jupyter_code>from joblib import Parallel, delayed
import multiprocessing
num_cores = multiprocessing.cpu_count()   

# ---- make a dictionary with results ----
resultDicts = []
makeDict = False
print("Creating results dictionary ...")
runResults = []
for rInd in range(len(dfResults)):
    result = trajLoaded.results[rInd].f_to_dict()
    runResults.append(result)
    if makeDict:
        resultDicts.append(result)
print("done.")

for measure in ['domfr_exc', 'max_exc', 'updowndiff_exc', 'domfr_power_exc', 'max_inh']:
    print("Analysing measure: {}".format(measure))
    parallelResults = Parallel(n_jobs=num_cores)(delayed(func.analyse_run)(measure = measure, result = i, dt = dt) for i in runResults)
    dfResults[measure] = parallelResults
print("done.")
import pickle
file_path_p = os.path.join(paths.PICKLE_DIR, 'pypet-aln-grid-bifurcation-adaptation-dfResults.p')
#pickle.dump(dfResults, file(file_path_p, "w"))<jupyter_output><empty_output><jupyter_text>Load processed results from file<jupyter_code>import pickle
file_path_p = os.path.join(paths.PICKLE_DIR, 'pypet-aln-grid-bifurcation-adaptation-dfResults.p')
#dfResults = pickle.load(file(file_path_p, "r"))<jupyter_output><empty_output><jupyter_text>Plot data<jupyter_code>def annotate_plot(point, model, ax):
    fs = 16
    fw = 'regular'
    C = 1 # Capacitance, convert mV/ms to pA
    pad = 3
    if point[0] == 'A':
        if model == 'brian':
            # attractor labels -------------
            ax.text(0.8*C,0.25*C, 'up', fontweight=fw, transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
            ax.text(0.35*C,0.06*C, 'LC$_{EI}$', fontweight=fw, transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
            ax.text(0.25*C,0.65*C, 'down', fontweight=fw, transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
            ax.text(0.8*C,0.65*C, 'bi', fontweight=fw, transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)

        else:
            # attractor labels -------------
            ax.text(3.5*C,0.75*C, 'up', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            ax.text(0.7*C,0.3*C, 'LC$_{EI}$', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            ax.text(0.4*C,1.9*C, 'down', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            ax.text(2.3*C,2.7*C, 'bi', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
    elif point[0] == 'B':
        if model == 'brian':
            # attractor labels -------------
            ax.text(0.88*C,0.25*C, 'up', fontweight=fw, transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
            ax.text(0.5*C,0.06*C, 'LC$_{EI}$', fontweight=fw,transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
            ax.text(0.3*C,0.7*C, 'down', fontweight=fw,transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
            ax.text(0.68*C,0.39*C, 'LC$_{aE}$', fontweight=fw,transform=ax.transAxes, #style='italic',
                    bbox={'facecolor':'white', 'alpha':0.85, 'pad':pad}, fontsize=fs)
        else:
            # attractor labels -------------
            ax.text(2.5*C, 0.3*C, u'LC$_{EI}$', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            ax.text(4.0*C, 2.8*C, u'LC$_{aE}$', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            ax.text(1.6*C, 2.75*C, 'down', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            ax.text(5.0*C, 1.0*C, 'up', fontweight=fw,#style='italic',
                    bbox={'facecolor':'white', 'alpha':0.9, 'pad':pad}, fontsize=fs)
            
up_down_threshold_hz = 10 # minimum difrernce of population rate between up and down state for bistable classification
domfr_threshold = 0.5 # minimum dominant frequency of an oscillation for classification

def plotResults(selectedResults, plot='max_exc', ax=None, pos=(0, 0)):
    plot_results = []
    amplitudes = []
    frequencies = []
    oscillation_powers = []
    max_inh_rates = []

    up_down_differences = []
    
    for i in range(len(selectedResults)):
        plot_results.append(selectedResults.iloc[i][plot])
        frequencies.append(selectedResults.iloc[i]['domfr_exc'])
        up_down_differences.append(selectedResults.iloc[i]['updowndiff_exc'])
    
    mues = np.unique(selectedResults['mue_ext_mean'])  
    muis = np.unique(selectedResults['mui_ext_mean'])    
    
    # prepare matrixes for plotting
    plot_results_matrix = np.reshape(plot_results, (len(mues), len(muis))).T
    frequencies_matrix = np.reshape(frequencies, (len(mues), len(muis))).T
    up_down_differences_matrix = np.reshape(up_down_differences, (len(mues), len(muis))).T
    
    # Fill in NaN's...
    mask = np.isnan(plot_results_matrix)
    plot_results_matrix[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), plot_results_matrix[~mask])
    mask = np.isnan(frequencies_matrix)
    frequencies_matrix[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), frequencies_matrix[~mask])
    mask = np.isnan(up_down_differences_matrix)
    up_down_differences_matrix[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), up_down_differences_matrix[~mask])  
    
    # --------------------------
    # -------  PLOTTING  -------
    # --------------------------
    
    C = 1 # Capacitance, convert mV/ms to nA (depricated!!)
    
    if ax==None:
        print('new figure')
        plt.figure(figsize=(3,3), dpi=300)
        ax = plt.gca()
        
    cmap = 'BuPu_r'
    if pos[0]==pos[1]:
        cmap = 'plasma' # change color on the diagonal
        
    im = ax.imshow(plot_results_matrix, origin='lower', aspect='auto', clim = (0, 80), cmap=cmap,\
                   extent = [mues[0], mues[-1], muis[0], muis[-1]])
    
    ax.tick_params(labelsize=20) 
    
    if pos==(0, 0):
        ax.set_xlabel("Input to E [nA]", fontsize=20)
        ax.set_ylabel("Input to I [nA]", fontsize=20)
        # force integer ticks
        plt.locator_params(nbins=4)
        
        # multiply axis ticks with capacitance C to get nA 
        ax.set_xticklabels(np.round(np.multiply(ax.get_xticks(),0.2),2)) 
        ax.set_yticklabels(np.round(np.multiply(ax.get_yticks(),0.2),2)) 
        ax.tick_params(labelsize=15) 
        
        annotate_plot('A', 'aln', ax) # annotate attractors
        
    elif pos == (2, 2):
        annotate_plot('B', 'aln', ax) # annotate attractors
        ax.set_xticks([])
        ax.set_yticks([])
    else:
        ax.set_xticks([])
        ax.set_yticks([])
        
    # ----- LC contour -----
    Z=frequencies_matrix.copy()
    blurredZ =scipy.ndimage.filters.gaussian_filter(frequencies_matrix,0.0) # blurred
    Z = blurredZ.copy()
    frequency_threshold = domfr_threshold
    Z[Z<frequency_threshold] = 0
    Z[Z>=frequency_threshold] = 1

    Xi,Yi = np.meshgrid(mues, muis)
    cset2 = ax.contour(Xi, Yi, Z, colors='white', linestyles='solid', levels=[0, 1.0001], linewidths=(7,),zorder=1)
    
    # ----- bistability contour -----
    Z=up_down_differences_matrix.copy()
    blurredZ =scipy.ndimage.filters.gaussian_filter(up_down_differences_matrix,0.0) # blurred
    Z = blurredZ.copy()
    up_down_threshold = up_down_threshold_hz
    Z[Z<up_down_threshold] = 0
    Z[Z>=up_down_threshold] = 1

    Xi,Yi = np.meshgrid(mues, muis)
    levels = np.arange(0,1,1)
    cset2 = ax.contour(Xi, Yi, Z, colors='springgreen', linestyles='dashed', levels=[0, 1.0001], linewidths=(7,),zorder=2)
    

    if ax==None:
        plt.show()
    return im<jupyter_output><empty_output><jupyter_text>Plot individually<jupyter_code># Plot individually
a_range = np.unique(trajLoaded.f_get('parameters.neuron.a').f_get_range())
b_range = np.unique(trajLoaded.f_get('parameters.neuron.b').f_get_range())

for a in a_range:
    for b in b_range:
        print(a, b)
        selector = (dfResults.a==a) & (dfResults.b==b)
        selectedResults = dfResults[selector]
        selectIndices = selectedResults.index
        #plotResults(selectedResults, 'max_exc')
        #plotResults(selectedResults, 'domfr_exc')
        #plotResults(selectedResults, 'domfr_power_exc')<jupyter_output><empty_output><jupyter_text>Plot grid<jupyter_code>import matplotlib.gridspec as gridspec
import matplotlib.patheffects as PathEffects
import string
labels = [l for l in string.ascii_lowercase]


a_range = np.unique(trajLoaded.f_get('parameters.neuron.a').f_get_range())
b_range = np.unique(trajLoaded.f_get('parameters.neuron.b').f_get_range())

fig = plt.figure(figsize=(14, 14), dpi=300)
plot_grid = gridspec.GridSpec(len(a_range),len(b_range))



for yg, a in enumerate(a_range): #zip(range(len(a_range))[::-1], a_range): 
    yg_flipped = len(a_range)-1-yg # flip for plotting, because gridspec (0,0) is upper left...
    for xg, b in enumerate(b_range):
        print(yg, xg, 'pars', a, b)
        selector = (dfResults.a==a) & (dfResults.b==b)
        selectedResults = dfResults[selector]
        selectIndices = selectedResults.index
        
        ax = plt.subplot(plot_grid[yg_flipped, xg]) # stupid: gridspec takes y first, then x ...
        
        imPlot = plotResults(selectedResults, 'max_exc', ax, (xg, yg))
        
        
        if (xg==0) and (yg==0):
            im = imPlot # save imPlot for colorbars later
            origin_ax = ax # save ax for grid arrows a and b
        if (xg==1) and (yg==0):
            im_grey = imPlot # save imPlot for colorbars later
            
        # a, b parameter legend upper left
        frametext = "a={}\nb={}".format(a, b)
        txt = ax.text(0.2 , 5.7, frametext, fontsize=19, color='white')
        txt.set_path_effects([PathEffects.withStroke(linewidth=2, foreground='k')])
        
# adjust plot
fig.tight_layout(pad=2.5, w_pad=0.5, h_pad=0.5)     


# colorbars
from mpl_toolkits.axes_grid1 import make_axes_locatable
cax=fig.add_axes([0.245,0.105,0.018,0.15])  ## the parameters are the specified position you set 
fig.colorbar(im, cax=cax)
fig.text(0.247, 0.162, 'max $r_E$', rotation=90, color='white', size=16)

cax=fig.add_axes([0.925,0.105,0.018,0.15])  ## the parameters are the specified position you set 
fig.colorbar(im_grey, cax=cax)
fig.text(0.927, 0.162, 'max $r_E$', rotation=90, color='white', size=16)


# outer grid a and b
bbox_props = dict(boxstyle="square,pad=0.1", fc="white", lw=1)
bbox_props = None

plt.annotate("b",size=28, bbox=bbox_props,
            xy=(0.8, 0.028), xycoords=fig.transFigure,
            xytext=(0.3, 0.02), textcoords=fig.transFigure,                 
            arrowprops=dict(arrowstyle="->", lw=2))
plt.annotate("Spike-triggered adaptation [pA]", size=18, xy=(0.4, 0.04), xycoords=fig.transFigure)

plt.annotate("a",size=28, bbox=bbox_props,
            xy=(0.046, 0.8), xycoords=fig.transFigure,
            xytext=(0.038, 0.3), textcoords=fig.transFigure,                 
            arrowprops=dict(arrowstyle="->", lw=2))
plt.annotate("Subthreshold adaptation [nS]", size=18, xy=(0.02, 0.66), rotation=90, xycoords=fig.transFigure)

fname = 'aln-grid-adaptation-bifurcation'
print("Saving {}".format(fname))
for extension in ['png', 'svg', 'jpg']:
    plt.savefig(os.path.join(paths.FIGURES_DIR, "bifurcation-plot-{}.{}".format(fname, extension)))
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16649752632246606, 'fraction_numerical': 0.022833946467081062, 'mean_word_length': 3.4678379144233493, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '467069', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 5946, 'n_tokens_neox': 5618, 'n_words': 939, 'path': 'notebooks/fig3/aln-grid-adaptation-bifurcation.ipynb', 'repo_name': 'caglorithm/stimulus_neural_populations', 'size': '15766'}