{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Monte Carlo Method\n",
        "\n",
        "As a simple example of a Monte Carlo method, we will approximate the value of $\\pi$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pyopencl as cl\n",
        "import pyopencl.array\n",
        "import pyopencl.clrandom"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "ctx = cl.create_some_context()\n",
        "queue = cl.CommandQueue(ctx)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Boilerplate for Random Number Generator"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "generator_preamble = \"\"\"\n",
        "#include <pyopencl-random123/philox.cl>\n",
        "\n",
        "typedef union {\n",
        "    uint4 v;\n",
        "    philox4x32_ctr_t c;\n",
        "} philox4x32_ctr_vec_union;\n",
        "\n",
        "\n",
        "uint4 philox4x32_bump(uint4 ctr)\n",
        "{\n",
        "    if (++ctr.x == 0)\n",
        "        if (++ctr.y == 0)\n",
        "            ++ctr.z;\n",
        "    return ctr;\n",
        "}\n",
        "\n",
        "uint4 philox4x32_gen(\n",
        "        uint4 ctr,\n",
        "        uint2 key,\n",
        "        uint4 *new_ctr)\n",
        "{\n",
        "    philox4x32_ctr_vec_union result;\n",
        "    result.c = philox4x32(\n",
        "        *(philox4x32_ctr_t *) &ctr,\n",
        "        *(philox4x32_key_t *) &key);\n",
        "    *new_ctr = philox4x32_bump(ctr);\n",
        "    return result.v;\n",
        "}\n",
        "\n",
        "float4 philox4x32_f32(\n",
        "        uint4 ctr,\n",
        "        uint2 key,\n",
        "        uint4 *new_ctr)\n",
        "{\n",
        "    *new_ctr = ctr;\n",
        "    return\n",
        "        convert_float4(philox4x32_gen(*new_ctr, key, new_ctr))\n",
        "        * 2.3283064365386963e-10f;\n",
        "}\n",
        "\"\"\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Reduction Code"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Complete the sampler code:\n",
        "\n",
        "```\n",
        "mc_preamble_src = \"\"\"\n",
        "\n",
        "#include <pyopencl-complex.h>\n",
        "\n",
        "float compute_sample(int i, unsigned int k1)\n",
        "{\n",
        "    uint4 ctr = { 0, 1, 2, 3 };\n",
        "    uint2 key2 = { i, k1 };\n",
        "    float4 rng_res = philox4x32_f32(ctr, key2, &(ctr));\n",
        "    ...\n",
        "}\n",
        "\"\"\"\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "mc_preamble_src = \"\"\"\n",
        "\n",
        "#include <pyopencl-complex.h>\n",
        "\n",
        "float compute_sample(int i, unsigned int k1)\n",
        "{\n",
        "    uint4 ctr = { 0, 1, 2, 3 };\n",
        "    uint2 key2 = { i, k1 };\n",
        "    float4 rng_res = philox4x32_f32(ctr, key2, &(ctr));\n",
        "    \n",
        "    cfloat_t samp0 = cfloat_new(rng_res.s0, rng_res.s1);\n",
        "    cfloat_t samp1 = cfloat_new(rng_res.s2, rng_res.s3);\n",
        "    \n",
        "    float result = 0;\n",
        "    if (cfloat_abs(samp0) <= 1)\n",
        "        result += 1;\n",
        "    if (cfloat_abs(samp1) <= 1)\n",
        "        result += 1;\n",
        "        \n",
        "    return result;\n",
        "}\n",
        "\"\"\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "from pyopencl.reduction import ReductionKernel"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Syntax:\n",
        "\n",
        "`ReductionKernel(context, dtype, netural, reduce_expr, map_expr, arguments)`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rknl = ReductionKernel(ctx, np.float32,\n",
        "        neutral=\"0\",\n",
        "        reduce_expr=\"a+b\", map_expr=\"compute_sample(i, k1)\",\n",
        "        arguments=\"unsigned int k1\", preamble=generator_preamble+mc_preamble_src)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "3.14154656\n"
          ]
        }
      ],
      "source": [
        "n = 100000000\n",
        "\n",
        "nsamples_accepted = rknl(15, range=slice(n), queue=queue).get()\n",
        "nsamples = 2*n\n",
        "approx_pi = 4 * nsamples_accepted/nsamples\n",
        "\n",
        "print(approx_pi)"
      ]
    }
  ],
  "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.5.1+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}