Python_LP_GIP.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "source": [
        "!pip install highspy"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cAXtEH8WP9cI",
        "outputId": "280d70ce-8a93-4ede-cf96-8400737acde6"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Requirement already satisfied: highspy in /usr/local/lib/python3.11/dist-packages (1.10.0)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.11/dist-packages (from highspy) (2.0.2)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "T_2hv8MhJ_4T"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import cvxpy as cp\n",
        "import time\n",
        "import networkx as nx\n",
        "import warnings\n",
        "\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "####################################################################################\n",
        "# Linear programming GIP proposed by Cajas (2025)\n",
        "####################################################################################\n",
        "\n",
        "def LP_GIP_Cajas(A, B, solver='HIGHS'):\n",
        "    n, _ = A.shape\n",
        "    vec_A = A.reshape((-1,1), order='F')\n",
        "    vec_B = B.reshape((-1,1), order='F')\n",
        "    J = np.ones((n,n))\n",
        "\n",
        "    X = cp.Variable((n,n))\n",
        "    Y = cp.Variable((n**2,n**2), symmetric=True)\n",
        "    vec_X = cp.reshape(cp.vec(X, order='F'), (-1,1))\n",
        "    opt_time = 0\n",
        "\n",
        "    for j in range(n):\n",
        "        Aj = cp.reshape(Y[j * n: j * n + 1].T, (n, n), order=\"F\")\n",
        "        for i in range(1,n):\n",
        "            Aij = cp.reshape(Y[i + j * n : i + j * n + 1].T, (n, n), order=\"F\")\n",
        "            Aj = cp.vstack([Aj, Aij])\n",
        "        if j == 0:\n",
        "            Z = Aj\n",
        "        else:\n",
        "            Z = cp.hstack([Z, Aj])\n",
        "\n",
        "    cost = cp.trace(J @ X)\n",
        "    objective = cp.Maximize(cost)\n",
        "\n",
        "    ones = np.ones((n,1))\n",
        "    ones2 = np.ones((n**2,1))\n",
        "\n",
        "    constraints = []\n",
        "    constraints += [Z @ ones2 == 1]\n",
        "    constraints += [Z.T @ ones2 == 1]\n",
        "    constraints += [Z @ vec_A - vec_B == 0]\n",
        "    constraints += [cp.reshape(cp.diag(Y), (-1,1), order='F') - vec_X == 0]\n",
        "    constraints += [Y @ ones2 - n * vec_X == 0]\n",
        "    constraints += [X @ ones == 1, X.T @ ones == 1]\n",
        "    constraints += [X @ A == B @ X]\n",
        "    constraints += [Y >= 0, Y <= 1]\n",
        "    constraints += [X >= 0, X <= 1]\n",
        "\n",
        "    # Solving the problem\n",
        "    prob = cp.Problem(objective, constraints)\n",
        "    start_time = time.time()\n",
        "    prob.solve(solver=solver,\n",
        "               verbose=False,\n",
        "               warm_start=False)\n",
        "    opt_time += (time.time() - start_time)\n",
        "\n",
        "    # Calculate the difference between calculated permuted matrix and original permuted matrix\n",
        "    B1 = (Z.value @ vec_A).reshape((n,n), order='F')\n",
        "    diff = np.abs(B1 - B).max()\n",
        "\n",
        "    return diff, opt_time"
      ],
      "metadata": {
        "id": "wdzrKlwKKQ_-"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "####################################################################################\n",
        "# Creating some instances of the GIP\n",
        "####################################################################################\n",
        "\n",
        "n = 13\n",
        "Gs = {}\n",
        "Gs[0] = nx.cycle_graph(n)\n",
        "Gs[1] = nx.complete_graph(n)\n",
        "Gs[2] = nx.paley_graph(n)\n",
        "Gs[3] = nx.turan_graph(n, 8)\n",
        "Gs[4] = nx.fast_gnp_random_graph(n, 0.2)\n",
        "Gs[5] = nx.random_regular_graph(d=4, n=n, seed=123)\n",
        "\n",
        "As = [(i, nx.adjacency_matrix(Gs[i]).toarray()) for i in Gs.keys()]\n",
        "As = dict(As)\n",
        "\n",
        "rs = np.random.RandomState(123)\n",
        "I = np.identity(n)\n",
        "a = rs.permutation(n)\n",
        "Bs = [(i, I[a] @ As[i] @ I[a].T) for i in As.keys()]\n",
        "Bs = dict(Bs)"
      ],
      "metadata": {
        "id": "08Rdm8ktKbyo"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "####################################################################################\n",
        "# Calculating the permutation matrix that transform A into B\n",
        "# I recommend MOSEK or other commercial solvers for larger instances (n > 13)\n",
        "####################################################################################\n",
        "\n",
        "for i in As.keys():\n",
        "  b, t = LP_GIP_Cajas(As[i], Bs[i], solver='HIGHS')\n",
        "  print('Optimization Time:', np.round(t, 2), ', Max difference:', np.round(b, 7))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bBxIFMpmL1pz",
        "outputId": "32a54f7e-74bf-4d6a-bc17-11acc1ae00d0"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Optimization Time: 78.41 , Max difference: 0.0\n",
            "Optimization Time: 6.42 , Max difference: 0.0\n",
            "Optimization Time: 48.75 , Max difference: 0.0\n",
            "Optimization Time: 9.49 , Max difference: 0.0\n",
            "Optimization Time: 7.87 , Max difference: 0.0\n",
            "Optimization Time: 52.83 , Max difference: 0.0\n"
          ]
        }
      ]
    }
  ]
}