introduction_pygeoda_0.0.3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to pygeoda 0.0.3\n",
    "\n",
    "Xun Li 17/10/2019\n",
    "\n",
    "pygeoda is a R package that wraps all core functions of spatial data\n",
    "analysis in GeoDa and libgeoda. Unlike the desktop software GeoDa,\n",
    "libgeoda is a non-UI and feature focused C++ library that is designed\n",
    "for programmers to do spatial data analysis using their favoriate\n",
    "programming languages (R, Python, Java etc.). It also aims to be easily\n",
    "integratd with other libraries, softwares or systems on different\n",
    "platforms.\n",
    "\n",
    "This book is used to introduce the pygeoda v0.0.3 library, includes all\n",
    "the functions that are currently provided in version 0.0.3.\n",
    "\n",
    "**Live Tutorials**\n",
    "\n",
    "There are python jupyter notebooks available to replicate all the content of this book: https://github.com/lixun910/pygeoda_tutorial/tree/v0.0.3\n",
    "\n",
    "You can try these python jupyter notebooks in your browser via MyBinder (no installation required): \n",
    "https://mybinder.org/v2/gh/lixun910/pygeoda/v0.0.3\n",
    "\n",
    "## 1. Install pygeoda v0.0.3\n",
    "\n",
    "Like GeoDa desktop software, `pygeoda` are avaiable to different platforms including: Mac, Linux and Windows. \n",
    "\n",
    "In current development stage, `pygeoda` 0.0.3 can be installed from source code or using `pip`.\n",
    "\n",
    "> Note: check pip page `pygeoda 0.0.3` at: https://pypi.org/project/pygeoda/\n",
    "\n",
    "\n",
    "### 1.1 Install pygeoda using pip\n",
    "\n",
    "You can use `pip` to install pygeoda 0.0.3:\n",
    "```Shell\n",
    "pip install pygeoda\n",
    "```\n",
    "\n",
    "### 1.2 Install from source\n",
    "\n",
    "If pip install doesn't work, you can also try to install pygeoda from its source code:\n",
    "```Shell\n",
    "    pip install git+https://github.com/lixun910/pygeoda\n",
    "```\n",
    "\n",
    "\n",
    "## 2. Using pygeoda\n",
    "\n",
    "If everything installed without error, you should be able to load pygeoda:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygeoda"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Load Spatial Data\n",
    "\n",
    "The data formats that pygeoda v0.0.3 can load directly includes:\n",
    "\n",
    "* ESRI Shapefile\n",
    "* MapInfo File\n",
    "* CSV\n",
    "* GML\n",
    "* GPX\n",
    "* KML\n",
    "* GeoJSON\n",
    "* TopoJSON\n",
    "* OpenFileGDB\n",
    "* GFT Google Fusion Tables\n",
    "* CouchDB\n",
    "\n",
    "In this tutorial, we only tested loading ESRI shapefiles\n",
    "    using pygeoda v0.0.3. Please create a ticket in pygeoda's\n",
    "    repo https://github.com/lixun910/pygeoda/issues if you\n",
    "    experience any issues when loading spatial data.\n",
    "    \n",
    "For example, to load the ESRI Shapefile **Guerry.shp**\n",
    "download from: https://geodacenter.github.io/data-and-lab/Guerry  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "guerry = pygeoda.open(\"./data/Guerry.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `pygeoda.open()` function returns a geoda object, which\n",
    "can be used to access the meta-data, fields, and columns of the input dataset.\n",
    "```\n",
    "    Usage\n",
    "    pygeoda.open(ds_path)\n",
    "\n",
    "    Arguments\n",
    "    ds_path\t(str) The path of the spatial dataset\n",
    "\n",
    "    Value\n",
    "    gda_obj An object of geoda class\n",
    "```\n",
    "\n",
    "#### Attributes of `geoda` object\n",
    "\n",
    "* n_cols\n",
    "* n_obs\n",
    "* field_names\n",
    "* field_type\n",
    "\n",
    "To access the meta-data of the loaded Guerry dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of columns: 26\n",
      "number of observations: 85\n",
      "field names: ('CODE_DE', 'COUNT', 'AVE_ID_', 'dept', 'Region', 'Dprtmnt', 'Crm_prs', 'Crm_prp', 'Litercy', 'Donatns', 'Infants', 'Suicids', 'MainCty', 'Wealth', 'Commerc', 'Clergy', 'Crm_prn', 'Infntcd', 'Dntn_cl', 'Lottery', 'Desertn', 'Instrct', 'Prsttts', 'Distanc', 'Area', 'Pop1831')\n",
      "field types: ('string', 'numeric', 'numeric', 'integer', 'string', 'string', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'integer', 'numeric', 'integer', 'numeric')\n"
     ]
    }
   ],
   "source": [
    "print(\"number of columns:\", guerry.num_cols)\n",
    "print(\"number of observations:\", guerry.num_obs)\n",
    "print(\"field names:\", guerry.field_names)\n",
    "print(\"field types:\", guerry.field_types)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Access Table Data\n",
    "\n",
    "The `geoda` instance also provide functions to get data from the dataset:\n",
    "\n",
    "* GetIntegerCol(col_name)\n",
    "* GetRealCol(col_name)\n",
    "* GetStringCol(col_name)\n",
    "\n",
    "For example, to get the integer values of the \"Crm_prp\" column:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "firs 10 values of Crm_prp: (15890, 5521, 7925, 7289, 8174, 10263, 8847, 9597, 4086, 10431)\n"
     ]
    }
   ],
   "source": [
    "crm_prp = guerry.GetIntegerCol(\"Crm_prp\")\n",
    "print(\"firs 10 values of Crm_prp:\", crm_prp[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Spatial Weights\n",
    "\n",
    "Spatial weights are central components in spatial data analysis. The spatial weights represents the possible spatial interaction between observations in space. Like GeoDa desktop software, `pygeoda` provides a rich variety of methods to create several different types of spatial weights:\n",
    "\n",
    "* Contiguity Based Weights: `queen_weights()`, `rook_weights()`\n",
    "* Distance Based Weights: `distance_weights()`\n",
    "* K-Nearest Neighbor Weights: `knn_weights()`\n",
    "* Kernel Weights: `kernel_weights()`\n",
    "\n",
    "### 2.3.1 Queen Contiguity Weights\n",
    "\n",
    "To create a Queen contiguity weights, we can call pygeoda's function \n",
    "```python\n",
    "pygeoda.weights.queen(gda, order=1,\n",
    "              include_lower_order = False, \n",
    "              precision_threshold = 0)\n",
    "``` \n",
    "by passing the GeoDa object `guerry` we just created: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      "\n",
      "is symmetric:True\n",
      "sparsity:0.0\n",
      "density:5.813148788927336\n",
      "min neighbors:2\n",
      "mean neighbors:4.9411764705882355\n",
      "median neighbors:5.0\n",
      "max neighbors:8\n"
     ]
    }
   ],
   "source": [
    "queen_w = pygeoda.weights.queen(guerry)\n",
    "print(queen_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The function `queen_weights()` returns an instance of \n",
    "`Weight` object. One can access the meta data of the spatial\n",
    "weights by accessing the attributes of `Weight` object:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Attributes of `Weight` object\n",
    "\n",
    "* num_obs\n",
    "* GetWeightsType() \n",
    "* IsSymmetric()\n",
    "* HasIsolates()\n",
    "* GetSparsity()\n",
    "* GetMinNbrs()\n",
    "* GetMedianNbrs()\n",
    "* GetMeanNbrs()\n",
    "* GetMaxNbrs()\n",
    "* GetDensity()\n",
    "* [] GetNeighbors(idx)\n",
    "* double SpatialLag(idx, [data])\n",
    "* SaveToFile()\n",
    "\n",
    "We can also access the details of the weights: e.g.\n",
    "list the neighbors of a specified observation, which\n",
    "is very helpful in exploratory spatial data analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neighbors of 0-st observation are: (35, 36, 66, 68)\n"
     ]
    }
   ],
   "source": [
    "nbrs = queen_w.GetNeighbors(0)\n",
    "print(\"Neighbors of 0-st observation are:\", nbrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute the spatial lag of a specified observation by passing the values of the selected variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spatial lag of 0-st observation is: 7899.25\n"
     ]
    }
   ],
   "source": [
    "lag0 = queen_w.SpatialLag(0, crm_prp)\n",
    "print(\"Spatial lag of 0-st observation is:\", lag0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "by passing the geoda object `guerry` we just created:\n",
    "### 2.3.2 Rook Contiguity Weights\n",
    "\n",
    "To create a Rook contiguity weights, we can call pygeoda's function \n",
    "```python\n",
    "pygeoda.weights.rook(gda, order=1, \n",
    "                     include_lower_order=False, \n",
    "                     precision_threshold = 0)\n",
    "``` \n",
    "by passing the GeoDa object `guerry` we just created: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      "\n",
      "is symmetric:True\n",
      "sparsity:0.0\n",
      "density:5.813148788927336\n",
      "min neighbors:2\n",
      "mean neighbors:4.9411764705882355\n",
      "median neighbors:5.0\n",
      "max neighbors:8\n"
     ]
    }
   ],
   "source": [
    "rook_w = pygeoda.weights.rook(guerry)\n",
    "print(rook_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The weights we created are in memory, which makes it straight forward for spatial data analysis and also are good for programming your application. To save the weights to a file, we need to call GeoDaWeight's function \n",
    "```python\n",
    "Weight.SaveToFile(ofname, layer_name, id_name, id_vec)\n",
    "```\n",
    "\n",
    "The `layer_name` is the layer name of loaded dataset. For a ESRI shapefile, the layer name is the file name without the suffix (e.g. Guerry). \n",
    "\n",
    "The `id_name` is a key (column name), which means the associated column contains unique values, that makes sure that the weights are connected to the correct observations in the data table. \n",
    "\n",
    "The `id_vec` is the actual column data of `id_name`, it could be a tuple of integer or string values.\n",
    "\n",
    "For example, in Guerry dataset, the column \"CODE_DE\" can be used as a key to save a weights file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rook_w.SaveToFile('./data/Guerry_r.gal', 'Guerry', 'CODE_DE', guerry.GetIntegerCol('CODE_DE'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we should find the file \"Guerry_r.gal\" in the same directory of Guerry shapefiles.\n",
    "\n",
    "### 2.3.3 Distance Based Weights\n",
    "\n",
    "To create a Distance based weights, we can call pygeoda's function \n",
    "```python\n",
    "pygeoda.weights.distance(gda, dist_thres, power=1.0, \n",
    "                         is_inverse=False, is_arc=False, \n",
    "                         is_mile=True)\n",
    "``` \n",
    "by passing the GeoDa object `guerry` we just created and the value of distance threshold. Like GeoDa, pygeoda provides a function to help you find a optimized distance threshold that guarantees that every observation has at least one neighbor:\n",
    "\n",
    "```python\n",
    "pygeoda.weights.min_threshold(GeoDa gda, bool is_arc = False, bool is_mile = True)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      "\n",
      "is symmetric:False\n",
      "sparsity:0.0\n",
      "density:4.346020761245675\n",
      "min neighbors:1\n",
      "mean neighbors:3.6941176470588237\n",
      "median neighbors:4.0\n",
      "max neighbors:7\n"
     ]
    }
   ],
   "source": [
    "dist_thres = pygeoda.weights.min_threshold(guerry)\n",
    "dist_w = pygeoda.weights.distance(guerry, dist_thres)\n",
    "print(dist_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3.4 K-Nearest Neighbor Weights\n",
    "\n",
    "A special case of distance based weights is K-Nearest neighbor weights, in which every obersvation will have exactly k neighbors. To create a KNN weights, we can call pygeoda's function:\n",
    "```python\n",
    "pygeoda.weights.knn(gda, k, power = 1.0, \n",
    "                    is_inverse = False,\n",
    "                    is_arc = False, is_mile = True)\n",
    "```\n",
    "\n",
    "For example, to create a 6-nearest neighbor weights using Guerry dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      "\n",
      "is symmetric:False\n",
      "sparsity:0.0\n",
      "density:7.0588235294117645\n",
      "min neighbors:6\n",
      "mean neighbors:6.0\n",
      "median neighbors:6.0\n",
      "max neighbors:6\n"
     ]
    }
   ],
   "source": [
    "knn6_w = pygeoda.weights.knn(guerry, 6)\n",
    "print(knn6_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3.5 Kernel Weights\n",
    "\n",
    "Kernel weights apply kernel function to determine the distance decay in the derived continuous weights kernel. The kernel weights are defined as a function K(z) of the ratio between the distance dij from i to j, and the bandwidth hi, with z=dij/hi. \n",
    "\n",
    "The kernl functions include:\n",
    "* triangular\n",
    "* uniform \n",
    "* quadratic\n",
    "* epanechnikov\n",
    "* quartic\n",
    "* gaussian\n",
    "\n",
    "Two functions are provided in pygeoda to create kernel weights:\n",
    "\n",
    "**Kernel Weights with fixed bandwidth**\n",
    "\n",
    "To create a kernel weights with fixed bandwith, using the following function\n",
    "```python\n",
    "pygeoda.weights.kernel_bandwidth(gda, bandwidth, kernel,\n",
    "                                 use_kernel_diagnals = False,\n",
    "                                 is_arc = False, is_mile = True)\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      "\n",
      "is symmetric:False\n",
      "sparsity:0.0\n",
      "density:4.346020761245675\n",
      "min neighbors:1\n",
      "mean neighbors:3.6941176470588237\n",
      "median neighbors:4.0\n",
      "max neighbors:7\n"
     ]
    }
   ],
   "source": [
    "kernel_w = pygeoda.weights.kernel_bandwidth(guerry, dist_thres, \"uniform\")\n",
    "print(kernel_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Besides the options `is_inverse`, `power`, `is_arc` and\n",
    "`is_mile` that are the same with the distance based weights,\n",
    "this kernel weights function has another option\n",
    "```\n",
    "    use_kernel_diagonals\n",
    "    (optional) FALSE (default) or TRUE, apply kernel on the\n",
    "    diagonal of weights matri\n",
    "```\n",
    "\n",
    "**Kernel Weights with adaptive bandwidth**\n",
    "\n",
    "To create a kernel weights with adaptive bandwidth or using max Knn distance as bandwidth, using the following function:\n",
    "```python\n",
    "pygeoda.weights.kernel(gda, k, kernel,\n",
    "                       adaptive_bandwidth = True,\n",
    "                       use_kernel_diagnals = False,\n",
    "                       is_arc = False, is_mile = True)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      "\n",
      "is symmetric:False\n",
      "sparsity:0.0\n",
      "density:7.0588235294117645\n",
      "min neighbors:6\n",
      "mean neighbors:6.0\n",
      "median neighbors:6.0\n",
      "max neighbors:6\n"
     ]
    }
   ],
   "source": [
    "adptkernel_w = pygeoda.weights.kernel(guerry, 6, \"uniform\")\n",
    "print(adptkernel_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This kernel weights function two more options:\n",
    "```\n",
    "\n",
    "    adaptive_bandwidth\n",
    "    (optional) TRUE (default) or FALSE: TRUE use adaptive bandwidth\n",
    "    calculated using distance of k-nearest neithbors, FALSE use max\n",
    "    distance of all observation to their k-nearest neighbors\n",
    "\n",
    "    use_kernel_diagonals\n",
    "    (optional) FALSE (default) or TRUE, apply kernel on the diagonal\n",
    "    of weights matrix\n",
    "```\n",
    "\n",
    "## 3 Spatial Data Analysis\n",
    "\n",
    "\n",
    "### 3.1 Local Spatial Autocorrelation\n",
    "\n",
    "pygeoda 0.0.3 provids following methods for univariate\n",
    "local spatial autocorrelation statistics:\n",
    "\n",
    "* Local Moran: local_moran()\n",
    "* Local Geary: local_geary()\n",
    "* Local Getis-Ord statistics: local_g() and local_gstar()\n",
    "* Local Join Count: local_joincount()\n",
    "\n",
    "\n",
    "Methods for bivariate and multivariate local spatial autocorrelation statistics, as well as global spatial autocorrelation satatistics, will be included in next release of pygeoda.\n",
    "\n",
    "In this tutorial, we will only introduce how to call these methods using pygeoda. For more information about the local spatial autocorrelation statisticis, please read: http://geodacenter.github.io/workbook/6a_local_auto/lab6a.html. \n",
    "\n",
    "#### 3.1.1 Local Moran\n",
    "\n",
    "The Local Moran statistic is a method to identify local clusters and local spatial outliers. For example, we can call  function `local_moran()` with the created Queen weights and the data “crm_prp” as input parameters:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "lisa = pygeoda.local_moran(queen_w, crm_prp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `local_moran()` function will return a `lisa` object,\n",
    "which we can call its functions to access the results of lisa\n",
    "computation. The functions include:\n",
    "\n",
    "* GetClusterIndicators(): Get the local cluster indicators returned from LISA computation.\n",
    "* GetColors(): Get the cluster colors of LISA computation.\n",
    "* GetLabels(): Get the cluster labels of LISA computation.\n",
    "* GetLISAValues(): Get the local spatial autocorrelation values returned from LISA computation.\n",
    "* GetNumNeighbors(): Get the number of neighbors of every observations in LISA computation.\n",
    "* GetPValues(): Get the local pseudo-p values of significance returned from LISA computation.\n",
    "* SetPermutations(num_perm): Set the number of permutations for the LISA computation\n",
    "* SetThreads(num_threads): Set the number of CPU threads for the LISA computation\n",
    "* Run(): Call to run LISA computation\n",
    "\n",
    "For example, we can call the function `GetLISAValues()`\n",
    "to get the values of local Moran:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.015431978309803657, 0.3270633223656033, 0.021295296214118884, 0.004610544790030418, -0.0028342407096540465, 0.41493771583040345, -0.13794630908086175, 0.09986576922564794, 0.2823176310018141, 0.1218745112146858, -0.09512054168698209, 0.032611193818022896, 0.3878324535340113, 1.1888723840162665, -0.6452792226077357, -0.30964927402750314, 0.3662775143008573, 2.0375343538940496, -0.005015479968822494, 0.06971105719113387)\n"
     ]
    }
   ],
   "source": [
    "lms = lisa.GetLISAValues()\n",
    "print(lms[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the pseudo-p values of significance of local Moran computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.414, 0.123, 0.001, 0.474, 0.452, 0.087, 0.243, 0.326, 0.299, 0.303, 0.237, 0.46, 0.258, 0.018, 0.199, 0.188, 0.131, 0.004, 0.456, 0.342)\n"
     ]
    }
   ],
   "source": [
    "pvals = lisa.GetPValues()\n",
    "print(pvals[:20])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the cluster indicators of local Moran computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0)\n"
     ]
    }
   ],
   "source": [
    "cats = lisa.GetClusterIndicators()\n",
    "print(cats[:20])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The predefined values of the indicators of LISA cluster are:\n",
    "```\n",
    "0 Not significant\n",
    "1 High-High\n",
    "2 Low-Low\n",
    "3 High-Low\n",
    "4 Low-High\n",
    "5 Neighborless\n",
    "6 Undefined\n",
    "```\n",
    "which can be accessed via function `GetLabels()`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('Not significant',\n",
       " 'High-High',\n",
       " 'Low-Low',\n",
       " 'High-Low',\n",
       " 'Low-High',\n",
       " 'Undefined',\n",
       " 'Isolated')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lbls = lisa.GetLabels()\n",
    "lbls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Different LISA objects (e.g. local_geary()) will return different\n",
    "    labels and colors.\n",
    "\n",
    "By default, the `local_moran()` function will run with some default parameters, e.g.:\n",
    "```\n",
    "permutation number: 999\n",
    "seed for random number generator: 123456789\n",
    "```\n",
    ", which are identical to GeoDa desktop software so that we can replicate the results as using  GeoDa software. It is also easy to change the paremter and re-run the LISA computation by calling Run() function. \n",
    "\n",
    "For example, re-run the above local Moran example using 9999 permutations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "lisa.SetPermutations(9999)\n",
    "lisa.Run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Then, we can use the same `lisa` object to get the new results after 9999 permutations:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.4187, 0.1265, 0.0004, 0.4679, 0.4545, 0.0728, 0.2312, 0.3071, 0.3115, 0.3088, 0.2187, 0.4803, 0.2623, 0.0113, 0.2, 0.1797, 0.1267, 0.0026, 0.4565, 0.3517)\n"
     ]
    }
   ],
   "source": [
    "pvals = lisa.GetPValues()\n",
    "print(pvals[:20])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "pygeoda uses GeoDa’s C++ code, in which multi-threading is used to accelerate the computation of LISA. We can specify how many threads to run the computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "lisa.SetThreads(4)\n",
    "lisa.Run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the new results after using different number of threads\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.4187, 0.1265, 0.0004, 0.4679, 0.4545, 0.0728, 0.2312, 0.3071, 0.3115, 0.3088, 0.2187, 0.4834, 0.2686, 0.0102, 0.2024, 0.1795, 0.1218, 0.0025, 0.459, 0.3588)\n"
     ]
    }
   ],
   "source": [
    "pvals = lisa.GetPValues()\n",
    "print(pvals[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.2 Local Geary\n",
    "\n",
    "Local Geary is a type of LISA that focuses on squared differences/dissimilarity. A small value of the local geary statistics suggest positive spatial autocorrelation, whereas large values suggest negative spatial autocorrelation. \n",
    "\n",
    "For example, we can call the function `local_geary()` with the created Queen weights and the data “crm_prp” as input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "geary_crmprp = pygeoda.local_geary(queen_w, crm_prp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the cluster indicators of the local Geary computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 2, 4, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2)\n"
     ]
    }
   ],
   "source": [
    "cats = geary_crmprp.GetClusterIndicators()\n",
    "print(cats[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the pseudo-p values of the local Geary computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.398, 0.027, 0.025, 0.126, 0.017, 0.314, 0.61, 0.141, 0.284, 0.11, 0.559, 0.462, 0.211, 0.236, 0.249, 0.229, 0.069, 0.041, 0.205, 0.02)\n"
     ]
    }
   ],
   "source": [
    "pvals = geary_crmprp.GetPValues()\n",
    "print(pvals[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.3 Local Getis-Ord Statistics\n",
    "\n",
    "There are two types of local Getis-Ord statistics: one is computing a ratio of the weighted average of the values in the neighboring locations, not including the value at the location; while another type of statistic includes the value at the location in both numerator and denominator.\n",
    "\n",
    "A value larger than the mean suggests a high-high cluster or hot spot, a value smaller than the mean indicates a low-low cluster or cold spot.\n",
    "\n",
    "For example, we can call the function `local_g()` with the created Queen weights and the data “crm_prp” as input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "localg_crmprp = pygeoda.local_g(queen_w, crm_prp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the cluster indicators of the local G computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0)\n"
     ]
    }
   ],
   "source": [
    "cats = localg_crmprp.GetClusterIndicators()\n",
    "print(cats[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the pseudo-p values of the local G computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.414, 0.123, 0.001, 0.474, 0.452, 0.087, 0.243, 0.326, 0.299, 0.303, 0.237, 0.46, 0.258, 0.018, 0.199, 0.188, 0.131, 0.004, 0.456, 0.342)\n"
     ]
    }
   ],
   "source": [
    "pvals = localg_crmprp.GetPValues()\n",
    "print(pvals[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the second type of local Getis-Ord statistics, we can call the function `local_gstar()` with the created Queen weights and the data “crm_prp” as input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "localgstar_crmprp = pygeoda.local_gstar(queen_w, crm_prp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.4 Local Join Count\n",
    "\n",
    "\n",
    "Local Join Count is a method to identify local clusters for binary data by using a local version of the so-called BB join count statistic. The statistic is only meaningful for those observations with value 1. \n",
    "\n",
    "For example, we can load the columbus dataset, and call the function `local_joincount()` with a Queen weights and the data “nsa”, which is a set of binary (0,1) values, as input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "columbus = pygeoda.open('./data/columbus.shp')\n",
    "columbus_w = pygeoda.weights.queen(columbus)\n",
    "localjc_nsa = pygeoda.local_joincount(columbus_w, columbus.GetIntegerCol('nsa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the cluster indicators of the local Join Count computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0)\n"
     ]
    }
   ],
   "source": [
    "cats = localjc_nsa.GetClusterIndicators()\n",
    "print(cats[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the pseudo-p values of the local Join Count computation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.213, 0.07, 0.017, 0.024, 0.001, 0.226, 0.055, 0.007, 0.006, 0.275, 0.017, 0.008, 0.043, 0.004, 0.002, 0.001, 0.147, 0.049, 0.087, 0.001)\n"
     ]
    }
   ],
   "source": [
    "pvals = localjc_nsa.GetPValues()\n",
    "print(pvals[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the number of neighbors of the local Join Count computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 3, 4, 4, 8, 2, 4, 6, 8, 4, 5, 6, 4, 6, 6, 8, 3, 4, 3, 10)\n"
     ]
    }
   ],
   "source": [
    "nn = localjc_nsa.GetNumNeighbors()\n",
    "print(nn[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### 3.2 Spatial Clustering\n",
    "\n",
    " \n",
    "\n",
    "Spatial clustering aims to group of a large number of\n",
    "geographic areas or points into a smaller number of regions\n",
    "based on similiarities in one or more variables.\n",
    "Spatially constrained clustering is needed when clusters are\n",
    "required to be spatially contiguous.\n",
    "\n",
    "In pygeoda v0.0.3, there are three different approaches\n",
    "explicitly incorporate the contiguity constraint in the\n",
    "optimization process: SKATER, Redcap and Max-p.\n",
    "More more details, please read the lab note that\n",
    "Dr. Luc Anselin wrote:\n",
    "http://geodacenter.github.io/workbook/8_spatial_clusters/lab8.html \n",
    "\n",
    "For example, to apply spatial clustering on the Guerry dataset,\n",
    "we use the queen weights to define the spatial contiguity\n",
    "and select 6 variables for similarity measure:\n",
    "\"Crm_prs\", \"Crm_prp\", \"Litercy\", \"Donatns\", \"Infants\", \"Suicids\".\n",
    "\n",
    "The following code is used to get a 2D data list for\n",
    "the selected variables:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "select_vars = [\"Crm_prs\", \"Crm_prp\", \"Litercy\", \"Donatns\", \"Infants\", \"Suicids\"]\n",
    "data = [guerry.GetRealCol(v) for v in select_vars]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.1 SKATER\n",
    "\n",
    "The Spatial C(K)luster Analysis by Tree Edge Removal(SKATER)\n",
    "algorithm introduced by Assuncao et al. (2006) is based on\n",
    "the optimal pruning of a minimum spanning tree that reflects\n",
    "the contiguity structure among the observations. It provides\n",
    "an optimized algorithm to prune to tree into several clusters\n",
    "that their values of selected variables are as similar as possible.\n",
    "\n",
    "The SKATER function in pygeoda:\n",
    "::\n",
    "\n",
    "    pygeoda.skater(k, w, data, distance_method='euclidean', bound_vals = [],  min_bound = 0, random_seed=123456789)\n",
    "\n",
    "\n",
    "Note: the parameters `distance_method`, `bound_vals`, `min_bound` and `random_seed` are optional.\n",
    "\n",
    "Note: See [Max-p] for the usage of  `bound_vals` and `min_bound`.\n",
    "\n",
    "For example, to create 4 spatially contiguous clusters using\n",
    "`skater()` with Guerry dataset, the queen weights and the\n",
    "values of the 6 selected variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((15, 74, 16, 55, 60, 39, 68, 33, 17, 82, 81, 0, 2, 40, 20, 80), (46, 50, 34, 38, 69, 47, 58, 19, 32, 41, 53, 26), (23, 79, 3, 29, 61, 21, 44, 11, 28, 13, 30, 35, 76, 77, 43, 9, 27, 45, 31, 78, 4, 10, 66, 37, 5, 14, 7, 63, 62), (49, 52, 72, 84, 8, 57, 56, 59, 42, 1, 25, 51, 48, 54, 64, 75, 18, 83, 73, 36, 24, 71, 6, 67, 65, 70, 22, 12))\n"
     ]
    }
   ],
   "source": [
    "skater_clusters = pygeoda.skater(4, queen_w, data)\n",
    "print(skater_clusters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `skater()` function returns a 2D tuple, which represents 4 clusters. Each cluster is composed by several contiguity areas, e.g. 15, 74, 16, 55, 60, 39, 68, 33, 17, 82, 81, 0, 2, 40, 20, 80\n",
    "\n",
    "pygeoda also provides utility functions to compute some descriptive statistics of the clustering results, e.g. to compute the ratio of between to total sum of squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The ratio of between to total sum of square: 0.3156446659311204\n"
     ]
    }
   ],
   "source": [
    "betweenss = pygeoda.between_sumofsquare(skater_clusters, data)\n",
    "totalss = pygeoda.total_sumofsquare( data)\n",
    "ratio =  betweenss / totalss\n",
    "print(\"The ratio of between to total sum of square:\", ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.2 REDCAP\n",
    "\n",
    "REDCAP (Regionalization with dynamically constrained agglomerative clustering and partitioning) is developed by D. Guo (2008). Like SKATER, REDCAP starts from building a spanning tree with 3 different ways (single-linkage, average-linkage, and the complete-linkage). The single-linkage way leads to build a minimum spanning tree. Then,REDCAP provides 2 different ways (first‐order and full-order constraining) to prune the tree to find clusters. The first-order approach with a minimum spanning tree is exactly the same with SKATER. In GeoDa and pygeoda, the following methods are provided:\n",
    "\n",
    "* First-order and Single-linkage\n",
    "* Full-order and Complete-linkage\n",
    "* Full-order and Average-linkage\n",
    "* Full-order and Single-linkage\n",
    "\n",
    "For example, to find same 4 clusters using the same dataset and weights as above using Full-order and Complete-linkage REDCAP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((13, 60, 76, 77, 43, 39, 10, 5, 2), (20,), (40,), (49, 52, 46, 50, 34, 38, 72, 84, 23, 79, 57, 56, 59, 29, 61, 1, 15, 74, 21, 44, 11, 3, 25, 42, 51, 8, 48, 54, 64, 16, 55, 18, 28, 35, 9, 27, 45, 69, 31, 47, 58, 36, 83, 19, 32, 24, 71, 6, 73, 41, 78, 75, 68, 33, 65, 70, 66, 53, 30, 37, 17, 82, 81, 0, 22, 14, 67, 63, 4, 26, 7, 62, 12, 80))\n",
      "The ratio of between to total sum of square: 0.1905641377254551\n"
     ]
    }
   ],
   "source": [
    "redcap_clusters = pygeoda.redcap(4, queen_w, data, \"fullorder-completelinkage\")\n",
    "print(redcap_clusters)\n",
    "\n",
    "betweenss = pygeoda.between_sumofsquare(redcap_clusters, data)\n",
    "totalss = pygeoda.total_sumofsquare( data)\n",
    "ratio =  betweenss / totalss\n",
    "print(\"The ratio of between to total sum of square:\", ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.3 Max-p\n",
    "\n",
    "The so-called max-p regions model (outlined in Duque, Anselin, and Rey 2012) uses a different approach and considers the regionalization problem as an application of integer programming. In addition, the number of regions is determined endogenously.\n",
    "\n",
    "The algorithm itself consists of a search process that starts with an initial feasible solution and iteratively improves upon it while maintaining contiguity among the elements of each cluster. Like Geoda, pygeoda provides three different heuristic algorithms to find an optimal solution for max-p:\n",
    "\n",
    "* greedy\n",
    "* Tabu Search\n",
    "* Simulated Annealing\n",
    "\n",
    "Unlike SKATER and REDCAP that one can specify the number of clusters as an input paramter, max-p doesn't take this parameter but the constrained variable and the minimum value that each cluster should reach. \n",
    "\n",
    "For example, we use `greedy` algorithm in maxp function with the same dataset and weights as above to find optimal clusters using max-p. \n",
    "\n",
    "First, we need to specify, for example, every cluster must have population >= 3236.67 thousands people:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(346.03, 513.0, 298.26, 155.9, 129.1, 340.73, 289.62, 253.12, 246.36, 270.13, 359.06, 359.47, 494.7, 258.59, 362.53, 445.25, 256.06, 294.83, 375.88, 598.87, 265.38, 482.75, 265.54, 299.56, 424.25, 278.82, 524.4, 357.38, 427.86, 312.16, 554.23, 346.3, 547.05, 245.29, 297.02, 550.26, 312.5, 281.5, 235.75, 391.22, 292.08, 470.09, 305.28, 283.83, 346.89, 140.35, 467.87, 591.28, 337.08, 249.83, 352.59, 415.57, 314.59, 433.52, 417.0, 282.52, 989.94, 397.73, 441.88, 655.22, 573.11, 428.4, 233.03, 157.05, 540.21, 424.26, 434.43, 338.91, 523.97, 457.37, 935.11, 693.68, 323.89, 448.18, 297.85, 543.7, 333.84, 242.51, 317.5, 239.11, 330.36, 282.73, 285.13, 397.99, 352.49)\n"
     ]
    }
   ],
   "source": [
    "bound_vals = guerry.GetRealCol(\"Pop1831\")\n",
    "min_bound = 3236.67 # 10% of Pop1831\n",
    "print(bound_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can call the max-p function with \"greedy\" algorithm,\n",
    "the bound values and minimum bound value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The ratio of between to total sum of square: 0.507018079733202\n"
     ]
    }
   ],
   "source": [
    "maxp_clusters = pygeoda.maxp(queen_w, data, bound_vals, min_bound, \"greedy\")\n",
    "\n",
    "betweenss = pygeoda.between_sumofsquare(maxp_clusters, data)\n",
    "ratio = betweenss / totalss\n",
    "print(\"The ratio of between to total sum of square:\", ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also specify using `tabu search` algorithm in maxp function\n",
    "with the parameter of tabu length:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The ratio of between to total sum of square: 0.5280299052291919\n"
     ]
    }
   ],
   "source": [
    "maxp_tabu_clusters = pygeoda.maxp(queen_w, data, bound_vals, min_bound, \"tabu\", tabu_length=95)\n",
    "\n",
    "betweenss = pygeoda.between_sumofsquare(maxp_tabu_clusters, data)\n",
    "ratio = betweenss / totalss\n",
    "print(\"The ratio of between to total sum of square:\", ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply `simulated annealing` algorithm in maxp function with\n",
    "the parameter of cooling rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The ratio of between to total sum of square: 0.5260248902417842\n"
     ]
    }
   ],
   "source": [
    "maxp_sa_clusters = pygeoda.maxp(queen_w, data, bound_vals, min_bound, \"sa\", cool_rate=0.75)\n",
    "\n",
    "betweenss = pygeoda.between_sumofsquare(maxp_sa_clusters, data)\n",
    "ratio = betweenss / totalss\n",
    "print(\"The ratio of between to total sum of square:\", ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also increase the number of iterations for local search process by specifying the parameter `initial` (default value is 99):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The ratio of between to total sum of square: 0.5260248902417843\n"
     ]
    }
   ],
   "source": [
    "maxp_clusters = pygeoda.maxp(queen_w, data, bound_vals, min_bound, \"greedy\", initial=1000)\n",
    "\n",
    "betweenss = pygeoda.between_sumofsquare(maxp_clusters, data)\n",
    "ratio = betweenss / totalss\n",
    "print(\"The ratio of between to total sum of square:\", ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "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.2"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 2
}