test_realworld.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Seed set to 42\n"
     ]
    }
   ],
   "source": [
    "from pytorch_lightning import seed_everything\n",
    "seed_everything(42)\n",
    "\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from networkx.drawing.nx_agraph import graphviz_layout\n",
    "from causaldag import partial_correlation_suffstat, partial_correlation_test, MemoizedCI_Tester, gsp, pcalg, hsic_test, kci_test\n",
    "\n",
    "from fci import fci\n",
    "from ccpg import ccpg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = './real-data/airfoil-self-noise.continuous.txt'\n",
    "samples = np.loadtxt(file, skiprows=1)\n",
    "nnodes = samples.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(file, 'r') as f:\n",
    "    names = f.readline().strip().split('\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "components [{'Chord'}, {'Velocity'}, {'Displacement', 'Pressure', 'Frequency', 'Attack'}]\n",
      "edges {(0, 2), (1, 2)}\n"
     ]
    }
   ],
   "source": [
    "# run ccpg\n",
    "suffstat = partial_correlation_suffstat(samples)\n",
    "ci_tester = MemoizedCI_Tester(partial_correlation_test, suffstat, alpha=1e-3)\n",
    "\n",
    "c, e = ccpg(set(range(nnodes)), ci_tester, verbose=False)\n",
    "print('components', [{names[j] for j in i} for i in c])\n",
    "print('edges', e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAAEDCAYAAAAx0WHLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAexAAAHsQEGxWGGAAAfMElEQVR4nO3de3RTdbrw8W9CW25TWqEoyCCXegovZYQBFS9jCwqIUNCBzgyICwvUkRmxYoEji0EoKHpGzgIV9eXoFDrwukYpeKs6Hh0F9IxnDssjra0WGW5alUur3O01+b1/JBtKoJCke2ffns8/2jZJf8luHr7Z2Uk8SimFEML1vGYvQAhhDTIMhBCADAMhRJAMAyEEIMNACBEkw0AIAcgwEEIEyTAQQgAyDIQQQTIMhBCADAMhRJAMAyEEIMNACBEkw0AIAcgwEEIEyTAQQgAyDIQQQXFmL8A+/MA+oAyoBhqBeKArMAjog8xWM8n2aS2PvO1ZOHYB24BTgAdofpNpX3cEMoG0mK9OyPbRgwyDi/oU+DCC02cAQwxaiziXbB+9SDdd0C4i+0MjePpdBqxFnEu2j55kGLTITyA9o7EteH5hHNk+epNh0KJ9BB6DRuNU8PwBR44cYe3atTQ1NemxMAHouX1EgOuHQVFREVlZWdx999088cQTAPj9fgJ7pT1RXqoHKOPIkSMsXLiQHj16MHPmTKqqqnRatXMVFRUxduxYZs2axWuvvQZo2yNU67ePOJvrdyAWFRWRkpJCVlYWXbp04bbbbuPqq68mPX0fb71VSm1tI5MmXUVNzSm2bt1NYmI7HntsHPffv5lOndqRnt6NESP+hWee+Yh///fbWbCghFmzbmTx4nfZsOF/zvpdiYmJdO/enfj4eJOurfUdOXKEuLg4EhMT2bt3L4mJiXTs2JHjx4+j/al269aNceN64PX62LnzMLm51/HJJ1VkZw9i4MDuTJ68npdemsby5e9RU3OSEyfqWbXqDqZM2cDw4VdSUXGA/Pxb6dVrHnl5eaSkpJCZmcnWrVuZO3cuKSkp5ObmsmHDBpNvjdiS4wyAF154gddff50lS5bQpUsXpk6dyvjxP2Po0B4kJ7dn+/avUQquuupybr99IG3bxnHw4Aluv30gt9ySxsGDJ865TK9XkZeXx/r16zl27BgAmZmZZGRkkJGREeuraBtvvvkm77//Pt26daNt27Y899xzfPzxx7z11lv07t2bgwcPkpuby+uv/5F1637DihUfnPdyKisP8eGHe7j++t7U1TVRWXkIgPz8TD7//CBvvllGUtKL3HXXXYwaNQqAAQMG8Pzzz9O3b18mT54cs+tsFTIMgHvuuYesrCyKiopISkoCwO+HRYtGERfX5vTpysq+Zf78Eh555DZefnkaW7bsJifnL6xadTtNTYGUPXWqIXhqL8uXL2fx4sWMHj0agA0bNpCcnBzLq2Y7lZWVDBs2jKysLLKzsxk2bBg1NTW0adOGvLw8AI4dO0ZJyQoAPJ7AQ4W2beNoavKjlKK2thG/3096ejcKCsacvux27eLwer3Ex7ehvj5wWq/3zCPlK6+8ku+++47KykpeeumlGF5ra5Bh0IK8vDvIzX2Zzp07cPXVPTl5sp5//rMGr9dDp07tmDv3Ddq3jyctrSvdunXi2LE6Vq7cymefHQheQjsAunTpwm9+8xuqqqpYtGgRM2bMYMgQeZ47ErfeeiuzZs1i/vz5HD16lNWrV3PZZd34138tYdeuw/Tvfym33fZ/ePrpD0lP74ZSivT07ni9HvLzX6O2tpGFC0c2u0QP0JGpU6cyZ84c3nvvPW666SbGjRvHLbfcwoEDB84aEm7h+n0GLdsDlLTi/OOBVJ3WIkLt3fs3+vatoKhoOykpHcnKSo/wEs7dPn//+9956qmnKCwsJDExUbe12oUMgxb5gUKiefqqurqW9PSVDBkylPT0dK655hpXPgY1QkNDA4899hjPPPM0X321mA4dwBPxkwodgZnIk2lnk2FwQbuAtyM+11//6mHs2Dmnvx4wYACff/65fstyqdLSUnJyckhOTqawsJDUVB/RbB8Yi7xG4VwyGi8ojcCx7Bd3ZqRmMGZMHoMHDz79s5tuuomGhobznU2EoaGhgYKCAkaMGEFubi4ffPABqampRLJ9zshABsH5yTC4qCEE/iXpGPw6tEkDX1dX/4jfPwYYgsfj4YUXXgBg1qxZVFZWcu2111JaWhqbJTtIaWkp1157LVu3buWTTz5h9uzZITv3wts+J08SPJ3svG2REmHyKaV2K6U2q/r6p9WpU08ov3+NUmqz8vn+qX72s4GquLj4rHNs27ZN1dfXK5/Pp55++mmVlJSkCgoKVH19vRlXwFbq6+vVkiVLVHJyslq9erXy+XwXOceZ7aPUGqXU6uB/N6vs7J+rNm28Kjc3V9XW1hq8cvuSYRCFG2+8UQGqqKjo9Pc2btyoBg4ceME/2t27d6uMjAw1aNAgtWPHjhis1J527NihBg0apDIzM9Xu3btbfXmdO3dWgEpISFD9+/eXgdACGQYRKi0tVfHx8QpQl19+uWpsbFRKKeXz+VR6evo5dRBKKqFlkddAeFJSUpTX61Ver1fNnj1bNTU16XK5TiPDIEITJkxQXq9XAcrj8agXX3zx9M/CqQONVMLZ9K6B5kaOHKmmT5+ufvKTn6idO3fqetlOIsMgQq+++qrKy8tTl1xyiVq+fLmqqKg4/bNw66D56d1eCUbVwPksWLBA3XXXXYZdvt3JMIhCeXm56tOnz3l/FkkdaNxaCUbWwPlUV1erTp06SR20QJ5a1NmkSZNQSvHKK6+EfZ7U1FS2bNnCzJkzGT58OEuXLnX0cQktHzdgrJSUFH7/+9/z6KOPGv67bMnsaWRHFyoDpaKrA43TKyHWNRBK6qBlUgYGiKYONE6tBLNqIJTUwQWYPY3s6GJloFTr6kDjlEowuwZCSR2cn5SBQVpTBxq7V4JVaiCU1EELzJ5GdhROGSilTx1o7FYJVquBUFIH55IyMJAedaCxSyVYtQZCSR2ch9nTyI7CLQOl9K0DjVUrweo1EErq4GxSBgbTsw40VqsEu9RAKKmDEGZPIzuKpAyUMqYONGZXgt1qIJTUwRlSBjFgRB1ozKoEu9ZAKKmDZsyeRnYUaRkoZWwdaGJVCXavgVBSBwFSBjFiZB1ojK4Ep9RAKKmDILOnkR1FUwZKxaYONHpXgtNqIJTUgZRBTMWiDjR6VYJTayCU1AFSBtGItgyUim0daKKtBKfXQCi314GUQYzFsg40kVaCW2oglOvrwOxpZEetKQOlzKkDzcUqwW01EMrNdSBlYAIz6kDTUiW4tQZCubkO5LMWo1BRUcGECRPYu3dv1JdRXFzMsmXLKCsrM+3jv/fs2cOMGTM4ePAgAN27dw9+hqH7hkBzNTU1pKamsn37dvr162f2cmJGysAkZtaBpmfPngwfPpyqqiqqqqoYPnw4PXv2NG09VuHWOpBhYBKv18uSJUtYunQpfr8/5r9f+wzDbdu2UV5eTnl5OVu2bJHPhAyaO3cub7zxBl9++aXZS4kZGQYmMqMOWto3YLVXQprNlXVg6u5Lm2rtswnNxfKZhXCfKTD7lZBW4bZnFqQMTBaLOoj0mQKphADX1YHZ08iO9CwDpYytg9YeN+D2SnBTHUgZWIARdaDXcQNurwRX1YHZ08iO9C4DpfStA6OOInRrJbilDqQMLEKPOjD6KEK3VoJr6sDsaWRHRpSBUq2rg1i/psBtleCGOpAysJBo6sCs1xS4rRJcUQdmTyM7MqoMlIqsDqzyCkO3VILT60DKwGLCqQOrvcLQLZXg+DowexrZkZFloNSF68AqNdASp1eCk+tAysCCzlcHVquBlji9EhxdB2ZPIzsyugyUOrsOrF4DLXFqJTi1DqQMLGrSpEn4fD4mT55s+RpoiVMrwbF1YPY0sqNYlMGOHTtUr169VIcOHdSuXbsM/V2x4LRKqK6uVklJSY6qAykDi2m+b2Du3Ln07t2bsrIys5fVak6rhJSUFH73u985qw7MnkZ2ZFQZnG/fgJnvpGwUp1SC0/YdSBlYwIWeKbDCeyXqzSmV4Lh9B2ZPIzvSswzCeabAiXWgsXslOKkOpAxMEslxA06sA43dK8FRdWD2NLKj1pZBNMcNOLkONHatBKfUgZRBDLXmKEIn14HGrpXgmDowexrZUTRloMdRhG6oA43dKsEJdSBlYDA9X1PghjrQ2K0SHFEHZk8jOwq3DIx4TYGb6kBjl0qwex1IGRjAyFcYuqkONHapBNvXgdnTyI4uVAaxeIWhG+tAY/VKsHMdSBnoJJbvN+DGOtBYvRJsXQdmTyM7Ci0DM95vwM11oLFqJdi1DqQMWsHMdx9ycx1orFoJtq0Ds6eRHZWXl6sePXqY/u5DUgdnWK0S7FgHUgYRamho4LnnnuPAgQOmv/uQ1MEZVqsEO9aBRymlzF6EXZSWlpKTk0N8fDyHDh3i66+/NntJFBcXs2zZMsrKyvB6ZbYD7NmzhxkzZnDs2DGKiooYPHiwKeuoqakhNTWV7du3069fP1PWEAn56wlD6L6BwsJC4uLizF4WIHVwPlapBNvVgbmPUqzvfM8UxOI9ECMh+w5aZva+BDvtO5AyaIFdPqcApA4uxOxKsFUdmD2NrOhixw1YrQyUkjoIh1mVYJc6kDJoxk41EErq4OLMqgTb1IHZ08gqIjmK0IploJTUQSRiXQl2qAPXl4GdayCU1EH4Yl0JtqgDs6eRmaJ9TYFVy0ApqYNoxKoSrF4HriwDJ9VAKKmDyMWqEixfB2ZPo1jT4xWGVi4DpaQOWsPoSrByHbimDJxcA6GkDqJndCVYug7MnkaxoPf7DVi9DJSSOtCDUZVg1TpwdBm4qQZCSR20nlGVYNk6MHsaGcXIdx+yQxkoJXWgJ70rwYp14LgycHMNhJI60I/elWDJOjB7GukpVu9FaJcyUErqwAh6VYLV6sARZSA10DKpA/3pVQmWqwOzp1FrmfHOxHYqA6WkDozU2kqwUh3YtgykBsIndWCc1laCperA7GkUDTNqoDm7lYFSUgexEG0lWKUObFUGUgPRkzowXrSVYJk6MHUURcDsGmjOjmWglNRBLEVaCVaoA8uXgdSAfqQOYifSSrBEHZg2hsJgpRpozq5loJTUgRnCrQSz68CSZSA1YBypg9gLtxLMrgODP1HJD+wDyoBqoBGIB7oCg4A+hH6Oi/apRcnJyRQWFlpyCFRUVDBhwgT27t1r9lKiIp/CZJ6LfdpTTU0NaWn/QkXFG1x+eTXh3m/0YOBfwi6gECgBqoBaoCn436rg9wuDp5MaiCWpA/NcrBJSUn7g22+Xcfnl/0s49xs9GVQGnwIfhn3qb79NZdy4hy1dA83ZvQxA6sAKzq0EP5HcbyADGKLbegz4K9hFZFcIevTYw7/92wypgRiSOjBf80p48sn7iPR+Ezi9foWg8zDwA9siPpdSMGZMe+QfqNjxer0sWbKEpUuX4vf7zV6Oa3m9Xu6//z7Wrp1GdI2+jcD9Toe16HIpp+0DTkV8Lo+H4Pn26bsccUFSB1axD6+3Nng/iJR+9xtdhsGnn35Kr169qK3dTnZ20envZ2evC+v8+/f/wLx5bxB41kHEitRB7E2cOJG6ujoAysvLeeihhwj83QcmQbj3mfnz3wievgi97je6DIP169ezbNkyiovf5YsvDlJQ8A7/9V97+eKLQzz66Lt8990x/vCHt5g1ayOvvVYOBK7Mgw++yqpVW4OXovjkk0/Iy8ujsbFRj2WJMEgdxNakSZPYvHkzAGvXrmXXrl3MmfMs99+/CZ/vzEBeuvQd5sx5lenT/0JDQxP/+Z87mT79L+Tnv8ahQyfYt+97vvnmaPD+tpZt27bxxBNPAJCfn09VVVXEa4tr7ZWrq6ujpqaGqVOnMmXKkwwYcBkFBWMAGDDgMhYtGs3hwyeor2/isssSefHF/yUtrSsJCW1YsWICECiDf/zjK5Ty8OSTf5O92zHUvA4mTpwot73BJk2axF133cWvfvUrioqKGD16NMnJbfnqq5N8++0xAI4fr2P//iOsWzeFP/95O++++yVr1/4PmzdPx9PsscRPf5rc7P6WyZ/+9CeOHj3KsWPH6NmzZ8Rra/WW37x5MwcPHmT27Nns2lWNz3dmL4i28A0bPmHChIEsXDiKEyfqUUrh9Z79AKl7904cOnSKH3/8sbVLEhGSOoiddu3a0b17d1atWsWTTz7JjTfeSEHB7axbN4UrrrjknNN7PB6UUmcNgdCfBw5ICjwEmTx5MnfffXdUa9NlGJSUlLBmzRo2bMgnLa0r8+a9TlXVEeLj2/Dww29zww19KCz8BytXbiUhoQ3p6d2prW1k/vw3eOaZjwDo1aszjzySw7333svJkydbuywRAdl3EFvTp09nxYoVTJw4kc8++4z589/mnntepq4u8PC4U6d29Op1CXPnvs6WLbu59db+5OZex29/u5EFC0o4fPjE6ctKT+/GvHlvU1VVxfjx4zlx4gQZGRlRrUvng472EDhCKlrjAesfZ+CEg45C+f1+rrrqKgoKCsjOzjZ7OS7T+vtNXV0P8vLyGDt2LHfccUdUl6LzA8Q+QMeIz+X3K06eBL+/l77LEWGTOjBTdPebgI5AH9q1a8fzzz8f9SAA3YeBF8iM/FxeD8uWbeXmm0eyZ88efZckwib7DswS3f0mIBO97sYG7DpOI3DMdCQyePTRjYwYMYKhQ4eyevVq+dfJBFIHZorufhM4nz4Meh5pCDCWM+lz9p5QpQJf+3ztg6cbQkJCAkuWLGHr1q0UFhZy8803SyWYQOrATBe+35z5uiPa/UZPBj6pnAbMJLBTsCfQnsBhDe3xeHqycuU/Wbz4K0In2+DBg9m+fbtUgkmkDszW8v0m8PX44M/1K4LTYv/mSgGVlZUqKSlJVVdXt3gaedszc/h8PpWenq6Ki4vNXoqIIdMON+vfvz9ZWVmsXLmyxdNIJZhD6sCdTD32dNGiRTz33HPU1NS0eBrZl2AO2XfgPqYOg3DqQCOVEFtSB+5j+qtSwqkDjVRCbEkduIvpwyCSOtBIJcSG1IG7mD4MILI60EglxIbUgXtYYhhEUwcaqQRjSR24hyWGAURXBxqpBGNJHbiDZYZBa+pAI5VgDKkDd7DMMIDW1YFGKsEYUgfOZ6lhoEcdaKQS9CV14HyWGgagTx1opBL0JXXgbJYbBnrWgUYqQR9SB85muWEA+taBRipBH1IHzmXJYWBEHWikElpH6sC5LDkMwJg60EgltI7UgTNZdhgYWQcaqYToSB04k2WHARhbBxqphOhIHTiPpYdBLOpAI5UQGakD57H0MIDY1IFGKiEyUgfOYvlhEMs60EglhEfqwFksPwwgtnWgkUoIj9SBc9hiGJhRBxqphAuTOnAOnT+F2Tg7d+7kuuuuY/fu3aSkpJiyhtLSUnJycoiPj+fQoUN8/fXXpqzDauQTnJ3BFmUA5taBRquEa665hm+++UYqIUjqwBlsUwZgjToAqKioYMyYMaSkpJCcnExhYSGpqammrccKpA7szzZlANaoA01CQoLsS2hG6sD+bDUMwJxnFloizzicTZ5ZsDfbDQMr1YFGnnEIkDqwN9sNA7BWHWikEgKkDuzLlsPAinWgcXslSB3Yly2HAVizDjRurwSpA3uy7TCwch1o3FoJUgf2ZNthANauA41bK0HqwH5sPQzsUAcat1WC1IH92HoYgD3qQOO2SpA6sBfbDwM71YHGLZUgdWAvth8GYK860LilEqQO7MMRw8COdaBxeiVIHdiHI4YB2LMONE6vBKkDe3DMMLBzHWicWglSBzahHGTnzp0qKSlJVVdXG/p7ysvLVZ8+fQz9HTt27FCDBg1SmZmZavfu3Yb+rljw+XwqPT1dFRcXm70U0QLHlAFAv379GD9+vK3rQOO0SpA6sD5HDQOw976DUE7blyD7DqzNccPASXWgcUolSB1Ym+OGATirDjROqQSpA+ty5DBwYh1o7F4JUgfW5chhAM6sA43dK0HqwJocOwycXAcau1aC1IE1OXYYgLPrQGPXSpA6sB5HDwM31IHGbpUgdWA9jh4G4I460NitEqQOrMXxw8BNdaCxSyVIHViL44cBuKsONHapBKkD63DFMHBjHWisXglSB9bhimEA7qwDjdUrQerAGlwzDNxcBxqrVoLUgTW4ZhiAu+tAY9VKkDown6uGgdTBGVarBKkD87lqGIDUQXNWqwSpA3O5bhhIHZzLKpUgdWAu1w0DkDo4H6tUgtSBeVw5DKQOWmZ2JUgdmMeVwwCkDi7E7EqQOjCHa4eB1MHFmVUJUgfmcO0wAKmDcJhVCVIHsefqYSB1EL5YV4LUQey5ehiA1EEkYl0JUgex5fphIHUQuVhVgtRBbLl+GIDUQTRiVQlSB7EjwwCpg9YwuhKkDmJHhkGQ1EH0jK4EqYPYkGEQJHXQekZVgtRBbMgwaEbqoPWMqgSpA+PJMGhG6kA/eleC1IHxZBiEkDrQj96VIHVgLBkGIaQO9KdXJUgdGEuGwXlIHehPr0qQOjCODIPzkDowTmsrQerAODIMWiB1YJzWVoLUgTFkGLRA6sB40VaC1IFBlGjRzp07VWJiosrJyVEjR448/f3y8nLVp08fE1fmPDt27FCDBg1SmZmZavfu3Rc9vc/nU+np6aq4uDgGq3MHKYMW+Hw+Vq9ezalTp1i/fj2lpaVmL8nRIq0EqQP9yTBoQWNjI++//z5xcXH4/X75g4uBSPclyL4DfckwaEG7du3YsWMH06ZNw+PxUFtba/aSXCPcSpA60JdHKaXMXoTVPfXUU/zHf/xfvvjiLaAMn+8gDQ0/0r59ItAVGAT0QWar/kpLS8nJySE5OZnCwkJSU1PP+rnf7+eqq65i6dIlTJo0BCgDqoFGIB7ZPuGTYRCWXcA24BTgAZrfZNrXHYFMIC3mq3O6hoYGHn/8cVatWsUjjzzCfffdh9d75o793//9Z1JTv+XSSzsg2yd6Mgwu6lPgwwhOnwEMMWgt7nb+SghsH6UUHo8njEuR7dMS6aYL2kVkg4Dg6XcZsBYRui/hr399Cm37hDcIQLZPy6QMWuQHCgk8NIhUR2AmMmuNU1r6KT16vENKSvsIBoFGts/5yK3Ron1ENwgInm+fjmsRoQYPTqJr1w5RDAKQ7XN+jh4GRUVFZGVl8cADDzBjxgw+/vhjVq1addHz7d+/n3nz8gnsfIqGh8Be7chs3bqVZ555Jsrf6UxHjx5l6NChvPDCCzQ2NtLQ0MD999/P7Nmz+e1vN9K79zJOnqwP67LOPP3oIScnl5MnTxq3cBuKM3sBRps1axZZWVl8//33PPjgg1x66aVUVlaydOlSrrjiCqZNm8amTZs4fvw4iYmJpKamkpGRAfzIjz/W89hjf+Po0VoGDbqce+65nhUrPuDrr4+QnNyeRx4Zy/Ll71FTc5ITJ+pZteoOpk79fwwd+lN27DhAWto/2L9/P127dmXatGmsW7eOY8eOcerUKebMmcPixYsZMmQIe/bsYcqUKXz00UeUl5dTU1PDbbfdZvZNZwnfffcdZWVl5OXl8dBDDzFs2DBmz57NuHHfALXceecG/vjH9/nyy8Pk5l7H8OFXcu+9xSQltaNz5w4sXnwrN974FOPHp3PDDX1Ys+Zj+vTpTHn5XrOvmuU4fhhounTpQmNjIwCHDx/mkksu4c4772TgwIFs2rSJiRMn8otf/IJf//rXwWHgx+OBpiYfnTt3YOPGUqZMGUJ5+QHWr58KQGXlIT78cA/XX9+buromKisP4fP5+cMfRvHKK+U88MAGunTpQlVVFe+++y4HDhygQ4cONDU1UVFRQXV1NUeOHKG+vp65c+fSoUMH6urq2LRpE5s2bTLx1rKOhoYG/H4/dXV11NXV8c477zBv3jwCxxFAQkIbZs26gY4d2/LQQyU0NPjIyOjL9OnDmDnzJY4fr6NDhwQWLBhJSUkFo0alMX36MD799Ftzr5gFuWYYfP/99yQkJACQmZlJ3759KSwspKwskPNNTU0ApwcGeHn77UoGDOjGtGnXMGLEs8Gnr85cpt/vJz29GwUFY05/r337eOLi2tCpUycef/xxcnJyyM7OZunSpRQWFp71Ksjs7Gw2bdrEzp07eemll7j55pv57LPPmD17trE3ho3s27ePvn370r59e/Lz80lKSgpuq3igiYYGH0lJ7YmL81Jf33TepxiTktqd/v+EhMCffNu2CTG8Fvbg+GGwZs0a3nvvPU6cOEFBQQHPPvssW7ZsoaSkhB9++IFRo0axb98+Nm7cyMsvv8wvf/nL4Dk78POf/5QFC0o4cOA4Pp+fxMR2DBjQjQcffJXOnTvy8MOj8Xo95Oe/Rm1tIwsXjgye1wMknbWO9PR0vF4v+fn51NbWsnDhwnPWmpaWxsqVK+nYsSPTp0839Haxi549e1JYWMjEiRNJTk6moaGB/Px8Skr24PMd58CB42edfvTo/syaVUx5+QF69kymU6czg2DUqH488MCrVFefZP/+Y7G+KpYnTy0CBQUFZGdnM3DgwGbf3QOUtOJSxwOpFz2ViJZsH73JMGiRHGdgbbJ99Ca3Rou8BI5lj0YmctMaTbaP3uQWuaA0AseyRyIDeTFMrMj20ZM8TAiLvGrR2mT76EGGQdj8BA5hldfLW5Nsn9aSYSCEAGRUCiGCZBgIIQAZBkKIIBkGQghAhoEQIkiGgRACkGEghAiSYSCEAGQYCCGCZBgIIQAZBkKIIBkGQghAhoEQIkiGgRACkGEghAiSYSCEAGQYCCGC/j9tyfO+bYDEbAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 250x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# run pc\n",
    "\n",
    "est_dag = pcalg(set(range(nnodes)), ci_tester)\n",
    "\n",
    "est_arcs = set()\n",
    "est_edges = set()\n",
    "for (i,j) in est_dag.arcs:\n",
    "    if (names[j], names[i]) in est_arcs:\n",
    "        est_edges.add(frozenset({names[j], names[i]}))\n",
    "        est_arcs.remove((names[j], names[i]))\n",
    "    else:\n",
    "        est_arcs.add((names[i], names[j]))\n",
    "\n",
    "G = nx.DiGraph(\n",
    "    est_arcs\n",
    ")\n",
    "\n",
    "E = nx.Graph(est_edges)\n",
    "\n",
    "plt.figure(figsize=(5, 5), dpi=50)\n",
    "pos = nx.circular_layout(G)\n",
    "nx.draw(\n",
    "    E, pos=pos, node_color='w',linewidths=2.0, width=1.5\n",
    ")\n",
    "nx.draw(\n",
    "    G,\n",
    "    arrowsize=12,\n",
    "    with_labels=True,\n",
    "    node_size=500,\n",
    "    node_color=\"#ffff8f\",\n",
    "    linewidths=2.0,\n",
    "    width=1.5,\n",
    "    font_size=10,\n",
    "    pos = pos\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dag",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}