{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0b85a598",
   "metadata": {},
   "outputs": [],
   "source": [
    "from brian2 import *\n",
    "\n",
    "from adex_sine import *\n",
    "\n",
    "defaultclock.dt = 10 * us\n",
    "\n",
    "\n",
    "class Model:\n",
    "    # Model type flags\n",
    "    SINE = 0\n",
    "    EXP2SYN = 1\n",
    "\n",
    "    # Noise flags\n",
    "    HIGH = 0\n",
    "    LOW = 1\n",
    "    OFF = -1\n",
    "\n",
    "    # Noise parameters\n",
    "    EXCITATORY_NOISE_VARIANCE = {HIGH: 0.5*2 * nS, LOW: 0.25 * nS, OFF: 0 * nS}\n",
    "    INHIBITORY_NOISE_VARIANCE = {HIGH: 1.25*2 * nS, LOW: 0.625 * nS, OFF: 0 * nS}\n",
    "\n",
    "    # Noise mean conductance\n",
    "    EXCITATORY_CONDUCTANCE = 1 * nS\n",
    "    INHIBITORY_CONDUCTANCE = 4 * nS\n",
    "\n",
    "    DEFAULT_PARAMETERS = {\n",
    "        \"sigma_flux\" : 6.75*pA,     \n",
    "        \"c\": 85 * pF,\n",
    "        \"tau_w\": 18 * ms,\n",
    "        \"b\": 0.25 * nA,\n",
    "        \"a\": 1.3 * nS,\n",
    "        \"v_T\": -45 * mV,\n",
    "        \"v_thresh\": 0 * mV,\n",
    "        \"DeltaT\": 0.2 * mV,\n",
    "        # EQUILIBRIUM POTENTIAL\n",
    "        \"e_l\": -65 * mV,\n",
    "        \"e_ex\": 0 * mV,\n",
    "        \"e_in\": -70 * mV,\n",
    "        # CONDUCTANCES\n",
    "        \"g_l\": 3 * nS,\n",
    "        \"mu_ex\": 0 * nS,\n",
    "        \"mu_in\": 0 * nS,\n",
    "        # EXCITATORY NOISE\n",
    "        \"sigma_ex\": 0 * nS,\n",
    "        \"tau_noise_ex\": 3 * ms,\n",
    "        # INHIBITORY NOISE\n",
    "        \"sigma_in\": 0 * nS,\n",
    "        \"tau_noise_in\": 10 * ms,\n",
    "        # SINE INPUT\n",
    "        \"f\": 100 * Hz,\n",
    "        \"A\": 0 * pA,\n",
    "        \"i_injected\": 0 * pA,\n",
    "        \"v_reset\": -70 * mV,\n",
    "        # m current\n",
    "        \"g_adapt\": 10 * nS,\n",
    "        \"e_k\": -90*mV,\n",
    "        \"beta_z\": -35*mV,\n",
    "        \"gamma_z\": 4*mV,     #5\n",
    "        \"tau_z\": 100*ms,\n",
    "    }\n",
    "\n",
    "    def __init__(\n",
    "        self, n, *, stim=None, noise=None, resistance=None, additional_vars=()\n",
    "    ):\n",
    "        if resistance is None:\n",
    "            raise ValueError(\"Resistance must be specified\")\n",
    "\n",
    "        if noise is None:\n",
    "            raise ValueError(\"Noise must be specified\")\n",
    "\n",
    "        self.stim_type = stim\n",
    "        self._input_resistance = None\n",
    "        self._noise_level = None\n",
    "        self._duration = 0\n",
    "        self.recorded_vars = (\"v\",) + additional_vars\n",
    "\n",
    "        self.neurons = self.set_default(n_neuron=n)\n",
    "        self.set_resistance(resistance)\n",
    "        self.set_noise(noise)\n",
    "\n",
    "        self.spikes = None\n",
    "        self.spiker = None\n",
    "        self.synapses = None\n",
    "        self.inhib_synapses = None\n",
    "        self.smon = None\n",
    "        self.network = None\n",
    "        self.build_network()\n",
    "\n",
    "    def create_model(self):\n",
    "        return ADEX_MODEL, self.DEFAULT_PARAMETERS\n",
    "\n",
    "    def set_default(self, n_neuron):\n",
    "        model, parameters = self.create_model()\n",
    "\n",
    "        neurons = NeuronGroup(\n",
    "            n_neuron,\n",
    "            model=model,\n",
    "            method=\"Euler\",\n",
    "            name=\"neurons\",\n",
    "            threshold=\"v > v_thresh\",\n",
    "            reset=\"v = v_reset; w += b\",\n",
    "        )\n",
    "\n",
    "        for parameter, value in parameters.items():\n",
    "            neurons.__setattr__(parameter, value)\n",
    "\n",
    "        neurons.v = neurons.e_l  # remove most of transient\n",
    "\n",
    "        return neurons\n",
    "\n",
    "    def set_resistance(self, level):\n",
    "        if level == self.LOW:\n",
    "            exc_conductance = self.EXCITATORY_CONDUCTANCE\n",
    "            inhib_conductance = self.INHIBITORY_CONDUCTANCE\n",
    "\n",
    "        else:\n",
    "            exc_conductance = inhib_conductance = 0\n",
    "\n",
    "        self._input_resistance = level\n",
    "        self._set_variable(\"mu_ex\", exc_conductance)\n",
    "        self._set_variable(\"mu_in\", inhib_conductance)\n",
    "\n",
    "    def set_noise(self, level):\n",
    "        if level == self.HIGH or level == self.LOW:\n",
    "            exc_noise = self.EXCITATORY_NOISE_VARIANCE[level]\n",
    "            inhib_noise = self.INHIBITORY_NOISE_VARIANCE[level]\n",
    "\n",
    "        else:\n",
    "            exc_noise = inhib_noise = 0\n",
    "\n",
    "        self._noise_level = level\n",
    "        self._set_variable(\"sigma_ex\", exc_noise)\n",
    "        self._set_variable(\"sigma_in\", inhib_noise)\n",
    "\n",
    "    def set_injected_current(self, amplitude):\n",
    "        self._set_variable(\"i_injected\", amplitude)\n",
    "        self._set_variable(\"A\", 0 * pA)\n",
    "\n",
    "    def set_stimulus_current(self, amplitude):\n",
    "        self._set_variable(\"A\", amplitude)\n",
    "        self._set_variable(\"i_injected\", 0 * pA)\n",
    "\n",
    "    @property\n",
    "    def f(self):\n",
    "        return self.neurons.f\n",
    "\n",
    "    @f.setter\n",
    "    def f(self, new_f):\n",
    "        self._set_variable(\"f\", new_f)  # this will reset smon\n",
    "        if self.stim_type == self.EXP2SYN:\n",
    "            self.spiker.T = 1 / new_f\n",
    "\n",
    "    def run(self, duration, report=\"stdout\"):\n",
    "        self._duration = duration\n",
    "        self.network.run(duration, report=report)\n",
    "\n",
    "    def build_network(self):\n",
    "        self.smon = StateMonitor(\n",
    "            self.neurons, self.recorded_vars, record=True, name=\"smon\"\n",
    "        )\n",
    "        self.spikes = SpikeMonitor(self.neurons, name=\"spikes\")\n",
    "\n",
    "        self.network = Network(self.neurons, self.smon, self.spikes)\n",
    "\n",
    "    def _set_variable(self, name, value):\n",
    "        self.neurons.__setattr__(name, value)\n",
    "        self.reset_recording()\n",
    "\n",
    "    def reset_recording(self):\n",
    "        try:\n",
    "            self.network\n",
    "        except AttributeError:\n",
    "            return  # network not yet initialized\n",
    "\n",
    "        self.network.remove(self.smon, self.spikes)\n",
    "\n",
    "        self.smon = StateMonitor(\n",
    "            self.neurons, self.recorded_vars, record=True, name=\"smon\"\n",
    "        )\n",
    "        self.spikes = SpikeMonitor(self.neurons, name=\"spikes\")\n",
    "\n",
    "        self.network.add(self.smon, self.spikes)\n",
    "\n",
    "    @property\n",
    "    def spike_train(self):\n",
    "        return self.spikes.spike_trains()\n",
    "\n",
    "    @property\n",
    "    def firing_rate(self):\n",
    "        return self.spikes.count / self.duration\n",
    "\n",
    "    @property\n",
    "    def duration(self):\n",
    "        return self._duration\n",
    "\n",
    "    @property\n",
    "    def input_resistance(self):\n",
    "        if self._input_resistance == self.HIGH:\n",
    "            return \"HIGH\"\n",
    "        else:\n",
    "            return \"LOW\"\n",
    "\n",
    "    @property\n",
    "    def noise_level(self):\n",
    "        if self._noise_level == self.HIGH:\n",
    "            return \"HIGH\"\n",
    "        elif self._noise_level == self.LOW:\n",
    "            return \"LOW\"\n",
    "        else:\n",
    "            return \"NO\"\n",
    "\n",
    "    def __repr__(self):\n",
    "        return f\"{self.neurons.N} Neurons with {self.input_resistance} input resistance and {self.noise_level} noise\"\n",
    "\n",
    "    def __str__(self):\n",
    "        return self.__repr__()\n",
    "\n",
    "    def store(self, name):\n",
    "        self.network.store(name)\n",
    "\n",
    "    def restore(self, name):\n",
    "        self.network.restore(name)\n",
    "\n",
    "    @property\n",
    "    def v(self):\n",
    "        return self.smon.v\n",
    "\n",
    "    @property\n",
    "    def t(self):\n",
    "        return self.smon.t\n",
    "\n",
    "    @property\n",
    "    def injected_current(self):\n",
    "        return self.neurons.i_injected\n",
    "\n",
    "    @property\n",
    "    def stimulus_amplitude(self):\n",
    "        return self.neurons.A\n",
    "\n",
    "\n",
    "class CurrentModel(Model):\n",
    "    def __init__(self, **kwargs):\n",
    "        super().__init__(stim=self.SINE, **kwargs)\n",
    "\n",
    "    def create_model(self):\n",
    "        model, parameters = super().create_model()\n",
    "        model += CURRENT_INPUT\n",
    "\n",
    "        return model, parameters\n",
    "\n",
    "\n",
    "class SineModel(CurrentModel):\n",
    "    def create_model(self):\n",
    "        model, parameters = super().create_model()\n",
    "        model += SINE_INPUT\n",
    "\n",
    "        return model, parameters\n",
    "\n",
    "\n",
    "class SawModel(CurrentModel):\n",
    "    def create_model(self):\n",
    "        model, parameters = super().create_model()\n",
    "        model += SAW_INPUT\n",
    "\n",
    "        return model, parameters\n",
    "\n",
    "\n",
    "class SynapticModel(Model):\n",
    "    def __init__(self, **kwargs):\n",
    "        super().__init__(stim=self.EXP2SYN, **kwargs)\n",
    "\n",
    "    SYNAPTIC_PARAMETERS = {\n",
    "        \"tau_input_1\": 0.4 * ms,\n",
    "        \"tau_input_2\": 4 * ms,\n",
    "        \"offset_A\": 1.48793507e-11,\n",
    "        \"offset_B\": -2.66359562e-08,\n",
    "        \"offset_C\": 1.77538800e-05,\n",
    "        \"offset_D\": -8.05925810e-04,\n",
    "        \"offset_E\": -3.51463644e-02,\n",
    "        \"offset_switch\": 0,\n",
    "    }\n",
    "\n",
    "    def create_model(self):\n",
    "        model, parameters = super().create_model()\n",
    "        model += EXP2SYN_WAVEFORM + SUMMATION_OFFSET\n",
    "        parameters = {**parameters, **self.SYNAPTIC_PARAMETERS}\n",
    "\n",
    "        return model, parameters\n",
    "\n",
    "    def build_network(self):\n",
    "        super().build_network()\n",
    "        self.spiker = NeuronGroup(\n",
    "            self.neurons.N,\n",
    "            \"\"\"T : second (constant)\n",
    "                                     lastspike : second\"\"\",\n",
    "            threshold=\"timestep(t-lastspike, dt)>=timestep(T, dt)\",\n",
    "            reset=\"lastspike=t\",\n",
    "        )\n",
    "        self.spiker.T = 1 / self.neurons.f\n",
    "        self.synapses = Synapses(\n",
    "            self.spiker, self.neurons, on_pre=\"input_aux += 1\"\n",
    "        )  # connect input to neurons\n",
    "        self.synapses.connect(\"i==j\")  # one synapse goes to one neuron\n",
    "\n",
    "        self.network.add(self.spiker, self.synapses)\n",
    "\n",
    "\n",
    "class SynapticCurrentModel(SynapticModel):\n",
    "    def __init__(self, offset=True, **kwargs):\n",
    "        self.offset = 1 if offset else 0\n",
    "        super().__init__(**kwargs)\n",
    "\n",
    "    def create_model(self):\n",
    "        model, parameters = super().create_model()\n",
    "        model += CURRENT_INPUT + SYNAPTIC_INPUT_CURRENT\n",
    "        parameters = {**parameters, **{\"offset_switch\": self.offset}}\n",
    "\n",
    "        return model, parameters\n",
    "\n",
    "\n",
    "class SynapticConductanceModel(SynapticModel):\n",
    "    FLAT = 0\n",
    "    ACTIVE = 1\n",
    "\n",
    "    CONDUCTANCE_PARAMETERS = {\n",
    "        \"A\": 0 * nS,  # overwrite A to be conductance\n",
    "        \"g_i\": 1 * nS,\n",
    "    }\n",
    "\n",
    "    INHIBITION_PARAMETERS = {\n",
    "        \"tau_inhibition_1\": 1 * ms,\n",
    "        \"tau_inhibition_2\": 10 * ms,\n",
    "    }\n",
    "\n",
    "    def __init__(self, offset=ACTIVE, **kwargs):\n",
    "        self.offset = offset\n",
    "        super().__init__(**kwargs)\n",
    "\n",
    "    def create_model(self):\n",
    "        model, parameters = super().create_model()\n",
    "        if self.offset == self.FLAT:\n",
    "            model += CONDUCTANCE_INPUT + SYNAPTIC_CONDUCTANCE_FLAT\n",
    "            parameters = {\n",
    "                **parameters,\n",
    "                **self.SYNAPTIC_PARAMETERS,\n",
    "                **self.CONDUCTANCE_PARAMETERS,\n",
    "                **{\"offset_switch\": 1},\n",
    "            }\n",
    "\n",
    "        elif self.offset == self.ACTIVE:\n",
    "            model += CONDUCTANCE_INPUT + SYNAPTIC_CONDUCTANCE_STIM\n",
    "            parameters = {\n",
    "                **parameters,\n",
    "                **self.SYNAPTIC_PARAMETERS,\n",
    "                **self.CONDUCTANCE_PARAMETERS,\n",
    "                **self.INHIBITION_PARAMETERS,\n",
    "            }\n",
    "\n",
    "        return model, parameters\n",
    "\n",
    "    def build_network(self):\n",
    "        super().build_network()\n",
    "        if self.offset != self.ACTIVE:\n",
    "            return\n",
    "\n",
    "        self.inhib_synapses = Synapses(\n",
    "            self.spiker, self.neurons, on_pre=\"input_inhib_aux += 1\", delay=2 * ms\n",
    "        )  # connect input to neurons\n",
    "        self.inhib_synapses.connect(\"i==j\")  # one synapse goes to one neuron\n",
    "\n",
    "        self.network.add(self.inhib_synapses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c82769f2",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING    Cannot use Cython, a test compilation failed: Microsoft Visual C++ 14.0 or greater is required. Get it with \"Microsoft C++ Build Tools\": https://visualstudio.microsoft.com/visual-cpp-build-tools/ (DistutilsPlatformError) [brian2.codegen.runtime.cython_rt.cython_rt.failed_compile_test]\n",
      "INFO       Cannot use compiled code, falling back to the numpy code generation target. Note that this will likely be slower than using compiled code. Set the code generation to numpy manually to avoid this message:\n",
      "prefs.codegen.target = \"numpy\" [brian2.devices.device.codegen_fallback]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting simulation at t=0. s for a duration of 12. s\n",
      "185.41 ms (1%) simulated in 10s, estimated 10m 37s remaining.\n",
      "0.36864 s (3%) simulated in 20s, estimated 10m 31s remaining.\n",
      "0.582 s (4%) simulated in 30s, estimated 9m 49s remaining.\n",
      "0.76694 s (6%) simulated in 40s, estimated 9m 46s remaining.\n",
      "0.95679 s (7%) simulated in 50s, estimated 9m 37s remaining.\n",
      "1.14061 s (9%) simulated in 1m 0s, estimated 9m 31s remaining.\n",
      "1.33212 s (11%) simulated in 1m 10s, estimated 9m 21s remaining.\n",
      "1.51794 s (12%) simulated in 1m 20s, estimated 9m 13s remaining.\n",
      "1.70214 s (14%) simulated in 1m 30s, estimated 9m 5s remaining.\n",
      "1.88638 s (15%) simulated in 1m 40s, estimated 8m 56s remaining.\n",
      "2.06909 s (17%) simulated in 1m 50s, estimated 8m 48s remaining.\n",
      "2.25583 s (18%) simulated in 2m 0s, estimated 8m 39s remaining.\n",
      "2.44241 s (20%) simulated in 2m 10s, estimated 8m 29s remaining.\n",
      "2.62839 s (21%) simulated in 2m 20s, estimated 8m 19s remaining.\n",
      "2.81246 s (23%) simulated in 2m 30s, estimated 8m 10s remaining.\n",
      "2.99237 s (24%) simulated in 2m 40s, estimated 8m 2s remaining.\n",
      "3.17645 s (26%) simulated in 2m 50s, estimated 7m 52s remaining.\n",
      "3.36452 s (28%) simulated in 3m 0s, estimated 7m 42s remaining.\n",
      "3.55011 s (29%) simulated in 3m 10s, estimated 7m 32s remaining.\n",
      "3.73966 s (31%) simulated in 3m 20s, estimated 7m 22s remaining.\n",
      "3.92144 s (32%) simulated in 3m 30s, estimated 7m 13s remaining.\n",
      "4.10594 s (34%) simulated in 3m 40s, estimated 7m 3s remaining.\n",
      "4.29375 s (35%) simulated in 3m 50s, estimated 6m 53s remaining.\n",
      "4.47866 s (37%) simulated in 4m 0s, estimated 6m 43s remaining.\n",
      "4.66617 s (38%) simulated in 4m 10s, estimated 6m 33s remaining.\n",
      "4.84779 s (40%) simulated in 4m 20s, estimated 6m 24s remaining.\n",
      "5.03443 s (41%) simulated in 4m 30s, estimated 6m 14s remaining.\n",
      "5.21888 s (43%) simulated in 4m 40s, estimated 6m 4s remaining.\n",
      "5.4002 s (45%) simulated in 4m 50s, estimated 5m 55s remaining.\n",
      "5.58733 s (46%) simulated in 5m 0s, estimated 5m 44s remaining.\n",
      "5.77236 s (48%) simulated in 5m 10s, estimated 5m 35s remaining.\n",
      "5.95877 s (49%) simulated in 5m 20s, estimated 5m 25s remaining.\n",
      "6.146 s (51%) simulated in 5m 30s, estimated 5m 14s remaining.\n",
      "6.3257 s (52%) simulated in 5m 40s, estimated 5m 5s remaining.\n",
      "6.51209 s (54%) simulated in 5m 50s, estimated 4m 55s remaining.\n",
      "6.69868 s (55%) simulated in 6m 0s, estimated 4m 45s remaining.\n",
      "6.88664 s (57%) simulated in 6m 10s, estimated 4m 35s remaining.\n",
      "7.07488 s (58%) simulated in 6m 20s, estimated 4m 25s remaining.\n",
      "7.26196 s (60%) simulated in 6m 30s, estimated 4m 15s remaining.\n",
      "7.44808 s (62%) simulated in 6m 40s, estimated 4m 5s remaining.\n",
      "7.63531 s (63%) simulated in 6m 50s, estimated 3m 54s remaining.\n",
      "7.82288 s (65%) simulated in 7m 0s, estimated 3m 44s remaining.\n",
      "8.02036 s (66%) simulated in 7m 10s, estimated 3m 33s remaining.\n",
      "8.20939 s (68%) simulated in 7m 20s, estimated 3m 23s remaining.\n",
      "8.39229 s (69%) simulated in 7m 30s, estimated 3m 14s remaining.\n",
      "8.5773 s (71%) simulated in 7m 40s, estimated 3m 4s remaining.\n",
      "8.76134 s (73%) simulated in 7m 50s, estimated 2m 54s remaining.\n",
      "8.95375 s (74%) simulated in 8m 0s, estimated 2m 43s remaining.\n",
      "9.147 s (76%) simulated in 8m 10s, estimated 2m 33s remaining.\n",
      "9.33577 s (77%) simulated in 8m 20s, estimated 2m 23s remaining.\n",
      "9.52439 s (79%) simulated in 8m 30s, estimated 2m 13s remaining.\n",
      "9.71086 s (80%) simulated in 8m 40s, estimated 2m 3s remaining.\n",
      "9.89765 s (82%) simulated in 8m 50s, estimated 1m 53s remaining.\n",
      "10.07797 s (83%) simulated in 9m 0s, estimated 1m 43s remaining.\n",
      "10.26109 s (85%) simulated in 9m 10s, estimated 1m 33s remaining.\n",
      "10.44746 s (87%) simulated in 9m 20s, estimated 1m 23s remaining.\n",
      "10.62559 s (88%) simulated in 9m 30s, estimated 1m 14s remaining.\n",
      "10.81089 s (90%) simulated in 9m 40s, estimated 1m 4s remaining.\n",
      "10.99173 s (91%) simulated in 9m 50s, estimated 54s remaining.\n",
      "11.18646 s (93%) simulated in 10m 0s, estimated 44s remaining.\n",
      "11.37715 s (94%) simulated in 10m 10s, estimated 33s remaining.\n",
      "11.61595 s (96%) simulated in 10m 20s, estimated 21s remaining.\n",
      "11.80217 s (98%) simulated in 10m 30s, estimated 11s remaining.\n",
      "11.98574 s (99%) simulated in 10m 40s, estimated 1s remaining.\n",
      "12. s (100%) simulated in 10m 40s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.collections.EventCollection at 0x2401d312160>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d31d670>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d32ecd0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d33f220>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d34e880>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d35ab80>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d36c2b0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d37a940>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d38e1c0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d39a850>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d3a8df0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d3b9280>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d3c9a90>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d3d5f40>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d3e7610>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d3f5dc0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d409400>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4168e0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d422ee0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4364f0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d445c10>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d456280>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d463760>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d471f40>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d483580>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d492880>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d49fca0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4b0220>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4bf910>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4ccf10>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4de5b0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4eca60>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d4fafa0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d50b640>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d519d00>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d52c100>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d53a700>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d547ca0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d558100>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d566550>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d573820>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d582d90>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5940a0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5a15e0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5afc40>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5c13d0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5d0910>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5dcfd0>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5ef610>,\n",
       " <matplotlib.collections.EventCollection at 0x2401d5fdd90>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD4CAYAAAAaT9YAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoBElEQVR4nO2db4xn13nXvw+eZNNNu8RL1tZgR0xSWRarrkjKKMoQhEJcd0KJYoMIm4jAIgweCZBSitTY5AXTV7gFVRVComu1gRVNkzFtgy23sLXcWqjSysmYphlnbdduu01Mhngao25aS6bbHl7MPeNnnnnOOc+9v5nZPdPvR1rNveeeP89z7rlH+/ue554rKSUQQgjpkz9zvQ0ghBAyHU7ihBDSMZzECSGkYziJE0JIx3ASJ4SQjpk7zMbe/va3p4WFhcNskhBCuueZZ575vZTSKe/aoU7iCwsLWF9fP8wmCSGke0Tkd0vXKKcQQkjHcBInhJCO4SROCCEdw0mcEEI6hpM4IYR0TCg6RUSuAPg2gD8GcC2ltCgiJwGsAVgAcAXA300p/d+DMZMQQojHmP+J//WU0rtTSovD+QMAnkwp3QHgyeGcEELIITKLnHIPgAvD8QUA985sDSGEkFFEJ/EE4JdF5BkRuX9IuzWltAkAw99bvIIicr+IrIvI+tbW1uwWE0II2SH6xub7U0rfEJFbADwhIs9HG0gpPQzgYQBYXFzkFygIIWQfCU3iKaVvDH9fEZEvAHgvgG+KyHxKaVNE5gG8coB27nD2/CVc3rwKADg9fwIAcHnz6q7jTE7TrK0s4czqxV3X11aWcPb8pZ08tTpyW2srSzv2RGnZqa97Zb38pbQxeQDs8kf7qM+t//aaV49tp+WbR609bZNnnyZyn0t9U+qHUpmSf5F+9HzIx+tXXsXxY3M4PX/CrSPiT9T36Bht2eGNLa98pF9b46Dlu7UxYp9XrtR3tj3bPwdBU04RkbeKyHflYwDfD+BZAI8BODdkOwfg0QOxkBBCSJHI/8RvBfAFEcn5fzal9D9E5EsAHhGR+wB8DcBHD85MQgghHs1JPKX02wD+kpP+LQB3HYRRhBBCYshhfu1+cXExTdmKNqoZ5rwZrR9mtBZmj6NaWEvDtrrdGJ3YanU1HU6XKWl7Nfttedt+zZ6axlrqo5qdpba9fsr5NlaXd9LsOkfJL5vH086trVPsjOaP9n3JXm89oaV/l+qo6fVevSVfZ1mfifZV7oszqxfx2uvXsLhwsmlvawwcpAafj6ciIs+od3R2wdfuCSGkYziJE0JIx3ASJ4SQjulmEvfiLK3Wq8/XVpawuHDS1atq+mUum9Oz1rqxurynrNZk11aWdsqdWb24097lzas4e/7Sji5Waj/if0l/1+meHq7LZ/1f91P2o6XxeTpqrkv3T8nOjdXlHZs8HVrbl8+1nZ5N2i47Rmxd1h7gDS1d+wK8cW+11l7SPD07dZ5anLzu69dev7bLFu+49F6C7Xttj73mrbNY/718pTI13/R48+yzfaPz6L712jl7/hJOz5/A4sLJnefMQ5e39QPAdz/4i1i/8ipOz5/Y8yxY9HNu061/ev44qBhxoKNJnBBCyF44iRNCSMdwEieEkI7pIk4c2KvP1eJS7XV9nonE19q8XoytzjOmzVbsb0u3r9lZ82lqTG9tr5haDH9rj4op/dmKXa+NCW8/EttmyfYxMfvA9nsKALC4cDLk29i9XWox7CXbav604qCjdllKcdiR2GudXiun27FlomM311HKW7KxVoe1YSqMEyeEkCMKJ3FCCOkYTuKEENIxXWjiJW01ci3j6VatvbNt2dI+CNG9nCN6n6WmG2paeqFXvqbDtnTLlo4Y2TejZFet/ej6gaallVsi6wQ5fexe2bZsbV+gTG1fD8CPZa/ZZ+9XbW/41j48Xp9FdO4xOrV3zVvbaK2z1OqPjpnoOkNtvpgCNXFCCDmicBInhJCO4SROCCEd04UmninF7mZqmqNXxtMdbR5dtqWlteLVI3qavlbSSSPfA63t2d2KEx+7h3Y0xjxT01/1cWsNwX4v1cZo673ko99ILMW72/5u6eQ1LTmqx3pt6m9tAm/EoEfeQYjE5Gtq+Ut2ev57+cf0Y6tdz49on1paazglG7z7Qk2cEEJIE07ihBDSMd3IKV4YX+RzZGPC83SeqMxg642G5kWOLa2fvZGfkJ7kAtQ/HVfr29bn60p9Eu2fVpu1uqe8Wh6RFiKhgdFw1VrIXIloiJ+VAvNnzEr2WhttHZ4kEAmtLNldGictO2p5av3glQfqz08k9Diaxs+zEUII2QMncUII6RhO4oQQ0jFdaOK1zy4BsRAgwA9da+mrupxuI6rFRjTgaAhWK6RvjH7o9Yttu9T+2Pq8/JE1C0+zrW134NlcaiM6Zlrabi0ksbW1Q01v9dryfKv5W1oDsX1ZCuOrafNj9PFS3bbsVMY+u9mGWtnoM1arS1/LtkyFmjghhBxROIkTQkjHcBInhJCO6XISz9rS6fkT2FhdBtB+JReIv7avNbvLm1d3Xu3O5/n62srSHk378uZVnD1/aedaxtPGctrp+RN77N5YXd5TvqSp5fI13a6k7Zb8s3ltn5ReW9b+WxusP9bOtZWlPX4Du+/b2fOXdmm6np01P7SdpXK5bq9+baN3f71Y743V5ZAe7vXrxurynvb0fc5p2i5bj9a+c5v5vnv3Pj9T9po+znXlOvI9189jqY899Ni248zeE22P9wxaO3W93n3LbbTGjn0G9XPn+WrnpoOiy0mcEELINpzECSGkY8KTuIjcJCK/LiKPD+cnReQJEXlx+HvzwZlJCCHEIxwnLiI/BGARwImU0odF5McAvJpSekhEHgBwc0rpU7U69mvvFK2BZd0p57GMjamuxQXX9kRo7avhtanx8ra20vVsqflQi2kulfN8mhL/XIqbHbtVbquOSJ2t/i7hxbBHPgNW25/Exr2X+qVmj6fbR7YrrtVh7bBYvyLjy6P17sZrr1/D8WNz2Fhdrn6CzZsb8nFrjOq8ub2SL7W1lFof5uOpzBwnLiK3A/ibAH5KJd8D4MJwfAHAvZMtJIQQMomonPITAH4YwJ+otFtTSpsAMPy9xSsoIveLyLqIrG9tbc1iKyGEEENzEheRDwN4JaX0zJQGUkoPp5QWU0qLp06dmlIFIYSQApH/ib8fwEdE5AqAzwP4oIj8DIBvisg8AAx/XzkoI7WWpXWu0/MndmJULV4srY3/1HHAOla0pV15ceHa1nwt48Wm6thf7xqwHadrY3dt3GvJroj+pmOGbfy7jse1/avjn/O57kt73dOhc5zyVD1cx43rvLr/a3H1ug3b3zYmGtg9VmzdXhyzbifHVHux87pPtU+2/zy77PsBWqPWbdlxZvvTxjp7fujxaseijQ/31ok02h9vjOtP7h0/NrcnrRT/ndHvGtj8pVjvXE7vt67HVW5fr8GV1ua0Xfv1ebYazUk8pfRgSun2lNICgI8B+JWU0icAPAbg3JDtHIBHD8xKQgghLrPEiT8E4G4ReRHA3cM5IYSQQ2RuTOaU0lMAnhqOvwXgrv03iRBCSJRu9xOvxXfXYr0zrZjvWoyuvVb61qe1vRVvPgueTa09zb14WkvND13WO275VtPwx95z2w+lb61G4809/0vtR78zGt23XdOKdZ+yH75Ni3zfs9Zmra/HPp+lZ8ezPVJnqd9b/VOqc4w91p9Z4H7ihBByROEkTgghHcNJnBBCOqYbTby2N0NEz67tmVE6r+lpnoZb26PCI/L9P11fbZ8La59X35g8JfsiWncuC/j9EtlTo6TX1zRmr46oNj1V+9TppfZLe+tke6KMvR9ePHprTEe/D+q1E9W/x/jW0u69NaiWPfqaba9kc2n9ItKfpTljDNTECSHkiMJJnBBCOoaTOCGEdEw3mnimpEPZ65moHlmL9fXayNe9/REibUTigi12r4ZIPLutsxVbPOse0DXsPhN2P4zWeoWX1roHY3RiS+te1Wz09gkv+dZ6HyGi4Y+Jw/bSrd2z7gMeid/2jrWNU3Rq75qt0zu3dUTHuO2v1rM1FWrihBByROEkTgghHdOdnJJpvbJc+ikTfb24VHft9enWa8+Rn9KRn70lO1p5s/+R16G9n4Bjw/Ja98he88q3+jsS3qf9LpWL9kmkXe1X9HX+mp+2n8aELpZkvdx2yX+vjlq9JdmpNbamhBV6+W1fltos9UnLh8ub259uy1vVantqktIN83k2QgghNyacxAkhpGM4iRNCSMd0o4m3tOKorjpVf2y1tbayhDOrF3d0s1rInQ21i4Rl2bZar/l6vtg6Sm3X+qGk8Y3ZNtSm1+zUbZfq9vKX6qzp/i291WuzpnV7ftb8aOnZY8aI1x9jy0Zfwx9r85htD6LhitF1p/ycAnvDdFtrXbZvgPoWzXztnhBCSBNO4oQQ0jGcxAkhpGO60cQzYzW2mj5c060jn9zKadGtSCN6Xc1fy/qVV3H82Nwuba/kq62rFWPr+eIRjQ326iv1TWv9omTbGF11TFx1pvZKeune1tZexsYtj43D1msvuj91e7VnxbaR+6B1b8auVeQY7OPH5nau1ba+tfXVnqdSjHZL+9baue4z216p/dKzPRVq4oQQckThJE4IIR3DSZwQQjqmC00caH92KRPReks6rNdeKy59zB4SrbjqqMbW2ttiamx2a6vcVvv53NPqS/1obdC25H6q+Vsqn/NF9seIaMKRrY/HxNtbPz1fW3pza+vX2rMQeVeh1QfRNZTaGkS0rtYzb+uMrEV4RPaVGbs2YeueAjVxQgg5onASJ4SQjuEkTgghHdOFJj5m3+OWxqvT87GOB7XXvTJe26VyJZtbZcdq0VNi0afE3I+JfS7ZqduvlRnzDkBtbWSMvjkmVjuq4Zbq1bbXxmppfHoafmT8e3bYvop+ps0yZr8Vz+4x96f1LkduU++Vkn0qxdXrvY80rXUI236pb6ZCTZwQQo4onMQJIaRjmpO4iLxFRL4oIr8hIl8VkR8Z0k+KyBMi8uLw9+aDN5cQQoimqYmLiAB4a0rpD0TkTQB+DcAnAfxtAK+mlB4SkQcA3JxS+lStrln3Ey9ph5F435ZGCezVAe0eE6VypfLedU+Lr8VT63Jj+iDTite1ZSJ7xkR1Sk9nHbMnivahpht7/k6JQ25p/a092KMx97X70NrTpZQ3Mr61X3bfexs/Xhvz0bWnSJnS9Sg1f8esebT2ECqNZ01r/6RZmEkTT9v8wXD6puFfAnAPgAtD+gUA985kJSGEkNGENHERuUlEvgzgFQBPpJSeBnBrSmkTAIa/txTK3i8i6yKyvrW1tU9mE0IIAYKTeErpj1NK7wZwO4D3isj3RBtIKT2cUlpMKS2eOnVqopmEEEI8RseJi8i/BvCHAP4JgA+klDZFZB7AUymlO2tlD2LvlExkjwydN7LnsW0rWm5s7HEpj9d+6Votb82HaCyvp8FGv1Na8jsa72/x9MvW3vCZ6F4jpXjmml/RWP2SH7r+SF37/c5BNNbe0hor9lqpPduuZ6cmshdRZO3La3Osft+KVZ+VmTRxETklIm8bjr8DwPcBeB7AYwDODdnOAXh0ZksJIYSMYi6QZx7ABRG5CduT/iMppcdF5BKAR0TkPgBfA/DRA7STEEKIQ3MSTyl9BcB7nPRvAbjrIIwihBASo6s3Nk/Pn9iJX9VaWd4TIaetrSw142RzjLYtt7G6vNOOrqOkAXuanZeu07S9tp2cz/O11Te6Tu2X7pPse6lOrRlbu3K5nK77UNugy2Y90OqDuu6N1eU99zCX9fTEs+cv7dyHXI++L6X4Xauha/SYsX1g+8LzS/dHzq/b0H7lY+2H5vLm1ZAGrceqbWNtZWnP2NH3pKY56zGY/aqNQ33f833UY0Vfz/Xb9nV/luq3nFm9uPPPq0OPKV2vV6aGvRdevdnPfK1m937T1SROCCFkN5zECSGkY7rYihYov8peek04ul2nrcNLt9eiW2p65zatFAIV3SrW86H1an2p7TGvPkfDPWuf5Wq9Yu6dj/k8mbb1tdev4fixuZ1xUgsnA/aGn43xS9c7ZivdyNaq1qa8ZUPEJo9ImTHbB5f6ydqtfbZ1a0r3ptSPut5aWOCU+9cKqbVtl57fKXArWkIIOaJwEieEkI7hJE4IIR3TvSY+9hVaXc77BFuprto2tV7es+cvYf3Kqzh+bG6Xbqkp6dg1zX6WV7u9/ilR07xr6fqataP0erm2b8qn31qvS0fWPLQ90W0V9LXW2kDLvpqu3PLZs7/lW2Qb1Zp/Y+7xfjBmTUifl+yIrDG06o08t56mPwVq4oQQckThJE4IIR3DSZwQQjqmG00caH+mLePFZEZ0q5L2N+ZTUl59NR2xZX9Jgy9tw1mL19ZE1xBsHH5Ud8/pXr+3PmPn+VvTF63uXdNNx+rZtfHmXbNl9bm1W9se0cS9ukt91IoVj6yptN6DGLOu4vVRLU0zZm0j6q9HNM5+7JpMPp4FauKEEHJE4SROCCEdw0mcEEI6pgtNPBIjPiYeXOPFckc020xrD4eardG9MiKxx609RWp9EY05t3nHxqyX6otooWPWKqJ7quQ2NBGNurYeU7NFMyWmvzbWvLqjcfPRMRRdv6mtV2jbvXUrz2+dpstFPhVYu2b7S/ttj70+LJW7vLm9V8/iwsl9+TQbQE2cEEKOLJzECSGkYziJE0JIx3SniUfjU1t6dCl/NB43kmes1huxOaJXazvstVZ8ayvWNbLHRnS/lprGaduPapytPS5aPtk1klnWAFrx47re3F+tOP/aGopXzto35j5Y+3T7rXWgsf3Wssurs9Qfnj4/dq1jjJ3WXovtwylQEyeEkCMKJ3FCCOkYTuKEENIxXWjiQD3eVDNFL7Z5WrGws2rowPjvFtbszvVFv8/otTElr243H59ZvbjTdlQ7zXlr5SL7WJSI6sXWH8/3Vh12XQDYu1d6RN+PxIa3dPiWX7oeL+46294aU5Hn0tobWfMZ+/5BZE/+kk01/d1ru9Wf+dx7/2QK1MQJIeSIwkmcEEI6hpM4IYR0TJeaeKamUWUi8ae6fDQWtcSYfTvG2FcqU7N3rPZv6x9TX60vM7W1Bt2OZapePzVfy16dr6Une+i1A6uZ1vYDsXVE4vBrunnNL8//sfV65XK9UY0+amM0BrxUX2tdJeJT65yaOCGEkD00J3EReYeI/KqIPCciXxWRTw7pJ0XkCRF5cfh788GbSwghRBP5n/g1AP8ypfQXAbwPwD8TkdMAHgDwZErpDgBPDueEEEIOkeYknlLaTCn9r+H42wCeA3AbgHsAXBiyXQBw7wHZuKNlnZ4/sfNPc3nzalOL87RWnabjm3W7tu1SeaAcpxrR29ZWlvZoZiW7a2hdWZfNx2srS8143o3VZQDbfXJm9aIb+6zR1+y90PdM59F9pPt5Y3V5J//G6vKOLbku65O1y17X8bo1TbK2TqDr0O1me2xsfO43z6bsbybnu7x5Fd/94C/uqce7h7lcSRvWftfGkI3N9/Lpe5bPdR/pcnZNo7YmYp9X2++t/b9t/jwOsn1efdoXbzx7vnrXtM92bGtbSuf7zShNXEQWALwHwNMAbk0pbQLbEz2AW/bdOkIIIVXCk7iIfCeAnwfwgyml8H8PReR+EVkXkfWtra0pNhJCCCkQCjEUkTcBeBzAxZTSjw9pLwD4QEppU0TmATyVUrqzVs9+f57N/gSOvMKsz6O0Xguu1Vl7JdmGp9m6vJ+R0a01I68sl0K+7M9s6180pG7K59VKNkdDuWp1l+qNvmau686MCWuthVh6cl6pD6f015QyrbE61sfMmC2KW7Zr+zx/auO25L+mFurZeja8a1OZKcRQRATATwN4Lk/gA48BODccnwPw6GQLCSGETGIukOf9AP4+gA0R+fKQ9q8APATgERG5D8DXAHz0QCwkhBBSpDmJp5R+DYAULt+1v+YQQggZQxev3U95vbilN7ZetW29Qqzz1LS3TE1b9bToVjuRNktt6/TaK+K1NYTWZ+dsu2Ne4R7Tbm2NYqwOa3VfuyYQ1WRbvrfWVUpjpKXfltptjfWobaW+H/MZujFbK0TWMUrbvJbWWiJjYuwaTq2PNHztnhBCyB44iRNCSMdwEieEkI7pQhMH2tuhWjwtq/Ya8BStzNPlxsbRTrnW2sbT+u/VVdtmtrWtbzRudmw8cc4bibWtxZeXKMX1R+OHSzZrovp31N4xWru1q6bx5y1wo/W27nnNn9YaB1B/7yG6xlRrq+Rnre6anZ4d3rOxfuVVHD82t2v7iClQEyeEkCMKJ3FCCOkYTuKEENIxXWnimalaq5e3prFmalroFP3XMlZ/jcQCX968itdev4bFhZNuPS1tcEzsdsZbH/DKeusPY3T2Vgy2tbMVz19rK+P1V7a/Fi9cijOvrU20Pr1W209lLLU49pp+XtOyo/1X86H2ucTcbsmH0jhs2VSyoUSpP2p2TYWaOCGEHFE4iRNCSMdwEieEkI7pQhMfo2lFdM9a3txeLY7Utq3TWnpcSw+s6e1e+17ddq+YMfsh12yo2ePZOyaOOdsJ1PcJae1pMSbmO3KvrN2tOHDv/kZsiqyJvPb6NRw/Nrdn7aOU3+uzaJvRfbFr4ziyfjKlH6L3PJcvrcHYtEi71gavD0ox59TECSGE7IGTOCGEdAwncUII6ZiuNPFaXG9E+/Zie1vadlSza2nm0b0casclWnr1GK26pq1G7fXaL8VL6/I21jaiU1s7I/cnondGjluxwa2Y/JIPuczY9RhLZC2kpLWXaO2rY49L1z0fo+tanh8lOyN6fWQf/KhG37JpKtTECSHkiMJJnBBCOoaTOCGEdEw3mngtljVCLZbU1hu5NlY39/Ree1461vk8pvRNS0Nu5avlLdkS6TcvBjlSZ83Hmr0ljbqld+d9om2ds8Y/W7uj+ne+bvc5adUR0dtzemtdI59H9yMqPQOte6rvAYCd/YGiWne+ltsu2dZ6l2Ps+gs1cUIIIXvgJE4IIR3DSZwQQjqmK028FG9rNShPA2vpZTr/lDjtXN7qXlYL1mkltA9Z49R4GrPXVmQvEG139Nui2k5bj1e/zRtdT6jFFbfsbNHKOzZGPRKLDLS/FRvRsmt7vtj2x3xDshWvbevN/tTssH556x61MrXxbcvPEjee8fT2MWWje86MhZo4IYQcUTiJE0JIx3Qjp1iiIVw1GWXsK7SlEKuaXNL6aVrK0wr9akkjtt8ivkXkqpZPmpJ8VLKx5E/EtkjIXkvCiI6nSHtTQlajttl6am3XpMWo5JZfy99YXd5J98ICW75FnokpoZulsWKJhAZG73vkebZ9QDmFEELIHjiJE0JIxzQncRH5jIi8IiLPqrSTIvKEiLw4/L35YM0khBDiEfmf+H8G8CGT9gCAJ1NKdwB4cjg/cLSmpLUwYFtLy9pT1qTWVpZ26XiAH5J0efMqLm9e3cmbP2+2sbq8R9/K2pZOz8dnz1/aqT/bkev2dLqN1WVsrC5jbWVply85f27f6nC5jMWzw/bX2fOXdpVdW1na8UeHNep+8HzVdWq79bntD41nYy6b7dF9bfvxzOpFN0zT6qN2TOg82g99Pfuvy9a0aV3G8zXryt41Ly3fEzvO7LlO89Zr8vj3xnAN+5xpjfrs+Us74wLYfR8ium9p/ORyx4/NAdjbv7W1mahv+dmx9dt7q/tNj2nvOcjXvXtzWDQn8ZTS/wTwqkm+B8CF4fgCgHv31yxCCCERpmrit6aUNgFg+HtLKaOI3C8i6yKyvrW1NbE5QgghHge+sJlSejiltJhSWjx16tRBN0cIIX+qCMWJi8gCgMdTSt8znL8A4AMppU0RmQfwVErpzlY9+xEnHokTLsWCl+JiW69ER2J1MznO1dsS1JZpxT97MarRz05FymS/9RpAbbveks+6rlo7pa1Jp8bpjinv9ZFnWyt2OvIOgc0bqavU3zaee+rWxFGi7wZ4GnwthlrbFhlbpTpb2wa3ngHto/U3Enfv2WZ9q9k4lYOIE38MwLnh+ByARyfWQwghZAYiIYafA3AJwJ0i8rKI3AfgIQB3i8iLAO4ezgkhhBwyc60MKaWPFy7dtc+2EEIIGUkXe6dkovt41DQ5S2T/BK89Xd7TnSNbdXpbckb2oqjZ47G2srRH9y7Zb9uJaOwR7Tf7WLo3rfa9eqw9pfY9ojpua8tizdh+bNkxZvx5dZVsKPnlldV5S21G1i28+r2xXuq/fL2mxev8dk0qqtdrO4B4fLpXj+fvVLh3CiGEHFE4iRNCSMdwEieEkI7pQhMfo5F5caW6jlJ+W38tflu3b8vWNPqahjdWY2zFytu6vPbGfOps7GfRatrv5c039qi2/gBlzdvT4Fu6aqvPx3ymr2SP7SPvXYGWHh3RoktlbR+24p69/F7ftNqMxnu3nqWI3bY/PO27tT42Rr9vxafX5iGdR4/rWaAmTgghRxRO4oQQ0jGcxAkhpGO60cRbccs6L1DXE73rLd0yqjdHteFMVIO29dXK13TvVt1R/XCsHl6LzS1p056fkX6InOe0SBx4bf2hND6snZE6S/564zvX5eXT+8BH+r50HBkntbRWH9i0Vj/X1pgi62GRPVNqbVqbPf+s/m33vJkKNXFCCDmicBInhJCO4SROCCEd040mXiLH5GqNDJi2j4WXLxK/Wkpr1acZExvrxaKW/Cvp/GO1b+uLl6+1p3gkNrrk59jY4toe02M00DFx9iX7I7r5mDUGm56padqRsRNdV/DanxqjXRoXY/TsltZdKmMZM1ZqfWCv7UesODVxQgg5onASJ4SQjuEkTgghHdPNJK61J32c9XBgW3vaWF12tce1laWda1qr8vJqHezM6kVXG7u8eXXneG1laU89aytLu3SwfF3H8Op/Oa9tR2uuOW9Ni/R0uZJ2CmzvPeHVNWbfB61Lnj1/aVd/eLHL1k5rY76PuQ8jenjOl/vnzOrF4lqK1ye63lJ/6D6+vHl1j96sy2n/SuNR33Ntv+43ALvGdE7XfaHHtbUl1+vl8e69tl/7rPPaMrVxouvw/Ml9WHpuvfrsOlcec7o+W1d+jnWZs+cv4czqxZ1/tbWJ0j20faDvjzdWDoJuJnFCCCF74SROCCEd002IoffzvhSqF9leVH+yzEosOp+l9rqutcO250kC2Y9MK6zP+qhtsvltiJVuy/pf8rHWJ6VwrlJe2771vdW2Z2st3QsDLNWrmRKOWLufXr0R36wvmTFbMrT60+azddXabd33ml+1vqqFSnp+lvyy10rlWz7UnmV732s+zAJDDAkh5IjCSZwQQjqGkzghhHRMl5p4RIdtvcpsy4zRmC9vbn9eDAAWF06GtrjVdbRe7c35PN16yqvzus5SWq47U3vdudSHtfZrOmdrncPaOYba2kXJjsh6gNdGpK9qZXPbJdtbuntpvcJe88pENH/Pbu841xN5jmr3tVY+8rr+2OfCUqpb1+/5WpoDZtHFqYkTQsgRhZM4IYR0DCdxQgjpmC40cU1LbwbKW1V6WmumFtOs02t6vFeXvRbN09IjS/lq8d+2jkgcfKkvWnXWYmit7bZdS2ktJBqnbOuI2K3Ll9YNMq13Elr1tuK7W+8ZeDZ49Ws/ajHcU+59rV6v30vbx1qmrO3U8kfHZ2QNIRpDT02cEEJIEU7ihBDSMTNN4iLyIRF5QUReEpEH9ssoQgghMSZr4iJyE4DfBHA3gJcBfAnAx1NKl0tlZtHESzHFOq2kV2dmiX0t5bU6l7YhojtmpsRg6zyZUnzuGL9K9WQ7z6xexGuvX8PxY3O72i3ZWPrUVuse2LJj7sms+8J4fWL7JXIPSzpoTWeOvkdQKlOr15aPrl3U9gyKridYGyy1vi3R2q/G2mfbKVHT4GfV36dwUJr4ewG8lFL67ZTS/wPweQD3zFAfIYSQkcwyid8G4Ovq/OUhbRcicr+IrIvI+tbW1gzNEUIIscwyiYuTtkebSSk9nFJaTCktnjp1aobmCCGEWGbRxJcArKaUlofzBwEgpfRvSmX2I06cEEL+tHFQmviXANwhIu8UkTcD+BiAx2aojxBCyEjm2ll8UkrXROSfA7gI4CYAn0kpfXXfLCOEENJk8iQOACmlXwLwS/tkCyGEkJHwjU1CCOkYTuKEENIxnMQJIaRjOIkTQkjHcBInhJCO4SROCCEdw0mcEEI6hpM4IYR0zKF+Y1NEtgD87sTibwfwe/tozvXkqPhyVPwA6MuNCn3Z5i+klNwdBA91Ep8FEVkvbQDTG0fFl6PiB0BfblToSxvKKYQQ0jGcxAkhpGN6msQfvt4G7CNHxZej4gdAX25U6EuDbjRxQgghe+npf+KEEEIMnMQJIaRjbvhJXEQ+JCIviMhLIvLA9bbHQ0TeISK/KiLPichXReSTQ/pJEXlCRF4c/t6syjw4+PSCiCyr9L8sIhvDtX8vIt4HqQ/an5tE5NdF5PHO/XibiPyciDw/3Juljn35F8PYelZEPicib+nFFxH5jIi8IiLPqrR9s11EjonI2pD+tIgsHLIv/3YYY18RkS+IyNsO1ZeU0g37D9ufffstAO8C8GYAvwHg9PW2y7FzHsD3DsffBeA3AZwG8GMAHhjSHwDwo8Px6cGXYwDeOfh403DtiwCWAAiA/w7gb1wHf34IwM8CeHw479WPCwD+8XD8ZgBv69EXALcB+B0A3zGcPwLgH/biC4C/BuB7ATyr0vbNdgD/FMBPDscfA7B2yL58P4C54fhHD9uXQ32oJnTYEoCL6vxBAA9eb7sCdj8K4G4ALwCYH9LmAbzg+YHt75QuDXmeV+kfB3D+kG2/HcCTAD6INybxHv04ge2JT0x6j77cBuDrAE5i+5OKjw8TRze+AFgwE9++2Z7zDMdz2H4rUg7LF3PtbwH47GH6cqPLKXnwZl4e0m5Yhp8/7wHwNIBbU0qbADD8vWXIVvLrtuHYph8mPwHghwH8iUrr0Y93AdgC8J8GaeinROSt6NCXlNL/BvDvAHwNwCaA308p/TI69EWxn7bvlEkpXQPw+wD+3IFZXucfYft/1rvsGjgQX270SdzT627YmEgR+U4APw/gB1NKV2tZnbRUST8UROTDAF5JKT0TLeKkXXc/Buaw/bP3P6aU3gPgD7H9s73EDevLoBffg+2f5H8ewFtF5BO1Ik7aDeFLgCm23xB+icinAVwD8Nmc5GTbd19u9En8ZQDvUOe3A/jGdbKlioi8CdsT+GdTSr8wJH9TROaH6/MAXhnSS369PBzb9MPi/QA+IiJXAHwewAdF5GfQnx8YbHg5pfT0cP5z2J7Ue/Tl+wD8TkppK6X0RwB+AcBfQZ++ZPbT9p0yIjIH4M8CePXALHcQkXMAPgzg76VBC8Eh+XKjT+JfAnCHiLxTRN6MbaH/sets0x6GleWfBvBcSunH1aXHAJwbjs9hWyvP6R8bVqLfCeAOAF8cflZ+W0TeN9T5D1SZAyel9GBK6faU0gK2+/pXUkqf6M2PwZf/A+DrInLnkHQXgMvo0BdsyyjvE5Hjgw13AXgOffqS2U/bdV1/B9vj9jB/wX4IwKcAfCSl9Jq6dDi+HMaixoyLCD+A7WiP3wLw6ettT8HGv4rtnzxfAfDl4d8PYFvLehLAi8Pfk6rMpwefXoCKEACwCODZ4dp/wAEu0DR8+gDeWNjs0g8A7wawPtyX/wbg5o59+REAzw92/BdsRzx04QuAz2Fby/8jbP9P8779tB3AWwD8VwAvYTvq412H7MtL2Nax87P/k4fpC1+7J4SQjrnR5RRCCCEVOIkTQkjHcBInhJCO4SROCCEdw0mcEEI6hpM4IYR0DCdxQgjpmP8P0FMLcKWbAr4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#from models import Model, SynapticConductanceModel\n",
    "from brian2 import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "model = SynapticConductanceModel(resistance=Model.LOW, # Model.LOW, Model.HIGH\n",
    "                                 noise=Model.HIGH, n=50, # Model.OFF, Model.LOW, Model,HIGH\n",
    "                                 offset=SynapticConductanceModel.ACTIVE)\n",
    "\n",
    "\n",
    "# model = SineModel(resistance=Model.LOW, # Model.LOW, Model.HIGH\n",
    "#                                   noise=Model.OFF, n=1, # Model.OFF, Model.LOW, Model,HIGH\n",
    "#                                 )\n",
    "\n",
    "\n",
    "model.f = 100 * Hz\n",
    "model.set_stimulus_current(400 * nS) # current should be scaled by 100x for Active Offset so (500nS is actually 5nS)\n",
    "model._set_variable(\"i_injected\", 65 * pA)\n",
    "\n",
    "\n",
    "model.run(12*second)\n",
    "\n",
    "spike_times = [s/ms for s in model.spike_train.values()]\n",
    "\n",
    "plt.eventplot(spike_times)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7f22e21e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(spike_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "20e3191c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.io\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "file_path = 'data0.mat'\n",
    "scipy.io.savemat(file_path, {'data0': spike_times[0]})\n",
    "\n",
    "file_path = 'data1.mat'\n",
    "scipy.io.savemat(file_path, {'data1': spike_times[1]})\n",
    "\n",
    "file_path = 'data2.mat'\n",
    "scipy.io.savemat(file_path, {'data2': spike_times[2]})\n",
    "\n",
    "file_path = 'data3.mat'\n",
    "scipy.io.savemat(file_path, {'data3': spike_times[3]})\n",
    "\n",
    "file_path = 'data4.mat'\n",
    "scipy.io.savemat(file_path, {'data4': spike_times[4]})\n",
    "\n",
    "file_path = 'data5.mat'\n",
    "scipy.io.savemat(file_path, {'data5': spike_times[5]})\n",
    "\n",
    "file_path = 'data6.mat'\n",
    "scipy.io.savemat(file_path, {'data6': spike_times[6]})\n",
    "\n",
    "file_path = 'data7.mat'\n",
    "scipy.io.savemat(file_path, {'data7': spike_times[7]})\n",
    "\n",
    "file_path = 'data8.mat'\n",
    "scipy.io.savemat(file_path, {'data8': spike_times[8]})\n",
    "\n",
    "file_path = 'data9.mat'\n",
    "scipy.io.savemat(file_path, {'data9': spike_times[9]})\n",
    "\n",
    "file_path = 'data10.mat'\n",
    "scipy.io.savemat(file_path, {'data10': spike_times[10]})\n",
    "\n",
    "file_path = 'data11.mat'\n",
    "scipy.io.savemat(file_path, {'data11': spike_times[11]})\n",
    "\n",
    "file_path = 'data12.mat'\n",
    "scipy.io.savemat(file_path, {'data12': spike_times[12]})\n",
    "\n",
    "file_path = 'data13.mat'\n",
    "scipy.io.savemat(file_path, {'data13': spike_times[13]})\n",
    "\n",
    "file_path = 'data14.mat'\n",
    "scipy.io.savemat(file_path, {'data14': spike_times[14]})\n",
    "\n",
    "file_path = 'data15.mat'\n",
    "scipy.io.savemat(file_path, {'data15': spike_times[15]})\n",
    "\n",
    "file_path = 'data16.mat'\n",
    "scipy.io.savemat(file_path, {'data16': spike_times[16]})\n",
    "\n",
    "file_path = 'data17.mat'\n",
    "scipy.io.savemat(file_path, {'data17': spike_times[17]})\n",
    "\n",
    "file_path = 'data18.mat'\n",
    "scipy.io.savemat(file_path, {'data18': spike_times[18]})\n",
    "\n",
    "file_path = 'data19.mat'\n",
    "scipy.io.savemat(file_path, {'data19': spike_times[19]})\n",
    "\n",
    "file_path = 'data20.mat'\n",
    "scipy.io.savemat(file_path, {'data20': spike_times[20]})\n",
    "\n",
    "file_path = 'data21.mat'\n",
    "scipy.io.savemat(file_path, {'data21': spike_times[21]})\n",
    "\n",
    "file_path = 'data22.mat'\n",
    "scipy.io.savemat(file_path, {'data22': spike_times[22]})\n",
    "\n",
    "file_path = 'data23.mat'\n",
    "scipy.io.savemat(file_path, {'data23': spike_times[23]})\n",
    "\n",
    "file_path = 'data24.mat'\n",
    "scipy.io.savemat(file_path, {'data24': spike_times[24]})\n",
    "\n",
    "file_path = 'data25.mat'\n",
    "scipy.io.savemat(file_path, {'data25': spike_times[25]})\n",
    "\n",
    "file_path = 'data26.mat'\n",
    "scipy.io.savemat(file_path, {'data26': spike_times[26]})\n",
    "\n",
    "file_path = 'data27.mat'\n",
    "scipy.io.savemat(file_path, {'data27': spike_times[27]})\n",
    "\n",
    "file_path = 'data28.mat'\n",
    "scipy.io.savemat(file_path, {'data28': spike_times[28]})\n",
    "\n",
    "file_path = 'data29.mat'\n",
    "scipy.io.savemat(file_path, {'data29': spike_times[29]})\n",
    "\n",
    "file_path = 'data30.mat'\n",
    "scipy.io.savemat(file_path, {'data30': spike_times[30]})\n",
    "\n",
    "file_path = 'data31.mat'\n",
    "scipy.io.savemat(file_path, {'data31': spike_times[31]})\n",
    "\n",
    "file_path = 'data32.mat'\n",
    "scipy.io.savemat(file_path, {'data32': spike_times[32]})\n",
    "\n",
    "file_path = 'data33.mat'\n",
    "scipy.io.savemat(file_path, {'data33': spike_times[33]})\n",
    "\n",
    "file_path = 'data34.mat'\n",
    "scipy.io.savemat(file_path, {'data34': spike_times[34]})\n",
    "\n",
    "file_path = 'data35.mat'\n",
    "scipy.io.savemat(file_path, {'data35': spike_times[35]})\n",
    "\n",
    "file_path = 'data36.mat'\n",
    "scipy.io.savemat(file_path, {'data36': spike_times[36]})\n",
    "\n",
    "file_path = 'data37.mat'\n",
    "scipy.io.savemat(file_path, {'data37': spike_times[37]})\n",
    "\n",
    "file_path = 'data38.mat'\n",
    "scipy.io.savemat(file_path, {'data38': spike_times[38]})\n",
    "\n",
    "file_path = 'data39.mat'\n",
    "scipy.io.savemat(file_path, {'data39': spike_times[39]})\n",
    "\n",
    "file_path = 'data40.mat'\n",
    "scipy.io.savemat(file_path, {'data40': spike_times[40]})\n",
    "\n",
    "file_path = 'data41.mat'\n",
    "scipy.io.savemat(file_path, {'data41': spike_times[41]})\n",
    "\n",
    "file_path = 'data42.mat'\n",
    "scipy.io.savemat(file_path, {'data42': spike_times[42]})\n",
    "\n",
    "file_path = 'data43.mat'\n",
    "scipy.io.savemat(file_path, {'data43': spike_times[43]})\n",
    "\n",
    "file_path = 'data44.mat'\n",
    "scipy.io.savemat(file_path, {'data44': spike_times[44]})\n",
    "\n",
    "file_path = 'data45.mat'\n",
    "scipy.io.savemat(file_path, {'data45': spike_times[45]})\n",
    "\n",
    "file_path = 'data46.mat'\n",
    "scipy.io.savemat(file_path, {'data46': spike_times[46]})\n",
    "\n",
    "file_path = 'data47.mat'\n",
    "scipy.io.savemat(file_path, {'data47': spike_times[47]})\n",
    "\n",
    "file_path = 'data48.mat'\n",
    "scipy.io.savemat(file_path, {'data48': spike_times[48]})\n",
    "\n",
    "file_path = 'data49.mat'\n",
    "scipy.io.savemat(file_path, {'data49': spike_times[49]})\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d353eff4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}