XOR model for computing the frustration index.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binary linear programming model (XOR model) for computing the frustration index of a signed graph\n",
    "\n",
    "Jupyter code written by Samin Aref in 2016\n",
    "\n",
    " Creative common licence: \n",
    " Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)\n",
    "\n",
    "Using this code for non-commercial purposes is permitted to all given that the two following publications are cited:\n",
    "\n",
    "    Aref, S., Mason, A. J., and Wilson, M. C. A Modeling and Computational Study of the Frustration Index in Signed Networks. Networks (2019). url: https://onlinelibrary.wiley.com/doi/full/10.1002/net.21907. doi: 10.1002/net.21907\n",
    "\n",
    "    Aref, S., Mason, A. J., and Wilson, M. C. Computing the line index of balance using integer programming optimisation. In Optimization Problems in Graph Theory, B. Goldengorin, Ed. Springer, 2018, pp. 65-84. url: https://www.springer.com/gp/book/9783319948294. doi: 10.1007/978-3-319-94830-0_3\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suggested articles:\n",
    "\n",
    "    Aref, S. Signed Network Structural Analysis and Applications with a Focus on Balance Theory. (Ph.D. thesis) University of Auckland (2019).\n",
    "    \n",
    "    Aref, S., and Wilson, M. C. Balance and frustration in signed networks. Journal of Complex Networks 7, 2 (2019), 163–189. doi: 10.1093/comnet/cny015.\n",
    "    \n",
    "    Aref, S., and Wilson, M. C. Measuring partial balance in signed networks. Journal of Complex Networks 6, 4 (2018), 566-595. doi: 10.1093/comnet/cnx044. \n",
    "    \n",
    "    Aref, S., and Neal, Z. Legislative effectiveness hangs in the balance: Studying balance and polarization through partitioning signed networks. arXiv:1906.01696 (2019). url: http://arxiv.org/pdf/1906.01696.\n",
    "    \n",
    "Related dataset: \n",
    "\n",
    "    Aref, S. Signed networks from sociology and political science, systems biology, international relations, finance, and computational chemistry. Figshare research data repository (2017). doi: 10.6084/m9.figshare.5700832.v2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "How many random graphs do you want to create?(1/2/...)2\n",
      "\n",
      " 1.  Erdős-Rényi binomial graph G_{n,p} \n",
      " 2.  Random graph of n nodes and m total edges G_{n,m} \n",
      " 3.  Newman-Watts-Strogatz small world graph (adding new edges to the ring) G_{n,k,p} \n",
      " 4.  Watts-Strogatz small world graph (rewiring the ring) G_{n,k,p} \n",
      " 5.  Random regular graph of n nodes each with degree d G_{d,n} \n",
      " 6.  Barabási-Albert preferential attachment model graph G_{n,m} \n",
      " 7.  Growing graph with powerlaw degree distribution and approximate average clustering G_{n,m,p}\n",
      " 8.  Two-dimensional grid n*n \n",
      "What type of random graph(s) do you want to create from the list above? (1/2/...)1\n",
      "What is the number of nodes?20\n",
      "What is the probability of having an edge for each two nodes? (0.1, 0.5, ...)0.5\n",
      "What is the number of negative edges for this instance? (0/1/2/.../m)5\n",
      "What is the number of nodes?40\n",
      "What is the probability of having an edge for each two nodes? (0.1, 0.5, ...)0.5\n",
      "What is the number of negative edges for this instance? (0/1/2/.../m)10\n"
     ]
    }
   ],
   "source": [
    "# This code is for producing synthetic data (generating random \"signed graphs\")\n",
    "# \"Signed graphs\" are graphs in which edges are postive or negative (signed)\n",
    "# This code asks you for type of random graph, several parameters for generating them such as size and order\n",
    "# and then it asks for \"number of negative edges\" so that it randomly makes some edges negative and produce \n",
    "# a signed graph out of an unsigned graph\n",
    "\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import random\n",
    "\n",
    "\n",
    "listOfGraphs=[]\n",
    "NumberOfNegative=[]\n",
    "run=int(input(\"How many random graphs do you want to create?(1/2/...)\"))\n",
    "\n",
    "print(\"\\n 1.  Erdős-Rényi binomial graph G_{n,p} \"\n",
    "      \"\\n 2.  Random graph of n nodes and m total edges G_{n,m} \"\n",
    "      \"\\n 3.  Newman-Watts-Strogatz small world graph (adding new edges to the ring) G_{n,k,p} \"\n",
    "      \"\\n 4.  Watts-Strogatz small world graph (rewiring the ring) G_{n,k,p} \"\n",
    "      \"\\n 5.  Random regular graph of n nodes each with degree d G_{d,n} \"\n",
    "      \"\\n 6.  Barabási-Albert preferential attachment model graph G_{n,m} \"\n",
    "      \"\\n 7.  Growing graph with powerlaw degree distribution and approximate average clustering G_{n,m,p}\"\n",
    "      \"\\n 8.  Two-dimensional grid n*n \")\n",
    "\n",
    "graph_type=int(input(\"What type of random graph(s) do you want to create from the list above? (1/2/...)\"))\n",
    "\n",
    "for i in range(run):\n",
    "\n",
    "    if graph_type == 1:\n",
    "        n=int(input(\"What is the number of nodes?\"))\n",
    "        p=float(input(\"What is the probability of having an edge for each two nodes? (0.1, 0.5, ...)\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.fast_gnp_random_graph(n,p)\n",
    "        listOfGraphs.append(G) \n",
    "    elif graph_type == 2:\n",
    "        n=int(input(\"What is the number of nodes?\"))\n",
    "        m=int(input(\"What is the total number of edges?\"))  \n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.gnm_random_graph(n,m)\n",
    "        listOfGraphs.append(G)\n",
    "    elif graph_type == 3:\n",
    "        n=int(input(\"What is the number of nodes?\"))\n",
    "        k=int(input(\"What is the number of nearest neighbors in ring topology?\"))\n",
    "        p=float(input(\"What is the probability of adding a new edge?\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.newman_watts_strogatz_graph(n,k,p)\n",
    "        listOfGraphs.append(G)\n",
    "    elif graph_type == 4:\n",
    "        n=int(input(\"What is the number of nodes?\"))\n",
    "        k=int(input(\"What is the number of nearest neighbors in ring topology?\"))\n",
    "        p=float(input(\"What is the probability of rewiring each edge?\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.watts_strogatz_graph(n,k,p)\n",
    "        listOfGraphs.append(G)\n",
    "    elif graph_type == 5:\n",
    "        n=int(input(\"What is the number of nodes?(n>d)\"))\n",
    "        d=int(input(\"What is the degree?(The value of n*d must be even)\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.random_regular_graph(d,n)\n",
    "        listOfGraphs.append(G)\n",
    "    elif graph_type == 6:\n",
    "        n=int(input(\"What is the number of nodes?\"))\n",
    "        m=int(input(\"What is the number of edges to attach from a new node to existing nodes?\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.barabasi_albert_graph(n,m)   \n",
    "        listOfGraphs.append(G)\n",
    "    elif graph_type == 7:\n",
    "        n=int(input(\"What is the number of nodes?\"))\n",
    "        m=int(input(\"What is the number of random edges to add for each new node?\"))\n",
    "        p=float(input(\"What is the probability of adding a triangle after adding a random edge?\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.powerlaw_cluster_graph(n,m,p)    \n",
    "        listOfGraphs.append(G)\n",
    "    elif graph_type == 8:\n",
    "        n=int(input(\"What is the size of the grid?(sqrt(n))\"))\n",
    "        NumberOfNegative.append(int(input(\"What is the number of negative edges for this instance? (0/1/2/.../m)\")))\n",
    "        G=nx.grid_2d_graph(n,n)\n",
    "        mapping=dict(zip(G.nodes(),range(len(G.nodes()))))\n",
    "        G=nx.relabel_nodes(G,mapping)\n",
    "        listOfGraphs.append(G)\n",
    "    else:\n",
    "        print(\"Try again and choose a graph type from the selection above!\")\n",
    "\n",
    "                    \n",
    "Graph=[]\n",
    "signedMatrix=[]\n",
    "unsignedMatrix=[]\n",
    "ShuffledEdges=[]\n",
    "sorted_weighted_edges=[]\n",
    "\n",
    "\n",
    "for i in range(run):\n",
    "        ShuffledEdges = list(listOfGraphs[i].edges())\n",
    "        random.shuffle(ShuffledEdges)\n",
    "        Graph.append(nx.Graph())\n",
    "        for u in listOfGraphs[i].nodes():\n",
    "            Graph[i].add_node(u)\n",
    "        for ind in range(0,NumberOfNegative[i]):\n",
    "            (u,v) = ShuffledEdges[ind]\n",
    "            Graph[i].add_edge(u,v,weight=-1)\n",
    "        for ind in range(NumberOfNegative[i],len(listOfGraphs[i].edges())):\n",
    "            (u,v) = ShuffledEdges[ind]\n",
    "            Graph[i].add_edge(u,v,weight=1)\n",
    "        signedMatrix.append(nx.to_numpy_matrix(Graph[i]))\n",
    "        unsignedMatrix.append(abs(signedMatrix[i]))\n",
    "        ShuffledEdges=[]\n",
    "\n",
    "        weighted_edges=nx.get_edge_attributes(Graph[i], 'weight') \n",
    "        sorted_weighted_edges.append({})\n",
    "        for (u,v) in weighted_edges:\n",
    "            if u<v:\n",
    "                (sorted_weighted_edges[i])[(u,v)] = weighted_edges[(u,v)]\n",
    "            if u>v:\n",
    "                (sorted_weighted_edges[i])[(v,u)] = weighted_edges[(u,v)]          \n",
    "          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter Method unchanged\n",
      "   Value: -1  Min: -1  Max: 5  Default: -1\n",
      "Changed value of parameter TimeLimit to 36000.0\n",
      "   Prev: 1e+100  Min: 0.0  Max: 1e+100  Default: 1e+100\n",
      "Parameter OutputFlag unchanged\n",
      "   Value: 1  Min: 0  Max: 1  Default: 1\n",
      "Changed value of parameter Threads to 4\n",
      "   Prev: 0  Min: 0  Max: 1024  Default: 0\n",
      "Optimize a model with 244 rows, 125 columns and 730 nonzeros\n",
      "Variable types: 0 continuous, 125 integer (125 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective 5.0000000\n",
      "Presolve removed 19 rows and 5 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 225 rows, 120 columns, 665 nonzeros\n",
      "Extracted 33 lazy constraints\n",
      "Variable types: 0 continuous, 120 integer (120 binary)\n",
      "\n",
      "Root relaxation: cutoff, 27 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0     cutoff    0         5.00000    5.00000  0.00%     -    0s\n",
      "\n",
      "Explored 0 nodes (27 simplex iterations) in 0.04 seconds\n",
      "Thread count was 4 (of 4 available processors)\n",
      "\n",
      "Solution count 1: 5 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 5.000000000000e+00, best bound 5.000000000000e+00, gap 0.0000%\n",
      "Instance 0  solution equals 5.0\n",
      "--------------------------------------------------------------------------------------------\n",
      "Parameter Method unchanged\n",
      "   Value: -1  Min: -1  Max: 5  Default: -1\n",
      "Changed value of parameter TimeLimit to 36000.0\n",
      "   Prev: 1e+100  Min: 0.0  Max: 1e+100  Default: 1e+100\n",
      "Parameter OutputFlag unchanged\n",
      "   Value: 1  Min: 0  Max: 1  Default: 1\n",
      "Changed value of parameter Threads to 4\n",
      "   Prev: 0  Min: 0  Max: 1024  Default: 0\n",
      "Optimize a model with 871 rows, 428 columns and 2611 nonzeros\n",
      "Variable types: 0 continuous, 428 integer (428 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective 10.0000000\n",
      "Presolve removed 37 rows and 12 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 834 rows, 416 columns, 2488 nonzeros\n",
      "Extracted 94 lazy constraints\n",
      "Variable types: 0 continuous, 416 integer (416 binary)\n",
      "\n",
      "Root relaxation: cutoff, 63 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0     cutoff    0        10.00000   10.00000  0.00%     -    0s\n",
      "\n",
      "Explored 0 nodes (63 simplex iterations) in 0.04 seconds\n",
      "Thread count was 4 (of 4 available processors)\n",
      "\n",
      "Solution count 1: 10 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.000000000000e+01, best bound 1.000000000000e+01, gap 0.0000%\n",
      "Instance 1  solution equals 10.0\n",
      "--------------------------------------------------------------------------------------------\n",
      "-------------------------------- ***  EXPERIMENT STATS  *** --------------------------------\n",
      "--------------------------------------------------------------------------------------------\n",
      "Frustration indices: [5.0, 10.0]\n",
      "Average frustrarion index 7.5\n",
      "\n",
      "Solve times: [0.06 0.05]\n",
      "Average solve time 0.05299997329711914\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# XOR model with L(G) as the optimal objective function\n",
    "# Based on the binary linear programming model discussed in the following publication\n",
    "# http://doi.org/10.1002/net.21907\n",
    "\n",
    "# This code solves graph optimization model(s) using \"Gurobi solver\"\n",
    "# to compute the measure called, frustration index, for the input signed graph(s)\n",
    "\n",
    "# Note that you must have installed Gurobi into Jupyter and registered a Gurobi license\n",
    "# in order to run this code\n",
    "\n",
    "\n",
    "#Setting parameters\n",
    "#lazyParam=int(input(\"What is the lazy parameter for unbalanced triangle lazy cuts? (0/1/2/3)\"))\n",
    "# See \"lazy\" as a tunable parameter in linear constraint attributes in Gurobi optimizer reference manual below:\n",
    "# https://www.gurobi.com/documentation/8.1/refman/lazy.html\n",
    "lazyParam=int(1)\n",
    "\n",
    "\n",
    "#methodParam=int(input(\"What method do you want to use?\"))\n",
    "#(-1=automatic, 0=primal simplex, 1=dual simplex, 2=barrier, 3=concurrent, 4=deterministic concurrent)\n",
    "methodParam=int(-1)\n",
    "\n",
    "#speedupParam=int(input(\"Do you want to use the speedups? (0=No, 1=Yes)\"))\n",
    "speedupParam=int(1)\n",
    "\n",
    "import time\n",
    "import numpy as np\n",
    "import multiprocessing\n",
    "from gurobipy import *\n",
    "objectivevalue=[]\n",
    "solveTime=[]\n",
    "effectiveBranchingFactors=[]\n",
    "\n",
    "for index in range(run):\n",
    "    \n",
    "    order=len(signedMatrix[index])\n",
    "    NumberOfNegative=((-1 == signedMatrix[index])).sum()/2\n",
    "    size=int(np.count_nonzero(signedMatrix[index])/2)\n",
    "\n",
    "    neighbors={}\n",
    "    Degree=[]\n",
    "    for u in sorted((Graph[index]).nodes()):\n",
    "        neighbors[u] = list((Graph[index])[u])\n",
    "        Degree.append(len(neighbors[u]))\n",
    "    unsignedDegree=Degree\n",
    "    \n",
    "    #Finding the node with the highest unsigned degree\n",
    "    maximum_degree = max(unsignedDegree)\n",
    "    [node_to_fix]=[([i for i, j in enumerate(unsignedDegree) if j == maximum_degree]).pop()]\n",
    "\n",
    "    # Model parameters\n",
    "    model = Model(\"XOR model for computing frustration index\")\n",
    "\n",
    "    # The method you have selected above will be used:\n",
    "    # (-1=automatic, 0=primal simplex, 1=dual simplex, 2=barrier, 3=concurrent, 4=deterministic concurrent)\n",
    "    model.setParam(GRB.param.Method, methodParam)\n",
    "    \n",
    "    # What is the time limit in second?\n",
    "    model.setParam('TimeLimit', 10*3600)\n",
    "    \n",
    "    # Do you want details of branching to be reported? (0=No, 1=Yes)\n",
    "    model.setParam(GRB.param.OutputFlag, 1) \n",
    "    \n",
    "    # Do you want a non-zero Mixed integer programming tolerance (MIP Gap)?\n",
    "    # Note that a non-zero MIP gap may prevent the model from computing the exact value of frustration index\n",
    "    # model.setParam('MIPGap', 0.0001)  \n",
    "    \n",
    "    # How many threads to be used for exploring the feasible space in parallel? \n",
    "    model.setParam(GRB.Param.Threads, multiprocessing.cpu_count())   \n",
    "    \n",
    "    #This chunk of code lists the graph triangles\n",
    "    GraphTriangles=[]\n",
    "    for n1 in sorted((Graph[index]).nodes()):\n",
    "        neighbors1 = set((Graph[index])[n1])\n",
    "        for n2 in filter(lambda x: x>n1, neighbors1):\n",
    "            neighbors2 = set((Graph[index])[n2])\n",
    "            common = neighbors1 & neighbors2\n",
    "            for n3 in filter(lambda x: x>n2, common):\n",
    "                GraphTriangles.append([n1,n2,n3])\n",
    "    len(GraphTriangles)\n",
    "    #print(\"--- %Listed\",len(GraphTriangles),\"triangles for the graph\")\n",
    "\n",
    "    #This chunk of code lists the balanced and unbalanced triangles\n",
    "    w=nx.get_edge_attributes(Graph[index], 'weight')  \n",
    "    unbalanced_triangles = []\n",
    "    balanced_triangles = []\n",
    "    for triad in GraphTriangles: \n",
    "        if  (sorted_weighted_edges[index])[(triad[0],triad[1])]*\\\n",
    "        (sorted_weighted_edges[index])[(triad[0],triad[2])]*(sorted_weighted_edges[index])[(triad[1],triad[2])] == -1:\n",
    "            unbalanced_triangles.append(triad)\n",
    "        elif (sorted_weighted_edges[index])[(triad[0],triad[1])]*\\\n",
    "        (sorted_weighted_edges[index])[(triad[0],triad[2])]*(sorted_weighted_edges[index])[(triad[1],triad[2])] == 1:\n",
    "            balanced_triangles.append(triad)  \n",
    "    \n",
    "\n",
    "    # Create decision variables and update model to integrate new variables\n",
    "    x=[]\n",
    "    for i in range(0,order):\n",
    "        x.append(model.addVar(vtype=GRB.BINARY, name='x'+str(i))) # arguments by name\n",
    "    model.update()\n",
    "    \n",
    "    f={}\n",
    "    for (i,j) in (sorted_weighted_edges[index]):\n",
    "        f[(i,j)]=model.addVar(lb=0.0, ub=1, vtype=GRB.BINARY, name='f'+str(i)+','+str(j))\n",
    "    model.update()\n",
    "\n",
    "    # Set the objective function\n",
    "    OFV=0\n",
    "    for (i,j) in (sorted_weighted_edges[index]):\n",
    "        OFV = OFV + f[(i,j)]                 \n",
    "    model.setObjective(OFV, GRB.MINIMIZE)\n",
    "\n",
    "    # Add constraints to the model and update model to integrate new constraints\n",
    "    \n",
    "    ## ADD CORE CONSTRAINTS ##\n",
    "\n",
    "    for (i,j) in (sorted_weighted_edges[index]):\n",
    "        model.addConstr( f[(i,j)] >= x[i] - ((sorted_weighted_edges[index])[(i,j)])*x[j] -\\\n",
    "                        (1-(sorted_weighted_edges[index])[(i,j)])/2\n",
    "                             , '1st Edge'+','+str(i)+','+str(j))\n",
    "        model.addConstr( f[(i,j)] >= -x[i] + ((sorted_weighted_edges[index])[(i,j)])*x[j] +\\\n",
    "                        (1-(sorted_weighted_edges[index])[(i,j)])/2\n",
    "                             , '2nd Edge'+','+str(i)+','+str(j))                  \n",
    "    model.update()  \n",
    "    \n",
    "    ## ADD ADDITIONAL CONSTRAINTS (speed-ups) ##\n",
    "    \n",
    "    if speedupParam==1:\n",
    "\n",
    "        # Triangle valid inequalities            \n",
    "        #lazyTriangles=int(input(\"What triangles do you want to implement the inequalities on? (1:all, 2:unbalanced)\"))\n",
    "        lazyTriangles=int(2)\n",
    "\n",
    "        if lazyTriangles == 1:      \n",
    "            triangleInequalityCount=len(GraphTriangles)*4     \n",
    "            for triangle in GraphTriangles:\n",
    "                [i,j,k]=triangle\n",
    "                model.addConstr( ((f[(i,j)] + ((sorted_weighted_edges[index])[(i,j)] -1)/2)/\\\n",
    "                                  (sorted_weighted_edges[index])[(i,j)])+\n",
    "                                ((f[(i,k)] + ((sorted_weighted_edges[index])[(i,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(i,k)]) >= \n",
    "                                ((f[(j,k)] + ((sorted_weighted_edges[index])[(j,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(j,k)]),\n",
    "                                'triangle1'+','+str(i)+','+str(j)+','+str(k))\n",
    "                model.addConstr( ((f[(i,j)] + ((sorted_weighted_edges[index])[(i,j)] -1)/2)/\\\n",
    "                                  (sorted_weighted_edges[index])[(i,j)])+\n",
    "                                ((f[(j,k)] + ((sorted_weighted_edges[index])[(j,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(j,k)]) >= \n",
    "                                ((f[(i,k)] + ((sorted_weighted_edges[index])[(i,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(i,k)]),\n",
    "                                'triangle2'+','+str(i)+','+str(j)+','+str(k))\n",
    "                model.addConstr( ((f[(j,k)] + ((sorted_weighted_edges[index])[(j,k)] -1)/2)/\\\n",
    "                                  (sorted_weighted_edges[index])[(j,k)])+\n",
    "                                ((f[(i,k)] + ((sorted_weighted_edges[index])[(i,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(i,k)]) >= \n",
    "                                ((f[(i,j)] + ((sorted_weighted_edges[index])[(i,j)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(i,j)]),\n",
    "                                'triangle3'+','+str(i)+','+str(j)+','+str(k))\n",
    "                model.addConstr( ((f[(i,j)] + ((sorted_weighted_edges[index])[(i,j)] -1)/2)/\\\n",
    "                                  (sorted_weighted_edges[index])[(i,j)])+\n",
    "                                ((f[(i,k)] + ((sorted_weighted_edges[index])[(i,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(i,k)]) + \n",
    "                                ((f[(j,k)] + ((sorted_weighted_edges[index])[(j,k)] -1)/2)/\\\n",
    "                                 (sorted_weighted_edges[index])[(j,k)]) <= 2,\n",
    "                                'triangle4'+','+str(i)+','+str(j)+','+str(k))\n",
    "\n",
    "        elif lazyTriangles == 2:\n",
    "            triangleInequalityCount=len(unbalanced_triangles)\n",
    "            for triangle in unbalanced_triangles:\n",
    "                [i,j,k]=triangle\n",
    "                model.addConstr(f[(i,j)] + f[(i,k)] + f[(j,k)] >= 1 ,'UnbalancedTriangle'+','+str(i)+','+str(j)+','+str(k))    \n",
    "\n",
    "        elif lazyTriangles == 3:\n",
    "            triangleInequalityCount=0\n",
    "\n",
    "        model.update()\n",
    "        model.setAttr('Lazy',model.getConstrs()[2*size:2*size+triangleInequalityCount],[lazyParam]*triangleInequalityCount)\n",
    "        model.update() \n",
    "\n",
    "        # Colour the node with the highest degree as 1\n",
    "        model.addConstr(x[node_to_fix]==1 , '1stnodecolour')   \n",
    "        model.update()\n",
    "\n",
    "        #branching priority is based on unsigned degree \n",
    "        model.setAttr('BranchPriority',model.getVars()[:order],unsignedDegree)\n",
    "        model.update()\n",
    "    \n",
    "\n",
    "    # Solve\n",
    "    start_time = time.time()\n",
    "    model.optimize()\n",
    "    solveTime.append(time.time() - start_time) \n",
    "    \n",
    "    \n",
    "    # Save optimal objective function values\n",
    "    obj = model.getObjective()\n",
    "    objectivevalue.append((obj.getValue()))\n",
    "    \n",
    "    # Compute the effective branching factors\n",
    "    if (model.NodeCount)**(1/(size+2*order)) >= 1:\n",
    "        effectiveBranchingFactors.append((model.NodeCount)**(1/(size+2*order)))\n",
    "        \n",
    "    # Report the optimal objective function value for each instance\n",
    "    print('Instance', index,' solution equals',objectivevalue[index]) \n",
    "    print(\"-\"*92)\n",
    "    \n",
    "    # Printing the solution (optional)\n",
    "    #print(\"Optimal values of the decision variables\")\n",
    "    #for v in model.getVars():\n",
    "    #    print (v.varName, v.x)\n",
    "    #print()    \n",
    "\n",
    "print(\"-\"*32,\"***  EXPERIMENT STATS  ***\",\"-\"*32)\n",
    "print(\"-\"*92)\n",
    "print(\"Frustration indices:\",objectivevalue)\n",
    "print(\"Average frustrarion index\",np.mean(objectivevalue))\n",
    "#print(\"Frustration index Standard Deviation\",np.std(objectivevalue))\n",
    "print()\n",
    "print(\"Solve times:\",np.around(solveTime, decimals=2))\n",
    "print(\"Average solve time\",np.mean(solveTime))\n",
    "#print(\"Solve time Standard Deviation\",np.std(solveTime))\n",
    "print()\n",
    "if (model.NodeCount)**(1/(size+2*order)) >= 1:\n",
    "    print(\"Effective branching factors (EBFs)\",effectiveBranchingFactors)\n",
    "    print(\"Average EBF\",np.mean(effectiveBranchingFactors))\n",
    "    #print(\"EBF Standard Deviation\",np.std(effectiveBranchingFactors))"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}