{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "constraint_based-0",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Constraint-based discovery (PC / FCI / GFCI / CD-NOD)\n",
    "\n",
    "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/Abel-ai-lab/andrey/blob/main/examples/constraint_based.ipynb)\n",
    "\n",
    "The constraint-based methods use conditional-independence tests instead of a score. `pc` returns a\n",
    "CPDAG; `fci` returns a PAG (partial ancestral graph) that stays sound under latent confounding.\n",
    "`gfci` builds a PAG from a GES graph, and `cdnod` learns from data pooled from several domains."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fb27b941602401d91542211134fc71a",
   "metadata": {},
   "source": [
    "## Start\n",
    "\n",
    "Choose a CPU runtime and run this cell first. In Colab it installs `andrey-core==0.9.0` from PyPI;\n",
    "local kernels use the installed environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "acae54e37e7d407bbb7b55eff062a284",
   "metadata": {
    "tags": [
     "colab-setup"
    ]
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "IN_COLAB = \"google.colab\" in sys.modules\n",
    "\n",
    "if IN_COLAB:\n",
    "    os.environ[\"ANDREY_DEVICE\"] = \"cpu\"\n",
    "    !pip install --quiet \"andrey-core[viz]==0.9.0\" pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "constraint_based-1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import andrey"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constraint_based-2",
   "metadata": {},
   "source": [
    "The same linear-Gaussian chain `x0 -> x1 -> x2`. Each example here lists its columns in an order\n",
    "drawn from the seed, so their order says nothing about the causal order; `shuffle=False` keeps the\n",
    "order the variables are written in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "constraint_based-3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "columns: ['x2', 'x1', 'x0']\n"
     ]
    }
   ],
   "source": [
    "CHAIN = [(\"x0\", \"x1\"), (\"x1\", \"x2\")]  # the true edges, cause first\n",
    "\n",
    "\n",
    "def sample_chain(seed: int, n: int = 1500, shuffle: bool = True) -> pd.DataFrame:\n",
    "    \"\"\"Draw n rows of the chain, its columns in an order drawn from the seed.\"\"\"\n",
    "    rng = np.random.default_rng(seed)\n",
    "    x0 = rng.standard_normal(n)\n",
    "    x1 = 0.8 * x0 + rng.standard_normal(n)\n",
    "    x2 = 0.8 * x1 + rng.standard_normal(n)\n",
    "    frame = pd.DataFrame({\"x0\": x0, \"x1\": x1, \"x2\": x2})\n",
    "    return frame[rng.permutation(frame.columns)] if shuffle else frame\n",
    "\n",
    "\n",
    "X = sample_chain(seed=2)\n",
    "print(\"columns:\", list(X.columns))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constraint_based-4",
   "metadata": {},
   "source": [
    "PC with a Fisher-Z independence test at alpha = 0.05."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "constraint_based-5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PC  -> cpdag [('x2', 'x1', 'undirected'), ('x1', 'x0', 'undirected')]\n"
     ]
    }
   ],
   "source": [
    "pc_out = andrey.pc(X, alpha=0.05)\n",
    "print(\"PC  ->\", pc_out.structure.kind, pc_out.structure.oriented_edges())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constraint_based-6",
   "metadata": {},
   "source": [
    "FCI returns a PAG, whose endpoint marks encode possible latent confounding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "constraint_based-7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FCI -> pag\n"
     ]
    }
   ],
   "source": [
    "fci_out = andrey.fci(X, alpha=0.05)\n",
    "print(\"FCI ->\", fci_out.structure.kind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "constraint_based-8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 360 240\" width=\"100%\" font-family=\"Inter, system-ui, sans-serif\"><rect width=\"360\" height=\"240\" fill=\"#fafbfc\"/><line x1=\"93.00\" y1=\"120.00\" x2=\"167.00\" y2=\"120.00\" stroke=\"#e0982b\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><line x1=\"193.00\" y1=\"120.00\" x2=\"267.00\" y2=\"120.00\" stroke=\"#e0982b\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><circle cx=\"80.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"80.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x2</text><circle cx=\"180.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"180.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x1</text><circle cx=\"280.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"280.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x0</text></svg>"
      ],
      "text/plain": [
       "'<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 360 240\" width=\"100%\" font-family=\"Inter, system-ui, sans-serif\"><rect width=\"360\" height=\"240\" fill=\"#fafbfc\"/><line x1=\"93.00\" y1=\"120.00\" x2=\"167.00\" y2=\"120.00\" stroke=\"#e0982b\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><line x1=\"193.00\" y1=\"120.00\" x2=\"267.00\" y2=\"120.00\" stroke=\"#e0982b\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><circle cx=\"80.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"80.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x2</text><circle cx=\"180.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"180.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x1</text><circle cx=\"280.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"280.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x0</text></svg>'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import andrey.viz\n",
    "\n",
    "andrey.viz.draw(pc_out.structure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constraint_based-10",
   "metadata": {},
   "source": [
    "## GFCI\n",
    "\n",
    "GFCI starts from the GES graph, removes the edges that independence tests reject, and\n",
    "orients the rest with FCI's rules, so it also returns a PAG. Here `c` and `d` share a hidden\n",
    "cause, and `a` and `b` each cause one of them; the PAG marks `c <-> d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "constraint_based-11",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('d', 'c', 'bidirected'), ('b', 'd', 'partially_directed'), ('a', 'c', 'partially_directed')]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 360 240\" width=\"100%\" font-family=\"Inter, system-ui, sans-serif\"><rect width=\"360\" height=\"240\" fill=\"#fafbfc\"/><line x1=\"148.76\" y1=\"142.17\" x2=\"211.24\" y2=\"97.83\" stroke=\"#d6495b\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"140.60,147.96 151.65,146.25 145.86,138.10\" fill=\"#d6495b\" stroke-linejoin=\"round\"/><polygon points=\"219.40,92.04 208.35,93.75 214.14,101.90\" fill=\"#d6495b\" stroke-linejoin=\"round\"/><line x1=\"153.00\" y1=\"155.48\" x2=\"208.60\" y2=\"155.48\" stroke=\"#7c5cd8\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"143.00,155.48 153.00,160.48 153.00,150.48\" fill=\"#7c5cd8\" stroke-linejoin=\"round\"/><circle cx=\"212.80\" cy=\"155.48\" r=\"4.20\" fill=\"#ffffff\" stroke=\"#7c5cd8\" stroke-width=\"2.0\"/><line x1=\"151.40\" y1=\"84.52\" x2=\"207.00\" y2=\"84.52\" stroke=\"#7c5cd8\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><circle cx=\"147.20\" cy=\"84.52\" r=\"4.20\" fill=\"#ffffff\" stroke=\"#7c5cd8\" stroke-width=\"2.0\"/><polygon points=\"217.00,84.52 207.00,79.52 207.00,89.52\" fill=\"#7c5cd8\" stroke-linejoin=\"round\"/><circle cx=\"130.00\" cy=\"155.48\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"130.00\" y=\"155.48\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">d</text><circle cx=\"130.00\" cy=\"84.52\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"130.00\" y=\"84.52\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">a</text><circle cx=\"230.00\" cy=\"84.52\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"230.00\" y=\"84.52\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">c</text><circle cx=\"230.00\" cy=\"155.48\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"230.00\" y=\"155.48\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">b</text></svg>"
      ],
      "text/plain": [
       "'<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 360 240\" width=\"100%\" font-family=\"Inter, system-ui, sans-serif\"><rect width=\"360\" height=\"240\" fill=\"#fafbfc\"/><line x1=\"148.76\" y1=\"142.17\" x2=\"211.24\" y2=\"97.83\" stroke=\"#d6495b\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"140.60,147.96 151.65,146.25 145.86,138.10\" fill=\"#d6495b\" stroke-linejoin=\"round\"/><polygon points=\"219.40,92.04 208.35,93.75 214.14,101.90\" fill=\"#d6495b\" stroke-linejoin=\"round\"/><line x1=\"153.00\" y1=\"155.48\" x2=\"208.60\" y2=\"155.48\" stroke=\"#7c5cd8\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"143.00,155.48 153.00,160.48 153.00,150.48\" fill=\"#7c5cd8\" stroke-linejoin=\"round\"/><circle cx=\"212.80\" cy=\"155.48\" r=\"4.20\" fill=\"#ffffff\" stroke=\"#7c5cd8\" stroke-width=\"2.0\"/><line x1=\"151.40\" y1=\"84.52\" x2=\"207.00\" y2=\"84.52\" stroke=\"#7c5cd8\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><circle cx=\"147.20\" cy=\"84.52\" r=\"4.20\" fill=\"#ffffff\" stroke=\"#7c5cd8\" stroke-width=\"2.0\"/><polygon points=\"217.00,84.52 207.00,79.52 207.00,89.52\" fill=\"#7c5cd8\" stroke-linejoin=\"round\"/><circle cx=\"130.00\" cy=\"155.48\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"130.00\" y=\"155.48\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">d</text><circle cx=\"130.00\" cy=\"84.52\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"130.00\" y=\"84.52\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">a</text><circle cx=\"230.00\" cy=\"84.52\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"230.00\" y=\"84.52\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">c</text><circle cx=\"230.00\" cy=\"155.48\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"230.00\" y=\"155.48\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">b</text></svg>'"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HIDDEN = [(\"h\", \"c\"), (\"h\", \"d\"), (\"a\", \"c\"), (\"b\", \"d\")]  # h has no column\n",
    "\n",
    "\n",
    "def sample_hidden(seed: int, n: int = 2000, shuffle: bool = True) -> pd.DataFrame:\n",
    "    \"\"\"Draw n rows of a, b, c, d, whose hidden common cause h is left out.\"\"\"\n",
    "    rng = np.random.default_rng(seed)\n",
    "    hidden = rng.standard_normal(n)\n",
    "    a = rng.standard_normal(n)\n",
    "    b = rng.standard_normal(n)\n",
    "    c = 0.8 * a + 0.8 * hidden + rng.standard_normal(n)\n",
    "    d = 0.8 * b + 0.8 * hidden + rng.standard_normal(n)\n",
    "    frame = pd.DataFrame({\"a\": a, \"b\": b, \"c\": c, \"d\": d})\n",
    "    return frame[rng.permutation(frame.columns)] if shuffle else frame\n",
    "\n",
    "\n",
    "gfci_out = andrey.gfci(sample_hidden(seed=3))\n",
    "print(gfci_out.structure.oriented_edges())\n",
    "andrey.viz.draw(gfci_out.structure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constraint_based-12",
   "metadata": {},
   "source": [
    "## CD-NOD\n",
    "\n",
    "CD-NOD takes rows pooled from several domains and a domain index, one value per row. Here the\n",
    "mean of `x` moves between two domains. PC leaves the chain `x - y - z` undirected; the index\n",
    "orients it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "constraint_based-13",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('x', 'y', 'directed'), ('y', 'z', 'directed')]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 360 240\" width=\"100%\" font-family=\"Inter, system-ui, sans-serif\"><rect width=\"360\" height=\"240\" fill=\"#fafbfc\"/><line x1=\"157.00\" y1=\"120.00\" x2=\"93.00\" y2=\"120.00\" stroke=\"#3a4658\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"167.00,120.00 157.00,115.00 157.00,125.00\" fill=\"#3a4658\" stroke-linejoin=\"round\"/><line x1=\"193.00\" y1=\"120.00\" x2=\"257.00\" y2=\"120.00\" stroke=\"#3a4658\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"267.00,120.00 257.00,115.00 257.00,125.00\" fill=\"#3a4658\" stroke-linejoin=\"round\"/><circle cx=\"180.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"180.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">y</text><circle cx=\"80.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"80.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x</text><circle cx=\"280.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"280.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">z</text></svg>"
      ],
      "text/plain": [
       "'<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 360 240\" width=\"100%\" font-family=\"Inter, system-ui, sans-serif\"><rect width=\"360\" height=\"240\" fill=\"#fafbfc\"/><line x1=\"157.00\" y1=\"120.00\" x2=\"93.00\" y2=\"120.00\" stroke=\"#3a4658\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"167.00,120.00 157.00,115.00 157.00,125.00\" fill=\"#3a4658\" stroke-linejoin=\"round\"/><line x1=\"193.00\" y1=\"120.00\" x2=\"257.00\" y2=\"120.00\" stroke=\"#3a4658\" stroke-width=\"2.0\" stroke-linecap=\"round\"/><polygon points=\"267.00,120.00 257.00,115.00 257.00,125.00\" fill=\"#3a4658\" stroke-linejoin=\"round\"/><circle cx=\"180.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"180.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">y</text><circle cx=\"80.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"80.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">x</text><circle cx=\"280.00\" cy=\"120.00\" r=\"13.00\" fill=\"#ffffff\" stroke=\"#182233\" stroke-width=\"1.8\"/><text x=\"280.00\" y=\"120.00\" text-anchor=\"middle\" font-size=\"11\" font-weight=\"600\" fill=\"#182233\" dominant-baseline=\"central\">z</text></svg>'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DOMAINS = [(\"x\", \"y\"), (\"y\", \"z\")]  # the true edges, cause first\n",
    "\n",
    "\n",
    "def sample_domains(\n",
    "    seed: int, n: int = 1000, shuffle: bool = True\n",
    ") -> tuple[pd.DataFrame, np.ndarray]:\n",
    "    \"\"\"Draw n rows from two domains, and the domain of each row; the mean of x moves with it.\"\"\"\n",
    "    rng = np.random.default_rng(seed)\n",
    "    domain = np.repeat([0.0, 1.0], n // 2)  # the first n / 2 rows come from domain 0\n",
    "    x = rng.standard_normal(n) + 1.5 * domain\n",
    "    y = 0.8 * x + rng.standard_normal(n)\n",
    "    z = 0.8 * y + rng.standard_normal(n)\n",
    "    frame = pd.DataFrame({\"x\": x, \"y\": y, \"z\": z})\n",
    "    return (frame[rng.permutation(frame.columns)] if shuffle else frame), domain\n",
    "\n",
    "\n",
    "cdnod_out = andrey.cdnod(*sample_domains(seed=4))\n",
    "print(cdnod_out.structure.oriented_edges())\n",
    "andrey.viz.draw(cdnod_out.structure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constraint_based-9",
   "metadata": {},
   "source": []
  }
 ],
 "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.14.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
