{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "a6f0700d",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "7e5da80c",
            "cell_type": "markdown",
            "source": "Halfspace intersection of planes forming some polygon\n",
            "metadata": {}
        },
        {
            "id": "2b4e66eb",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.spatial import HalfspaceIntersection\nimport numpy as np\nhalfspaces = np.array([[-1, 0., 0.],\n                       [0., -1., 0.],\n                       [2., 1., -4.],\n                       [-0.5, 1., -2.]])\nfeasible_point = np.array([0.5, 0.5])\nhs = HalfspaceIntersection(halfspaces, feasible_point)",
            "outputs": []
        },
        {
            "id": "8038eca3",
            "cell_type": "markdown",
            "source": "Plot halfspaces as filled regions and intersection points:\n",
            "metadata": {}
        },
        {
            "id": "2a667d28",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nfig = plt.figure()\nax = fig.add_subplot(1, 1, 1, aspect='equal')\nxlim, ylim = (-1, 3), (-1, 3)\nax.set_xlim(xlim)\nax.set_ylim(ylim)\nx = np.linspace(-1, 3, 100)\nsymbols = ['-', '+', 'x', '*']\nsigns = [0, 0, -1, -1]\nfmt = {\"color\": None, \"edgecolor\": \"b\", \"alpha\": 0.5}\nfor h, sym, sign in zip(halfspaces, symbols, signs):\n    hlist = h.tolist()\n    fmt[\"hatch\"] = sym\n    if h[1]== 0:\n        ax.axvline(-h[2]/h[0], label='{}x+{}y+{}=0'.format(*hlist))\n        xi = np.linspace(xlim[sign], -h[2]/h[0], 100)\n        ax.fill_between(xi, ylim[0], ylim[1], **fmt)\n    else:\n        ax.plot(x, (-h[2]-h[0]*x)/h[1], label='{}x+{}y+{}=0'.format(*hlist))\n        ax.fill_between(x, (-h[2]-h[0]*x)/h[1], ylim[sign], **fmt)\nx, y = zip(*hs.intersections)\nax.plot(x, y, 'o', markersize=8)",
            "outputs": []
        },
        {
            "id": "6d4ea39a",
            "cell_type": "markdown",
            "source": "By default, qhull does not provide with a way to compute an interior point.\nThis can easily be computed using linear programming. Considering halfspaces\nof the form $Ax + b \\leq 0$, solving the linear program:\n\n$$ max \\: y $$\n$$ s.t. Ax + y ||A_i|| \\leq -b $$\nWith $A_i$ being the rows of A, i.e. the normals to each plane.\n\nWill yield a point x that is furthest inside the convex polyhedron. To\nbe precise, it is the center of the largest hypersphere of radius y\ninscribed in the polyhedron. This point is called the Chebyshev center\nof the polyhedron (see [1] 4.3.1, pp148-149). The\nequations outputted by Qhull are always normalized.\n",
            "metadata": {}
        },
        {
            "id": "a049b92f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.optimize import linprog\nfrom matplotlib.patches import Circle\nnorm_vector = np.reshape(np.linalg.norm(halfspaces[:, :-1], axis=1),\n    (halfspaces.shape[0], 1))\nc = np.zeros((halfspaces.shape[1],))\nc[-1] = -1\nA = np.hstack((halfspaces[:, :-1], norm_vector))\nb = - halfspaces[:, -1:]\nres = linprog(c, A_ub=A, b_ub=b, bounds=(None, None))\nx = res.x[:-1]\ny = res.x[-1]\ncircle = Circle(x, radius=y, alpha=0.3)\nax.add_patch(circle)\nplt.legend(bbox_to_anchor=(1.6, 1.0))\nplt.show()",
            "outputs": []
        }
    ]
}