1b_constrained_em_reduced-withD.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy import stats\n",
    "import collections\n",
    "from collections import defaultdict\n",
    "import random\n",
    "import copy\n",
    "import math\n",
    "import itertools\n",
    "from scipy.optimize import linprog\n",
    "from pulp import LpMaximize, LpProblem, LpStatus, lpSum, LpVariable\n",
    "import time\n",
    "import pulp\n",
    "# add paths\n",
    "# path_to_cplex = r'/opt/anaconda3/lib/python3.8/site-packages/pulp/solverdir/cbc/osx/64/cbc'\n",
    "# solver = pulp.CPLEX_CMD(path=path_to_cplex)\n",
    "from collections import Counter\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams.update({'font.size': 14})\n",
    "mpl.rc('text', usetex = False)\n",
    "mpl.rc('font', family = 'serif')\n",
    "import networkx as nx\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Load input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load episode set $E$ with the users that retweeted each original tweet in the trace. \n",
    "\n",
    "Each episode $E_{s}$ includes the users that retweeted s, ordered chronologically, as they appear in the trace. The first user in each episode is the user that originally tweeted the tweet, and is denoted by $r_{s}$. Subsequent users in $E_{s}$ are users that retweeted s, either directly from user $r_{s}$ or from another user that appears in Es before them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "lines = 224667\n",
    "dataset = 'covid'\n",
    "# random.seed(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = pickle.load(open(\"./Extracted/\"+dataset+\"/E\"+ str(lines) + \".p\", \"rb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "for tweet in E:\n",
    "    E[tweet] = list(dict.fromkeys(E[tweet]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load episode set $D$ with the users that retweeted each original tweet in the trace. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = pickle.load(open(\"./Extracted/\"+dataset+\"/D\"+ str(lines) + \".p\", \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the set of original tweets denoted by $S$. \n",
    "\n",
    "The set of original tweets is denoted by $S$, where $S$ = S is the total number of original tweets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = pickle.load(open(\"./Extracted/\"+dataset+\"/S\"+ str(lines) + \".p\", \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load $U$ set xith unique users"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = pickle.load(open(\"./Extracted/\"+dataset+\"/U\"+ str(lines) + \".p\", \"rb\"))\n",
    "U = list(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load $M_{ij}$ variables that count number of episodes where the ordered pair (i,j) appears"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = pickle.load(open(\"./Extracted/\"+dataset+\"/M_d\"+ str(lines) + \".p\", \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Find important quantities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Number of unique users N\n",
    "- Number of episodes / original tweets S\n",
    "- Number of active pairs (i,j)\n",
    "\n",
    "For every active pair we have to count the $\\sigma_{ij}$'s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique users N = 44602\n"
     ]
    }
   ],
   "source": [
    "N = len(U)\n",
    "print('Number of unique users N =', N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Episodes (original tweets) S = 9722\n"
     ]
    }
   ],
   "source": [
    "print('Number of Episodes (original tweets) S =',len(E))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "active_pairs_n = 0 \n",
    "for i in M:\n",
    "    active_pairs_n+=len(M[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of active user pairs in the trace: 41020780 out of the 1989293802 possible pairs\n"
     ]
    }
   ],
   "source": [
    "print('Number of active user pairs in the trace:', active_pairs_n, 'out of the', N*(N-1), 'possible pairs')\n",
    "\n",
    "# 50000 lines: Number of active user pairs in the trace: 415297 out of the 47837972 possible pairs\n",
    "# 70000 lines: Number of active user pairs in the trace: 818995 out of the 90069590 possible pairs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. CONSTRAINED-EM Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define important functions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_check(E, Q):\n",
    "    '''\n",
    "    Function that checks feasibility of results\n",
    "\n",
    "    '''\n",
    "    retweets = 0 # minimum existing edges\n",
    "    total_feasible_edges = []\n",
    "    total_inf = 0 \n",
    "    for s in E:\n",
    "        feasible_edges = 0 \n",
    "        for j in E[s]:\n",
    "            indx = E[s].index(j)\n",
    "            if indx!=0:\n",
    "                u_before = E[s][:indx]\n",
    "                for i in u_before: \n",
    "                    if j in Q[i] and Q[i][j] > 0.5:\n",
    "                        feasible_edges +=1\n",
    "                        total_feasible_edges.append((i,j))\n",
    "\n",
    "        infeasible = (len(E[s]) - 1) - feasible_edges\n",
    "        if infeasible > 0:\n",
    "            print('Tweet', s, 'retweeted by', len(E[s])-1, 'users in total. But, we only found:', feasible_edges, 'feasible edges, so the infeasible ones are:', infeasible)            \n",
    "            total_inf+=infeasible\n",
    "            \n",
    "        retweets += len(E[s])-1\n",
    "        total_feasible_edges = list(set(total_feasible_edges))\n",
    "        \n",
    "    print('Total feasbile edges:', len(total_feasible_edges), 'Number of retweets', retweets)\n",
    "    print('Total infeasible edges:', total_inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flatten(obj):\n",
    "    if type(obj) == list:\n",
    "        return [l for L in obj for l in L]\n",
    "    if type(obj) == dict:\n",
    "        return [l for i in obj for l in obj[i].values()]\n",
    "    if type(obj) == defaultdict:\n",
    "        return [l for i in obj for l in obj[i].values()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_Q(M, Q, s, a, b, r):\n",
    "    for i in M:\n",
    "        for j in M[i]:\n",
    "            e = M[i][j] * s[i][j]\n",
    "            numer = r * (a**e) * ((1-a)**(M[i][j]-e))\n",
    "            denom = r * (a**e) * ((1-a)**(M[i][j]-e)) + (1-r) * (b**e) * ((1-b)**(M[i][j]-e))\n",
    "            Q[i][j] = numer/denom   \n",
    "    return Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_a(M, Q, s):\n",
    "    a = ((sum([M[i][j] * Q[i][j] * s[i][j]  for i in M for j in M[i]]))\n",
    "        /((sum([M[i][j] * Q[i][j] for i in M for j in M[i]])))) \n",
    "    return(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_b(M, Q, s):\n",
    "    b = ((sum([M[i][j] * (1-Q[i][j]) * s[i][j]  for i in M for j in M[i]]))\n",
    "        /((sum([M[i][j] * (1-Q[i][j]) for i in M for j in M[i]]))))\n",
    "    return(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_r(Q, N, active_pairs_n):\n",
    "#     r = (sum([Q[i][j] for i in Q for j in Q[i]]))/(N*(N-1)) \n",
    "    r = (sum([Q[i][j] for i in Q for j in Q[i]]))/(active_pairs_n)      \n",
    "    return(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_s(model, M, Q, x, a, b, s, active_pairs, lam):\n",
    "    \"\"\" \n",
    "    This function updates the s_{ij} parameters of the optimization problemm\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    M : dict\n",
    "        Dictionary with M_{ij} values tha show how many times an ordered pair (i,j) appears in the trace.\n",
    "    Q : dict\n",
    "        Dictionary with the Q_{ij} values\n",
    "    W : dict\n",
    "        Dictionary with the coefficients of the problem \n",
    "    x : dict\n",
    "        Dictionary with decision variables\n",
    "    active_pairs : list\n",
    "        Active pairs of the problem\n",
    "    s : dict\n",
    "        Dictionary with the updated s_{ij} parameters\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    s : dict\n",
    "        Dictionary with the updated s_{ij} parameters\n",
    "    \"\"\"\n",
    "    \n",
    "    W = defaultdict(dict)\n",
    "#     W_ones = defaultdict(dict)\n",
    "\n",
    "#     start = time.time()\n",
    "    for pair in active_pairs:\n",
    "        i = pair[0]\n",
    "        j = pair[1]\n",
    "        W[i][j] = M[i][j]*((Q[i][j]*math.log(a/(1-a)))+ (1-Q[i][j])*math.log(b/(1-b))) + random.uniform(0,0.001)\n",
    "\n",
    "#     minv = min((set(min(list((i.values())) for i in W.values()))))\n",
    "#     print('min W:', minv)\n",
    "    \n",
    "        \n",
    "#     end = time.time()\n",
    "#     print('(s) W:', end-start)\n",
    "#     for i in M:\n",
    "#         for j in M[i]:\n",
    "#             if (i,j) in active_pairs:\n",
    "#                 W[i][j] = M[i][j]*((Q[i][j]*math.log(a/(1-a)))+ (1-Q[i][j])*math.log(b/(1-b))) + random.uniform(0,0.001)\n",
    "#             else: # s[i][j] = 1\n",
    "#                 W_ones[i][j] = M[i][j]*((Q[i][j]*math.log(a/(1-a)))+ (1-Q[i][j])*math.log(b/(1-b))) + random.uniform(0,0.001)\n",
    "    \n",
    "#     start = time.time()\n",
    "    s = pulp_solve(model, active_pairs, W, x, s, lam)\n",
    "#     end = time.time()\n",
    "#     print('(s) pulp_solve:', end-start)\n",
    "    return(s, W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pulp_create_reduct(E, M):\n",
    "    \"\"\" \n",
    "    This function initializes the optimization problem with reduced constraints. Runs the first time\n",
    "\n",
    "    Returns\n",
    "    ----------\n",
    "    model : LpProblem\n",
    "        Pulp model of the problem \n",
    "    x : Dict\n",
    "        Dictionary with decision variables\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    # Initialize the maximization problem\n",
    "    model = LpProblem(name=\"constr-newman\", sense=LpMaximize)\n",
    "    \n",
    "    # For each episode in $E$, each line in constraints_list includes:\n",
    "    # - for each user j in $E$, the users before them in pos 0 and \n",
    "    # - the user itself in pos 1\n",
    "    \n",
    "    sij = defaultdict(dict)\n",
    "    constraints_list = []\n",
    "    \n",
    "    for s in D:\n",
    "        for users_list in list(D[s].values()):\n",
    "            index_now = list(D[s].values()).index(users_list)\n",
    "            if index_now == 0:\n",
    "                for u in list(D[s].values())[1]:\n",
    "                    sij[D[s][0][0]][u] = 1\n",
    "            else:\n",
    "                for u in users_list:\n",
    "                    u_before_l = list(D[s].values())[:index_now]\n",
    "                    u_before = [item for sublist in u_before_l for item in sublist]\n",
    "                    if u in u_before: u_before.remove(u)\n",
    "                    constraints_list.append([u_before, u])\n",
    "\n",
    "    len_b = len(constraints_list)\n",
    "    print('constraints before', len(constraints_list))\n",
    "\n",
    "    # phase I: delete all constraints that include a pair with sij = 1\n",
    "    \n",
    "    for c in list(constraints_list):\n",
    "        j = c[1]\n",
    "        for i in c[0]:\n",
    "            if i in sij and j in sij[i] and sij[i][j]==1:\n",
    "                constraints_list.remove(c)\n",
    "                break\n",
    "                \n",
    "    # For each episode in $E$, the constraints_dictionary includes:\n",
    "    # - the user j as keys\n",
    "    # - the users that come before them for each constraint that they appear in pos j as values\n",
    "    \n",
    "    constraints_list.sort()\n",
    "    constraints_list = list(constraints_list for constraints_list,_ in itertools.groupby(constraints_list))\n",
    "    \n",
    "#     print('3rd:', [[924604094999908352, 952286580890365952, 367664497],787929838774517760] in constraints_list)\n",
    "    \n",
    "    constraints_dict = dict()\n",
    "    ind = 0 \n",
    "    for c in constraints_list:\n",
    "            j = c[1]\n",
    "            if j in constraints_dict:\n",
    "                constraints_dict[j].append([c[0],ind])\n",
    "            else:\n",
    "                constraints_dict[j] = []\n",
    "                constraints_dict[j].append([c[0],ind])\n",
    "            ind+=1\n",
    "\n",
    "    # 2nd phase ----------------------------\n",
    "#     print('3rdb:', [[924604094999908352, 952286580890365952, 367664497],787929838774517760] in constraints_list)\n",
    "\n",
    "    \n",
    "    for j in constraints_dict:\n",
    "        for constraint1 in constraints_dict[j]:\n",
    "            for constraint2 in constraints_dict[j]:\n",
    "                if constraint1[1]!=constraint2[1]:\n",
    "                    if set(constraint1[0]).issubset(constraint2[0]):\n",
    "                        if constraint1[0]!=constraint2[0] and collections.Counter(constraint1[0]) != collections.Counter(constraint2[0]): \n",
    "                            if [constraint2[0],j] in constraints_list: # we may have already deleted it\n",
    "                                constraints_list.remove([constraint2[0],j])\n",
    "                    \n",
    "#     print('4th:', [[924604094999908352, 952286580890365952, 367664497],787929838774517760] in constraints_list)\n",
    "                            \n",
    "\n",
    "    len_a = len(constraints_list)          \n",
    "    \n",
    "    print('constraints after', len(constraints_list))\n",
    "    print('change of constraints:', abs(len_a - len_b))\n",
    "    \n",
    "    # Add decision variables \n",
    "    \n",
    "    # Create active pairs and sij's randomly \n",
    "    active_pairs = []\n",
    "    for i in M:\n",
    "        for j in M[i]:\n",
    "            if i not in sij or j not in sij[i]:\n",
    "                active_pairs.append((i,j))\n",
    "                sij[i][j] = random.uniform(0, 1)\n",
    "\n",
    "    active_pairs = list(set(active_pairs))\n",
    "        \n",
    "    x = dict()\n",
    "    x = {f\"{pair[0]}-{pair[1]}\":LpVariable(name=f\"x{pair[0]}-{pair[1]}\", lowBound=0, upBound=0.9999999999999) for pair in active_pairs}\n",
    "    \n",
    "    # Add constraints in model\n",
    "    for c in constraints_list: \n",
    "        j = c[1]\n",
    "        tmp = 0\n",
    "        for i in c[0]:\n",
    "            # x[k] = LpVariable(name=f\"x{u}-{j}\", lowBound=0, upBound=1)\n",
    "            tmp+=x[f\"{i}-{j}\"]\n",
    "        model += (tmp >= 1)\n",
    "        \n",
    "    print('Active pairs:', len(active_pairs))\n",
    "        \n",
    "    return model, x, active_pairs, sij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pulp_solve(model, active_pairs, W, x, s, lam):\n",
    "    \"\"\" \n",
    "    This function solves the optimization problem with PULP.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    model : PulpModel\n",
    "        Model initialized from pulp\n",
    "    active_pairs : list\n",
    "        Active pairs of the problem\n",
    "    W : dict\n",
    "        Dictionary with the coefficients of the problem \n",
    "    x : dict\n",
    "        Dictionary with decision variables\n",
    "    lam : float\n",
    "        Lambda value by which we decreace the W's\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    s : dict\n",
    "        Dictionary with found s_{ij} parameters\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    maxv = max((set(max(list(i.values()) for i in W.values()))))\n",
    "    c = maxv + 0.0001\n",
    "    \n",
    "    # Update the objective function\n",
    "    model += lpSum((W[i][j]-c)*x[f\"{i}-{j}\"] for i in W for j in W[i])\n",
    "#     print(model)\n",
    "#     start = time.time()\n",
    "    # And solve it \n",
    "    status = model.solve()\n",
    "#     end = time.time()\n",
    "#     print('(s) model.solve:', end-start)\n",
    "    # print(f\"status: {model.status}, {LpStatus[model.status]}\")\n",
    "\n",
    "#     start = time.time()\n",
    "    for var in x.values():\n",
    "        i = var.name.split(\"_\")[0]\n",
    "        i = int(i.split(\"x\")[1])\n",
    "        j = int(var.name.split(\"_\")[1])\n",
    "        s[i][j] = var.value()\n",
    "#     end = time.time()\n",
    "#     print('(s) assign to sij:', end-start)        \n",
    "#     print('status:', LpStatus[model.status])\n",
    "\n",
    "#     costp = pulp.value(model.objective)\n",
    "#     print('Max objective function:', costp)\n",
    "    \n",
    "    return s\n",
    "                    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newman(eps, N, M, active_pairs_n, lam):\n",
    "        \"\"\" \n",
    "        This function is the main algorithm for path inference with constraints. \n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        rep : int \n",
    "            Number of times to run the algorithm.\n",
    "        N : int\n",
    "            Number of users in the trace.\n",
    "        M : dict\n",
    "            Dictionary with M_{ij} values tha show how many times an ordered pair (i,j) appears in the trace.\n",
    "\n",
    "        Returns\n",
    "        ----------\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        iterat = 1\n",
    "        # ======================== INITIALIZE ========================\n",
    "        print('Initialize...')\n",
    "\n",
    "        # 0.1 -- Initialize Q dictionary \n",
    "\n",
    "        Q = defaultdict(dict)\n",
    "\n",
    "        # 0.2 -- Initialize a, b, r $ parameters\n",
    "        a = random.uniform(0.5, 1)\n",
    "        b = random.uniform(0, 0.5)\n",
    "        r = random.uniform(0, 1)\n",
    "\n",
    "#         a= 0.9999 \n",
    "#         b= 0.0001 \n",
    "#         r= 0.003\n",
    "\n",
    "        print('Initial a=', a,'b=', b, 'r=', r)\n",
    "\n",
    "        print('Create pulp model with reduced constraints and variables...')\n",
    "        # Create pulp model, one time\n",
    "        model, x, active_pairs, s = pulp_create_reduct(E, M)\n",
    "        \n",
    "        # ======================== START algorithm ========================\n",
    "\n",
    "        while True: # Repeat until convergence\n",
    "            # Step 1 ==== UPDATE VALUES ====\n",
    "\n",
    "            # 1.1 -- Update Q dictionary\n",
    "            \n",
    "            # print('Update Q...')\n",
    "#             start = time.time()\n",
    "            Q = update_Q(M, Q, s, a, b, r)\n",
    "#             end = time.time()\n",
    "#             print('Q', start-end)\n",
    "            # 1.2 -- Update a,b,r,s_{ij} parameters according to Equations\n",
    "            \n",
    "            # print('Update a...')\n",
    "#             start = time.time()\n",
    "            a = update_a(M, Q, s)\n",
    "#             end = time.time()\n",
    "#             print('a', start-end)\n",
    "            \n",
    "            # print('Update b...')\n",
    "#             start = time.time()\n",
    "            b = update_b(M, Q, s)\n",
    "#             end = time.time()\n",
    "#             print('b', start-end)\n",
    "            \n",
    "            # print('Update r...')\n",
    "#             start = time.time()\n",
    "            r = update_r(Q, N, active_pairs_n)\n",
    "#             end = time.time()\n",
    "#             print('r:', start-end)            \n",
    "            if a ==1: \n",
    "                a = 0.9999999\n",
    "                flag = True\n",
    "            if b ==0: \n",
    "                b = 0.0001\n",
    "                flag = True\n",
    "            \n",
    "            # print('Update s...')\n",
    "#             start = time.time()\n",
    "            s, W  = update_s(model, M, Q, x, a, b, s, active_pairs, lam)\n",
    "#             end = time.time()\n",
    "#             print('s:', start-end)\n",
    "            # Step 2 ==== CHECK CONVERGENCE ====\n",
    "            \n",
    "            # print('Check convergence...')\n",
    "\n",
    "            if iterat > 1:\n",
    "                new_q = np.array(flatten(Q))\n",
    "                new_a = a\n",
    "                new_b = b\n",
    "                \n",
    "                cost = sum(s[i][j]*M[i][j]*(Q[i][j]*math.log(a/(1-a))+(1-Q[i][j])*math.log(b/(1-b))) for i in s for j in s[i])\n",
    "                print('Objective cost:', cost)\n",
    "                new_cost = cost\n",
    "                \n",
    "                changeq = np.linalg.norm(new_q - old_q)  \n",
    "                changea = abs(new_a-old_a)\n",
    "                changeb = abs(new_b-old_b)\n",
    "                changecost = abs(new_cost - old_cost)\n",
    "\n",
    "#                 if changeq < eps and changea < eps and changeb < eps or flag == True: \n",
    "#                 if changeq < eps or flag == True:\n",
    "                if changeq<eps:\n",
    "                    print('-----')\n",
    "#                     print('flag', flag)\n",
    "                    print('- Change q:', changeq, 'Change a:', changea, 'Change b:', changeb)\n",
    "                    print('- Changecost:', changecost)\n",
    "                    print('- a=', a,'b=',b,'r=', r)\n",
    "                    print('cost:', cost)\n",
    "                \n",
    "                    Q_old = copy.deepcopy(Q)\n",
    "                    old_q = np.array(flatten(Q_old))\n",
    "                    old_a = a\n",
    "                    old_b = b\n",
    "                    old_cost = cost        \n",
    "                    break\n",
    "                else: \n",
    "                    print('-----')\n",
    "                    print('- Change q:', changeq, 'Change a:', changea, 'Change b:', changeb)\n",
    "                    print('- Changecost:', changecost)\n",
    "                    print('- a=', a,'b=',b,'r=', r)\n",
    "                    print('cost:', cost)\n",
    "                \n",
    "                    Q_old = copy.deepcopy(Q)\n",
    "                    old_q = np.array(flatten(Q_old))\n",
    "                    old_a = a\n",
    "                    old_b = b\n",
    "                    old_cost = cost\n",
    "                    \n",
    "                    iterat += 1\n",
    "                    \n",
    "            if iterat == 1:\n",
    "                flag = False\n",
    "                old_q = np.array(flatten(Q))\n",
    "                old_a = a\n",
    "                old_b = b\n",
    "                \n",
    "                cost = sum(s[i][j]*M[i][j]*(Q[i][j]*math.log(a/(1-a))+(1-Q[i][j])*math.log(b/(1-b))) for i in s for j in s[i])\n",
    "#                 print('Objective cost:', cost)\n",
    "                old_cost = cost\n",
    "                \n",
    "                iterat+=1\n",
    "                \n",
    "        # ======================== END algorithm ========================\n",
    "    \n",
    "        return iterat, W, a, b, r, s, Q, cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------SEED: 1197612772160344064\n",
      "Initialize...\n",
      "Initial a= 0.7857012973449568 b= 0.2144445273375573 r= 0.5780913011344704\n",
      "Create pulp model with reduced constraints and variables...\n",
      "constraints before 60675\n",
      "constraints after 49916\n",
      "change of constraints: 10759\n",
      "Active pairs: 41011483\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "random.seed(10)\n",
    "print('------------------SEED:', i)\n",
    "eps = 0.001\n",
    "lam = 1\n",
    "iterat, W, a, b, r, s, Q, cost = newman(eps, N, M, active_pairs_n, lam)\n",
    "end = time.time()\n",
    "print('Time:', end-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "496"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iterat "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'Q_30' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-103-9c6017e8af3e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0msame\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mn_same\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mQ_30\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mQ_30\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mQ_40\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mQ_40\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'Q_30' is not defined"
     ]
    }
   ],
   "source": [
    "# np.linalg.norm(np.array(flatten(Q_10))-np.array(flatten(Q_20)))\n",
    "exists = 0\n",
    "n_exists = 0\n",
    "same = 0 \n",
    "n_same = 0\n",
    "for i in Q_30:\n",
    "    for j in Q_30[i]:\n",
    "        if i in Q_40 and j in Q_40[i]: \n",
    "            exists+=1\n",
    "            if Q_30[i][j] > 0.9:\n",
    "                if Q_40[i][j] > 0.9:\n",
    "                    same +=1\n",
    "                else:\n",
    "                    n_same +=1\n",
    "            else:\n",
    "                if Q_40[i][j] > 0.9:\n",
    "                    n_same +=1\n",
    "                else:\n",
    "                    same +=1\n",
    "        else: \n",
    "            n_exists+=1\n",
    "exists, n_exists, same, n_same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- a= 0.9996750478190183 1-a 0.0003249521809817191\n",
      "- b= 3.7354632604484254e-06 1-b: 0.9999962645367395\n",
      "- r= 0.02690363162162662\n",
      "0.9996750478190183 3.7354632604484254e-06 0.02690363162162662\n",
      "Feasibility check:\n",
      "Total feasbile edges: 11140 Number of retweets 11669\n",
      "Total infeasible edges: 0\n"
     ]
    }
   ],
   "source": [
    "print('- a=', a, '1-a', 1-a)\n",
    "print('- b=', b, '1-b:', 1-b)\n",
    "print('- r=', r)\n",
    "print(a,b,r)\n",
    "print('Feasibility check:')\n",
    "f_check(E,Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "try: \n",
    "    pickle.dump(Q, open(\"./Extracted/AlgorithmResults/Q_constrained_\" + str(lines) + \".p\", \"wb\"))\n",
    "except: \n",
    "    print(\"Unable to write to file\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "try: \n",
    "    pickle.dump(s, open(\"./Extracted/AlgorithmResults/s_constrained_\" + str(lines) + \".p\", \"wb\"))\n",
    "except: \n",
    "    print(\"Unable to write to file\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list((set(max(list(i.values()) for i in W.values()))))[0]\n",
    "# max((set(max(list(i.values()) for i in W.values()))))\n",
    "\n",
    "# maxv = max((set(max(list(i.values()) for i in M.values()))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}