phase_angle_correlation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/YeoniH/geometry-driven-hyperuniformity/blob/main/phase_angle_correlation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "In this Notebook, we will see how the phase angle correlations develop along the evolution of a point pattern via Lloyd iterations.\n",
        "\n",
        "The *hexatic order parameter* $\\psi_6$ of a point $\\mathbf{x}$ can be defined as follows:\n",
        "\\begin{equation}\n",
        "  \\psi_6(\\mathbf{x}) = \\frac{1}{n_\\mathbf{x}} \\sum_{\\mathbf{y}}e^{6i\\theta_{xy}} = |\\psi_6(\\mathbf{x})|e^{i\\Theta(\\mathbf{x})},\n",
        "\\end{equation}\n",
        "which is a complex number with magnitude $|\\psi_6(\\mathbf{x})|\\leq 1$ and phase angle $\\Theta(\\mathbf{x}) \\in [-\\pi,\\pi]$, which indicates the average bond orientation.\n",
        "The hexatic order of a point configuration can be assessed by the average hexatic order of the points contained.\n",
        "\n",
        "To proceed with the followings, please copy the Notebook into your local/Google Drive and run the cells below."
      ],
      "metadata": {
        "id": "qDckEW6DyfLs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, we create a folder named `my_projects` and change directory (via shell command `cd`) to the newly created folder."
      ],
      "metadata": {
        "id": "eUJxwdZc0e4S"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!mkdir my_projects\n",
        "%cd my_projects"
      ],
      "metadata": {
        "id": "I3uF3z_1A93x"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We will need to clone the following GitHub repository, where all the codes we need are stored."
      ],
      "metadata": {
        "id": "KjLW2cHm034s"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oMlWYgEA6Eo4"
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/YeoniH/geometry-driven-hyperuniformity.git"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Whenever you run this Notebook afresh, please make sure you run the following line to pull all the new updates in the GitHub repository."
      ],
      "metadata": {
        "id": "w_WZzsjy1DYr"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!git pull"
      ],
      "metadata": {
        "id": "Po9AgKea9uNe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, we will install one of the fundamental libraries, called `freud-analysis`, for our computational analyses of an evolving point configuration."
      ],
      "metadata": {
        "id": "16PqJOgR1SQY"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install freud-analysis"
      ],
      "metadata": {
        "id": "EmBNRsHb9wKT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Just to double check if we have all we need, we use `ls` (list command).\n",
        "\n",
        "(If you see `geometry-driven-hyperuniformity` and `/content/my_projects/geometry-driven-hyperuniformity` as output, then you're good to go.)"
      ],
      "metadata": {
        "id": "m3h1FpY_1fbQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!ls\n",
        "%cd geometry-driven-hyperuniformity"
      ],
      "metadata": {
        "id": "TedC9PjmKAqs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "In the following cell, we define a new function to generate a realisation of Poisson point process with `density` value and side length `L` of the\n",
        "(square) simulation box as input."
      ],
      "metadata": {
        "id": "_VP94jOw27Bv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from os import system\n",
        "import numpy as np\n",
        "\n",
        "def generate_2d_poisson_pp(density, L):\n",
        "  \"\"\"\n",
        "  Args:\n",
        "    density: rate, i.e., number of points per unit volume\n",
        "    L: side length of the square simulation box\n",
        "  \"\"\"\n",
        "  mean = density * (L**2)\n",
        "  N = np.random.poisson(mean)\n",
        "\n",
        "  X = np.random.uniform(-L/2, L/2, N)\n",
        "  Y = np.random.uniform(-L/2, L/2, N)\n",
        "\n",
        "  coord = np.zeros((N, 3))\n",
        "  for i in range(N):\n",
        "    coord[i] = [X[i], Y[i], 0.0]\n",
        "  return coord"
      ],
      "metadata": {
        "id": "OyAKYn3U-cNk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, let us generate a new point pattern using the function `generate_2d_poisson_pp()` that we just defined.\n",
        "\n",
        "(You can play around with different values for `L`, but it is recommended to use `L` larger than 50 to obtain a nice result at the end, which is otherwise make it hard to see the develpment of hexagonal domains and the correesponding correlation length.)"
      ],
      "metadata": {
        "id": "VUv1DSVA3LEr"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "L = 20\n",
        "init_points = generate_2d_poisson_pp(density=1, L=L)\n",
        "num_points = len(init_points)"
      ],
      "metadata": {
        "id": "g5VxxJ0SDkSr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let us print out how many numbers the newly generated system (point pattern/configuration) has and also see how it looks like."
      ],
      "metadata": {
        "id": "YP4VlzKR3w-0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "print(\"There are\", num_points, \"points\")\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "ax.tick_params(axis='both', labelsize=11)\n",
        "ax.set_aspect('equal', 'box')\n",
        "ax.scatter(init_points[:, 0], init_points[:, 1], s=2, color=\"k\")\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "Wc9d6YJGD0o4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Given the initial point pattern, we run the Lloyd's algorithm to drive the system into a frozen converged state (effectively hyperuniform state) and plot the temporal evolution of the total quantizer energy of the system along the Lloyd dynamic.\n",
        "\n",
        "You can choose the variable `num_iter` as small or large as you want, but it is highly recommendable to run it sufficiently (say $\\geq$ 10000 iterations for systems with more than 5000 points) to attain hyperuniformity."
      ],
      "metadata": {
        "id": "oOcO-pJZ38Gg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from Analysis import PeriodicVoro as PV\n",
        "from tqdm import tqdm\n",
        "\n",
        "dim = 2\n",
        "num_iter = 1000   # number of iterations steps of the Lloyd's algorithm\n",
        "\n",
        "init_config = PV(dim=dim, points=init_points, L=L, step=0)\n",
        "total_qe = []\n",
        "\n",
        "points = init_points\n",
        "for t in tqdm(range(num_iter)):\n",
        "  config = PV(dim=dim, points=points, L=L, step=t)\n",
        "  total_qe.append(config.compute_total_quantizer_energy())\n",
        "  points = config.update_by_lloyd_algorithm()\n",
        "final_points = points\n",
        "\n",
        "plt.plot(list(range(1, num_iter)), total_qe[1:])\n",
        "plt.xlabel('Lloyd step', fontsize=12)\n",
        "plt.ylabel('Total quantizer energy', fontsize=12)\n",
        "plt.xscale('log')\n",
        "plt.yscale('log')\n",
        "plt.show()\n",
        "\n",
        "print('Total quantizer energy at step', num_iter, ':', total_qe[-1])"
      ],
      "metadata": {
        "id": "SRbw5l0jECxe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can now plot the Voronoi landscape of the final point pattern obtained above (at Lloyd step `num_iter`), highlighting topological defects, namely pentagons and heptagons in blue and red, respectively."
      ],
      "metadata": {
        "id": "HR1sMW9z4fI3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig, ax = plt.subplots(figsize=(8, 8))\n",
        "ax.tick_params(axis='both', labelsize=11)\n",
        "ax.set_aspect('equal', 'box')\n",
        "\n",
        "final_config = PV(dim, final_points, L, num_iter)\n",
        "for idx, polytope in enumerate(final_config.voro.polytopes):\n",
        "  facecolor = {5: \"b\", 6: \"w\", 7: \"r\"}\n",
        "  ax.fill(*zip(*polytope[:, :2]), facecolor=facecolor[len(polytope)], edgecolor=\"k\", linewidth=1.)\n",
        "if num_points <= 2000:\n",
        "  ax.scatter(final_points[:, 0], final_points[:, 1], s=2, color=\"k\")\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "GixwVuXvju2D"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let us compute the correlation length of the final point configuration, which represents the average hexagonal domain size."
      ],
      "metadata": {
        "id": "lkpm0jVB43b0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "final_config.correlation_length()"
      ],
      "metadata": {
        "id": "qtnFssirIcIp"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}