M2E_statistical_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bc96f39-b94f-40d5-be5c-a0d9ebdf8180",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "import json\n",
    "\n",
    "sys.path.append('./figures/MAM2EBRAINS')\n",
    "\n",
    "from multiarea_model import MultiAreaModel\n",
    "from multiarea_model import Analysis\n",
    "from M2E_compute_pop_rates import compute_pop_rates\n",
    "\n",
    "base_path = os.path.abspath(\".\")\n",
    "data_path = os.path.abspath(\"simulations\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c907631c-4611-469a-941c-ddb31ef753c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the threshold of matching rate of mean firing rate of all populations\n",
    "# if smaller than this value, raise error so that the test won't pass and something is wrong and needs to be checked manually\n",
    "match_rate_threshold = 0.95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97fbab67-7282-4800-a6d4-78bb3a36640a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Down-scaled model.\n",
    "Neurons and indegrees are both scaled down to 10 %.\n",
    "Can usually be simulated on a local machine.\n",
    "\n",
    "Warning: This will not yield reasonable dynamical results from the\n",
    "network and is only meant to demonstrate the simulation workflow.\n",
    "\"\"\"\n",
    "d = {}\n",
    "conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',\n",
    "               'cc_weights_factor': 1.9, # run model in Ground State\n",
    "               'cc_weights_I_factor': 2.0}\n",
    "network_params = {'N_scaling': 0.006,\n",
    "                  'K_scaling': 0.006,\n",
    "                  'fullscale_rates': os.path.join(base_path, 'tests/fullscale_rates.json')}\n",
    "\n",
    "sim_params = {'t_sim': 2000.,\n",
    "              'num_processes': 1,\n",
    "              'local_num_threads': 1}\n",
    "\n",
    "M = MultiAreaModel(network_params, simulation=True,\n",
    "                   sim_spec=sim_params,\n",
    "                   theory=True)\n",
    "\n",
    "M.simulation.simulate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6730fe7b-e8bf-4aff-9497-bc614f44febd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an instance of Analysis to load data\n",
    "A = Analysis(M, M.simulation, data_list=['spikes'], load_areas=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c907303-d2ad-4f17-83df-8b0e21749a1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create pop_rates, load stationary firing rates, and calculate mean firing rate for all populations\n",
    "label = M.simulation.label\n",
    "complete_area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd',\n",
    "                      'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd',\n",
    "                      'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp',\n",
    "                      'STPa', '46', 'AITd', 'TH']\n",
    "\n",
    "area_list = complete_area_list\n",
    "\n",
    "compute_pop_rates(M, data_path, label) # Compute pop_rates\n",
    "\n",
    "fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json') # Load stationary firing rates\n",
    "with open(fn, 'r') as f:\n",
    "    pop_rates = json.load(f)\n",
    "\n",
    "# Process spike data and calculate mean firing rate for all populations, save them in a dict\n",
    "rates = np.zeros((len(area_list), 8))\n",
    "for i, area in enumerate(area_list):\n",
    "    for j, pop in enumerate(M.structure[area][::-1]):\n",
    "        rate = pop_rates[area][pop]\n",
    "        if rate == 0.0:\n",
    "            rate = 1e-5\n",
    "        if area == 'TH' and j > 3:  # To account for missing layer 4 in TH\n",
    "            rates[i][j + 2] = rate\n",
    "        else:\n",
    "            rates[i][j] = rate\n",
    "\n",
    "rates = np.transpose(rates)\n",
    "rates = np.array(rates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34d46019-38f2-4b5d-af09-575906b1d4d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the calculated mean firing rate with the reference values saved in M2E_reference_MFRs.json and calculate the percentage of matching\n",
    "# Load the array\n",
    "pop_rates_ref = np.load('M2E_pop_rates_ref.npy')\n",
    "\n",
    "# Count the percentage of cells that values match (the same)\n",
    "comparison_result = pop_rates_ref == rates\n",
    "# print(comparison_result)\n",
    "same_cells = np.sum(comparison_result)\n",
    "# print(same_cells)\n",
    "total_cells = np.prod(pop_rates_ref.shape)\n",
    "# print(total_cells)\n",
    "match_rate = same_cells / total_cells\n",
    "# print(match_rate)\n",
    "\n",
    "# Judge if the matching rate is lower than threshold set and raise error if yes\n",
    "if match_rate < match_rate_threshold:\n",
    "    raise TestError(\"Mismatch of mean firing rates over populations between the latest simulation and saved reference values, there may be problems of recent updates of dependencies, please take a check!\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EBRAINS-23.09",
   "language": "python",
   "name": "ebrains-23.09"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}