{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "\n",
    "# Ring of polynomials for the weights of the base\n",
    "R.<k> = PolynomialRing(QQ)\n",
    "\n",
    "# Ring of rational functions to express PageRank analytically\n",
    "S.<a> = PolynomialRing(R)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sturm(p):\n",
    "    \"\"\"\n",
    "        Computes the Sturm polynomials.\n",
    "        \n",
    "        Args:\n",
    "            p: a polynomial in a with rational coefficients in k\n",
    "            \n",
    "        Returns:\n",
    "            the Sturm polynomials of p.\n",
    "    \"\"\"\n",
    "    S = []\n",
    "    S.append(p)\n",
    "    S.append(diff(p, a))\n",
    "    while True:\n",
    "        next_p = -S[-2].mod(S[-1])\n",
    "        if next_p == 0:\n",
    "            break\n",
    "        S.append(next_p)\n",
    "    return S\n",
    "\n",
    "# Maybe there's a better way, but I couldn't find it\n",
    "\n",
    "ge_op = (x>=0).operator()\n",
    "gt_op = (x>0).operator()\n",
    "\n",
    "def changes(S, *points):\n",
    "    \"\"\"\n",
    "        Computes changes in the Sturm polynomials at given points.\n",
    "        \n",
    "        Args:\n",
    "            S: the Sturm polynomials of a polynomial.\n",
    "            \n",
    "        Returns:\n",
    "            a dictionary mapping each point to the number of sign changes and the first valid value k for such changes.\n",
    "    \"\"\"\n",
    "\n",
    "    changes={}\n",
    "    for point in points:\n",
    "        prev_op = None\n",
    "        c = 0\n",
    "        m = 0\n",
    "        for i in range(len(S)):\n",
    "            # This should be fixed to handle the case of integer points, which corresponds to zeros\n",
    "            sol = solve(S[i](a=point)(k=x)>0,x)\n",
    "            op = '+'\n",
    "            if len(sol) > 0:\n",
    "                ineq = sol[-1][-1]\n",
    "                op = '+' if ineq.operator() == ge_op or ineq.operator() == gt_op else '-'\n",
    "                m = max(m, ineq.rhs())\n",
    "            if i > 0 and op != prev_op:\n",
    "                c += 1\n",
    "            prev_op = op\n",
    "        \n",
    "        # Now we try to improve m manually\n",
    "            \n",
    "        for r in range(math.ceil(m) - 1, 0, -1):\n",
    "            prev_op = '0'\n",
    "            d = 0\n",
    "            for i in range(len(S)):\n",
    "                if S[i].denominator()(a=point)(k=r) == 0:\n",
    "                    d = -1\n",
    "                    break\n",
    "                v = S[i](a=point)(k=r)\n",
    "                if v == 0:\n",
    "                    op = '0'\n",
    "                elif v > 0:\n",
    "                    op = '+'\n",
    "                else:\n",
    "                    op = '-'\n",
    "                if op != prev_op and op != '0' and prev_op != '0':\n",
    "                    d += 1\n",
    "                prev_op = op\n",
    "                \n",
    "            if d == c:\n",
    "                m = m -1\n",
    "            else:\n",
    "                break\n",
    "        \n",
    "        changes[point] = [c, m]\n",
    "        \n",
    "    return changes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Score and rank monotonicity counterexample\n",
    "# Node names correspond to the labeling in the paper\n",
    "\n",
    "# 0: k-clique A (k-1 elements)\n",
    "# 1: special node of the k-clique A attached to two paths (nodes 2 and 4)\n",
    "# 2: first node of the left path [*]\n",
    "# 3: second node of the left path\n",
    "# 4: first node of the right path\n",
    "# 5: second node of the right path [*]\n",
    "# 6: special node of the k-clique B attached to 3\n",
    "# 7: special node of the k-clique B attached to 5 ONLY AFTER\n",
    "# 8: k-clique B (k-2 elements)\n",
    "\n",
    "A_pre  = matrix(9, 9, [\n",
    "    (k-2)/(k-1), (k - 1)/(k-1), 0, 0, 0, 0, 0, 0, 0, #0: 0+1\n",
    "    1/(k+1), 0, 1/(k+1), 0, 1/(k+1), 0, 0, 0, 0, #1: 0+2+4\n",
    "    0,   1/2, 0, 1/2, 0, 0, 0, 0, 0, #2: 1+3\n",
    "    0,   0, 1/2, 0, 0, 0, 1/2, 0, 0, #3: 2+6\n",
    "    0,   1/2, 0, 0, 0, 1/2, 0, 0, 0, #4: 1+5\n",
    "    0,   0, 0, 0, 1, 0, 0, 0, 0, #5: 4\n",
    "    0,   0, 0, 1/k, 0, 0, 0, 1/k, 1/k, #6: 3+7+8\n",
    "    0,   0, 0, 0, 0, 0, 1/(k-1), 0, 1/(k-1), #7: 6+8\n",
    "    0,   0, 0, 0, 0, 0, (k - 2)/(k-1), (k -2)/(k-1), (k-3)/(k-1), #8: 6+7+8\n",
    "])\n",
    "\n",
    "A_post  = matrix(9, 9, [\n",
    "    (k-2)/(k-1), (k - 1)/(k-1), 0, 0, 0, 0, 0, 0, 0, #0: 0+1\n",
    "    1/(k+1), 0, 1/(k+1), 0, 1/(k+1), 0, 0, 0, 0, #1: 0+2+4\n",
    "    0,   1/2, 0, 1/2, 0, 0, 0, 0, 0, #2: 1+3\n",
    "    0,   0, 1/2, 0, 0, 0, 1/2, 0, 0, #3: 2+6\n",
    "    0,   1/2, 0, 0, 0, 1/2, 0, 0, 0, #4: 1+5\n",
    "    0,   0, 0, 0, 1/2, 0, 0, 1/2, 0, #5: 4+7\n",
    "    0,   0, 0, 1/k, 0, 0, 0, 1/k, 1/k, #6: 3+7+8\n",
    "    0,   0, 0, 0, 0, 1/k, 1/k, 0, 1/k, #7: 5+6+8\n",
    "    0,   0, 0, 0, 0, 0, (k - 2)/(k-1), (k -2)/(k-1), (k-3)/(k-1), #8: 6+7+8\n",
    "])\n",
    "\n",
    "\n",
    "# Initialization vector from the total graph\n",
    "\n",
    "r_pre = (1-a) * vector([1/(k - 1 + 7 + k - 2)]*9) * ~(identity_matrix(9) - a * A_pre)\n",
    "r_post = (1-a) * vector([1/(k - 1 + 7 + k - 2)]*9) * ~(identity_matrix(9) - a * A_post)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Score Counterexample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, we show that the score delta of node 5 at ⍺ = 2/3 is negative for a large enough $k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[x > -2, x < 0.08162966219681982], [x > 11.97222985977454]]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve((r_post[5]-r_pre[5])(a=2/3)(k=x)<0,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we show that its numerator has exactly two roots in (0..1] using Sturm's theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of roots in [0..1] 2 from 11.99235993208829\n"
     ]
    }
   ],
   "source": [
    "p = (r_post[5]-r_pre[5]).numerator()\n",
    "\n",
    "St = sturm(p)\n",
    "result = changes(St, 0, 1)\n",
    "print(\"Number of roots in [0..1]\", result[0][0] - result[1][0], \"from\", max(result[0][1], result[1][1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "Finally, we sandwich 2/3 between two points between which there are no sign changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of roots in [0..375/2/(k + 250)]: 1 from 12.9156626506024\n",
      "375/2/(k + 250) goes to 0 for k → ∞\n",
      "Number of roots in [(k + 250)/(k + 500)..1]: 1 from 12.928571428565192\n",
      "(k + 250)/(k + 500) goes to 1 for k → ∞\n"
     ]
    }
   ],
   "source": [
    "e = 3/4 - 3*k/(4*k+ 1000)\n",
    "result = changes(St, 0, e)\n",
    "print(\"Number of roots in [0..\", e, \"]: \", result[0][0] - result[e][0], \" from \", max(result[0][1], result[e][1]), sep = \"\")\n",
    "print(e, \"goes to\", limit(e, k=infinity), \"for k → ∞\")\n",
    "\n",
    "\n",
    "e = 1/2 + k/(2*k+1000)\n",
    "result = changes(St, e, 1)\n",
    "print(\"Number of roots in [\", e, \"..1]: \", result[e][0] - result[1][0], \" from \", max(result[e][1], result[1][1]), sep = \"\")\n",
    "print(e, \"goes to\", limit(e, k=infinity), \"for k → ∞\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rank Counterexample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PRE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, we show that pre delta value of nodes 5 and 2  at ⍺ = 2/3 is positive for a large enough $k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[x > -2, x < 0.1411419796174787], [x > 13.26981707317073]]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve((r_pre[5]-r_pre[2])(a=2/3)(k=x)>0,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we show that its numerator has exactly two roots in (0..1] using Sturm's theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of roots in [0..1]: 2 from 13.26593406593407\n"
     ]
    }
   ],
   "source": [
    "p = (r_pre[5]-r_pre[2]).numerator()\n",
    "\n",
    "St = sturm(p)\n",
    "result = changes(St, 0, 1)\n",
    "print(\"Number of roots in [0..1]:\", result[0][0] - result[1][0], \"from\", max(result[0][1], result[1][1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we sandwich 2/3 between two points between which there are no sign changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of roots in [0..75/2/(k + 50)]: 1 from 13.26593406593407\n",
      "75/2/(k + 50) goes to 0 for k → ∞\n",
      "\n",
      "Number of roots in [(k + 50)/(k + 100)..1]: 1 from 13.26593406593407\n",
      "(k + 50)/(k + 100) goes to 1 for k → ∞\n"
     ]
    }
   ],
   "source": [
    "e = 3/4 - 3*k/(4*k+200)\n",
    "result = changes(St, 0, e)\n",
    "print(\"Number of roots in [0..\", e, \"]: \", result[0][0] - result[e][0], \" from \", max(result[0][1], result[e][1]), sep = \"\")\n",
    "print(e, \"goes to\", limit(e, k=infinity), \"for k → ∞\")\n",
    "\n",
    "print()\n",
    "\n",
    "e = 1/2 + k/(2*k+200)\n",
    "result = changes(St, e, 1)\n",
    "print(\"Number of roots in [\", e, \"..1]: \", result[e][0] - result[1][0], \" from \", max(result[e][1], result[1][1]), sep = \"\")\n",
    "print(e, \"goes to\", limit(e, k=infinity), \"for k → ∞\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## POST\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, we check that the post delta value of nodes 5 and 2 at ⍺ = 2/3 is negative for a large enough $k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[x < -2], [x > (1/3), x < 1], [x > 5]]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve((r_post[5]-r_post[2])(a=2/3)(k=x)<0, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we show that its numerator has exactly three roots in (0..1] using Sturm's theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of roots in [0..1]: 3 from sqrt(5) + 4\n"
     ]
    }
   ],
   "source": [
    "p = (r_post[5]-r_post[2]).numerator()\n",
    "\n",
    "St = sturm(p)\n",
    "result = changes(St, 0, 1)\n",
    "print(\"Number of roots in [0..1]:\", result[0][0] - result[1][0], \"from\", max(result[0][1], result[1][1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we sandwich 2/3 between two points between which there are no sign changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of roots in [0..20/(k + 200)]: 2 from 220/9\n",
      "20/(k + 200) goes to 0 for k → ∞\n",
      "\n",
      "Number of roots in [(k + 50)/(k + 100)..1]: 1 from 5.07389428263215\n",
      "(k + 50)/(k + 100) goes to 1 for k → ∞\n"
     ]
    }
   ],
   "source": [
    "e = 1/10 - 1*k/(10*k+2000)\n",
    "result = changes(St, 0, e)\n",
    "print(\"Number of roots in [0..\", e, \"]: \", result[0][0] - result[e][0], \" from \", max(result[0][1], result[e][1]), sep = \"\")\n",
    "print(e, \"goes to\", limit(e, k=infinity), \"for k → ∞\")\n",
    "\n",
    "print()\n",
    "\n",
    "e = 1/2 + k/(2*k+200)\n",
    "result = changes(St, e, 1)\n",
    "print(\"Number of roots in [\", e, \"..1]: \", result[e][0] - result[1][0], \" from \", max(result[e][1], result[1][1]), sep = \"\")\n",
    "print(e, \"goes to\", limit(e, k=infinity), \"for k → ∞\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.3",
   "language": "sage",
   "name": "sagemath"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}