diff --git a/figures/SchueckerSchmidt2017/Snakefile b/figures/SchueckerSchmidt2017/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..38ccf53e82c0b0262efbd23af0e96029e4a6cf36
--- /dev/null
+++ b/figures/SchueckerSchmidt2017/Snakefile
@@ -0,0 +1,75 @@
+rule all:
+ input:
+ 'iteration_1/K_prime.npy',
+ 'iteration_2/K_prime.npy',
+ 'iteration_3/K_prime.npy',
+ 'iteration_4/K_prime.npy',
+ 'Fig1_scheme.eps',
+# 'Fig2_EE_example.eps',
+ 'Fig3_bistability.eps',
+ 'Fig4_meanfield_mam.eps',
+ 'Fig5_eigenspace.eps',
+ 'Fig6_unstable_FP.eps',
+ 'Fig7_stabilization_analysis.eps',
+ 'Fig8_conn_analysis.eps'
+# 'Fig9_simulation.eps'
+
+rule Stabilization:
+ output:
+ 'iteration_1/K_prime.npy',
+ 'iteration_2/K_prime.npy',
+ 'iteration_3/K_prime.npy',
+ 'iteration_4/K_prime.npy'
+ shell:
+ 'python3 stabilization.py'
+
+rule Fig2_EE_example:
+ output:
+ 'Fig2_EE_example.eps'
+ shell:
+ 'python3 Fig2_EE_example.py'
+
+rule Fig3_bistability:
+ output:
+ 'Fig3_bistability.eps'
+ shell:
+ 'python3 Fig3_bistability.py'
+
+rule Fig4_meanfield_mam:
+ output:
+ 'Fig4_meanfield_mam.eps'
+ shell:
+ 'python3 Fig4_meanfield_mam.py'
+
+rule Fig5_eigenspace:
+ output:
+ 'Fig5_eigenspace.eps'
+ shell:
+ 'python3 Fig5_eigenspace.py'
+
+rule Fig6_unstable_FP:
+ output:
+ 'Fig6_unstable_FP.eps'
+ shell:
+ 'python3 Fig6_unstable_FP.py'
+
+rule Fig7_stabilization_analysis:
+ output:
+ 'Fig7_stabilization_analysis.eps'
+ shell:
+ 'python3 Fig7_stabilization_analysis.py'
+
+rule Fig8_conn_analysis:
+ output:
+ 'Fig8_conn_analysis.eps'
+ shell:
+ 'python3 Fig8_conn_analysis.py'
+
+
+rule Fig9_simulation:
+ output:
+ 'Fig9_simulation.eps'
+ shell:
+ 'python3 Fig9_simulation.py'
+
+
diff --git a/figures/SchueckerSchmidt2017/utils.py b/figures/SchueckerSchmidt2017/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..92f0933fb543c1264b779cc210c76152e4088c4b
--- /dev/null
+++ b/figures/SchueckerSchmidt2017/utils.py
@@ -0,0 +1,147 @@
+import copy
+import data_wrapper.data_wrapper as dw
+import multiprocessing as mp
+import numpy as np
+import os
+import pylab as pl
+from functools import partial
+
+from multiarea_model.multiarea_helpers import create_mask
+from scipy.signal import find_peaks_cwt
+
+
+def integrate(f5, M_base):
+ mask5 = create_mask(M_base.structure, target_pops=['5E'],
+ source_areas=[], external=True)
+ mask6 = create_mask(M_base.structure, target_pops=['6E'],
+ source_areas=[], external=True)
+
+ f6 = 10/3.*f5-7/3.
+ conn_params = copy.deepcopy(M_base.params['connection_params'])
+ conn_params.update({'fac_nu_ext_5E': f5,
+ 'fac_nu_ext_6E': f6})
+ M = copy.deepcopy(M_base)
+ M.K_matrix[mask5] *= f5
+ M.K_matrix[mask6] *= f6
+ p, r = M.theory.integrate_siegert()
+ return r
+
+
+def iteration_results(fac_nu_ext_5E_list, M_base, threads=None):
+ if threads is None:
+ results = np.array([integrate(f5, M_base) for f5 in fac_nu_ext_5E_list])
+ else:
+ integrate_part = partial(integrate,
+ M_base=M_base)
+ pool = mp.Pool(processes=threads)
+ results = np.array(pool.map(integrate_part, fac_nu_ext_5E_list))
+ pool.close()
+ pool.join()
+ return results
+
+
+def velocity_peaks(time, result, threshold=0.05):
+ d_nu = np.abs(np.diff(np.mean(result, axis=0)) / np.diff(time)[0])
+ ind = np.where(d_nu < threshold)
+ minima = find_peaks_cwt(-1. * np.log(d_nu[ind]), np.array([0.1]))
+ if len(minima) > 0:
+ t_min = time[ind][minima]
+ else:
+ t_min = []
+ min_full = [np.argmin(np.abs(time - t)) for t in t_min]
+ return d_nu, min_full
+
+
+def plot_iteration(results, theory_params, threshold=0.05, full=True):
+ import seaborn as sns
+
+ traj = np.mean(results, axis=1)
+ if full:
+ ind = list(range(0, len(traj)))
+ else:
+ i = np.argmax(np.diff(traj[:, -1]))
+ ind = [i, i+1]
+ cp = sns.color_palette(n_colors=len(ind))
+
+ time = np.arange(0., theory_params['T'], theory_params['dt'])
+ fig = pl.figure()
+ ax = fig.add_subplot(121)
+ [ax.plot(time,
+ traj[i]) for i in ind]
+
+ ax.set_yscale('Log')
+ ax.set_xlabel('Time (a.u.)')
+ ax.set_ylabel(r'$\langle \nu \rangle$')
+
+ ax = fig.add_subplot(122)
+ for n, i in enumerate(ind):
+ d_nu, minima = velocity_peaks(time, results[i], threshold=threshold)
+ ax.plot(time[:-1], d_nu, color=cp[n])
+ if not full:
+ ax.vlines(time[minima],
+ 1e-6,
+ 1e0,
+ linestyles='dashed',
+ color=cp[n])
+ ax.set_yscale('Log')
+ pl.show()
+
+
+def save_iteration(step, data):
+ data_dir = 'iteration_{}'.format(step)
+ try:
+ os.mkdir(data_dir)
+ except FileExistsError:
+ pass
+ dw.save('iteration_{}'.format(step), data)
+
+
+def load_iteration(step):
+ data_dir = 'iteration_{}'.format(step)
+ data = {}
+ files = os.listdir(data_dir)
+ for f in files:
+ data[os.path.splitext(f)[0]] = np.load(os.path.join(data_dir, f))
+ return data
+
+
+def compute_iteration(max_iter, fac_nu_ext_5E_list, theory_params, M_base, threads=None):
+ par_list = np.zeros(0)
+ results = np.zeros((0, 254, int(theory_params['T'] / theory_params['dt'])))
+
+ i = 0
+ while i < max_iter:
+ print("Iteration: {}".format(i))
+ print(fac_nu_ext_5E_list)
+ r = iteration_results(fac_nu_ext_5E_list, M_base, threads=threads)
+ results = np.vstack((results, r))
+ par_list = np.append(par_list, fac_nu_ext_5E_list)
+ i = np.argmax(np.diff(np.mean(r, axis=1)[:, -1]))
+ i += 1
+ fac_nu_ext_5E_list = np.arange(fac_nu_ext_5E_list[i],
+ # to ensure that the array includes the last value, we add a small epsilon
+ fac_nu_ext_5E_list[i+1] + 1.e-10,
+ 10**(-(i+2)))
+
+ ind = np.argsort(par_list)
+ par_list = par_list[ind]
+ results = results[ind]
+
+ return {'results': results, 'parameters': par_list}
+
+
+def determine_velocity_minima(time, data, threshold=0.05):
+
+ traj = np.mean(data['results'], axis=1)
+ i = np.argmax(np.diff(traj[:, -1]))
+ r_low = data['results'][i]
+ r_high = data['results'][i+1]
+
+ dnu_low, minima_low = velocity_peaks(time, r_low, threshold=threshold)
+ dnu_high, minima_high = velocity_peaks(time, r_high, threshold=threshold)
+
+ par_transition = data['parameters'][i]
+ return par_transition, r_low, r_high, minima_low, minima_high
+
+
+