pycuda-stack-gp-5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline \n",
    "\n",
    "import numpy as np\n",
    "import numpy\n",
    "import pandas as pd\n",
    "import re\n",
    "import random as rnd\n",
    "from sklearn.metrics import accuracy_score\n",
    "import math\n",
    "import operator\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "import random\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pycuda.driver as cuda\n",
    "import pycuda.autoinit\n",
    "from pycuda.compiler import SourceModule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the arrays\n",
    "a = [float(x) for x in open('titanic_cpp_x_old.txt').read().split()]\n",
    "cols = int(a[0])\n",
    "bdx = array(a[1:]).reshape(-1, cols)\n",
    "bdx.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = [float(x) for x in open('titanic_cpp_test_x_old.txt').read().split()]\n",
    "cols = int(a[0])\n",
    "bdx_test = array(a[1:]).reshape(-1, cols)\n",
    "bdx_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdy = array([float(x) for x in open('titanic_cpp_y.txt').read().split()])\n",
    "bdy.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdy_test = array([float(x) for x in open('titanic_cpp_test_y.txt').read().split()])\n",
    "bdy_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdx.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data organization\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_insts = 19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_rows = 1024\n",
    "num_features = bdx.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_threads_x = 16384*4\n",
    "num_threads_y = num_rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prog_len = 32\n",
    "gprogs = np.zeros((num_threads_x, prog_len*2), dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. The data is a single array of 8 features, many rows, and each thread on the Y axis takes an individual row\n",
    "gdata = np.zeros((num_rows, num_features), dtype=np.float32)\n",
    "gdatay = np.zeros((num_rows, 1), dtype=np.float32)\n",
    "gdatatest = np.zeros((num_rows, num_features), dtype=np.float32)\n",
    "gdatatesty = np.zeros((num_rows, 1), dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3. The results - for each thread, single column of results\n",
    "gresult = np.zeros((num_rows, 1*num_threads_x), dtype=np.float32)\n",
    "gscores = np.zeros((1, 1*num_threads_x), dtype=np.float32)\n",
    "gresulttest = np.zeros((num_rows, 1*num_threads_x), dtype=np.float32)\n",
    "gscorestest = np.zeros((1, 1*num_threads_x), dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_idx = np.zeros((1,num_threads_x), dtype=np.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Allocate the buffers in GPU memory\n",
    "gprogs_gpu = cuda.mem_alloc(gprogs.nbytes)\n",
    "gprogsnew_gpu = cuda.mem_alloc(gprogs.nbytes) # to hold next generation\n",
    "sorted_idx_gpu = cuda.mem_alloc(sorted_idx.nbytes)\n",
    "\n",
    "gdata_gpu = cuda.mem_alloc(gdata.nbytes)\n",
    "gdatay_gpu = cuda.mem_alloc(gdatay.nbytes)\n",
    "gresult_gpu = cuda.mem_alloc(gresult.nbytes)\n",
    "gscores_gpu = cuda.mem_alloc(gscores.nbytes)\n",
    "\n",
    "gdatatest_gpu = cuda.mem_alloc(gdatatest.nbytes)\n",
    "gdatatesty_gpu = cuda.mem_alloc(gdatatesty.nbytes)\n",
    "gresulttest_gpu = cuda.mem_alloc(gresulttest.nbytes)\n",
    "gscorestest_gpu = cuda.mem_alloc(gscorestest.nbytes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threads_per_block = 16\n",
    "blocks_x = num_threads_x // threads_per_block\n",
    "blocks_y = num_threads_y // threads_per_block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mutation_rate = 0.1\n",
    "mutation_prob = 0.5\n",
    "crossover_prob = 0.5\n",
    "survival_rate = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The GP kernel\n",
    "s = f\"\"\"\n",
    "#define PROGLEN {prog_len}\n",
    "#define NUM_FEATURES {num_features}\n",
    "#define NUM_THREADS_X {num_threads_x}\n",
    "\"\"\"\n",
    "\n",
    "sf = f\"\"\"\n",
    "#define MUTATION_RATE {mutation_rate}\n",
    "#define MUTATION_PROB {mutation_prob}\n",
    "#define CROSSOVER_PROB {crossover_prob}\n",
    "#define SURVIVAL_RATE {survival_rate}\n",
    "\"\"\"\n",
    "\n",
    "mod = SourceModule(s + sf + \"\"\"\n",
    "\n",
    "#include <curand_kernel.h>\n",
    "\n",
    "#define NOP 0\n",
    "#define PUSHV 1\n",
    "#define PUSHC 2\n",
    "#define ADD 3\n",
    "#define MUL 4\n",
    "#define DIV 5\n",
    "#define NEG 6\n",
    "#define MIN 7\n",
    "#define MAX 8\n",
    "#define GREATER 9\n",
    "#define LESS 10\n",
    "#define EQUAL 11\n",
    "#define SIN 12\n",
    "#define COS 13\n",
    "#define EXP 14\n",
    "#define LOG 15\n",
    "#define SQR 16\n",
    "#define SQRT 17\n",
    "#define TANH 18\n",
    "\n",
    "\n",
    "\n",
    "const int nstates = NUM_THREADS_X;\n",
    "__device__ curandState_t* states[nstates];\n",
    "\n",
    "extern \"C\" {\n",
    "\n",
    "__global__ void initkernel(int seed)\n",
    "{\n",
    "    int tidx = threadIdx.x + blockIdx.x * blockDim.x;\n",
    "\n",
    "    if (tidx < nstates) \n",
    "    {\n",
    "        curandState_t* s = new curandState_t;\n",
    "        if (s != 0) \n",
    "        {\n",
    "            curand_init(seed, tidx, 0, s);\n",
    "        }\n",
    "\n",
    "        states[tidx] = s;\n",
    "    }\n",
    "}\n",
    "\n",
    "__device__ float randfloat(int tidx)\n",
    "{\n",
    "    curandState_t s = *states[tidx];\n",
    "    float x = curand_uniform(&s);\n",
    "    *states[tidx] = s;\n",
    "    return x;\n",
    "}\n",
    "\n",
    "__device__ unsigned int randint(int tidx, int aX, int aY)\n",
    "{\n",
    "    curandState_t s = *states[tidx];\n",
    "    unsigned int x = curand(&s);\n",
    "    *states[tidx] = s;\n",
    "    \n",
    "    if (aX == aY)\n",
    "    {\n",
    "        return aX;\n",
    "    }\n",
    "    if (aX == (aY-1))\n",
    "    {\n",
    "        // for two consecutives, pick either with equal probability\n",
    "        if (randfloat(tidx) < 0.5)\n",
    "        {\n",
    "            return aX;\n",
    "        }\n",
    "        else\n",
    "        {\n",
    "            return aY;\n",
    "        }\n",
    "    }\n",
    "    return aX + (x % (aY - aX + 1));\n",
    "}\n",
    "\n",
    "\n",
    "__global__ void init_progs(float *gprogs)\n",
    "{\n",
    "    int tidx = (blockIdx.x * blockDim.x + threadIdx.x);\n",
    "    int prog_idx =  tidx * PROGLEN * 2;\n",
    "    float* start_idx = (gprogs + prog_idx);\n",
    "    \n",
    "    for(int ip=0; ip<(PROGLEN*2); ip += 2)\n",
    "    {\n",
    "        // random instruction\n",
    "        float inst = (float)(randint(tidx, 0, TANH));\n",
    "        float operand = 0;\n",
    "        if (inst == PUSHV)\n",
    "        {\n",
    "            operand = (float)(randint(tidx, 0, NUM_FEATURES-1));\n",
    "        }\n",
    "        if (inst == PUSHC)\n",
    "        {\n",
    "            operand = (randfloat(tidx) - 0.5) * 5;\n",
    "        }\n",
    "        \n",
    "        \n",
    "        *(start_idx + ip) = inst; // instruction\n",
    "        *(start_idx + ip + 1) = operand; // operand\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "// crossover\n",
    "__device__ void crossover(int tidx, float* mom, float* dad, float* baby)\n",
    "{\n",
    "    unsigned int p1 = randint(tidx, 0, PROGLEN - 1) * 2;\n",
    "\tunsigned int p2 = randint(tidx, 0, PROGLEN - 1) * 2;\n",
    "    \n",
    "\twhile (p1 == p2)\n",
    "\t{\n",
    "        p1 = randint(tidx, 0, PROGLEN - 1) * 2;\n",
    "        p2 = randint(tidx, 0, PROGLEN - 1) * 2;\n",
    "\t}\n",
    "\n",
    "\tif (p1 > p2)\n",
    "\t{\n",
    "\t\tint tmp = p1;\n",
    "\t\tp1 = p2;\n",
    "\t\tp2 = tmp;\n",
    "\t}\n",
    "    \n",
    "    float x;\n",
    "    \n",
    "    for (unsigned int i = 0; i < p1; i += 2)\n",
    "\t{\n",
    "        x = mom[i];\n",
    "        baby[i] = x;\n",
    "        x =  mom[i+1];\n",
    "        baby[i+1] = x;\n",
    "\t}\n",
    "    \n",
    "\tfor (unsigned int i = p1; i < p2; i += 2)\n",
    "\t{\n",
    "        x = dad[i];\n",
    "        baby[i] = x;\n",
    "        x =  dad[i+1];\n",
    "        baby[i+1] = x;\n",
    "\t}\n",
    "    \n",
    "\tfor (unsigned int i = p2; i < PROGLEN*2; i += 2)\n",
    "\t{\n",
    "        x = mom[i];\n",
    "        baby[i] = x;\n",
    "        x =  mom[i+1];\n",
    "        baby[i+1] = x;\n",
    "\t}\n",
    "}\n",
    "\n",
    "\n",
    "// mutation\n",
    "__device__ void mutate(int tidx, float* baby)\n",
    "{\n",
    "    for (int i = 0; i < PROGLEN*2; i += 2)\n",
    "\t{\n",
    "\t\tif (randfloat(tidx) < MUTATION_RATE)\n",
    "\t\t{\n",
    "\t\t\tif ((baby[i] == PUSHV) && (randfloat(tidx) < 0.5))\n",
    "\t\t\t{\n",
    "\t\t\t\tbaby[i + 1] = (float)(randint(tidx, 0, NUM_FEATURES - 1));\n",
    "\t\t\t}\n",
    "\t\t\telse if ((baby[i] == PUSHC) && (randfloat(tidx) < 0.5))\n",
    "\t\t\t{\n",
    "\t\t\t\tbaby[i + 1] = (randfloat(tidx) - 0.5) * 5;\n",
    "\t\t\t}\n",
    "\t\t\telse\n",
    "\t\t\t{\n",
    "                float inst = (float)(randint(tidx, 0, TANH));\n",
    "                float operand = 0;\n",
    "                if (inst == PUSHV)\n",
    "                {\n",
    "                    operand = (float)(randint(tidx, 0, NUM_FEATURES-1));\n",
    "                }\n",
    "                if (inst == PUSHC)\n",
    "                {\n",
    "                    operand = (randfloat(tidx) - 0.5) * 5;\n",
    "                }\n",
    "                baby[i] = inst;\n",
    "                baby[i+1] = operand;\n",
    "\t\t\t}\n",
    "\t\t}\n",
    "\t}\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "// will produce one baby per thread and store it in gprogsnew\n",
    "__global__ void reproduce(int* sorted_idx, float *gprogs, float *gprogsnew)\n",
    "{\n",
    "    int tidx = (blockIdx.x * blockDim.x + threadIdx.x);\n",
    "    int prog_idx =  tidx * PROGLEN * 2;\n",
    "    float* baby_start_idx = (gprogsnew + prog_idx);\n",
    "    \n",
    "    // elitism\n",
    "    if (tidx == 0)\n",
    "    {\n",
    "        int idx = sorted_idx[0];\n",
    "        float* mom = gprogs + idx*PROGLEN*2;\n",
    "        \n",
    "        for(int i=0; i<PROGLEN*2; i++)\n",
    "        {\n",
    "            baby_start_idx[i] = mom[i];\n",
    "        }\n",
    "        \n",
    "        return;\n",
    "    }\n",
    "    \n",
    "    // decide whether to do crossover or not\n",
    "    if (randfloat(tidx) < CROSSOVER_PROB)\n",
    "    {\n",
    "        // will do crossover\n",
    "        int idx1 = sorted_idx[randint(tidx, 0, (int)((float)(NUM_THREADS_X) * SURVIVAL_RATE))];\n",
    "        int idx2 = sorted_idx[randint(tidx, 0, (int)((float)(NUM_THREADS_X) * SURVIVAL_RATE))];\n",
    "        \n",
    "        float* mom = gprogs + idx1*PROGLEN*2;\n",
    "        float* dad = gprogs + idx2*PROGLEN*2;\n",
    "        \n",
    "        crossover(tidx, mom, dad, baby_start_idx);\n",
    "        mutate(tidx, baby_start_idx);\n",
    "    }\n",
    "    else \n",
    "    {\n",
    "        // will only mutate, but copy the baby first\n",
    "        int idx = sorted_idx[randint(tidx, 0, (int)((float)(NUM_THREADS_X) * SURVIVAL_RATE))];\n",
    "        float* mom = gprogs + idx*PROGLEN*2;\n",
    "        \n",
    "        for(int i=0; i<PROGLEN*2; i++)\n",
    "        {\n",
    "            baby_start_idx[i] = mom[i];\n",
    "        }\n",
    "        \n",
    "        if (randfloat(tidx) < MUTATION_PROB)\n",
    "        {\n",
    "            mutate(tidx, baby_start_idx);\n",
    "        }\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "__global__ void compute(int num_real_rows, float *gprogs, float *gdata, float *gresult)\n",
    "{\n",
    "    int prog_idx = (blockIdx.x * blockDim.x + threadIdx.x) * PROGLEN * 2;\n",
    "    int data_row_idx = (blockIdx.y * blockDim.y + threadIdx.y) * NUM_FEATURES;\n",
    "\n",
    "    int th_y = blockIdx.y * blockDim.y + threadIdx.y;\n",
    "    int th_x = blockIdx.x * blockDim.x + threadIdx.x;\n",
    "\n",
    "    int result_idx = NUM_THREADS_X * th_y + th_x; \n",
    "\n",
    "    if (th_y >= num_real_rows)\n",
    "        return; // don't compute where there's no data\n",
    "\n",
    "    float stack[PROGLEN+4]; // a little bigger just in case\n",
    "\n",
    "    int sp = 0;\n",
    "\n",
    "    // push a 0 to the stack\n",
    "    stack[sp] = 0;\n",
    "    sp++;\n",
    "\n",
    "    float x = 0;\n",
    "    float y = 0;\n",
    "\n",
    "    float inst = 0;\n",
    "    float operand = 0;\n",
    "\n",
    "    /////////////////////////////////\n",
    "    // Evaluate program on data row\n",
    "    float* start_idx = (gprogs + prog_idx);\n",
    "\n",
    "    for(int ip=0; ip<(PROGLEN*2); ip += 2)\n",
    "    {\n",
    "        // fetch instruction and operand\n",
    "        inst = *(start_idx + ip);\n",
    "        operand = *(start_idx + ip + 1);\n",
    "\n",
    "        // execute instruction\n",
    "        if (inst == NOP)\n",
    "            continue; \n",
    "\n",
    "        if (inst == PUSHV)\n",
    "        {\n",
    "            // fetch variable from data\n",
    "            x = *(gdata + data_row_idx + (int)(operand));\n",
    "\n",
    "            // push to stack\n",
    "            stack[sp] = x;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if (inst == PUSHC)\n",
    "        {\n",
    "            // push constant to stack\n",
    "            stack[sp] = operand;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        // math/logic instructions go below\n",
    "        if ((inst == ADD) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = x + y;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == MUL) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = x * y;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == DIV) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = x / y;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == NEG) && (sp >= 1))\n",
    "        {\n",
    "            /*\n",
    "            // pop one value from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = -x;\n",
    "            sp++;*/\n",
    "\n",
    "            stack[sp] = -stack[sp];\n",
    "        }\n",
    "\n",
    "        if ((inst == MIN) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = (x > y)?y:x;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == MAX) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = (x > y)?x:y;\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == GREATER) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = (float)(x > y);\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == LESS) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = (float)(x < y);\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == EQUAL) && (sp >= 2))\n",
    "        {\n",
    "            // pop two values from stack\n",
    "            x = stack[sp];\n",
    "            sp--;\n",
    "            y = stack[sp];\n",
    "            sp--;\n",
    "\n",
    "            // push result to stack\n",
    "            stack[sp] = (float)(x == y);\n",
    "            sp++;\n",
    "        }\n",
    "\n",
    "        if ((inst == SIN) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = sin(stack[sp]);\n",
    "        }\n",
    "\n",
    "        if ((inst == COS) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = cos(stack[sp]);\n",
    "        }\n",
    "\n",
    "        if ((inst == EXP) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = exp(stack[sp]);\n",
    "        }\n",
    "\n",
    "        if ((inst == LOG) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = log(stack[sp]);\n",
    "        }\n",
    "\n",
    "        if ((inst == SQR) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = stack[sp]*stack[sp];\n",
    "        }\n",
    "\n",
    "        if ((inst == SQRT) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = sqrt(stack[sp]);\n",
    "        }\n",
    "\n",
    "        if ((inst == TANH) && (sp >= 1))\n",
    "        {\n",
    "            stack[sp] = tanh(stack[sp]);\n",
    "        }\n",
    "    }\n",
    "\n",
    "    // store result\n",
    "    float f = 1/(1+exp(-stack[sp]));\n",
    "    if (f > 0.5) \n",
    "    {\n",
    "        f = 1.0;\n",
    "    }\n",
    "    else \n",
    "    {\n",
    "        f = 0.0;\n",
    "    }\n",
    "    gresult[result_idx] = f;    \n",
    "}\n",
    "\n",
    "}\n",
    "  \"\"\", no_extern_c=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "func = mod.get_function(\"compute\")\n",
    "init_func = mod.get_function(\"initkernel\")\n",
    "initprog_func = mod.get_function(\"init_progs\")\n",
    "reproduce_func = mod.get_function(\"reproduce\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modsc = SourceModule(\"\"\"\n",
    "\n",
    "#define NUM_THREADS_X %d\n",
    "\n",
    "  __global__ void scores(int num_real_rows, float *gdatay, float *gresult, float *gscores)\n",
    "  {\n",
    "      int th_x = blockIdx.x * blockDim.x + threadIdx.x;\n",
    "      float sum = 0;\n",
    "      \n",
    "      for(int i=0; i<num_real_rows; i++)\n",
    "      {\n",
    "          float y1 = gdatay[i];\n",
    "          float y2 = *(gresult + i*NUM_THREADS_X + th_x);\n",
    "          if (y1 == y2)\n",
    "          {\n",
    "              sum ++;\n",
    "          }\n",
    "      }\n",
    "      \n",
    "      sum /= num_real_rows;\n",
    "      gscores[th_x] = sum;\n",
    "  }\n",
    "  \"\"\" % (num_threads_x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "funcsc = modsc.get_function(\"scores\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# call this only once\n",
    "seed = rnd.randint(0, 10000)\n",
    "init_func(np.int32(seed), block=(threads_per_block, 1, 1), grid=(blocks_x, 1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init first population\n",
    "initprog_func(gprogs_gpu, block=(threads_per_block, 1, 1), grid=(blocks_x, 1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# store data\n",
    "gdata[0 : bdx.shape[0], :] = bdx\n",
    "gdatay[0 : bdx.shape[0], :] = bdy.reshape(-1, 1)\n",
    "\n",
    "gdatatest[0 : bdx_test.shape[0], :] = bdx_test\n",
    "gdatatesty[0 : bdx_test.shape[0], :] = bdy_test.reshape(-1, 1)\n",
    "\n",
    "num_real_rows = np.int32(bdx.shape[0])\n",
    "num_real_rows_test = np.int32(bdx_test.shape[0])\n",
    "\n",
    "cuda.memcpy_htod(gdata_gpu, gdata)\n",
    "cuda.memcpy_htod(gdatay_gpu, gdatay)\n",
    "cuda.memcpy_htod(gdatatest_gpu, gdatatest)\n",
    "cuda.memcpy_htod(gdatatesty_gpu, gdatatesty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evolution loop\n",
    "best_ever = 0\n",
    "for i in tqdm(range(10000000)):\n",
    "    # evaluate\n",
    "    func(num_real_rows, gprogs_gpu, gdata_gpu, gresult_gpu, \n",
    "         block=(threads_per_block, threads_per_block, 1), grid=(blocks_x, blocks_y, 1))\n",
    "    # compute scores\n",
    "    funcsc(num_real_rows, gdatay_gpu, gresult_gpu, gscores_gpu, \n",
    "            block=(threads_per_block, 1, 1), grid=(blocks_x, 1, 1))\n",
    "    # get scores\n",
    "    cuda.memcpy_dtoh(gscores, gscores_gpu)\n",
    "    accs = gscores.reshape(-1)\n",
    "\n",
    "    best = np.max(accs)\n",
    "    best_idx = np.argmax(accs)\n",
    "    #print('best', best)\n",
    "    if best > best_ever:\n",
    "        best_ever = best\n",
    "        # eval on test data\n",
    "        # evaluate\n",
    "        func(num_real_rows_test, gprogs_gpu, gdatatest_gpu, gresulttest_gpu, \n",
    "             block=(threads_per_block, threads_per_block, 1), grid=(blocks_x, blocks_y, 1))\n",
    "        # compute scores\n",
    "        funcsc(num_real_rows_test, gdatatesty_gpu, gresulttest_gpu, gscorestest_gpu, \n",
    "                block=(threads_per_block, 1, 1), grid=(blocks_x, 1, 1))\n",
    "        # get scores\n",
    "        cuda.memcpy_dtoh(gscorestest, gscorestest_gpu)\n",
    "        taccs = gscorestest.reshape(-1)\n",
    "        print('new best', best, 'test:',taccs[best_idx])\n",
    "        \n",
    "    # sort scores, put the sorted_idx in the GPU, and reproduce next generation\n",
    "    sorted_idx = np.argsort(accs)[::-1].astype(np.int32)\n",
    "    cuda.memcpy_htod(sorted_idx_gpu, sorted_idx)\n",
    "    \n",
    "    reproduce_func(sorted_idx_gpu, gprogs_gpu, gprogsnew_gpu, block=(threads_per_block, 1, 1), grid=(blocks_x, 1, 1))\n",
    "    \n",
    "    # replace the old pop\n",
    "    cuda.memcpy_dtod(gprogs_gpu, gprogsnew_gpu, gprogs.nbytes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": []
  }
 ],
 "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.6.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}