{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Check Random Sampling\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom matplotlib import pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\n\nimport gstools as gs\n\n\ndef norm_rad(vec):\n    \"\"\"Direction on the unit sphere.\"\"\"\n    vec = np.array(vec, ndmin=2)\n    norm = np.zeros(vec.shape[1])\n    for i in range(vec.shape[0]):\n        norm += vec[i] ** 2\n    norm = np.sqrt(norm)\n    return np.einsum(\"j,ij->ij\", 1 / norm, vec), norm\n\n\ndef plot_rand_meth_samples(generator):\n    \"\"\"Plot the samples of the rand meth class.\"\"\"\n    norm, rad = norm_rad(generator._cov_sample)\n\n    fig = plt.figure(figsize=(10, 4))\n\n    if generator.model.dim == 3:\n        ax = fig.add_subplot(121, projection=Axes3D.name)\n        u = np.linspace(0, 2 * np.pi, 100)\n        v = np.linspace(0, np.pi, 100)\n        x = np.outer(np.cos(u), np.sin(v))\n        y = np.outer(np.sin(u), np.sin(v))\n        z = np.outer(np.ones(np.size(u)), np.cos(v))\n        ax.plot_surface(x, y, z, rstride=4, cstride=4, color=\"b\", alpha=0.1)\n        ax.scatter(norm[0], norm[1], norm[2])\n    elif generator.model.dim == 2:\n        ax = fig.add_subplot(121)\n        u = np.linspace(0, 2 * np.pi, 100)\n        x = np.cos(u)\n        y = np.sin(u)\n        ax.plot(x, y, color=\"b\", alpha=0.1)\n        ax.scatter(norm[0], norm[1])\n        ax.set_aspect(\"equal\")\n    else:\n        ax = fig.add_subplot(121)\n        ax.bar(-1, np.sum(np.isclose(norm, -1)), color=\"C0\")\n        ax.bar(1, np.sum(np.isclose(norm, 1)), color=\"C0\")\n        ax.set_xticks([-1, 1])\n        ax.set_xticklabels((\"-1\", \"1\"))\n    ax.set_title(\"Direction sampling\")\n\n    ax = fig.add_subplot(122)\n    x = np.linspace(0, 10 / generator.model.integral_scale)\n    y = generator.model.spectral_rad_pdf(x)\n    ax.plot(x, y, label=\"radial spectral density\")\n    sample_in = np.sum(rad <= np.max(x))\n    ax.hist(rad[rad <= np.max(x)], bins=sample_in // 50, density=True)\n    ax.set_xlim([0, np.max(x)])\n    ax.set_title(\"Radius samples shown {}/{}\".format(sample_in, len(rad)))\n    ax.legend()\n    plt.show()\n\n\nmodel = gs.Stable(dim=3, alpha=1.5)\nsrf = gs.SRF(model, seed=2020)\nplot_rand_meth_samples(srf.generator)"
      ]
    }
  ],
  "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}