Theta_Critical_Gradient_Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-Spike Tempotron Learning Algorithm for Computing the Theta Critical Gradient "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Execute all cells containing functions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import optimize\n",
    "import multiprocessing\n",
    "rtol = 1e-5 # rtol from numpy.isclose\n",
    "from zipfile import BadZipFile\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gen_omega(n, omega_coefficient):\n",
    "    omega = np.random.random(n)*omega_coefficient\n",
    "    return omega"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def desired_fea(n_fea, desired_n, cond = 0, multi = 0):\n",
    "    fea_num = np.zeros(n_fea)\n",
    "    for i in range(desired_n):\n",
    "        if cond:\n",
    "            if multi:\n",
    "                fea_num[i] = multi\n",
    "            else:\n",
    "                fea_num[i] = i + 1\n",
    "        else:\n",
    "            fea_num[i] = 1\n",
    "    return fea_num"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def STS(unweighted_input, theta, omega):\n",
    "    spike_time = []\n",
    "    unresetted_V = []\n",
    "    datalen = unweighted_input.shape[1]\n",
    "    V = (omega[:,np.newaxis] * unweighted_input).sum(axis=0)\n",
    "    done = False\n",
    "    while not done:\n",
    "        spike_idx = find_first_spike(V, theta)\n",
    "        if spike_idx == -1:\n",
    "            done = True\n",
    "        else:\n",
    "            spike_time.append(spike_idx)\n",
    "            unresetted_V.append(V[spike_idx])\n",
    "            mem_len = min(ref_memory_len, datalen - spike_idx)\n",
    "            V[spike_idx:spike_idx+mem_len] -= theta * ref_kernel[:mem_len]\n",
    "    return spike_time, unresetted_V, max(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bisect(unweighted_input, theta_spikes, theta_range, smallest_range, omega):\n",
    "    theta_list = [theta_range[0], sum(theta_range)/2, theta_range[1]]\n",
    "    high = STS(unweighted_input, theta_range[0], omega)\n",
    "    medium = STS(unweighted_input, sum(theta_range)/2, omega)\n",
    "    low = STS(unweighted_input, theta_range[1], omega)\n",
    "    spike_list = [len(high[0]), len(medium[0]), len(low[0])]\n",
    "    if theta_list[2] - theta_list[0] >= smallest_range:\n",
    "        if theta_spikes <= spike_list[0] and theta_spikes > spike_list[1]:\n",
    "            theta_range = [theta_list[0], theta_list[1]]\n",
    "        elif theta_spikes <= spike_list[1] and theta_spikes > spike_list[2]:\n",
    "            theta_range = [theta_list[1], theta_list[2]]\n",
    "        mid_theta, mid_spike = bisect(unweighted_input, theta_spikes, theta_range, smallest_range, omega)\n",
    "        return sorted(np.append(theta_list, mid_theta[1:-1]), reverse=True), sorted(np.append(spike_list, mid_spike[1:-1])) \n",
    "    else:\n",
    "        if max(spike_list) - min(spike_list) == 1:\n",
    "            return theta_list, spike_list\n",
    "        else:\n",
    "            new_smallest_range = smallest_range/10\n",
    "            theta_range = [theta_list[0], theta_list[2]]\n",
    "            mid_theta, mid_spike = bisect(unweighted_input, theta_spikes, theta_range, new_smallest_range, omega)\n",
    "            return sorted(np.append(theta_list, mid_theta[1:-1]), reverse=True), sorted(np.append(spike_list, mid_spike[1:-1])) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def root_fun(theta_range, unweighted_input, omega):\n",
    "    V_max = STS(unweighted_input, theta_range[0], omega)[2]\n",
    "    return [theta_range[0] - V_max, V_max - theta_range[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def theta_critical(unweighted_input, num_spikes, desired_spikes, theta_range, smallest_range, omega):\n",
    "    if num_spikes > desired_spikes:\n",
    "        theta_spikes = num_spikes\n",
    "    elif num_spikes < desired_spikes:\n",
    "        theta_spikes = num_spikes + 1\n",
    "    theta_list, spike_list = bisect(unweighted_input, theta_spikes, theta_range, smallest_range, omega)\n",
    "    theta_critical_range = optimize.root(root_fun, [theta_list[spike_list.index(theta_spikes)-1], theta_list[spike_list.index(theta_spikes)]], args=(unweighted_input, omega)).x \n",
    "    spike_time, unresetted_V, max_V = STS(unweighted_input, theta_critical_range[0] * (1 - rtol), omega)\n",
    "    t_star_list = spike_time[:unresetted_V.index(min(unresetted_V))+1]\n",
    "    return theta_critical_range[1], t_star_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Optimized voltage calculation code by Alex ###\n",
    "\n",
    "def get_memory_len(kernel_array, ratio):\n",
    "    \"\"\"\n",
    "    Return the number of time bins until kernel_array has decreased by the factor 'ratio'\n",
    "    'kernel_array' may initially rise, but must be monotonically decreasing from \n",
    "    the maximum.\n",
    "    \"\"\"\n",
    "    arr = (kernel_array - ratio*kernel_array.max())[::-1]\n",
    "        # The point where this array reaches zero is the desired memory time\n",
    "        # We flip the order with [::-1] because np.searchsorted expects increasing order\n",
    "    memory_len = len(kernel_array) - np.searchsorted(arr, 0)\n",
    "    return memory_len\n",
    "\n",
    "def find_first_spike(V, threshold):\n",
    "    # Based on the equivalent code to py_find_1st (https://github.com/roebel/py_find_1st)\n",
    "    ind = np.flatnonzero(V >= threshold)\n",
    "    if len(ind):\n",
    "        return ind[0]\n",
    "    else:\n",
    "        return -1\n",
    "\n",
    "def calculate_Vt(unweighted_input, theta, omega):\n",
    "    spike_time = []\n",
    "    datalen = unweighted_input.shape[1]\n",
    "    V = (omega[:,np.newaxis] * unweighted_input).sum(axis=0)\n",
    "    done = False\n",
    "    while not done:\n",
    "        spike_idx = find_first_spike(V, theta)\n",
    "        if spike_idx == -1:\n",
    "            done = True\n",
    "        else:\n",
    "            spike_time.append(spike_idx)\n",
    "            mem_len = min(ref_memory_len, datalen - spike_idx)\n",
    "            V[spike_idx:spike_idx+mem_len] -= theta * ref_kernel[:mem_len]\n",
    "    return V, spike_time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def C_tx(spike_time):    #eqn 29\n",
    "    C = np.zeros(len(spike_time))\n",
    "    count = 0\n",
    "    for spike in spike_time:\n",
    "        C[count] = 1\n",
    "        for idx in range(count):\n",
    "            if spike_time[idx] > spike - ref_memory_len:\n",
    "                C[count] += ref_kernel[spike - spike_time[idx]] \n",
    "        count += 1\n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dVtx_dwi(unweighted_input, spike_time, C):    #eqn 30\n",
    "    kernel_t = unweighted_input[:, spike_time]\n",
    "    return kernel_t/C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dVtx_dtsk(tot_input, C, spike_time, x, k):    #eqn 31\n",
    "    return (-tot_input[spike_time[x]]/C[x]**2)*(np.exp(-(spike_time[x]-spike_time[k])/tau_mem)/tau_mem)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dif_kernel(length, tau_mem, tau_syn, time_ij):\n",
    "    time = np.arange(0., length, 1.) #ms\n",
    "    kernel = np.zeros(length)\n",
    "    eta = tau_mem/tau_syn\n",
    "    V_norm = eta**(eta/(eta-1))/(eta-1)\n",
    "    for count in range(length):\n",
    "        kernel[count] = V_norm*(-(np.exp(-(time[count]-time_ij)/tau_mem))/tau_mem+(np.exp(-(time[count]-time_ij)/tau_syn))/tau_syn)  \n",
    "    return kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dV0t_dt(data, diff_kernel, omega):    #eqn 32 first part\n",
    "    diff_kernel_array = np.zeros((n, data.shape[1]))\n",
    "    for neuron, ith_bin in zip(*np.where(data)):\n",
    "        kernel_len = min(len(diff_kernel), data.shape[1] - ith_bin)\n",
    "        diff_kernel_array[neuron, ith_bin:ith_bin+kernel_len] += np.multiply(diff_kernel, omega[neuron])[:kernel_len]\n",
    "    return diff_kernel_array.sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def t_derivative(sum_diff_kernel, tot_input, C, spike_time, x):    #eqn 32\n",
    "    return sum_diff_kernel[spike_time[x]]/C[x] + tot_input[spike_time[x]]*sum(np.exp(-(spike_time[x]-spike_time[j])/tau_mem) for j in range(x))/(tau_mem*C[x]**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def A_cache(n_spikes, sum_diff_kernel, tot_input, C, spike_time):\n",
    "    A_cache = np.zeros(n_spikes)\n",
    "    for k in range(n_spikes):\n",
    "        A_cache[k] = fn_A(sum_diff_kernel, tot_input, A_cache, C, spike_time, k)\n",
    "    return A_cache"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fn_A(sum_diff_kernel, tot_input, A, C, spike_time, k):    #eqn 23\n",
    "    return 1 - sum(A[j]/t_derivative(sum_diff_kernel, tot_input, C, spike_time, j)*dVtx_dtsk(tot_input, C, spike_time, k, j) for j in range(k))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def B_cache(n_spikes, sum_diff_kernel, tot_input, dVt_dwi, C, spike_time):\n",
    "    B_cache = np.zeros(n_spikes)\n",
    "    for k in range(n_spikes):\n",
    "        B_cache[k] = fn_B(sum_diff_kernel, tot_input, B_cache, C, spike_time, dVt_dwi, k)\n",
    "    return B_cache"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fn_B(sum_diff_kernel, tot_input, B, C, spike_time, dVt_dwi, k):    #eqn 24\n",
    "    if len(dVt_dwi) != 0:\n",
    "        return -dVt_dwi[k] - sum(B[j]/t_derivative(sum_diff_kernel, tot_input, C, spike_time, j)*dVtx_dtsk(tot_input, C, spike_time, k, j) for j in range(k))\n",
    "    else:\n",
    "        return 0 - sum(B[j]/t_derivative(sum_diff_kernel, tot_input, C, spike_time, j)*dVtx_dtsk(tot_input, C, spike_time, k, j) for j in range(k))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def theta_grad_i(A_cache, B_cache, star):   #Eqn 27\n",
    "    return -B_cache[star]/A_cache[star]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def theta_critical_grad(n, sum_diff_kernel, tot_input, dVt_dw, A, C, spike_time, n_spikes, star):\n",
    "    grad = np.zeros(n)\n",
    "    for count in range(n):\n",
    "        dVt_dwi = dVt_dw[count]\n",
    "        B = B_cache(n_spikes, sum_diff_kernel, tot_input, dVt_dwi, C, spike_time)\n",
    "        grad[count] += theta_grad_i(A, B, star)\n",
    "    return grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def update_w(omega, learning_rate, desired_spikes, num_spikes, sum_diff_kernel, tot_input, dVt_dw, A, C, spike_time, n_spikes, star):\n",
    "    if num_spikes > desired_spikes:\n",
    "        return omega - learning_rate*theta_critical_grad(n, sum_diff_kernel, tot_input, dVt_dw, A, C, spike_time, n_spikes, star)\n",
    "    elif num_spikes < desired_spikes:\n",
    "        return omega + learning_rate*theta_critical_grad(n, sum_diff_kernel, tot_input, dVt_dw, A, C, spike_time, n_spikes, star)\n",
    "    else:\n",
    "        return omega"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def multispike_training(current_omega, n_cycles, data_sets, fea_idx, learning_rate):\n",
    "    cur_omega_list = []\n",
    "    for cycle in range(n_cycles):\n",
    "        for trial in range(cycle*100, cycle*100+100):\n",
    "            clear_output(wait=True)\n",
    "            print(\"Cycle {}, trial {}.\".format(cycle, trial))\n",
    "            input_data = None\n",
    "            while input_data is None:\n",
    "                idx = int(np.random.random()*data_sets)\n",
    "                try:\n",
    "                    input_data = np.load(\"data/data_\"+str(idx)+\".npz\")\n",
    "                    data, presyn_input, markers, n_fea_occur, fea_time, fea_order = input_data['arr_0'], input_data['arr_1'], input_data['arr_2'], input_data['arr_3'], input_data['arr_4'], input_data['arr_5']\n",
    "                except (FileNotFoundError, BadZipFile, KeyError):\n",
    "                    pass\n",
    "            desired_spikes = sum(n_fea_occur * fea_idx)\n",
    "            V_t, spike_time = calculate_Vt(presyn_input, theta, current_omega)\n",
    "            tot_input = (current_omega[:,np.newaxis] * presyn_input).sum(axis=0)\n",
    "            num_spikes = len(spike_time)\n",
    "            if num_spikes != desired_spikes:\n",
    "                theta_star, tx_list = theta_critical(presyn_input, num_spikes, desired_spikes, theta_range, 0.00001, current_omega)\n",
    "                C = C_tx(tx_list)\n",
    "                dVt_dw = dVtx_dwi(presyn_input, tx_list, C)\n",
    "                sum_diff_kernel = dV0t_dt(data, diff_kernel, current_omega)\n",
    "                A = A_cache(len(tx_list), sum_diff_kernel, tot_input, C, tx_list)\n",
    "                t_star = len(tx_list)-1\n",
    "                current_omega = update_w(current_omega, learning_rate, desired_spikes, num_spikes, sum_diff_kernel, tot_input, dVt_dw, A, C, tx_list, len(tx_list), t_star)\n",
    "        cur_omega_list.append(current_omega)\n",
    "    return cur_omega_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Training Procedures for Multi-Spike Tempotron"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Pre-set Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Points to consider before executing the cell below:\n",
    "\n",
    "a. Omega_coefficient of 0.022 was chosen based on the average firing rate of approximately 5 Hz for the postsynaptic neuron (multi-spike tempotron) with initial synaptic efficacies (omega).\n",
    "\n",
    "b. Theta_range was determined based on the omega_coefficient, it is recommended to reset the theta_range if the omega_coefficient is changed to avoid potential \"out-of-range\" errors.\n",
    "\n",
    "c. The effective dif_kernel_len was selected based on the membrane integration time constant (tau_mem) and synaptic decay time constant (tau_syn). It is highly recommended to determine the effective dif_kernel_len again if tau_mem and tau_syn are changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 500 # number of input neurons\n",
    "omega_coefficient = 0.022\n",
    "theta = 1 # threshold\n",
    "theta_range = [0.5, 1.5] # initial threshold range for bisection\n",
    "tau_mem = 20 # ms\n",
    "tau_syn = 5 # ms\n",
    "time_ij = 0 # ms\n",
    "dif_kernel_len = 150 # length of differentiated kernel (postsynaptic potential --> filter), negligible tail is truncated\n",
    "\n",
    "# Generate the initial omega for each input neuron\n",
    "np.random.seed(100000)\n",
    "omega = gen_omega(n, omega_coefficient)\n",
    "\n",
    "# Create the refractory kernel, and then remove the negligible tail\n",
    "ref_kernel = np.exp(- np.arange(1000) / tau_mem)\n",
    "ref_memory_len = get_memory_len(ref_kernel, ratio=0.001)\n",
    "# synaptic memory length\n",
    "ref_kernel = ref_kernel[:ref_memory_len]\n",
    "\n",
    "# Generate the differentiated kernel\n",
    "diff_kernel = dif_kernel(dif_kernel_len, tau_mem, tau_syn, time_ij)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Training "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multispike_training function takes 5 arguments, which are current_omega, n_cycles, data_sets, fea_idx, and learning_rate, as inputs and return a list of updated omegas. The updated omega list stores the updated omega after every cycle of training (100 trials each). \n",
    "\n",
    "a. current_omega: intial omega\n",
    "\n",
    "b. n_cycles: number of training cycles required\n",
    "\n",
    "c. data_sets: available number of precomputed data sets in the pool\n",
    "\n",
    "d. fea_idx: indices of features that the multispike tempotron has to identify\n",
    "\n",
    "e. learning_rate: the size of update after each error trial\n",
    "\n",
    "\n",
    "Note: The learning_rate of 0.00001 was chosen according to the original publication (Robert Gütig, 2016)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.1 Identification of one feature with one spike "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The desired_fea function is used to generate a list of desired feature numbers.\n",
    "\n",
    "To generate one spike to identify only one feature, set desired_n to 1 and cond to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cycles = 1000\n",
    "data_sets = 19999\n",
    "learning_rate = 0.00001\n",
    "n_fea = 10 # Total number of features and distractors\n",
    "desired_n = 1 # Number of desired features\n",
    "\n",
    "# To generate a list of desired feature indices\n",
    "fea_idx = desired_fea(n_fea, desired_n, cond = 0, multi = 0) \n",
    "\n",
    "np.random.seed(int(1e9))\n",
    "latest_omega1F1S_list = multispike_training(omega, n_cycles, data_sets, fea_idx, learning_rate)\n",
    "np.save(\"theta_omega1F1S_list\", latest_omega1F1S_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.2 Identification of one feature with multiple spikes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate multiple spikes to identify one feature, set desired_n to 1, cond to 1, and multi to the desired number of spikes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cycles = 1000\n",
    "data_sets = 19999\n",
    "learning_rate = 0.00001\n",
    "n_fea = 10 # Total number of features and distractors\n",
    "desired_n = 1 # Number of desired features\n",
    "\n",
    "# To generate a list of desired feature indices\n",
    "fea_idx = desired_fea(n_fea, desired_n, cond = 1, multi = 5) \n",
    "\n",
    "np.random.seed(int(1e9))\n",
    "latest_omega1F5S_list = multispike_training(omega, n_cycles, data_sets, fea_idx, learning_rate)\n",
    "np.save(\"theta_omega1F5S_list\", latest_omega1F5S_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.3 Identification of multiple features with one spike per feature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate one spike to identify multiple features, set desired_n to the number of desired features and cond to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cycles = 1000\n",
    "data_sets = 19999\n",
    "learning_rate = 0.00001\n",
    "n_fea = 10 # Total number of features and distractors\n",
    "desired_n = 5 # Number of desired features\n",
    "\n",
    "# To generate a list of desired feature indices\n",
    "fea_idx = desired_fea(n_fea, desired_n, cond = 0, multi = 0) \n",
    "\n",
    "np.random.seed(int(1e9))\n",
    "latest_omega5F1S_list = multispike_training(omega, n_cycles, data_sets, fea_idx, learning_rate)\n",
    "np.save(\"theta_omega5F1S_list\", latest_omega5F1S_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.4 Identification of multiple features with different number of spikes for different features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate different number of spikes to identify different features, set desired_n to the number of desired features, cond to 1, and multi to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cycles = 1000\n",
    "data_sets = 19999\n",
    "learning_rate = 0.00001\n",
    "n_fea = 10 # Total number of features and distractors\n",
    "desired_n = 5 # Number of desired features\n",
    "\n",
    "# To generate a list of desired feature indices\n",
    "fea_idx = desired_fea(n_fea, desired_n, cond = 1, multi = 0) \n",
    "\n",
    "np.random.seed(int(1e9))\n",
    "latest_omega5FmS_list = multispike_training(omega, n_cycles, data_sets, fea_idx, learning_rate)\n",
    "np.save(\"theta_omega5FmS_list\", latest_omega5FmS_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.5 Identification of multiple features with a fixed number of spikes (> 1) for each feature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate a fixed number of spikes that is greater than 1 for multiple features, set desired_n to the number of desired features, cond to 1, and multi to desired number of spikes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cycles = 1000\n",
    "data_sets = 19999\n",
    "learning_rate = 0.00001\n",
    "n_fea = 10 # Total number of features and distractors\n",
    "desired_n = 5 # Number of desired features\n",
    "\n",
    "# To generate a list of desired feature indices\n",
    "fea_idx = desired_fea(n_fea, desired_n, cond = 1, multi = 5) \n",
    "\n",
    "np.random.seed(int(1e9))\n",
    "latest_omega5F5S_list = multispike_training(omega, n_cycles, data_sets, fea_idx, learning_rate)\n",
    "np.save(\"theta_omega5F5S_list\", latest_omega5F5S_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Multiple runs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cycles = 1000\n",
    "data_sets = 19999\n",
    "learning_rate = 0.00001\n",
    "n_fea = 10 # Total number of features and distractors\n",
    "desired_n = 5 # Number of desired features\n",
    "\n",
    "fea_idx_list = [desired_fea(n_fea, desired_n=1, cond=0, multi=0),\n",
    "                desired_fea(n_fea, desired_n=1, cond=1, multi=5),\n",
    "                desired_fea(n_fea, desired_n=5, cond=0, multi=0),\n",
    "                desired_fea(n_fea, desired_n=5, cond=1, multi=0),\n",
    "                desired_fea(n_fea, desired_n=5, cond=1, multi=5)\n",
    "               ]\n",
    "filenames = [\"theta_omega1F1S_list_new1\",\n",
    "             \"theta_omega1F5S_list_new1\",\n",
    "             \"theta_omega5F1S_list_new1\",\n",
    "             \"theta_omega5FmS_list_new1\",\n",
    "             \"theta_omega5F5S_list_new1\"\n",
    "            ]\n",
    "\n",
    "def run_and_save(filename, fea_idx):\n",
    "    res = multispike_training(omega, n_cycles, data_sets, fea_idx, learning_rate)\n",
    "    np.save(filename, res)\n",
    "\n",
    "with multiprocessing.Pool(5) as pool:\n",
    "    pool.starmap(run_and_save, zip(filenames, fea_idx_list))"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}