2DofArm_simulation-Main.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulates 20 trials for each one of the trajectories\n",
    "# These trials will only be used to train the readout.\n",
    "## The results are saved into the base_dir to files following this template: \n",
    "\n",
    "### \"States\\_Torques\\_movement\"+str(pos\\_i)+\"\\_LSM_\"+str(run_i)+\".gzpickle\") or\n",
    "### \"States\\_Torques\\_movement\"+str(pos\\_i)+\"\\_LSM_\"+str(run_i)+\"_STP_ON.gzpickle\")\n",
    "\n",
    "\n",
    "## Which things you can change:\n",
    "1)Liquid's parameters.  \n",
    "2)The amount of noise injected into the inputs.  \n",
    "3)STP: ON and OFF.  \n",
    "4)Stratification or simply 30% random excitatory neurons instead.  \n",
    "5)Width of the gaussian used to distribute the inputs into the liquid.  \n",
    "6)Lambda value.  \n",
    "7)Noisy: different starting and background noise levels.  \n",
    "8)Initial seed for the liquid's structure. If my_seed==0 the system uses a random liquid.  \n",
    "9)Number of trials.  \n",
    "10)Number of different trajectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Makes possible to show the output from matplotlib inline\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Makes the figures in the PNG format:\n",
    "# For more information see %config InlineBackend\n",
    "%config InlineBackend.figure_formats=set([u'png'])\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 20, 10\n",
    "\n",
    "import numpy\n",
    "import sys\n",
    "import save_load_file as slf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loads the modules and starts the object to be used with the parallel processing iPython stuff...\n",
    "\n",
    "# Remember to execute at the shell: ipcluster start -n 4\n",
    "# or from the iPython notebook interface!\n",
    "# from IPython.parallel import Client\n",
    "from ipyparallel import Client\n",
    "\n",
    "# When using the ipython in my desktop, launch the cluster in the right profile :)\n",
    "cli = Client()\n",
    "\n",
    "# lbview = cli.load_balanced_view()\n",
    "dview = cli[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Brian Step-By-Step class\n",
    "## This class alows the user to do step-by-step simulations using Brian WITHOUT having to generate the input spikes before the simulation starts. It's also useful because the simulation can last, idealy, forever without using a giant input spike list.\n",
    "\n",
    "https://github.com/ricardodeazambuja/BrianStep-By-Step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing step_by_step_brian.py\n"
     ]
    }
   ],
   "source": [
    "%%file step_by_step_brian.py\n",
    "\n",
    "class step_by_step_brian_sim(object):\n",
    "    '''\n",
    "    Step-by-Step Brian Simulation (Nov/2014 - ricardo.deazambuja@plymouth.ac.uk)\n",
    "    \n",
    "    This class was created to make it easier to run a Brian simulation step-by-step, passing input spikes without \n",
    "    running out of memory or having to create the input spikes beforehand. It also makes the code more clear because\n",
    "    you generate a separated function with your Brian simulation code. This function (here called simulation) receives:\n",
    "    simulation(brian.defaultclock, brian)\n",
    "    brian.defaultclock: Brian defaultclock to be used\n",
    "    brian: the result of the command \"import brian\", so the user doesn't need to import Brian and have to use \"brian.\".\n",
    "    And must return a tuple:\n",
    "    (Input_layer, Output_layer, pop_objects, syn_objects, monitors_objects)\n",
    "    Input_layer: is a SpikeGenerator\n",
    "    Output_layer: the layer the user wants to output the spikes\n",
    "    pop_objects: a list with all the NeuronGroups used\n",
    "    syn_objects: a list with all the Synapses objects used\n",
    "    monitors_objects: a list with all the Monitors or functions used with the @network_operation decorator\n",
    "    \n",
    "    At initialization the simulation step size (in ms) can be passed (default is 2).\n",
    "    After the creation of the instance, calling the method \"run_step(input_spike_index_list)\" sends the spikes \n",
    "    to the simulation and simulates one step (according to the initialization).\n",
    "    The method run_step returns returns a tuple:\n",
    "    int(number_of_the_run), \n",
    "    float(current_simulation_time), \n",
    "    numpy.array(tuple(processed_received_spikes)),\n",
    "    list(list(output_spikes)), \n",
    "    list(float(output_spikes_times))\n",
    "    '''\n",
    "    \n",
    "    def __init__(self, simulation, init_step_size=2):\n",
    "        print \"Initializing the simulation...\"\n",
    "        self.step_size = init_step_size\n",
    "        self._generator = self._return_generator(simulation)\n",
    "        print \"Initializing the simulation...Done\"\n",
    "        print \"Call .run_step(input_spikes_list) to run one step of the simulation!\"\n",
    "    \n",
    "    def run_step(self,input_spikes=None):\n",
    "        '''\n",
    "        Calls the generator .next and send methods and returns the spikes and times generated.\n",
    "        '''\n",
    "        self._generator.next() # Runs up to the first yield (where the generator waits for the .send method)\n",
    "        ans = self._generator.send(input_spikes) # Sends the spikes and runs to the second yield \n",
    "                                                 #(where the generator returns the result of the simulation)\n",
    "        return ans\n",
    "    \n",
    "    def _return_generator(self, simulation):\n",
    "        '''\n",
    "        Defines a simulation using a python generator.\n",
    "        '''\n",
    "\n",
    "        import brian\n",
    "        import numpy\n",
    "\n",
    "        print \"Starting the simulation!\"\n",
    "\n",
    "        print \"Reseting the Brian Simulation object...\",\n",
    "        brian.reinit() # This is only necessary when using the same enviroment over and over (like with iPython).\n",
    "        print \"Done!\"\n",
    "\n",
    "        clock_mult = self.step_size\n",
    "        brian.defaultclock.dt = clock_mult*brian.ms\n",
    "        \n",
    "        print \"Initial simulation time:\", brian.defaultclock.t\n",
    "        print \"Simulation step:\", brian.defaultclock.dt\n",
    "        \n",
    "        # Calls the user function with the Brian objects to be used in the simulation\n",
    "        Input_layer, Output_layer, pop_objects, syn_objects, monitors_objects = simulation(brian.defaultclock, brian)\n",
    "\n",
    "        output_spikes = []\n",
    "        output_spikes_time = []\n",
    "\n",
    "        # Every time spikes occur at the SpikeMonitor related to the output neuron group, this function is called\n",
    "        def output_spikes_proc(spikes):\n",
    "            if len(spikes):\n",
    "                output_spikes.append(spikes.tolist()) # Saves the indexes of the neurons who generated spikes\n",
    "                output_spikes_time.append(1000*float(brian.defaultclock.t)) # Converts and save the actual time in ms\n",
    "                # The spike monitor and all this code could be replaced by the .get_spikes() method of neurongroups.\n",
    "                # I need to check what is fastest way!\n",
    "\n",
    "        OutputMonitor=brian.SpikeMonitor(Output_layer, record=False, function=output_spikes_proc)\n",
    "        # Because it is not saving, the system is not going to run out of memory after a long simulation.\n",
    "        \n",
    "        net = brian.Network(pop_objects + syn_objects + monitors_objects + [OutputMonitor])\n",
    "\n",
    "        r=0\n",
    "        while True:\n",
    "            spiketimes = yield # Receives the content from the Python generator method .send()          \n",
    "            if spiketimes:\n",
    "                spiketimes = [(i,brian.defaultclock.t) for i in spiketimes] # The spikes received are inserted as the last simulated time\n",
    "                Input_layer.set_spiketimes(spiketimes)\n",
    "            net.run(clock_mult*brian.ms) # I'm running one step each time this function is called\n",
    "            r+=1\n",
    "            yield (\n",
    "                  r,\n",
    "                  float(brian.defaultclock.t)*1000, \n",
    "                  numpy.array(Input_layer.get_spiketimes()).astype(dtype=numpy.float), # I'm doing this way to prove the spikes were received\n",
    "                  output_spikes,\n",
    "                  output_spikes_time \n",
    "                  )# After the .send method, the generator executes this line and stops here\n",
    "\n",
    "            output_spikes=[] # Cleans the output_spikes list so only the last spikes generated are sent\n",
    "            output_spikes_time=[] # Cleans the output_spikes list so only the last spikes generated are sent            \n",
    "            \n",
    "            # I'm using the .astype(numpy.float) because the arrays have Brian objects (units),\n",
    "            # and I think using only floats the memory footprint can be smaller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import step_by_step_brian\n",
    "reload(sys.modules['step_by_step_brian']) # Makes sure the interpreter is going to reload the module\n",
    "step_by_step_brian_sim = step_by_step_brian.step_by_step_brian_sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Membrane low-pass filter class\n",
    "## This class automates the process to read and generate the equivalent output of a low-pass membrane behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing membrane_lowpass_md.py\n"
     ]
    }
   ],
   "source": [
    "%%file membrane_lowpass_md.py\n",
    "\n",
    "import numpy\n",
    "\n",
    "class membrane_lowpass(object):\n",
    "    \n",
    "    def __init__(self, Number_of_Neurons, tau):\n",
    "        '''\n",
    "        Initializes the neuron membranes.\n",
    "        Number_of_Neurons: total number of neurons to be simulated\n",
    "        tau: time constant (in seconds)\n",
    "        '''\n",
    "        self.neurons = numpy.zeros(Number_of_Neurons)\n",
    "        self.times = numpy.zeros(Number_of_Neurons)\n",
    "        self.tau = tau\n",
    "    \n",
    "    def process_spikes(self, spikes, current_time):\n",
    "        '''\n",
    "        Processes the received spikes at the current time updating their membrane values.\n",
    "        spikes: list with the indexes of the neurons who spiked.\n",
    "        current_time: the time the neurons spiked (float)\n",
    "        '''\n",
    "        delta_t = current_time-self.times[spikes] # Calculates the difference between the last time they spiked\n",
    "        current_values = self.neurons[spikes]*numpy.exp(-delta_t/self.tau) # Calculates the current values\n",
    "        self.times[spikes]=numpy.ones(len(spikes))*current_time # Updates the last time they spiked\n",
    "        self.neurons[spikes] = current_values + numpy.ones(len(spikes)) # Updates the neuron membrane values\n",
    "\n",
    "    def check_values(self, current_time):\n",
    "        '''\n",
    "        Returns the current membrane values at the specified time.\n",
    "        current_time: time used to calculate the membrane values.\n",
    "        '''\n",
    "        delta_t = current_time-self.times # Calculates the time since last spike\n",
    "        return self.neurons*numpy.exp(-delta_t/self.tau) # Calculates the current values\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import membrane_lowpass_md\n",
    "reload(sys.modules['membrane_lowpass_md']) # Makes sure the interpreter is going to reload the module\n",
    "membrane_lowpass = membrane_lowpass_md.membrane_lowpass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting simulation_2DoF_Arm_md.py\n"
     ]
    }
   ],
   "source": [
    "%%file simulation_2DoF_Arm_md.py\n",
    "# These variables I'm expecting to receive from outsite the namespace of this function\n",
    "# using the func_globals dictionary. It's kind of an ugly hack...\n",
    "# So it's necessary to 'create' them, or an error message is generated.\n",
    "input_gain=0\n",
    "STP_OFF=0\n",
    "noisy=0\n",
    "stratification=0\n",
    "gaussian_pop=0\n",
    "w_SD=0\n",
    "i_noiseA=0\n",
    "offset_noise=0\n",
    "initial_volt=0\n",
    "random_reset=0\n",
    "Net_shape=0\n",
    "Number_of_input_layers=0\n",
    "Number_of_neurons_inputs=0\n",
    "lbd_value=0\n",
    "noisy_currents=0\n",
    "noisy_input_weights=0\n",
    "noisy_input=0\n",
    "input_noise_level=0\n",
    "AII=0\n",
    "AIE=0\n",
    "AEI=0\n",
    "AEE=0\n",
    "\n",
    "killed_neurons=[] # indices of the liquid's neurons that should not spike\n",
    "\n",
    "def simulation_2DoF_Arm(brian_clock, brian):\n",
    "    \"\"\"    \n",
    "    #############################################################################\n",
    "    The variables:\n",
    "    my_seed: value (integer) used as the random seed to keep the liquid structure the same (or not)\n",
    "    input_gain: the multiplier (float) used within the input layer\n",
    "    noisy: 0=>always uses my_seed to start the random generator; 1=>uses the default random state every run\n",
    "    gaussian_pop: 1=>gaussian distribution as the input layers gains; 0=>random weights (just like in Maass's papers)\n",
    "                 -1=>gaussian distribution with an offset noise\n",
    "    w_SD: standard deviation of the gaussian used to distribute the input weights\n",
    "    stratification: 0=>no stratification; 1=>stratified inputs\n",
    "    STP_OFF: True=>No STP; False=>STP\n",
    "    Must be set from outside using the .func_globals dictionary. Example:\n",
    "    simulation_2DoF_Arm.func_globals['my_seed']=93200\n",
    "    #############################################################################\n",
    "    \"\"\"\n",
    "    \n",
    "    import numpy # I could use \"brian.\" because Brian imports numpy, but I prefer not to.\n",
    "    import time\n",
    "    import sys\n",
    "\n",
    "#\n",
    "#     THE VARIABLES BELOW MUST BE SET OUTSIDE OR THE SIMULATION IS NOT GOING TO WORK!\n",
    "#\n",
    "#     Example:\n",
    "#     simulation_2DoF_Arm.func_globals['my_seed']=93200         # This seed is important to always repeat the liquid's structure\n",
    "#     simulation_2DoF_Arm.func_globals['input_gain']=70.0    # Base value used within the input weights\n",
    "#     simulation_2DoF_Arm.func_globals['noisy']=1            # Controls if the liquid is going to use the noisy currents/input weights\n",
    "#     simulation_2DoF_Arm.func_globals['gaussian_pop']=1         # Controls which type of input weight configuration is used (1=>gaussian distributed)\n",
    "#     simulation_2DoF_Arm.func_globals['w_SD']=3.0           # The width of the gaussian used above.\n",
    "#     simulation_2DoF_Arm.func_globals['stratification']=1   # Type of connections: stratified or 30% random\n",
    "#     simulation_2DoF_Arm.func_globals['STP_OFF']=False      # Uses or not STP (False means it uses STP)\n",
    "\n",
    "#     simulation_2DoF_Arm.func_globals['i_noiseA'] = 1.0             # amplitude (dimensionless) of the random currents (updated every time step).\n",
    "#     simulation_2DoF_Arm.func_globals['offset_noise'] = (13.5,14.5) \n",
    "                              # Range of the uniform distribution used to generate the constant offset current.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [14.975nA,15.025nA]\n",
    "                              # Joshi/Maass 2005 does [13.5nA,14.5nA]. Maybe it is to avoid too many spikes without inputs...\n",
    "\n",
    "#     simulation_2DoF_Arm.func_globals['initial_volt'] = (13.5,14.9) \n",
    "                              # Range of the uniform distribution used to generate the initial membrane voltage.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [13.5mV,15.0mV]\n",
    "                              # Joshi/Maass 2005 does [13.5mV,14.9mV], this way it's impossible to generate spikes (14.9mV is below threshold)\n",
    "\n",
    "#     simulation_2DoF_Arm.func_globals['random_reset'] = (13.8,14.5)    # Range of the uniform distribution used to generate the constant reset voltages.\n",
    "#     simulation_2DoF_Arm.func_globals['Net_shape'] = (20,5,6)          # Defines the shape of the network (liquid state machine)\n",
    "#     simulation_2DoF_Arm.func_globals['Number_of_input_layers'] = 6    # It means I will have 6 input layers...\n",
    "#     simulation_2DoF_Arm.func_globals['Number_of_neurons_inputs'] = 50 # ...with 50 neurons each.\n",
    "#     simulation_2DoF_Arm.func_globals['lbd_value'] = 1.2               # lbd controls the connection probability\n",
    "\n",
    "#     simulation_2DoF_Arm.func_globals['noisy_currents'] = True # Injects noisy currents the neurons\n",
    "#     simulation_2DoF_Arm.func_globals['noisy_input_weights'] = True # Injects noise into the input weights\n",
    "#     simulation_2DoF_Arm.func_globals['noisy_input'] = True # Ignores the random seed for the input noise\n",
    "#     simulation_2DoF_Arm.func_globals['input_noise_level'] = 0.1 # 1/SNR\n",
    "\n",
    "      # STP parameters\n",
    "#     simulation_2DoF_Arm.func_globals['AII']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>2.8\n",
    "#     simulation_2DoF_Arm.func_globals['AIE']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>3.0\n",
    "#     simulation_2DoF_Arm.func_globals['AEI']=150.0 # Joshi/Maass2006=>150.0 / Maass2002=>1.6\n",
    "#     simulation_2DoF_Arm.func_globals['AEE']=70.0  # Joshi/Maass2006=>70.0 / Maass2002=>1.2\n",
    "\n",
    "    Number_of_neurons_lsm=Net_shape[0]*Net_shape[1]*Net_shape[2] # Just to avoid calculating this everytime it's used\n",
    "\n",
    "    \n",
    "    # These lines make easier to use the Brian objects without the \"brian.\" at the beginning\n",
    "    ms = brian.ms\n",
    "    mV = brian.mV\n",
    "    nA = brian.nA\n",
    "    nF = brian.nF\n",
    "    NeuronGroup = brian.NeuronGroup\n",
    "    SpikeGeneratorGroup = brian.SpikeGeneratorGroup\n",
    "    Synapses = brian.Synapses\n",
    "    SpikeMonitor = brian.SpikeMonitor\n",
    "    network_operation = brian.network_operation \n",
    "    defaultclock = brian.defaultclock\n",
    "\n",
    "    if my_seed!=0:\n",
    "        # This seed can guarantee that the liquid is going to have the same structure because this object \n",
    "        # is passed to the lsm_connections_probability module.\n",
    "        liq_rs = numpy.random.RandomState(my_seed) # this is the random state used within the liquid structure generation\n",
    "    else:\n",
    "        liq_rs = numpy.random.RandomState() #random liquid\n",
    "    \n",
    "    import lsm_connections_probability # Creates the 3D Grid and the connections according to Maass 2002\n",
    "    reload(sys.modules['lsm_connections_probability']) # Makes sure the interpreter is going to reload the module\n",
    "    lm = lsm_connections_probability\n",
    "\n",
    "    import lsm_dynamical_synapses_v1 # Creates the dynamical synapses according to Maass 2002 and using the output \n",
    "                                     # from the lsm_connections_probability    \n",
    "    reload(sys.modules['lsm_dynamical_synapses_v1']) # Makes sure the interpreter is going to reload the module\n",
    "    ls = lsm_dynamical_synapses_v1\n",
    "\n",
    "    \n",
    "    # Here I'm creating the random state for the rest of the simulation\n",
    "    # This way I can have a fixed liquid structure, but with all the other noisy sources.\n",
    "    if noisy or (my_seed==0):\n",
    "        sim_rs = numpy.random.RandomState()\n",
    "    else:\n",
    "        sim_rs = numpy.random.RandomState(my_seed)\n",
    "\n",
    "    if noisy_input or (my_seed==0):\n",
    "        input_rs = numpy.random.RandomState()\n",
    "    else:\n",
    "        input_rs = numpy.random.RandomState(my_seed)\n",
    "        \n",
    "    initial_time = time.time()\n",
    "\n",
    "    print \"Initial time:\",initial_time\n",
    "\n",
    "    print \"#\"*78\n",
    "    print \"#\"*78\n",
    "    print \"Liquid State Machine - 2 DoF arm experiments!\"\n",
    "    print \"#\"*78\n",
    "    print \"#\"*78\n",
    "\n",
    "\n",
    "    defaultclock = brian_clock  # Receives the clock from the step-by-step simulator\n",
    "                                # I'm setting it to the defaultclock because all connected NeuronGroups \n",
    "                                # must use the same clock.\n",
    "\n",
    "\n",
    "    lsm_3dGrid_flat = numpy.zeros(Number_of_neurons_lsm) \n",
    "    # This creates a numpy 1D array with 'Number_of_neurons_lsm' positions\n",
    "    # I'm using a numpy array to be able to use the reshape method to \n",
    "    # change from 1D (vector) to 3D (matrix)\n",
    "\n",
    "\n",
    "    def randon_connections_gen(Number_of_neurons, ratio, number=None):\n",
    "        '''\n",
    "        Generate the random neuron indexes list according to the number of neurons and the ratio\n",
    "        '''\n",
    "        # List used to generate the randoms 'ratio' indexes for the Liquid\n",
    "        l_range = range(Number_of_neurons)\n",
    "\n",
    "        # Generates the random \"indexes\" of the flattened version of the 3DGrid.\n",
    "        # At each iteration one random item is extracted from l_range and inserted in connection_index.\n",
    "        # This is the way I've found to sample without repetitions.\n",
    "        # - Another way is using the shuffle from numpy and then grabbing the first N values!      \n",
    "        if number==None:\n",
    "            connection_index = [l_range.pop(liq_rs.randint(0,len(l_range))) for i in range(int(Number_of_neurons*ratio))] \n",
    "        else:\n",
    "            connection_index = [l_range.pop(liq_rs.randint(0,len(l_range))) for i in range(int(number))] \n",
    "        connection_index.sort() #This is only useful to make easier to human beings to read the list :)\n",
    "\n",
    "        return connection_index\n",
    "\n",
    "\n",
    "    #\n",
    "    # Number of Inhibitory and Excitatory neurons - LIQUID - 20% of the total neurons\n",
    "    inhibitory_index_L = randon_connections_gen(Number_of_neurons_lsm, 0.2)\n",
    "\n",
    "\n",
    "    # This is the dictionary that has all the connections parameters according to Maass 2002.\n",
    "    # It is necessary to create the 3D connections and the STP configuration matrices\n",
    "    # E=>1 (excitatory) and I=>0 (inhibitory)\n",
    "    # Ex.: (0,0) => II\n",
    "    # Dynamical Synapses Parameters (STP):\n",
    "    Connections_Parameters={\n",
    "                  (0,0):[ # II\n",
    "                          0.1,       # CGupta=0.1        # Parameter used at the connection probability - from Maass2002 paper\n",
    "                          0.32,      # UMarkram=0.32     # Use (U) - Parameter used at the Dynamic Synapse - from Maass2002 paper\n",
    "                          0.144,     # DMarkram=0.144    # Time constant for Depression (tau_rec) - used at the Dynamic Synapse - from Maass2002 paper                    \n",
    "                          0.06,      # FMarkram=0.06     # Time constant for Facilitation (tau_facil) - used at the Dynamic Synapse - from Maass2002 paper\n",
    "                          AII,       # AMaass=2.8        # (nA) In the Maass2002 paper the value is negative, but because I need a positive scale (random.normal parameter) and there is a negative sign in front of the abs function I changed this to positive\n",
    "                          0.8        # Delay_trans = 0.8 # In Maass paper the transmission delay is 0.8 to II, IE and EI        \n",
    "                      ],\n",
    "                  (0,1):[ # IE\n",
    "                          0.4,    # CGupta=0.4\n",
    "                          0.25,   # UMarkram=0.25\n",
    "                          0.7,    # DMarkram=0.7\n",
    "                          0.02,   # FMarkram=0.02\n",
    "                          AIE,    # AMaass=3.0 #in the Maass2002 paper the value is negative, but because I need a positive scale (random.normal parameter) and there is a negative sign in front of the abs function I changed this to positive\n",
    "                          0.8     # Delay_trans = 0.8 # in Maass paper the transmission delay is 0.8 to II, IE and EI\n",
    "                      ],\n",
    "                  (1,0):[ # EI\n",
    "                          0.2,    # CGupta=0.2\n",
    "                          0.05,   # UMarkram=0.05\n",
    "                          0.125,  # DMarkram=0.125\n",
    "                          1.2,    # FMarkram=1.2\n",
    "                          AEI,    # AMaass=1.6\n",
    "                          0.8     # Delay_trans = 0.8 # in Maass paper the transmission delay is 0.8 to II, IE and EI\n",
    "                      ],\n",
    "                  (1,1):[ # EE\n",
    "                          0.3,    # CGupta=0.3 \n",
    "                          0.5,    # UMarkram=0.5\n",
    "                          1.1,    # DMarkram=1.1\n",
    "                          0.05,   # FMarkram=0.05\n",
    "                          AEE,    # AMaass=1.2 #scaling parameter or absolute synaptic efficacy or weight - from Maass2002 paper\n",
    "                          1.5     # Delay_trans = 1.5 # in Maass paper the transmission delay is 1.5 to EE connection\n",
    "                      ]\n",
    "                  }\n",
    "\n",
    "\n",
    "    # Utilizes the functions in the lsm_connections_probability.py\n",
    "    # =>output = {'exc':connections_list_exc,'inh':connections_list_inh, '3Dplot_a':positions_list_a, '3Dplot_b':positions_list_b}\n",
    "    # connections_list_exc= OR connections_list_inh=\n",
    "        # ((i,j), # PRE and POS synaptic neuron indexes\n",
    "        # pconnection, # probability of the connection\n",
    "        # (W_n, U_ds, D_ds, F_ds), # parameters according to Maass2002\n",
    "        # Delay_trans, # parameters according to Maass2002\n",
    "        # connection_type)\n",
    "\n",
    "    # Generate the connections matrix inside the Liquid (Liquid->Liquid) - according to Maass2002\n",
    "    #\n",
    "    print \"Liquid->Liquid connections...\"\n",
    "\n",
    "    print \"Generating the Liquid->Liquid connections...\"\n",
    "    output_L_L = lm.generate_connections(lsm_3dGrid_flat, inhibitory_index_L, Net_shape, \n",
    "                                      CParameters=Connections_Parameters, lbd=lbd_value, randomstate=liq_rs) \n",
    "                                    # lbd controls the connections\n",
    "                                    # randomstate receives a numpy.random.RandomState object, or automatically sets a random one.\n",
    "                                                                    \n",
    "    \n",
    "\n",
    "    print \"Liquid->Liquid connections...Done!\"\n",
    "\n",
    "    \n",
    "    #\n",
    "    # These are the cell (neuron) parameters according to Maass 2002\n",
    "    #\n",
    "    cell_params_lsm = {'cm'          : 30.0,    # [nF]  Capacitance of the membrane \n",
    "                                                # =>>>> MAASS PAPER DOESN'T MENTION THIS PARAMETER DIRECTLY\n",
    "                                                #       but the paper mention a INPUT RESISTANCE OF 1MEGA Ohms and tau_m=RC=30ms, so cm=30nF\n",
    "                       'tau_m'       : 30.0,    # [ms] Membrane time constant => Maass2002\n",
    "                       'tau_refrac_E': 3.0,     # [ms] Duration of refractory period - 3mS for EXCITATORY => Maass2002\n",
    "                       'tau_refrac_I': 2.0,     # [ms] Duration of refractory period - 2mS for INHIBITORY => Maass2002\n",
    "                       'tau_syn_E'   : 3.0,     # [ms] Decay time of excitatory synaptic current => Maass2002\n",
    "                       'tau_syn_I'   : 6.0,     # [ms] Decay time of inhibitory synaptic current => Maass 2002\n",
    "                       'v_reset'     : 13.5,    # [mV] Reset potential after a spike => Maass 2002\n",
    "                       'v_rest'      : 0.0,     # [mV] Resting membrane potential => Maass 2002\n",
    "                       'v_thresh'    : 15.0,    # [mV] Spike threshold => Maass 2002\n",
    "                       'i_noise'     : i_noiseA # [nA] Used in Joshi/Maass 2005: mean 0 and SD=1nA; Maass 2002=>SD=0.2nA\n",
    "                    }\n",
    "\n",
    "    offset_cur = (offset_noise[0],offset_noise[1])  # Values for the uniform distribution of constant offset currents\n",
    "    init_volt = (initial_volt[0],initial_volt[1])   # Values for the uniform distribution of initial membrane voltages\n",
    "    rand_res = (random_reset[0],random_reset[1])    # Values for the uniform distribution of constant random resets\n",
    "\n",
    "\n",
    "    # IF_curr_exp - MODEL EXPLAINED\n",
    "    # Leaky integrate and fire model with fixed threshold and\n",
    "    # decaying-exponential post-synaptic current. \n",
    "    # (Separate synaptic currents for excitatory and inhibitory synapses)\n",
    "    lsm_neuron_eqs='''\n",
    "      dv/dt  = (ie + ii + i_offset + i_noise)/c_m + (v_rest-v)/tau_m : mV\n",
    "      die/dt = -ie/tau_syn_E                : nA\n",
    "      dii/dt = -ii/tau_syn_I                : nA\n",
    "      tau_syn_E                             : ms\n",
    "      tau_syn_I                             : ms\n",
    "      tau_m                                 : ms\n",
    "      c_m                                   : nF\n",
    "      v_rest                                : mV\n",
    "      i_offset                              : nA\n",
    "      i_noise                               : nA\n",
    "      '''\n",
    "\n",
    "\n",
    "\n",
    "    ########################################################################################################################\n",
    "    #\n",
    "    # LIQUID - Setup\n",
    "    #\n",
    "    print \"LIQUID - Setup...\"\n",
    "\n",
    "    # Creates a vector with the corresponding refractory period according to the type of neuron (inhibitory or excitatory)\n",
    "    # IT MUST BE A NUMPY ARRAY OR BRIAN GIVES CRAZY ERRORS!!!!!\n",
    "    refractory_vector = [ cell_params_lsm['tau_refrac_E'] ]*Number_of_neurons_lsm # fills the list with the value corresponding to excitatory neurons\n",
    "    for i in range(Number_of_neurons_lsm):\n",
    "        if i in inhibitory_index_L:\n",
    "            refractory_vector[i]=cell_params_lsm['tau_refrac_I'] # only if the neuron is inibitory, changes the refractory period value!\n",
    "    refractory_vector=numpy.array(refractory_vector)*ms # Here it is converted to a NUMPY ARRAY\n",
    "    \n",
    "\n",
    "    \n",
    "    # This is the population (neurons) used exclusively to the Liquid (pop_lsm).\n",
    "    pop_lsm = NeuronGroup(Number_of_neurons_lsm, model=lsm_neuron_eqs, \n",
    "                                                 threshold=cell_params_lsm['v_thresh']*mV, \n",
    "                                                 reset='v=sim_rs.uniform(rand_res[0],rand_res[1])*mV', \n",
    "                                                 refractory=refractory_vector, \n",
    "                                                 max_refractory=max(cell_params_lsm['tau_refrac_E']*ms, \n",
    "                                                                    cell_params_lsm['tau_refrac_I']*ms))\n",
    "#     else:\n",
    "#         print \"Noisy resets OFF!\"\n",
    "#         # This is the population (neurons) used exclusively to the Liquid (pop_lsm).\n",
    "#         pop_lsm = NeuronGroup(Number_of_neurons_lsm, model=lsm_neuron_eqs, \n",
    "#                                                      threshold=cell_params_lsm['v_thresh']*mV, \n",
    "#                                                      reset=cell_params_lsm['v_reset']*mV,\n",
    "#                                                      refractory=refractory_vector, \n",
    "#                                                      max_refractory=max(cell_params_lsm['tau_refrac_E']*ms, \n",
    "#                                                                         cell_params_lsm['tau_refrac_I']*ms))\n",
    "\n",
    "\n",
    "    # Here I'm mixing numpy.fill with the access of the state variable \"c_m\" in Brian (because Brian is using a numpy.array)\n",
    "    # Sets the value of the capacitance according to the cell_params_lsm (same value to all the neurons)\n",
    "    pop_lsm.c_m.fill(cell_params_lsm['cm']*nF)\n",
    "\n",
    "\n",
    "    # Sets the value of the time constant RC (or membrane constant) according to the cell_params_lsm (same value to all the neurons)\n",
    "    pop_lsm.tau_m.fill(cell_params_lsm['tau_m']*ms)\n",
    "\n",
    "    # Sets the fixed i_offset current.\n",
    "    # The i_offset current is random, but never changes during the same simulation.\n",
    "    # This current, according to Maass 2002, should be drawn from a uniform distr [14.975nA,15.025nA]\n",
    "    # Joshi/Maass 2005 does [13.5nA,14.5nA]. Maybe is to avoid too many spikes without inputs...\n",
    "    pop_lsm.i_offset=sim_rs.uniform(offset_cur[0],offset_cur[1], Number_of_neurons_lsm)*nA\n",
    "\n",
    "    pop_lsm.tau_syn_E.fill(cell_params_lsm['tau_syn_E']*ms) # (same value to all the neurons)\n",
    "    pop_lsm.tau_syn_I.fill(cell_params_lsm['tau_syn_I']*ms) # (same value to all the neurons)\n",
    "\n",
    "    pop_lsm.v_rest.fill(cell_params_lsm['v_rest']*mV) # (same value to all the neurons)\n",
    "\n",
    "    # This current changes (randomly) at each time step (here is only the initialization).\n",
    "    pop_lsm.i_noise=sim_rs.normal(loc=0, scale=cell_params_lsm['i_noise'],size=Number_of_neurons_lsm)*nA\n",
    "\n",
    "    \n",
    "    # Sets the initial membrane voltage. Doesn't change during the simulation.\n",
    "    # According to Maass 2002, this voltage should be drawn from a uniform distr [13.5mV,15.0mV]\n",
    "    # Joshi/Maass 2005 does [13.5mV,14.9mV], this way it's impossible to generate spikes (14.9mV is below threshold)\n",
    "    pop_lsm.v=sim_rs.uniform(init_volt[0],init_volt[1], Number_of_neurons_lsm)*mV\n",
    "\n",
    "    #Kills neurons by setting their initial voltage to a crazily LOW value\n",
    "    pop_lsm.v[killed_neurons]=-1E12*mV\n",
    "    \n",
    "    #\n",
    "    # Loading or creating the Synapses objects used within the Liquid\n",
    "    print \"Liquid->Liquid connections...\"\n",
    "\n",
    "    syn_lsm_obj = ls.LsmConnections(pop_lsm, pop_lsm, output_L_L, nostp=STP_OFF)\n",
    "\n",
    "    # Generates the Liquid->Liquid - EXCITATORY synapses\n",
    "    syn_lsm_exc = syn_lsm_obj.create_synapses('exc')\n",
    "\n",
    "    # Generates the Liquid->Liquid - INHIBITORY synapses\n",
    "    syn_lsm_inh = syn_lsm_obj.create_synapses('inh')\n",
    "    \n",
    "\n",
    "    print \"Liquid->Liquid connections...Done!\"\n",
    "\n",
    "    total_number_of_connections_liquid = len(syn_lsm_exc) + len(syn_lsm_inh)\n",
    "\n",
    "    print \"Number of excitatory synapses in the Liquid: \" + str(len(syn_lsm_exc)) # DEBUG to verify if it is working\n",
    "    print \"Number of inhibitory synapses in the Liquid: \" + str(len(syn_lsm_inh)) # DEBUG to verify if it is working\n",
    "\n",
    "    # To understand what is being returned:\n",
    "    # pop_lsm: it is necessary to connect the neuron network with the rest of the world\n",
    "    # [syn_lsm_obj, syn_lsm_exc, syn_lsm_inh]: to include these objects at the simulation (net=Net(...); net.run(total_sim_time*ms)); \n",
    "    # It is a list because is easy to concatenate lists :D\n",
    "\n",
    "    print \"LIQUID - Setup...Done!\"\n",
    "\n",
    "    #\n",
    "    # End of the LIQUID - Setup\n",
    "    ########################################################################################################################\n",
    "\n",
    "\n",
    "\n",
    "    ########################################################################################################################\n",
    "    #\n",
    "    # INPUT - Setup\n",
    "    #\n",
    "    print \"INPUT - Setup...\"\n",
    "  \n",
    "    spiketimes = [] # The spikes are going to be received during the simulation, \n",
    "                    # so this is always an empty list when using the step_by_step_brian_sim!\n",
    "    \n",
    "    # I'm using only one big input layer because Brian docs say it is better for the performance\n",
    "    SpikeInputs = SpikeGeneratorGroup(Number_of_input_layers*Number_of_neurons_inputs, spiketimes)\n",
    "    \n",
    "\n",
    "    #\n",
    "    #\n",
    "    # Here the synapses are created. The synapses created are ALWAYS excitatory because it is \n",
    "    # connecting through 'ie' in the neuron model!\n",
    "\n",
    "    syn_world_Input = Synapses(SpikeInputs, pop_lsm,\n",
    "                                         model='''w : 1''',\n",
    "                                         pre='''ie+=w''')\n",
    "\n",
    "\n",
    "    weights_input_liquid = [] # remember that the weights must follow the same order of the creation of synapses\n",
    "\n",
    "\n",
    "    def gaussian(lamb,n,nt):\n",
    "        '''\n",
    "        Generates a gaussian centered at 'n'\n",
    "        '''\n",
    "        # return input_gain*(1/(lamb*numpy.sqrt(2*numpy.pi)))*numpy.exp(-((nt-n)**2)/(2*(lamb)**2))  # Energy normalized version\n",
    "        return input_gain*numpy.exp(-((nt-n)**2)/(2*(lamb)**2))  # No energy normalization\n",
    "  \n",
    "    def simple_inputs(*args):\n",
    "        '''\n",
    "        This function ignores the arguments and return a random value from a gaussian distribution with mean input_gain and SD=input_gain/2.0\n",
    "        '''\n",
    "        return (abs(sim_rs.normal(loc=input_gain, scale=input_gain/2.0))) # The \"abs\" function is to guarantee all inputs are excitatories!    \n",
    "\n",
    "    def gaussian_noise(lamb,n,nt):\n",
    "        '''\n",
    "        Generates a gaussian centered at 'n' with a background noise (1/3 of the amplitude)\n",
    "        '''\n",
    "        return 3*input_gain*numpy.exp(-((nt-n)**2)/(2*(lamb)**2)) + (abs(sim_rs.normal(loc=input_gain, scale=input_gain/2.0)))\n",
    "\n",
    "\n",
    "    # Verifies the user selection and sets the proper weight generation function.\n",
    "    if gaussian_pop==1:\n",
    "        weight_func = gaussian\n",
    "    elif gaussian_pop==0:\n",
    "        weight_func = simple_inputs\n",
    "    else:\n",
    "        weight_func = gaussian_noise\n",
    "        \n",
    "\n",
    "\n",
    "    if stratification==0:\n",
    "        # List with the indexes of all the excitatory neurons in the liquid\n",
    "        excitatory_index_L = [i for i in range(Number_of_neurons_lsm) if i not in inhibitory_index_L]\n",
    "\n",
    "        # Here I connect the input neurons only to excitatories neurons in the liquid:\n",
    "        sim_rs.shuffle(excitatory_index_L) # Shuffles the excitatory index vector.\n",
    "        rand_connections = excitatory_index_L[:int(len(excitatory_index_L)*0.3)] # Gets randomly 30% of the excitatory connections\n",
    "\n",
    "        # Here I connect the input neurons to any type of neurons in the liquid:\n",
    "        # rand_connections = randon_connections_gen(Number_of_neurons_lsm, 0, number=int(Number_of_neurons_lsm*0.3)) # generates random connections to 30% of the neurons in the liquid!\n",
    "\n",
    "        # Goes through the liquid to generate the proper connections\n",
    "        for inp in range(Number_of_input_layers):\n",
    "            for i in range(inp*Number_of_neurons_inputs,Number_of_neurons_inputs*(inp+1)):\n",
    "                for j,ji in zip(rand_connections,range(len(rand_connections))):\n",
    "                    syn_world_Input[i,j] = True # So it is one-to-random NoIN neurons, each input neuron is connect to all the \"input layer\" of the liquid.\n",
    "                    # All inputs have the same connections to the Liquid (I would say that they are all connect to the input layer of the liquid)\n",
    "                    # If they are all connected to the same neurons, seems to me that is less probable that the readout is going to learn\n",
    "                    # only to filter the input...\n",
    "                    centre_position=(i-(inp*Number_of_neurons_inputs))*(len(rand_connections)-1)/float(Number_of_neurons_inputs)\n",
    "                    weights_input_liquid.append(weight_func(w_SD,centre_position,ji)*nA)\n",
    "    else:\n",
    "\n",
    "        # Goes through the liquid to generate the proper connections, but dividing the liquid into the same number of input layers\n",
    "        liquid_input_layer_size = int(Number_of_neurons_lsm/float(Number_of_input_layers))\n",
    "        for inp in range(Number_of_input_layers):\n",
    "            for i in range(inp*Number_of_neurons_inputs,Number_of_neurons_inputs*(inp+1)):\n",
    "                for j,ji in zip(range(inp*liquid_input_layer_size,liquid_input_layer_size*(inp+1)),range(liquid_input_layer_size)):\n",
    "                    if j not in inhibitory_index_L:\n",
    "                        syn_world_Input[i,j] = True\n",
    "                        centre_position=(i-(inp*Number_of_neurons_inputs))*(liquid_input_layer_size-1)/float(Number_of_neurons_inputs)\n",
    "                        weights_input_liquid.append(weight_func(w_SD,centre_position,ji)*nA)\n",
    "                        \n",
    "\n",
    "\n",
    "\n",
    "  \n",
    "    weights_input_liquid = numpy.array(weights_input_liquid)\n",
    "    \n",
    "    \n",
    "    syn_world_Input.w = weights_input_liquid\n",
    "    syn_world_Input.delay=0*ms\n",
    "\n",
    "    print \"INPUT - Setup...Done!\"\n",
    "\n",
    "    #\n",
    "    # End of the INPUT - Setup (creation of the connections between the Poisson input and the Liquid!)\n",
    "    #\n",
    "    ########################################################################################################################\n",
    "\n",
    "\n",
    "    # Generates the noisy current at each time step (as seen in Joshi/Maass 2005)\n",
    "    @network_operation(clock=defaultclock)\n",
    "    def generate_i_noise():\n",
    "        if noisy_currents:\n",
    "            # These are the noise currents inside each liquid's neuron\n",
    "            pop_lsm.i_noise=cell_params_lsm['i_noise']*sim_rs.normal(loc=0,scale=1,size=Number_of_neurons_lsm)*nA\n",
    "        \n",
    "        if noisy_input_weights:\n",
    "            # This is my new version where the noise level can be controlled\n",
    "            syn_world_Input.w[:]=weights_input_liquid+((input_noise_level*input_gain)*input_rs.normal(loc=0,scale=1,size=len(weights_input_liquid))*nA)        \n",
    "\n",
    "\n",
    "    populations_sim = [pop_lsm, SpikeInputs]\n",
    "\n",
    "    synapses_sim = [syn_lsm_exc, syn_lsm_inh, syn_world_Input]\n",
    "\n",
    "    monitors_sim = [generate_i_noise] \n",
    "\n",
    "    Input_layer, Output_layer, pop_objects, syn_objects, monitors_objects = SpikeInputs, pop_lsm, populations_sim, synapses_sim, monitors_sim\n",
    "\n",
    "    print \"Setup time:\", time.time()-initial_time\n",
    "    \n",
    "    # Returns these objects to be used with the step_by_step_brian_sim class\n",
    "    return Input_layer, Output_layer, pop_objects, syn_objects, monitors_objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main functions to the 2 joint arm simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'simulation_2DoF_Arm_physics' from 'simulation_2DoF_Arm_physics.pyc'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from simulation_2DoF_Arm_physics import *\n",
    "reload(sys.modules['simulation_2DoF_Arm_physics']) # Makes sure the interpreter is going to reload the module"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Start of the simulation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def generates_input_spikes_lsm(delta_t = 200):\n",
    "    # First the arm parameters according to:\n",
    "    # Joshi, Prashant, and Wolfgang Maass. “Movement Generation with Circuits of Spiking Neurons.” \n",
    "    # Neural Computation 17, no. 8 (2005): 1715–1738.\n",
    "\n",
    "    aparams = {\n",
    "        'l1' : 0.5, # metres\n",
    "        'l2' : 0.5,\n",
    "        'lc1' : 0.25,\n",
    "        'lc2' : 0.25,\n",
    "        'm1' : 1.0, # kg\n",
    "        'm2' : 1.0,\n",
    "        'i1' : 0.03, # kg*m*m\n",
    "        'i2' : 0.03\n",
    "    }\n",
    "\n",
    "    # Total time spent during the movement (in seconds)\n",
    "    MT = 0.5 # 500ms\n",
    "\n",
    "    # Simulation time step (in seconds)\n",
    "    time_step = 2/1000.0 # 2ms\n",
    "\n",
    "    # Number of neurons at each input layer (defines the resolution of the system)\n",
    "    Ninput = 50\n",
    "\n",
    "    t_mov=numpy.arange(0, MT, time_step) # t starts in 0s and steps time_step(s) until reaches MT(s)\n",
    "\n",
    "    # Reads all the experiments (generated using the 2DofArm_simulation_data_generator.ipynb) \n",
    "    # to check which range of torques is necessary to control the arm.\n",
    "\n",
    "    # These values are based ONLY in the \"number_of_trajectories\" trajectories used!\n",
    "    torques_all_1 = numpy.array([]);\n",
    "    torques_all_2 = numpy.array([]);\n",
    "    for i in xrange(1,number_of_trajectories+1):\n",
    "        states,torques = slf.load_from_file_gz(\"./\"+base_dir+\"/\"+sim_set+\"/States_Torques_movement\"+str(i)+\".gzpickle\")\n",
    "        torques_all_1 = numpy.concatenate((torques_all_1,torques[:,0])) # torques joint 1\n",
    "        torques_all_2 = numpy.concatenate((torques_all_2,torques[:,1])) # torques joint 2\n",
    "\n",
    "    # Creates the 50 possible discrete values of torques for each joint (probably they will be different)\n",
    "    torques_values_1=numpy.linspace(torques_all_1.min(),torques_all_1.max(),num=50)\n",
    "    torques_values_2=numpy.linspace(torques_all_2.min(),torques_all_2.max(),num=50)\n",
    "\n",
    "    # Creates the 50 possible discrete values of angles for each joint\n",
    "    teta1_values=numpy.linspace(-numpy.pi/6,numpy.pi,num=50)\n",
    "    teta2_values=numpy.linspace(0,numpy.pi,num=50)\n",
    "\n",
    "    # These are the 50 possible values that the input neurons will represent (x,y) \n",
    "    x_values=numpy.linspace(-1,1,num=50)\n",
    "    y_values=numpy.linspace(-1,1,num=50)\n",
    "\n",
    "    # These are the normalized starting positions\n",
    "    # xstart_normalized =  x_values[abs(x_values-xstart).argmin()]\n",
    "    # ystart_normalized =  y_values[abs(y_values-ystart).argmin()]\n",
    "\n",
    "    # The data file below was created using the iPython notebook: \"2DofArm_simulation_data_generator.ipynb\"\n",
    "    states,torques = slf.load_from_file_gz(\"./\"+base_dir+\"/\"+sim_set+\"/States_Torques_movement\"+str(tji)+\".gzpickle\")\n",
    "\n",
    "    # These are the normalized torques according to the minimum and maximum torques available in the set (experiments)\n",
    "    torques_normalized = numpy.array([[torques_values_1[abs(torques_values_1-xy[0]).argmin()], torques_values_2[abs(torques_values_2-xy[1]).argmin()]] for xy in torques])\n",
    "\n",
    "    # These are the normalized teta1 and teta2 angles in a matrix format (numpy.array)\n",
    "    states_normalized = numpy.array([[teta1_values[abs(teta1_values-st[0]).argmin()], teta2_values[abs(teta2_values-st[1]).argmin()]] for st in states])\n",
    "\n",
    "    # HERE SEEMS TO BE THE ONLY PLACE I COULD INJECT AN \"ANALOG\" NOISE INTO THE INPUTS.\n",
    "\n",
    "    # Indexes of the normalized torque values to use in the LSM simulation \n",
    "    # Remember: the LSM has only 50 neurons to encode each variable and the LSM receives only indexes!\n",
    "    torques_normalized_idx = numpy.array([[abs(torques_values_1-xy[0]).argmin(), abs(torques_values_2-xy[1]).argmin()] for xy in torques])\n",
    "\n",
    "    # Indexes of the normalized joint angle values to use in the LSM simulation \n",
    "    joint_states_idx = numpy.array([[abs(teta1_values-xy[0]).argmin(), abs(teta2_values-xy[1]).argmin()] for xy in states[:,[0,1]]])\n",
    "    \n",
    "    # The input (one big neurongroup with 300 neurons) will be divided like this:\n",
    "    # 6 groups of 50 neurons.\n",
    "    # - Group 1: xdest => offset:0\n",
    "    # - Group 2: ydest => offset:50\n",
    "    # - Group 3: teta1 => offset:100\n",
    "    # - Group 4: teta2 => offset:150\n",
    "    # - Group 5: tau1  => offset:200\n",
    "    # - Group 6: tau2  => offset:250\n",
    "\n",
    "    #\n",
    "    # Generates the input values used in the LSM simulation.\n",
    "    #\n",
    "    \n",
    "    ############# Using ONLY the initial (start) position\n",
    "    # Indexes of the normalized INITIAL position values to use in the LSM simulation\n",
    "    # xstart_idx =  abs(x_values-xstart).argmin()\n",
    "    # ystart_idx =  abs(y_values-ystart).argmin()\n",
    "    \n",
    "    # x_lsm = xstart_idx + 0\n",
    "    # y_lsm = ystart_idx + 50    \n",
    "\n",
    "    ############# Using ONLY the final (dest) position\n",
    "    # Indexes of the normalized FINAL position values to use in the LSM simulation    \n",
    "    xdest_idx =  abs(x_values-xdest).argmin()\n",
    "    ydest_idx =  abs(y_values-ydest).argmin()    \n",
    "    \n",
    "    x_lsm = xdest_idx + 0\n",
    "    y_lsm = ydest_idx + 50\n",
    "    \n",
    "    \n",
    "    teta1_lsm = joint_states_idx[:,0] + 100\n",
    "    teta2_lsm = joint_states_idx[:,1] + 150\n",
    "    tau1_lsm = torques_normalized_idx[:,0] + 200\n",
    "    tau2_lsm = torques_normalized_idx[:,1] + 250\n",
    "\n",
    "#     delta_t = 200 # Time delay imposed to the joint angles - proprioceptive (in ms)\n",
    "    delta_t_idx = int(delta_t/(1000*time_step)) # Number of steps the delay (delta_t) takes\n",
    "\n",
    "    input_spikes_pos = [[x_lsm, y_lsm]]\n",
    "\n",
    "    if delta_t_idx != 0:\n",
    "        input_spikes_teta = [[]]\n",
    "    else:\n",
    "        input_spikes_teta = [[teta1_lsm[0], teta2_lsm[0]]]\n",
    "        \n",
    "    input_spikes_tau = [[tau1_lsm[0], tau2_lsm[0]]]\n",
    "    \n",
    "    if delta_t_idx!=0:\n",
    "        delta_t_idx_counter = 1 # It's 1 because input_spikes already has one set of values\n",
    "    else:\n",
    "        delta_t_idx_counter = 0\n",
    "\n",
    "    for idx in range(1, len(tau1_lsm)):\n",
    "        if delta_t_idx_counter==delta_t_idx or delta_t_idx==0:\n",
    "            # Here is where I insert the proprioceptive information about the joint angles after the delta_t time\n",
    "            input_spikes_pos.append([x_lsm, y_lsm])\n",
    "            input_spikes_teta.append([teta1_lsm[idx-delta_t_idx_counter], teta2_lsm[idx-delta_t_idx_counter]])\n",
    "            input_spikes_tau.append([tau1_lsm[idx], tau2_lsm[idx]])\n",
    "        else:\n",
    "            input_spikes_pos.append([x_lsm, y_lsm])\n",
    "            input_spikes_teta.append([])\n",
    "            input_spikes_tau.append([tau1_lsm[idx], tau2_lsm[idx]])        \n",
    "            delta_t_idx_counter+=1\n",
    "\n",
    "    input_spikes = (input_spikes_pos,input_spikes_teta,input_spikes_tau)\n",
    "    return input_spikes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set C - WITH STP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@dview.parallel(block=True)\n",
    "def run_simulations_2(sim_num): \n",
    "    import numpy\n",
    "    import time\n",
    "    import sys\n",
    "    import save_load_file as slf\n",
    "        \n",
    "    import step_by_step_brian\n",
    "    reload(sys.modules['step_by_step_brian']) # Makes sure the interpreter is going to reload the module\n",
    "    step_by_step_brian_sim = step_by_step_brian.step_by_step_brian_sim\n",
    "    \n",
    "    import simulation_2DoF_Arm_md  \n",
    "    reload(sys.modules['simulation_2DoF_Arm_md']) # Makes sure the interpreter is going to reload the module\n",
    "    simulation_2DoF_Arm = simulation_2DoF_Arm_md.simulation_2DoF_Arm    \n",
    "    \n",
    "    pos_i, run_i, input_spikes, folder_names = sim_num\n",
    "\n",
    "    sim_set,base_dir = folder_names\n",
    "    \n",
    "    input_spikes_pos,input_spikes_teta,input_spikes_tau = input_spikes\n",
    "    \n",
    "    # Initializing the simulation...\n",
    "    \n",
    "    # Simulation's variables:\n",
    "    simulation_2DoF_Arm.func_globals['my_seed']=93200      # This seed is important to always repeat the liquid's structure\n",
    "    simulation_2DoF_Arm.func_globals['input_gain']=70.0    # Base value used within the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy']=1            # Controls if the liquid is going to use the noisy currents/input weights\n",
    "    simulation_2DoF_Arm.func_globals['gaussian_pop']=1     # Controls which type of input weight configuration is used (1=>gaussian distributed)\n",
    "    simulation_2DoF_Arm.func_globals['w_SD']=3.0           # The width of the gaussian used above.\n",
    "    simulation_2DoF_Arm.func_globals['stratification']=1   # Type of connections: stratified or 30% random\n",
    "    simulation_2DoF_Arm.func_globals['STP_OFF']=False      # Uses or not STP (False means it uses STP)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['i_noiseA'] = 1.0/100     # amplitude (dimensionless) of the random currents (updated every time step).\n",
    "    simulation_2DoF_Arm.func_globals['offset_noise'] = (13.5/100,14.5/100) \n",
    "                              # Range of the uniform distribution used to generate the constant offset current.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [14.975nA,15.025nA]\n",
    "                              # Joshi/Maass 2005 does [13.5nA,14.5nA]. Maybe it is to avoid too many spikes without inputs...\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['initial_volt'] = (13.5,14.9) \n",
    "                              # Range of the uniform distribution used to generate the initial membrane voltage.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [13.5mV,15.0mV]\n",
    "                              # Joshi/Maass 2005 does [13.5mV,14.9mV], this way it's impossible to generate spikes (14.9mV is below threshold)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['random_reset'] = (13.8,14.5)    # Range of the uniform distribution used to generate the constant reset voltages.\n",
    "    simulation_2DoF_Arm.func_globals['Net_shape'] = (20,5,6)          # Defines the shape of the network (liquid state machine)\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_input_layers'] = 6    # It means I will have 6 input layers...\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_neurons_inputs'] = 50 # ...with 50 neurons each.\n",
    "    simulation_2DoF_Arm.func_globals['lbd_value'] = 1.2               # lbd controls the connection probability\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['noisy_currents'] = True # Injects noisy currents the neurons\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input_weights'] = False # Injects noise into the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input'] = True # Ignores the random seed for the input noise\n",
    "    simulation_2DoF_Arm.func_globals['input_noise_level'] = 0 # 1/SNR\n",
    "    \n",
    "    # STP parameters\n",
    "    simulation_2DoF_Arm.func_globals['AII']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>2.8\n",
    "    simulation_2DoF_Arm.func_globals['AIE']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>3.0\n",
    "    simulation_2DoF_Arm.func_globals['AEI']=150.0 # Joshi/Maass2006=>150.0 / Maass2002=>1.6\n",
    "    simulation_2DoF_Arm.func_globals['AEE']=70.0  # Joshi/Maass2006=>70.0 / Maass2002=>1.2\n",
    "    \n",
    "    \n",
    "    maximum_input_index_noise=1 # This is the amount of noise inserted in the feedback loop as variations in the input neurons that spike\n",
    "    my_noisy_input=1            # Controls if the input spike noise is going to use my_seed(0) or just be randomly noisy(1)\n",
    "                                # This way I can have a simulation where all the other noise sources are kept the same, but this.\n",
    "    \n",
    "    if my_noisy_input:\n",
    "        input_rs = numpy.random.RandomState()\n",
    "    else:\n",
    "        input_rs = numpy.random.RandomState(simulation_2DoF_Arm.func_globals['my_seed'])\n",
    "    \n",
    "    step_size = 2\n",
    "    \n",
    "    # Starting the simulation object!\n",
    "    s = step_by_step_brian_sim(simulation_2DoF_Arm, step_size)\n",
    "\n",
    "    init_t = time.time()\n",
    "    simulated_values = []\n",
    "#####DEBUG\n",
    "    inputs = []\n",
    "#####END_DEBUG    \n",
    "\n",
    "    # Runs the simulation going step-by-step into the whole input_spikes list!        \n",
    "    for ipos,iteta,itau in zip(input_spikes_pos,input_spikes_teta,input_spikes_tau):\n",
    "        if maximum_input_index_noise>0:\n",
    "            if len(iteta):\n",
    "                simulated_values.append(s.run_step([ipos[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),ipos[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)]+[iteta[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),iteta[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)]+[itau[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),itau[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)])) # Sends the spikes and runs to the second yield \n",
    "            else:\n",
    "                simulated_values.append(s.run_step([ipos[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),ipos[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)]+iteta+[itau[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),itau[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)])) # Sends the spikes and runs to the second yield \n",
    "        else:\n",
    "            simulated_values.append(s.run_step(ipos+iteta+itau)) # Sends the spikes and runs to the second yield \n",
    "#####DEBUG\n",
    "            inputs.append(ipos+iteta+itau)\n",
    "#####END_DEBUG                \n",
    "\n",
    "    # This last step is to make sure the result from the last input was processed and then the simulation\n",
    "    # can reach the final time step:\n",
    "    simulated_values.append(s.run_step([])) # Sends the spikes and runs to the second yield \n",
    "    \n",
    "    print \"Total simulation time (seconds):\",time.time()-init_t\n",
    "\n",
    "    if simulation_2DoF_Arm.func_globals['STP_OFF']:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\".gzpickle\"\n",
    "    else:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\"_STP_ON.gzpickle\"\n",
    "\n",
    "    slf.save_to_file_gz(simulated_values,\"./\"+base_dir+\"/\"+sim_set+filename)\n",
    "    \n",
    "#####DEBUG\n",
    "    slf.save_to_file_gz(inputs,\"./\"+base_dir+\"/\"+sim_set+\"/States_Torques_movement\"+str(pos_i)+\"_INPUTS_LSM_\"+str(run_i)+\".gzpickle\")\n",
    "#####END_DEBUG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "base_dir = \"2DofArm_simulation_data\"\n",
    "sim_set = \"set_C\"\n",
    "\n",
    "\n",
    "total_trials = 20 # total number of trials\n",
    "\n",
    "trajectories = [[[0.75,0.25],[0.0,0.5]], [[0.25,0.60],[-0.25,0.60]], [[-0.10,0.75],[-0.10,0.25]],[[-0.75,0.50],[-0.40,0.00]]]\n",
    "\n",
    "number_of_trajectories = len(trajectories)\n",
    "\n",
    "for tji, positions in [(tji,positions) for tji,positions in zip(range(1,len(trajectories)+1),trajectories)]:\n",
    "\n",
    "    xstart,ystart = positions[0]\n",
    "    xdest,ydest = positions[1]\n",
    "\n",
    "    # The default value to delta_t is 200, so when simulating the same as Joshi/Maass 2006 it's \n",
    "    # not necessary to especify the delta_t.\n",
    "    input_spikes = generates_input_spikes_lsm(delta_t=200)\n",
    "\n",
    "    print \"Trajectory Number:\",tji\n",
    "    %time results = run_simulations_2.map([(tji,i,input_spikes,(sim_set,base_dir)) for i in range(1,total_trials+1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set D - WITHOUT STP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@dview.parallel(block=True)\n",
    "def run_simulations_2(sim_num): \n",
    "    import numpy\n",
    "    import time\n",
    "    import sys\n",
    "    import save_load_file as slf\n",
    "        \n",
    "    import step_by_step_brian\n",
    "    reload(sys.modules['step_by_step_brian']) # Makes sure the interpreter is going to reload the module\n",
    "    step_by_step_brian_sim = step_by_step_brian.step_by_step_brian_sim\n",
    "    \n",
    "    import simulation_2DoF_Arm_md  \n",
    "    reload(sys.modules['simulation_2DoF_Arm_md']) # Makes sure the interpreter is going to reload the module\n",
    "    simulation_2DoF_Arm = simulation_2DoF_Arm_md.simulation_2DoF_Arm    \n",
    "    \n",
    "    pos_i, run_i, input_spikes, folder_names = sim_num\n",
    "\n",
    "    sim_set,base_dir = folder_names\n",
    "    \n",
    "    input_spikes_pos,input_spikes_teta,input_spikes_tau = input_spikes\n",
    "    \n",
    "    # Initializing the simulation...\n",
    "    \n",
    "    # Simulation's variables:\n",
    "    simulation_2DoF_Arm.func_globals['my_seed']=93200      # This seed is important to always repeat the liquid's structure\n",
    "    simulation_2DoF_Arm.func_globals['input_gain']=70.0    # Base value used within the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy']=1            # Controls if the liquid is going to use the noisy currents/input weights\n",
    "    simulation_2DoF_Arm.func_globals['gaussian_pop']=1     # Controls which type of input weight configuration is used (1=>gaussian distributed)\n",
    "    simulation_2DoF_Arm.func_globals['w_SD']=3.0           # The width of the gaussian used above.\n",
    "    simulation_2DoF_Arm.func_globals['stratification']=1   # Type of connections: stratified or 30% random\n",
    "    simulation_2DoF_Arm.func_globals['STP_OFF']=True       # Uses or not STP (False means it uses STP)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['i_noiseA'] = 1.0/100     # amplitude (dimensionless) of the random currents (updated every time step).\n",
    "    simulation_2DoF_Arm.func_globals['offset_noise'] = (13.5/100,14.5/100) \n",
    "                              # Range of the uniform distribution used to generate the constant offset current.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [14.975nA,15.025nA]\n",
    "                              # Joshi/Maass 2005 does [13.5nA,14.5nA]. Maybe it is to avoid too many spikes without inputs...\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['initial_volt'] = (13.5,14.9) \n",
    "                              # Range of the uniform distribution used to generate the initial membrane voltage.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [13.5mV,15.0mV]\n",
    "                              # Joshi/Maass 2005 does [13.5mV,14.9mV], this way it's impossible to generate spikes (14.9mV is below threshold)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['random_reset'] = (13.8,14.5)    # Range of the uniform distribution used to generate the constant reset voltages.\n",
    "    simulation_2DoF_Arm.func_globals['Net_shape'] = (20,5,6)          # Defines the shape of the network (liquid state machine)\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_input_layers'] = 6    # It means I will have 6 input layers...\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_neurons_inputs'] = 50 # ...with 50 neurons each.\n",
    "    simulation_2DoF_Arm.func_globals['lbd_value'] = 1.2               # lbd controls the connection probability\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['noisy_currents'] = True # Injects noisy currents the neurons\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input_weights'] = False # Injects noise into the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input'] = True # Ignores the random seed for the input noise\n",
    "    simulation_2DoF_Arm.func_globals['input_noise_level'] = 0 # 1/SNR\n",
    "    \n",
    "    # STP parameters\n",
    "    simulation_2DoF_Arm.func_globals['AII']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>2.8\n",
    "    simulation_2DoF_Arm.func_globals['AIE']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>3.0\n",
    "    simulation_2DoF_Arm.func_globals['AEI']=150.0 # Joshi/Maass2006=>150.0 / Maass2002=>1.6\n",
    "    simulation_2DoF_Arm.func_globals['AEE']=70.0  # Joshi/Maass2006=>70.0 / Maass2002=>1.2\n",
    "    \n",
    "    \n",
    "    maximum_input_index_noise=1 # This is the amount of noise inserted in the feedback loop as variations in the input neurons that spike\n",
    "    my_noisy_input=1            # Controls if the input spike noise is going to use my_seed(0) or just be randomly noisy(1)\n",
    "                                # This way I can have a simulation where all the other noise sources are kept the same, but this.\n",
    "    \n",
    "    if my_noisy_input:\n",
    "        input_rs = numpy.random.RandomState()\n",
    "    else:\n",
    "        input_rs = numpy.random.RandomState(simulation_2DoF_Arm.func_globals['my_seed'])\n",
    "    \n",
    "    step_size = 2\n",
    "    \n",
    "    # Starting the simulation object!\n",
    "    s = step_by_step_brian_sim(simulation_2DoF_Arm, step_size)\n",
    "\n",
    "    init_t = time.time()\n",
    "    simulated_values = []\n",
    "#####DEBUG\n",
    "    inputs = []\n",
    "#####END_DEBUG    \n",
    "\n",
    "    # Runs the simulation going step-by-step into the whole input_spikes list!        \n",
    "    for ipos,iteta,itau in zip(input_spikes_pos,input_spikes_teta,input_spikes_tau):\n",
    "        if maximum_input_index_noise>0:\n",
    "            if len(iteta):\n",
    "                simulated_values.append(s.run_step([ipos[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),ipos[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)]+[iteta[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),iteta[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)]+[itau[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),itau[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)])) # Sends the spikes and runs to the second yield \n",
    "            else:\n",
    "                simulated_values.append(s.run_step([ipos[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),ipos[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)]+iteta+[itau[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),itau[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)])) # Sends the spikes and runs to the second yield \n",
    "        else:\n",
    "            simulated_values.append(s.run_step(ipos+iteta+itau)) # Sends the spikes and runs to the second yield \n",
    "#####DEBUG\n",
    "            inputs.append(ipos+iteta+itau)\n",
    "#####END_DEBUG                \n",
    "\n",
    "    # This last step is to make sure the result from the last input was processed and then the simulation\n",
    "    # can reach the final time step:\n",
    "    simulated_values.append(s.run_step([])) # Sends the spikes and runs to the second yield \n",
    "    \n",
    "    print \"Total simulation time (seconds):\",time.time()-init_t\n",
    "\n",
    "    if simulation_2DoF_Arm.func_globals['STP_OFF']:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\".gzpickle\"\n",
    "    else:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\"_STP_ON.gzpickle\"\n",
    "\n",
    "    slf.save_to_file_gz(simulated_values,\"./\"+base_dir+\"/\"+sim_set+filename)\n",
    "    \n",
    "#####DEBUG\n",
    "    slf.save_to_file_gz(inputs,\"./\"+base_dir+\"/\"+sim_set+\"/States_Torques_movement\"+str(pos_i)+\"_INPUTS_LSM_\"+str(run_i)+\".gzpickle\")\n",
    "#####END_DEBUG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "base_dir = \"2DofArm_simulation_data\"\n",
    "sim_set = \"set_D\"\n",
    "\n",
    "\n",
    "total_trials = 20 # total number of trials\n",
    "\n",
    "trajectories = [[[0.75,0.25],[0.0,0.5]], [[0.25,0.60],[-0.25,0.60]], [[-0.10,0.75],[-0.10,0.25]],[[-0.75,0.50],[-0.40,0.00]]]\n",
    "\n",
    "number_of_trajectories = len(trajectories)\n",
    "\n",
    "for tji, positions in [(tji,positions) for tji,positions in zip(range(1,len(trajectories)+1),trajectories)]:\n",
    "\n",
    "    xstart,ystart = positions[0]\n",
    "    xdest,ydest = positions[1]\n",
    "\n",
    "    # The default value to delta_t is 200, so when simulating the same as Joshi/Maass 2006 it's \n",
    "    # not necessary to especify the delta_t.\n",
    "    input_spikes = generates_input_spikes_lsm(delta_t=200)\n",
    "\n",
    "    print \"Trajectory Number:\",tji\n",
    "    %time results = run_simulations_2.map([(tji,i,input_spikes,(sim_set,base_dir)) for i in range(1,total_trials+1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set A - WITH STP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@dview.parallel(block=True)\n",
    "def run_simulations_2(sim_num): \n",
    "    import numpy\n",
    "    import time\n",
    "    import sys\n",
    "    import save_load_file as slf\n",
    "        \n",
    "    import step_by_step_brian\n",
    "    reload(sys.modules['step_by_step_brian']) # Makes sure the interpreter is going to reload the module\n",
    "    step_by_step_brian_sim = step_by_step_brian.step_by_step_brian_sim\n",
    "    \n",
    "    import simulation_2DoF_Arm_md  \n",
    "    reload(sys.modules['simulation_2DoF_Arm_md']) # Makes sure the interpreter is going to reload the module\n",
    "    simulation_2DoF_Arm = simulation_2DoF_Arm_md.simulation_2DoF_Arm    \n",
    "    \n",
    "    pos_i, run_i, input_spikes, folder_names = sim_num\n",
    "\n",
    "    sim_set,base_dir = folder_names\n",
    "    \n",
    "    input_spikes_pos,input_spikes_teta,input_spikes_tau = input_spikes\n",
    "    \n",
    "    # Initializing the simulation...\n",
    "    \n",
    "    # Simulation's variables:\n",
    "    simulation_2DoF_Arm.func_globals['my_seed']=93200      # This seed is important to always repeat the liquid's structure\n",
    "    simulation_2DoF_Arm.func_globals['input_gain']=70.0    # Base value used within the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy']=1            # Controls if the liquid is going to use the noisy currents/input weights\n",
    "    simulation_2DoF_Arm.func_globals['gaussian_pop']=1     # Controls which type of input weight configuration is used (1=>gaussian distributed)\n",
    "    simulation_2DoF_Arm.func_globals['w_SD']=3.0           # The width of the gaussian used above.\n",
    "    simulation_2DoF_Arm.func_globals['stratification']=1   # Type of connections: stratified or 30% random\n",
    "    simulation_2DoF_Arm.func_globals['STP_OFF']=False      # Uses or not STP (False means it uses STP)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['i_noiseA'] = 1.0             # amplitude (dimensionless) of the random currents (updated every time step).\n",
    "    simulation_2DoF_Arm.func_globals['offset_noise'] = (13.5,14.5) \n",
    "                              # Range of the uniform distribution used to generate the constant offset current.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [14.975nA,15.025nA]\n",
    "                              # Joshi/Maass 2005 does [13.5nA,14.5nA]. Maybe it is to avoid too many spikes without inputs...\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['initial_volt'] = (13.5,14.9) \n",
    "                              # Range of the uniform distribution used to generate the initial membrane voltage.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [13.5mV,15.0mV]\n",
    "                              # Joshi/Maass 2005 does [13.5mV,14.9mV], this way it's impossible to generate spikes (14.9mV is below threshold)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['random_reset'] = (13.8,14.5)    # Range of the uniform distribution used to generate the constant reset voltages.\n",
    "    simulation_2DoF_Arm.func_globals['Net_shape'] = (20,5,6)          # Defines the shape of the network (liquid state machine)\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_input_layers'] = 6    # It means I will have 6 input layers...\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_neurons_inputs'] = 50 # ...with 50 neurons each.\n",
    "    simulation_2DoF_Arm.func_globals['lbd_value'] = 1.2               # lbd controls the connection probability\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['noisy_currents'] = True # Injects noisy currents the neurons\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input_weights'] = False # Injects noise into the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input'] = True # Ignores the random seed for the input noise\n",
    "    simulation_2DoF_Arm.func_globals['input_noise_level'] = 10E-5 # 1/SNR\n",
    "\n",
    "    # STP parameters\n",
    "    simulation_2DoF_Arm.func_globals['AII']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>2.8\n",
    "    simulation_2DoF_Arm.func_globals['AIE']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>3.0\n",
    "    simulation_2DoF_Arm.func_globals['AEI']=150.0 # Joshi/Maass2006=>150.0 / Maass2002=>1.6\n",
    "    simulation_2DoF_Arm.func_globals['AEE']=70.0  # Joshi/Maass2006=>70.0 / Maass2002=>1.2    \n",
    "    \n",
    "    maximum_input_index_noise=1 # This is the amount of noise inserted in the feedback loop as variations in the input neurons that spike\n",
    "    my_noisy_input=1            # Controls if the input spike noise is going to use my_seed(0) or just be randomly noisy(1)\n",
    "                                # This way I can have a simulation where all the other noise sources are kept the same, but this.        \n",
    "    \n",
    "    if my_noisy_input:\n",
    "        input_rs = numpy.random.RandomState()\n",
    "    else:\n",
    "        input_rs = numpy.random.RandomState(simulation_2DoF_Arm.func_globals['my_seed'])\n",
    "    \n",
    "    step_size = 2\n",
    "    \n",
    "    # Starting the simulation object!\n",
    "    s = step_by_step_brian_sim(simulation_2DoF_Arm, step_size)\n",
    "\n",
    "    init_t = time.time()\n",
    "    simulated_values = []\n",
    "#####DEBUG\n",
    "    inputs = []\n",
    "#####END_DEBUG    \n",
    "    \n",
    "    # Runs the simulation going step-by-step into the whole input_spikes list!        \n",
    "    for ipos,iteta,itau in zip(input_spikes_pos,input_spikes_teta,input_spikes_tau):\n",
    "        if maximum_input_index_noise>0:\n",
    "            simulated_values.append(s.run_step(ipos+iteta+[itau[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),itau[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)])) # Sends the spikes and runs to the second yield \n",
    "        else:\n",
    "            simulated_values.append(s.run_step(ipos+iteta+itau)) # Sends the spikes and runs to the second yield \n",
    "#####DEBUG\n",
    "            inputs.append(ipos+iteta+itau)\n",
    "#####END_DEBUG                \n",
    "\n",
    "    # This last step is to make sure the result from the last input was processed and then the simulation\n",
    "    # can reach the final time step:\n",
    "    simulated_values.append(s.run_step([])) # Sends the spikes and runs to the second yield \n",
    "    \n",
    "    print \"Total simulation time (seconds):\",time.time()-init_t\n",
    "\n",
    "    if simulation_2DoF_Arm.func_globals['STP_OFF']:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\".gzpickle\"\n",
    "    else:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\"_STP_ON.gzpickle\"\n",
    "\n",
    "    slf.save_to_file_gz(simulated_values,\"./\"+base_dir+\"/\"+sim_set+filename)\n",
    "    \n",
    "#####DEBUG\n",
    "    slf.save_to_file_gz(inputs,\"./\"+base_dir+\"/\"+sim_set+\"/States_Torques_movement\"+str(pos_i)+\"_INPUTS_LSM_\"+str(run_i)+\".gzpickle\")\n",
    "#####END_DEBUG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "base_dir = \"2DofArm_simulation_data\"\n",
    "sim_set = \"set_A\"\n",
    "\n",
    "\n",
    "total_trials = 20 # total number of trials\n",
    "\n",
    "trajectories = [[[0.75,0.25],[0.0,0.5]], [[0.25,0.60],[-0.25,0.60]], [[-0.10,0.75],[-0.10,0.25]],[[-0.75,0.50],[-0.40,0.00]]]\n",
    "\n",
    "number_of_trajectories = len(trajectories)\n",
    "\n",
    "for tji, positions in [(tji,positions) for tji,positions in zip(range(1,len(trajectories)+1),trajectories)]:\n",
    "\n",
    "    xstart,ystart = positions[0]\n",
    "    xdest,ydest = positions[1]\n",
    "\n",
    "    # The default value to delta_t is 200, so when simulating the same as Joshi/Maass 2006 it's \n",
    "    # not necessary to especify the delta_t.\n",
    "    input_spikes = generates_input_spikes_lsm(delta_t=200)\n",
    "\n",
    "    print \"Trajectory Number:\",tji\n",
    "    %time results = run_simulations_2.map([(tji,i,input_spikes,(sim_set,base_dir)) for i in range(1,total_trials+1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set B - WITHOUT STP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@dview.parallel(block=True)\n",
    "def run_simulations_2(sim_num): \n",
    "    import numpy\n",
    "    import time\n",
    "    import sys\n",
    "    import save_load_file as slf\n",
    "        \n",
    "    import step_by_step_brian\n",
    "    reload(sys.modules['step_by_step_brian']) # Makes sure the interpreter is going to reload the module\n",
    "    step_by_step_brian_sim = step_by_step_brian.step_by_step_brian_sim\n",
    "    \n",
    "    import simulation_2DoF_Arm_md  \n",
    "    reload(sys.modules['simulation_2DoF_Arm_md']) # Makes sure the interpreter is going to reload the module\n",
    "    simulation_2DoF_Arm = simulation_2DoF_Arm_md.simulation_2DoF_Arm    \n",
    "    \n",
    "    pos_i, run_i, input_spikes, folder_names = sim_num\n",
    "\n",
    "    sim_set,base_dir = folder_names\n",
    "    \n",
    "    input_spikes_pos,input_spikes_teta,input_spikes_tau = input_spikes\n",
    "    \n",
    "    # Initializing the simulation...\n",
    "    \n",
    "    # Simulation's variables:\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['my_seed']=93200      # This seed is important to always repeat the liquid's structure\n",
    "    simulation_2DoF_Arm.func_globals['input_gain']=70.0    # Base value used within the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy']=1            # Controls if the liquid is going to use the noisy currents/input weights\n",
    "    simulation_2DoF_Arm.func_globals['gaussian_pop']=1     # Controls which type of input weight configuration is used (1=>gaussian distributed)\n",
    "    simulation_2DoF_Arm.func_globals['w_SD']=3.0           # The width of the gaussian used above.\n",
    "    simulation_2DoF_Arm.func_globals['stratification']=1   # Type of connections: stratified or 30% random\n",
    "    simulation_2DoF_Arm.func_globals['STP_OFF']=True       # Uses or not STP (False means it uses STP)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['i_noiseA'] = 1.0     # amplitude (dimensionless) of the random currents (updated every time step).\n",
    "    simulation_2DoF_Arm.func_globals['offset_noise'] = (13.5,14.5) \n",
    "                              # Range of the uniform distribution used to generate the constant offset current.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [14.975nA,15.025nA]\n",
    "                              # Joshi/Maass 2005 does [13.5nA,14.5nA]. Maybe it is to avoid too many spikes without inputs...\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['initial_volt'] = (13.5,14.9) \n",
    "                              # Range of the uniform distribution used to generate the initial membrane voltage.\n",
    "                              # Maass 2002, should be drawn from a uniform distr [13.5mV,15.0mV]\n",
    "                              # Joshi/Maass 2005 does [13.5mV,14.9mV], this way it's impossible to generate spikes (14.9mV is below threshold)\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['random_reset'] = (13.8,14.5)    # Range of the uniform distribution used to generate the constant reset voltages.\n",
    "    simulation_2DoF_Arm.func_globals['Net_shape'] = (20,5,6)          # Defines the shape of the network (liquid state machine)\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_input_layers'] = 6    # It means I will have 6 input layers...\n",
    "    simulation_2DoF_Arm.func_globals['Number_of_neurons_inputs'] = 50 # ...with 50 neurons each.\n",
    "    simulation_2DoF_Arm.func_globals['lbd_value'] = 1.2               # lbd controls the connection probability\n",
    "\n",
    "    simulation_2DoF_Arm.func_globals['noisy_currents'] = True # Injects noisy currents the neurons\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input_weights'] = False # Injects noise into the input weights\n",
    "    simulation_2DoF_Arm.func_globals['noisy_input'] = True # Ignores the random seed for the input noise\n",
    "    simulation_2DoF_Arm.func_globals['input_noise_level'] = 0 # 1/SNR\n",
    "    \n",
    "    # STP parameters\n",
    "    simulation_2DoF_Arm.func_globals['AII']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>2.8\n",
    "    simulation_2DoF_Arm.func_globals['AIE']=47.0  # Joshi/Maass2006=>47.0 / Maass2002=>3.0\n",
    "    simulation_2DoF_Arm.func_globals['AEI']=150.0 # Joshi/Maass2006=>150.0 / Maass2002=>1.6\n",
    "    simulation_2DoF_Arm.func_globals['AEE']=70.0  # Joshi/Maass2006=>70.0 / Maass2002=>1.2\n",
    "    \n",
    "    \n",
    "    maximum_input_index_noise=1 # This is the amount of noise inserted in the feedback loop as variations in the input neurons that spike\n",
    "    my_noisy_input=1            # Controls if the input spike noise is going to use my_seed(0) or just be randomly noisy(1)\n",
    "                                # This way I can have a simulation where all the other noise sources are kept the same, but this.    \n",
    "    if my_noisy_input:\n",
    "        input_rs = numpy.random.RandomState()\n",
    "    else:\n",
    "        input_rs = numpy.random.RandomState(simulation_2DoF_Arm.func_globals['my_seed'])\n",
    "    \n",
    "    step_size = 2\n",
    "    \n",
    "    # Starting the simulation object!\n",
    "    s = step_by_step_brian_sim(simulation_2DoF_Arm, step_size)\n",
    "\n",
    "    init_t = time.time()\n",
    "    simulated_values = []\n",
    "#####DEBUG\n",
    "    inputs = []\n",
    "#####END_DEBUG    \n",
    "\n",
    "    # Runs the simulation going step-by-step into the whole input_spikes list!        \n",
    "    for ipos,iteta,itau in zip(input_spikes_pos,input_spikes_teta,input_spikes_tau):\n",
    "        if maximum_input_index_noise>0:\n",
    "            simulated_values.append(s.run_step(ipos+iteta+[itau[0]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1),itau[1]+input_rs.randint(-maximum_input_index_noise,maximum_input_index_noise+1)])) # Sends the spikes and runs to the second yield \n",
    "        else:\n",
    "            simulated_values.append(s.run_step(ipos+iteta+itau)) # Sends the spikes and runs to the second yield \n",
    "#####DEBUG\n",
    "            inputs.append(ipos+iteta+itau)\n",
    "#####END_DEBUG                \n",
    "\n",
    "    # This last step is to make sure the result from the last input was processed and then the simulation\n",
    "    # can reach the final time step:\n",
    "    simulated_values.append(s.run_step([])) # Sends the spikes and runs to the second yield \n",
    "    \n",
    "    print \"Total simulation time (seconds):\",time.time()-init_t\n",
    "\n",
    "    if simulation_2DoF_Arm.func_globals['STP_OFF']:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\".gzpickle\"\n",
    "    else:\n",
    "        filename = \"/States_Torques_movement\"+str(pos_i)+\"_LSM_\"+str(run_i)+\"_STP_ON.gzpickle\"\n",
    "\n",
    "    slf.save_to_file_gz(simulated_values,\"./\"+base_dir+\"/\"+sim_set+filename)\n",
    "    \n",
    "#####DEBUG\n",
    "    slf.save_to_file_gz(inputs,\"./\"+base_dir+\"/\"+sim_set+\"/States_Torques_movement\"+str(pos_i)+\"_INPUTS_LSM_\"+str(run_i)+\".gzpickle\")\n",
    "#####END_DEBUG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "base_dir = \"2DofArm_simulation_data\"\n",
    "sim_set = \"set_B\"\n",
    "\n",
    "\n",
    "total_trials = 20 # total number of trials\n",
    "\n",
    "trajectories = [[[0.75,0.25],[0.0,0.5]], [[0.25,0.60],[-0.25,0.60]], [[-0.10,0.75],[-0.10,0.25]],[[-0.75,0.50],[-0.40,0.00]]]\n",
    "\n",
    "number_of_trajectories = len(trajectories)\n",
    "\n",
    "for tji, positions in [(tji,positions) for tji,positions in zip(range(1,len(trajectories)+1),trajectories)]:\n",
    "\n",
    "    xstart,ystart = positions[0]\n",
    "    xdest,ydest = positions[1]\n",
    "\n",
    "    # The default value to delta_t is 200, so when simulating the same as Joshi/Maass 2006 it's \n",
    "    # not necessary to especify the delta_t.\n",
    "    input_spikes = generates_input_spikes_lsm(delta_t=200)\n",
    "\n",
    "    print \"Trajectory Number:\",tji\n",
    "    %time results = run_simulations_2.map([(tji,i,input_spikes,(sim_set,base_dir)) for i in range(1,total_trials+1)])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}