draft_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d76569e-6b75-4bd1-9efc-2d41bc2fd2d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gglasso\n",
    "import pickle\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "from numpy.linalg import matrix_rank\n",
    "from matplotlib.pyplot import figure\n",
    "from scipy import stats\n",
    "from scipy.linalg import eigh\n",
    "from numpy import genfromtxt\n",
    "\n",
    "from gglasso.solver.admm_solver import ADMM_MGL\n",
    "from gglasso.problem import glasso_problem\n",
    "\n",
    "from gglasso.helper.data_generation import generate_precision_matrix, group_power_network, sample_covariance_matrix\n",
    "from gglasso.helper.basic_linalg import adjacency_matrix\n",
    "from gglasso.helper.data_generation import time_varying_power_network, sample_covariance_matrix\n",
    "from gglasso.helper.experiment_helper import lambda_grid, discovery_rate, error\n",
    "from gglasso.helper.utils import get_K_identity\n",
    "from gglasso.helper.experiment_helper import plot_evolution, plot_deviation, surface_plot, single_heatmap_animation\n",
    "from gglasso.helper.model_selection import aic, ebic, K_single_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22acf91b-9772-49b5-a1af-03cd317c48a3",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86ad6809-d3d3-4736-912d-6b3bd188da14",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "start = 0\n",
    "stop = 10\n",
    "\n",
    "sub_corr = []\n",
    "\n",
    "### on cloud\n",
    "# for i in range(start, stop):\n",
    "#     sub_corr.append(genfromtxt(\"/storage/groups/bds01/datasets/brains/corr_matrices/corr{0}.csv\".format(i), delimiter=','))\n",
    "\n",
    "### on premisis\n",
    "for i in range(start, stop):\n",
    "    sub_corr.append(genfromtxt(\"../data/sub_corr50/sub_corr{0}.csv\".format(i), delimiter=','))\n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "sub_corr = np.array(sub_corr)\n",
    "sub_corr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64256d68-7625-4424-b2b6-6cd28a196550",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "lambda1_range = np.logspace(-0.9, -1.5, 2)\n",
    "\n",
    "N = [sub_corr.shape[1]]*sub_corr.shape[0]\n",
    "\n",
    "\n",
    "est_uniform, est_indv, statistics = K_single_grid(sub_corr, lambda1_range, N, \n",
    "                                                  method = 'eBIC', gamma = 0.3, \n",
    "                                                  latent = False, use_block = True)\n",
    "\n",
    "\n",
    "K = \"test\"\n",
    "\n",
    "if not os.path.exists(\"../data/est_uniform{0}/\".format(K)):\n",
    "    os.mkdir(\"../data/est_uniform{0}/\".format(K))\n",
    "    \n",
    "if not os.path.exists(\"../data/est_individ{0}/\".format(K)):\n",
    "    os.mkdir(\"../data/est_individ{0}/\".format(K))\n",
    "\n",
    "# dump matrices into csv\n",
    "for i in range(start, stop):\n",
    "    np.savetxt(\"../data/est_uniform{0}/est_uniform{1}.csv\".format(K, i), est_uniform[\"Theta\"][i], \n",
    "               delimiter=\",\", header='')\n",
    "    np.savetxt(\"../data/est_individ{0}/est_individ{1}.csv\".format(K, i), est_indv[\"Theta\"][i], \n",
    "               delimiter=\",\", header='')\n",
    "    \n",
    "with open(\"../data/statistics{0}.txt\".format(K), 'w') as f:\n",
    "    print(statistics, file=f)\n",
    "\n",
    "\n",
    "### cloud\n",
    "# os.mkdir(\"/storage/groups/bds01/datasets/brains/est_uniform{0}/\".format(K))\n",
    "# os.mkdir(\"/storage/groups/bds01/datasets/brains/est_individ{0}/\".format(K))\n",
    "\n",
    "# # dump matrices into csv\n",
    "# for i in range(start, stop):\n",
    "#     np.savetxt(\"/storage/groups/bds01/datasets/brains/est_uniform{0}/est_uniform{1}.csv\".format(K, i), est_uniform[\"Theta\"][i], \n",
    "#                delimiter=\",\", header='')\n",
    "#     np.savetxt(\"/storage/groups/bds01/datasets/brains/est_individ{0}/est_individ{1}.csv\".format(K, i), est_indv[\"Theta\"][i], \n",
    "#                delimiter=\",\", header='')\n",
    "    \n",
    "# with open(\"statistics{0}.txt\".format(K), 'w') as f:\n",
    "#     print(statistics, file=f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8728fa3-e458-41d1-9c55-e338a88b59d9",
   "metadata": {
    "tags": []
   },
   "source": [
    "### SGL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1466aba-f9a9-41d5-a20f-e246bc1ad854",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "lambda1_range = np.logspace(-0.9, -1.5, 10)\n",
    "lambda1_range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a61b5da-86f9-4c61-9916-adbee95074d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = [sub_corr.shape[1]]*sub_corr.shape[0]\n",
    "\n",
    "est_uniform, est_indv, statistics = K_single_grid(sub_corr, lambda1_range, N, \n",
    "                                                  method = 'eBIC', gamma = 0.3, \n",
    "                                                  latent = False, use_block = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1886a751-7e81-4269-9398-f93fc291d4b4",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## GLasso problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac391209-aaaf-4eb3-af46-acb95d5fb8be",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "N = sub_corr.shape[1]\n",
    "\n",
    "P = glasso_problem(sub_corr[:10,:], N, reg_params = {'lambda1': lambda1_range}, latent = False, do_scaling = False)\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf33a867-9fa5-4b00-829c-9d99558b1f11",
   "metadata": {},
   "outputs": [],
   "source": [
    "modelselect_params = {'lambda1_range': lambda1_range}\n",
    "gamma = 0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93f901d7-0332-406c-be63-8603203ef8f5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "P.model_selection(modelselect_params = modelselect_params, method = 'eBIC', gamma = gamma)\n",
    "\n",
    "# regularization parameters are set to the best ones found during model selection\n",
    "print(P.reg_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cff86b7f-d6f0-4e69-8d94-88677c2e1b81",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "P.solution.precision_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f83c67bf-99e0-43cf-b04b-5da0d73519d5",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sgl_stats = pd.DataFrame({'SP': list(P.modelselect_stats[\"SP\"]), 'lambda': list(P.modelselect_stats[\"LAMBDA\"])})\n",
    "sgl_stats[sgl_stats[\"lambda\"] == 0.14237976633439411]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b7224e4-00bf-434e-aa21-1b615b9bd6d0",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "Opt lambda is 0.142 and corresponding sparsity level is 8.5%"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3caf018b-eaf0-43ac-8101-21fd0791fd88",
   "metadata": {
    "tags": []
   },
   "source": [
    "### SGL + low rank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da3c4ef9-18bb-4826-877d-fdd0e3e4a09b",
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_lambda = [0.142379]\n",
    "# mu1_range = np.logspace(0.8, 0.9, 30)\n",
    "opt_mu = [7.94]\n",
    "\n",
    "N = sub_corr.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c53f7d46-eddf-45a0-84ca-e85aafcc7d7d",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "P_low = glasso_problem(sub_corr[0,:], N, reg_params = {'lambda1': opt_lambda, 'mu1_range': opt_mu},\n",
    "                       latent = True, do_scaling = False)\n",
    "print(P_low)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a8af825-e228-42b3-ba86-751e4791ffd5",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "modelselect_params = {'lambda1_range': opt_lambda, 'mu1_range': opt_mu}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7524abce-8419-4f8b-ad21-8083d2975fd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_low.model_selection(modelselect_params = modelselect_params, method = 'eBIC', gamma = 0.3)\n",
    "\n",
    "print(P_low.reg_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "220e73c3-6dcd-4d4c-89bb-c001b1f2c20a",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_low.modelselect_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "874fa151-c1ed-4d56-855a-1453e4bec236",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "# 10 Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88221a53-9a63-4e3c-81c9-7a043eea86cb",
   "metadata": {},
   "source": [
    "### SGL over separate 10 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a198e7a9-25fa-4eb6-b424-9e07bb8f8c33",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "lambda1_range = np.logspace(-0.3, -1, 10)\n",
    "N = sub_corr.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4aa8992-f2f9-47ac-9ef0-44702a153705",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "est_uniform, est_indv, statistics = K_single_grid(sub_corr[:9, :], lambda1_range, N,\n",
    "                                                  method = 'eBIC', gamma = gamma, \n",
    "                                                  latent = False, use_block = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e5c5476-3d92-4187-8aee-c243c71c3a4b",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3829230f-2452-4849-86fa-bd71323f1f86",
   "metadata": {},
   "outputs": [],
   "source": [
    "Theta = est_uniform[\"Theta\"]\n",
    "\n",
    "Theta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f95e1d7d-ae75-417d-a78f-39b12ad1854d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "K = Theta.shape[0]\n",
    "\n",
    "\n",
    "for i in range(0, K):\n",
    "    fig, axes = plt.subplots(1, 2, sharex=True, figsize=(30,30))\n",
    "\n",
    "    ax = sns.heatmap(sub_corr[i, :], ax=axes[0], center=0, vmin = -0.5, vmax = 0.5, \n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Covariance')\n",
    "    \n",
    "    ax = sns.heatmap(Theta[i, :], ax=axes[1], center=0, vmin = -0.5, vmax = 0.5,\n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Inverse Covariance')\n",
    "    \n",
    "    fig.show()\n",
    "    \n",
    "    # fig.savefig(\"plots/SGL_plots/SGL_lowrank_heatmap{0}.png\".format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a17cf2f2-e812-4db7-b863-9645c2ac433d",
   "metadata": {},
   "source": [
    "stopped here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a5fce08-272b-466e-961e-255ccee6beed",
   "metadata": {},
   "source": [
    "## The patterns are somewhat similar, it means we can solve group graphical lasso"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca25c3cb-ab34-4868-8f59-160dcdc3f7a5",
   "metadata": {},
   "source": [
    "### SGL + low-rank over separate 10 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3b7b723-eeed-4c60-b6e2-d31709904146",
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda1_range = np.logspace(-0.9, -1.5, 4)\n",
    "# mu1_range = np.linspace(8.75, 6.25, 5)\n",
    "mu1_range = np.arange(6.25, 11,  0.5)[::-1]\n",
    "gamma = 0.3\n",
    "\n",
    "N = sub_corr.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a362a15-77d8-43dc-bc43-4321958ad7f6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "low_est_uniform, low_est_indv, low_statistics = K_single_grid(sub_corr, lambda1_range, N,\n",
    "                                                  mu_range = mu1_range,\n",
    "                                                  method = 'eBIC', gamma = gamma, \n",
    "                                                  latent = True, use_block = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1d83cc2-442e-48e1-8601-8e4be657102e",
   "metadata": {},
   "outputs": [],
   "source": [
    "low_est_indv[\"Theta\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49de6fe3-5a83-4815-9363-8b347464ef2f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "low_statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b34433a3-a07c-4976-a9e9-6034041c1c47",
   "metadata": {},
   "source": [
    "The solver always picks the largest mu until rank becomes 0, so we make a conclusion there are no latent variables?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afd637aa-9c01-4743-b173-de0b0792c4af",
   "metadata": {},
   "outputs": [],
   "source": [
    "Theta = low_est_uniform[\"Theta\"]\n",
    "L = low_est_uniform[\"L\"]\n",
    "\n",
    "Theta.shape, L.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8667ba6f-cdc4-4ea4-92b1-cb2459804805",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "K = Theta.shape[0]\n",
    "\n",
    "\n",
    "for i in range(0, K):\n",
    "    fig, axes = plt.subplots(1, 3, sharex=True, figsize=(30,30))\n",
    "\n",
    "    ax = sns.heatmap(sub_corr[i, :], ax=axes[0], center=0, vmin = -0.5, vmax = 0.5, \n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Covariance')\n",
    "    \n",
    "    ax = sns.heatmap(Theta[i, :], ax=axes[1], center=0, vmin = -0.5, vmax = 0.5,\n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Inverse Covariance')\n",
    "    \n",
    "    ax = sns.heatmap(L[i, :], ax=axes[2], center=0, vmin = -0.5, vmax = 0.5,\n",
    "                 square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Low-rank')\n",
    "    \n",
    "    fig.show()\n",
    "    \n",
    "    # fig.savefig(\"plots/SGL_plots/SGL_lowrank_heatmap{0}.png\".format(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ddf9135-929b-430c-865d-0a18c643c2d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "124bf8b7-18b2-4cd3-bc6e-6d1b533b1a7d",
   "metadata": {},
   "source": [
    "## Low-rank matrix is empty, so no latent confouners!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2276ff46-0beb-40e4-93f6-aa65efaf2dea",
   "metadata": {},
   "source": [
    "### MGL (Group Graphical Lasso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "604fabfb-1e84-4fec-a66e-d7c062fdfeaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lambda1_range = [0.17, 0.14237977, 0.1]\n",
    "lambda1_range = np.logspace(-0.8, -1.1, 5)\n",
    "lambda2_range = np.logspace(-0.8, -1.1, 5)\n",
    "N = sub_corr.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30c37966-e28b-45ad-87bf-92fb5862fc57",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_10 = glasso_problem(sub_corr[:9,:], N, reg_params = {'lambda1_range': lambda1_range, 'lambda2_range': lambda2_range}, \n",
    "                      latent = False, do_scaling = False)\n",
    "print(P_10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6463177d-d0d1-48d3-aa90-761056316947",
   "metadata": {},
   "outputs": [],
   "source": [
    "modelselect_params = {'lambda1_range': lambda1_range, 'lambda2_range': lambda2_range}\n",
    "gamma = 0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9036d9e-746f-4329-84ad-9cd3d78a0431",
   "metadata": {},
   "outputs": [],
   "source": [
    "modelselect_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f2646b-9e6f-4055-a214-22bfcabcf494",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "P_10.model_selection(modelselect_params = modelselect_params, method = 'eBIC', gamma = gamma)\n",
    "\n",
    "# regularization parameters are set to the best ones found during model selection\n",
    "print(P_10.reg_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e6c665d-6b49-4fdd-8211-24c3a8d47932",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.logspace(-0.8, -1.1, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a4a8b4d-f94d-4558-a05b-cc058640a08f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "statistics_10 = P_10.modelselect_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76db3d69-5640-4762-a069-f0129d53c07c",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_10.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "349de764-9abc-43b6-9dc9-70453156b3cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_10.solution.precision_.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "502bc778-4d55-4a42-a99c-0d7bf9a1e0d7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "K = P_10.solution.precision_.shape[0]\n",
    "\n",
    "\n",
    "for i in range(0, K):\n",
    "    fig, axes = plt.subplots(1, 2, sharex=True, figsize=(30,30))\n",
    "\n",
    "    ax = sns.heatmap(sub_corr[i, :], ax=axes[0], center=0, vmin = -0.5, vmax = 0.5, \n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Covariance')\n",
    "    \n",
    "    ax = sns.heatmap(P_10.solution.precision_[i, :], ax=axes[1], center=0, vmin = -0.5, vmax = 0.5,\n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Inverse Covariance')\n",
    "    \n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2b241f-2336-43ea-8979-4996a4f50623",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "statistics_10['L1'].shape, statistics_10['L2'].shape, statistics_10['BIC'][0.3].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "488b5274-1df4-42f0-b460-a6f29c03e7da",
   "metadata": {},
   "outputs": [],
   "source": [
    "statistics_10['BIC'][0.3].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "223b13c8-9e65-4762-bb22-33c611fbb424",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot for GGL\n",
    "fig = surface_plot(statistics_10['L1'], statistics_10['L2'], statistics_10['BIC'][0.7], name = 'eBIC')\n",
    "fig.savefig(\"plots/GGL_l1_l2_10.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6303338-7e9f-4306-b0b4-b1891ff56968",
   "metadata": {},
   "source": [
    "### Read solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "803a9ff4-ee10-48a0-854f-4ba70a84d4d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = []\n",
    "sol = []\n",
    "\n",
    "for i in range(0, 951):\n",
    "    corr.append(genfromtxt(\"/storage/groups/bds01/datasets/brains/corr_matrices/corr{0}.csv\".format(i), delimiter=','))\n",
    "    sol.append(genfromtxt(\"/storage/groups/bds01/datasets/brains/est_uniform/est_uniform{0}.csv\".format(i), delimiter=','))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8211cef5-fbbe-4ba0-a1d9-01f6cd121bee",
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = np.array(sol)\n",
    "corr = np.array(corr)\n",
    "corr.shape, sol.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "836de53d-47cc-44bf-a30c-664446bfa567",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = pd.read_csv(\"/storage/groups/bds01/datasets/brains/train.csv\", sep=',')\n",
    "df_test = pd.read_csv(\"/storage/groups/bds01/datasets/brains/submission_valid.csv\", sep=',')\n",
    "\n",
    "df = pd.concat([df_train, df_test], ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e40687a0-7831-4f31-a03f-8a687251c669",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ca6c7dc-1f38-4d74-b519-948d8c41f822",
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"sex_f0_m1\"].isna().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "137817a2-954b-4a42-b2fc-88e450643363",
   "metadata": {},
   "outputs": [],
   "source": [
    "sex = df[\"sex_f0_m1\"].reset_index()\n",
    "\n",
    "male_ixs = np.array(sex[sex[\"sex_f0_m1\"] == 1.0][\"index\"])\n",
    "female_ixs = np.array(sex[sex[\"sex_f0_m1\"] == 0.0][\"index\"])\n",
    "\n",
    "male_ixs.shape, female_ixs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59a0e61e-b8ac-4a04-ab8f-d9424384ea20",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_edge_probablity(data=list, indices=list):\n",
    "    \n",
    "    P = np.zeros((436, 436))\n",
    "    \n",
    "    for i in indices:\n",
    "        x = data[i, :]\n",
    "        x = (x != 0).astype(int)\n",
    "\n",
    "        P = P + x\n",
    "    \n",
    "    P = P / len(indices)\n",
    "    \n",
    "    return pd.DataFrame(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a66dcb2-ca49-4169-84fa-5038c51ea14d",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_male = calculate_edge_probablity(data=sol, indices=male_ixs)\n",
    "P_female = calculate_edge_probablity(data=sol, indices=female_ixs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23bcb860-1e64-4357-8811-5ce6abd970fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_male[P_male < 0.8] = 0\n",
    "P_female[P_female < 0.8] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9dc4b90-bdaf-4798-9ea6-bcfffb9ed898",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# plt.figure(figsize=(16,16))\n",
    "# ax = sns.heatmap(P_male, center=0, vmin = -0.5, vmax = 0.5, square = True, cbar = False, \n",
    "#                  cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "fig, axes = plt.subplots(1, 2, sharex=True, figsize=(30,30))\n",
    "\n",
    "ax = sns.heatmap(P_male, ax=axes[0], center=0, vmin = -0.5, vmax = 0.5, \n",
    "                 square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "ax.set_title(\"Male's connectivity network\")\n",
    "\n",
    "ax = sns.heatmap(P_female, ax=axes[1], center=0, vmin = -0.5, vmax = 0.5,\n",
    "                 square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "ax.set_title(\"Female's connectivity network\")\n",
    "\n",
    "fig.show()\n",
    "\n",
    "fig.savefig(\"/mnt/home/icb/oleg.vlasovetc/brain_challenge/f-threshold-select/plots/connectivity_heatmap.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea99ce3d-960c-481d-b750-bdb0cf897ac7",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = P_male.astype(bool)\n",
    "b = P_female.astype(bool)\n",
    "\n",
    "a == b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ea59753-5cc0-4ed9-955b-68d121b6b79b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,16))\n",
    "ax = sns.heatmap(a == b, center=0, vmin = -0.5, vmax = 0.5, square = True, cbar = False, \n",
    "                 cmap = \"coolwarm\", xticklabels=False, yticklabels=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b29cbcd-7c3a-4e32-af0e-98784a48fc1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,16))\n",
    "ax = sns.heatmap(P_df, center=0, vmin = -0.5, vmax = 0.5, \n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4583333b-3eee-4c73-9c69-6768e9cb7ec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "test = sol[:5, :]\n",
    "test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e2bb9f5-d9de-46c5-bfbe-7b6d98eb44a6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "empty = np.zeros((436, 436))\n",
    "\n",
    "K = test.shape[0]\n",
    "\n",
    "for i in range(0, K):\n",
    "    x = test[i, :]\n",
    "    x = (x != 0).astype(int)\n",
    "    \n",
    "    empty = empty + x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20a8ef56-ab0c-43f0-8919-b60c691aed3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "empty = empty / K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45def222-64db-4e44-a71d-14d950c911d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "empty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c2f5e11-8c52-4cbb-9704-c99b57abc7b1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85c17ba0-460d-4603-8afd-cbe8ce21cded",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "K = sol.shape[0]\n",
    "\n",
    "\n",
    "for i in range(0, K):\n",
    "    fig, axes = plt.subplots(1, 2, sharex=True, figsize=(30,30))\n",
    "    \n",
    "    ax = sns.heatmap(corr[i, :], ax=axes[0], center=0, vmin = -0.5, vmax = 0.5, \n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Covariance')\n",
    "    \n",
    "    ax = sns.heatmap(sol[i, :], ax=axes[1], center=0, vmin = -0.5, vmax = 0.5,\n",
    "                     square = True, cbar = False, cmap = \"coolwarm\", xticklabels=False, yticklabels=False)\n",
    "    ax.set_title('Inverse Covariance')\n",
    "    \n",
    "    fig.savefig(\"/storage/groups/bds01/datasets/brains/plots/SGL_plots/SGL_heatmap{0}.png\".format(i))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}