{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Genetic Algorithm on GPU + Constraint Programming on CPU "},{"metadata":{"trusted":false},"cell_type":"code","source":"import multiprocessing\nimport os\nimport random\nimport time\nimport numpy as np\nimport pandas as pd\nimport torch\n\nfrom datetime import datetime\nfrom joblib import Parallel, delayed\nfrom ortools.sat.python import cp_model\nfrom scipy.signal import convolve2d\nfrom sklearn.metrics import mean_absolute_error","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Constants"},{"metadata":{"trusted":false},"cell_type":"code","source":"N = 25  # grid dimension\ndevice = 'cuda'\nTEST_CSV = '../input/conways-reverse-game-of-life-2020/test.csv'\nOUTPUT_CSV = 'submission.csv'\nGA_OUTPUT_CSV = 'submission-ga.csv'\nOR_OUTPUT_CSV = 'submission-or.csv'\nTIME_BUDGET = 1.85*3600\nOR_MAX_BUDGET = 10","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"start_time = time.time()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Genetic algorithm\n\nUseful reading: https://www.mathworks.com/help/gads/how-the-genetic-algorithm-works.html"},{"metadata":{"trusted":false},"cell_type":"code","source":"torch.backends.cudnn.benchmark = True\ntorch.backends.cudnn.deterministic = False","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef loss(input, target):\n    return torch.sum(input ^ target, dim=(-1,-2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"cv = torch.nn.Conv2d(1, 1, kernel_size=3, padding=1, padding_mode='circular', bias=False)\ncv.requires_grad=False\ncv.weight = torch.nn.Parameter(\n    torch.tensor(\n        [[[[ 1., 1., 1.],\n           [ 1., 0., 1.],\n           [ 1., 1., 1.]]]],\n        device=device,\n        dtype=torch.float16\n    ),\n    requires_grad=False,\n)\n\n\n@torch.jit.script\ndef forward(grid, delta: int):\n    N=25\n    g = grid.reshape(-1, 1, N, N)\n    for _ in torch.arange(delta):\n        g = g.to(torch.float16)\n        neighbor_sum = cv(g)\n        g = ((neighbor_sum == 3) | ((g == 1) & (neighbor_sum == 2)))\n    return g.reshape(-1, N, N)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef random_parents(data, n_parents: int, device: str):\n    N = 25\n    RANDOM_ALIVE = .25\n    return torch.rand((n_parents, N, N), device=device) > (1-RANDOM_ALIVE)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef select_best(parents, delta: int, target, n_best: int):\n    scores = loss(forward(parents, delta), target)\n    best_values, best_indices = torch.topk(scores, n_best, dim=0, largest=False, sorted=True)\n    new_parents = parents[best_indices, ...]\n    return new_parents, best_values[0], new_parents[0, ...]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def precomputes_masks():\n    N = 25\n    BLOCK_SIZE = 17\n\n    block = torch.nn.Conv2d(1, 1, kernel_size=BLOCK_SIZE, padding=BLOCK_SIZE//2,\n                            padding_mode='circular', bias=False)\n    block.requires_grad=False\n    block.weight = torch.nn.Parameter(\n        torch.ones((1, 1, BLOCK_SIZE, BLOCK_SIZE),\n            device=device,\n            dtype=torch.float16\n        ),\n        requires_grad=False,\n    )\n\n    masks = torch.zeros((N * N, 1, N, N), device=device, dtype=torch.float16)\n    \n    for x in range(N):\n        for y in range(N):\n            masks[x * N + y, 0, x, y] = 1.\n    masks = block(masks)\n    \n    return masks[:, 0, ...] > .5","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef random_combine(parents, n_offsprings: int, device: str, pre_masks):\n    N = 25\n    \n    dads = torch.randint(low=0, high=parents.shape[0], size=(n_offsprings,),\n                         device=device, dtype=torch.long)\n    dads = parents[dads, ...]\n    \n    moms = torch.randint(low=0, high=parents.shape[0], size=(n_offsprings,),\n                         device=device, dtype=torch.long)\n    moms = parents[moms, ...]\n    \n    masks = pre_masks[torch.randint(low=0, high=pre_masks.shape[0], size=(n_offsprings,),\n                                    device=device, dtype=torch.long)]\n\n    return torch.where(masks, dads, moms)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef mutate(parents, device: str):\n    MUTATION = .0016  # .005 \n    mutations = torch.rand(parents.shape, device=device) < MUTATION\n    return parents ^ mutations","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef optimize_one_puzzle(delta: int, data, device: str, pre_masks):\n    N = 25\n    N_GENERATION = 30  # Number of generations\n    P = 4_150  # population\n    N_BEST = P // 30  # best to keep as new parents\n    N_ELITES = 8  # parents unchanged for next generation\n    \n    best_score = torch.tensor([N*N], device=device)\n    best = torch.zeros((N,N), device=device).to(torch.bool)\n    parents = random_parents(data, P, device)\n\n    elites = torch.empty((1, N, N), dtype=torch.bool, device=device)\n    elites[0, ...] = data  # set target as potential dad ;)\n\n    for i in range(N_GENERATION):\n        parents = random_combine(parents, P, device, pre_masks)\n        parents = mutate(parents, device)\n        parents[:N_ELITES, ...] = elites\n        parents, best_score, best = select_best(parents, delta, data, N_BEST)\n        elites = parents[:N_ELITES, ...]\n        if best_score == 0:  # early stopping\n            break\n\n    return best_score, best","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":false},"cell_type":"code","source":"@torch.jit.script\ndef optimize_all_puzzles(deltas, df, device: str, pre_masks):\n    sub = df.clone()\n    \n    for n in torch.arange(df.shape[0]):\n        delta = deltas[n]\n        data = df[n, ...]\n        _, sub[n, ...] = optimize_one_puzzle(delta, data, device, pre_masks)\n\n    return sub","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def leaderboard_score(deltas, df, sub, device: str):\n    result = torch.empty(sub.shape[0], device=device, dtype=torch.long)\n    for delta in range(1, 6):\n        start = sub[deltas == delta]\n        end   = df[deltas == delta]\n        result[deltas == delta] = loss(forward(start, delta), end)\n    print('Leaderboard score (TORCH):', torch.sum(result).item() / (result.shape[0]*N*N))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def genetic():\n    df = pd.read_csv(TEST_CSV, index_col='id')\n#     df = df.iloc[:1000, :]\n    submission = df.copy()\n    submission.drop(['delta'], inplace=True, axis=1)\n    indexes = df.index\n    deltas = torch.from_numpy(df.delta.values).to(device)\n    df = torch.BoolTensor(df.values[:, 1:].reshape((-1, N, N))).to(device)\n    start_time = time.time()\n    pre_masks = precomputes_masks()\n    sub = optimize_all_puzzles(deltas, df, device, pre_masks)\n    print(f'Processed {sub.shape[0]:,} puzzles in {time.time() - start_time:.2f} seconds 🔥🔥🔥')\n    submission.iloc[:sub.shape[0], :] = sub.reshape((-1, N*N)).cpu().numpy().astype(int)\n    submission.to_csv(GA_OUTPUT_CSV)\n    leaderboard_score(deltas, df, sub, device)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Ortools"},{"metadata":{"trusted":false},"cell_type":"code","source":"def life_step(X):\n    nbrs_count = convolve2d(X, np.ones((3, 3)), mode='same', boundary='wrap') - X\n    return (nbrs_count == 3) | (X & (nbrs_count == 2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def individual_scores(df, pred):\n    scores = []\n    for i in pred.index:\n        delta = df.loc[i][0]\n        start = np.asarray(pred.loc[i]).reshape(N, N)\n        end   = np.asarray(df.loc[i][1:]).reshape(N, N)\n        x = start\n        for _ in range(delta):\n            x = life_step(x)\n        scores.append(mean_absolute_error(x, end))\n    return scores","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def move_forward(start, delta):\n    x = start\n    for _ in range(delta):\n        x = life_step(x)\n    return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def pd_leaderboard_score(df, pred):\n    scores = individual_scores(df, pred)\n    print('Leaderboard score (PD):', sum(scores) / len(scores))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"class Grid:\n    def __init__(self, model, dim):\n        self._cells = {}\n        self._dim = dim\n        for x in range(dim):\n            for y in range(dim):\n                self._cells[x,y] = model.NewBoolVar(f'cell {x}, {y}')\n\n    def __getitem__(self, pos):\n        x, y = pos\n        n = self._dim\n        return self._cells[(x+n) % n, (y+n) % n]\n    \n    def solution(self, solver):\n        n = self._dim\n        s = np.ones((n, n), dtype=np.uint8)\n        for x in range(n):\n            for y in range(n):\n                s[x, y] = solver.Value(self._cells[x, y])\n        return s        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def solve_one_step(end, use_zero_point):\n    budget = TIME_BUDGET - (time.time() - start_time)\n    if budget <= 0:\n        return None\n\n    # Creates the model.\n    model = cp_model.CpModel()\n\n    grid = Grid(model, N)\n\n    # Creates the constraints with final state\n    for x in range(N):\n        for y in range(N):\n            crown = (grid[x-1,y-1] + grid[x-1,y] + grid[x-1,y+1] + \n                     grid[x  ,y-1] +               grid[x  ,y+1] + \n                     grid[x+1,y-1] + grid[x+1,y] + grid[x+1,y+1])\n            if end[x,y] == 1:\n                # (crown == 3) or (crown == 2 and grid[x,y] == 1)\n                b = model.NewBoolVar(\"\")\n                model.Add(crown == 3).OnlyEnforceIf(b)\n                model.Add(crown == 2).OnlyEnforceIf(b.Not())\n                model.Add(grid[x,y] == 1).OnlyEnforceIf(b.Not())\n            # zero point idea stolen from\n            # https://www.kaggle.com/jamesmcguigan/game-of-life-z3-constraint-satisfaction\n            elif(use_zero_point and\n                 (end[(x-1)%N,(y-1)%N] + end[(x-1)%N,y  ] + end[(x-1)%N,(y+1)%N] + \n                  end[x      ,(y-1)%N] +                    end[x      ,(y+1)%N] + \n                  end[(x+1)%N,(y-1)%N] + end[(x+1)%N,y  ] + end[(x+1)%N,(y+1)%N]) == 0):\n#                  and delta != 1):\n                model.Add(grid[x,y] == 0)\n                model.Add(crown != 3)\n            else:\n                # (crown != 3) and (crown != 2 or grid[x,y] == 0)\n                b = model.NewBoolVar(\"\")\n                model.Add(crown != 3)\n                model.Add(grid[x,y] == 0).OnlyEnforceIf(b)\n                model.Add(crown != 2).OnlyEnforceIf(b.Not())\n            \n    # Creates a solver and solves the model.\n    solver = cp_model.CpSolver()\n\n    # Sets a time limit.\n    solver.parameters.max_time_in_seconds = min(budget, OR_MAX_BUDGET)\n\n    status = solver.Solve(model)\n\n    if status in (cp_model.FEASIBLE, cp_model.OPTIMAL):\n        return grid.solution(solver)\n    return None","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def solve(idx, delta, end, max_err):\n    start = end\n    \n    for d in range(delta):\n        start = solve_one_step(start, True)\n        if start is None:\n            break\n\n    return start","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def load_solve(df, i, err):\n    budget = TIME_BUDGET - (time.time() - start_time)\n    if budget <= 0:\n        return None\n    \n    start = time.time()\n    delta = df.loc[i][0]\n    end = np.asarray(df.loc[i][1:]).reshape(N, N)\n    ret = solve(i, delta, end, err)\n    if ret is None:\n        print(f'Not solved {i}: delta={delta} in {time.time() - start}s')\n#     else:\n#         print(f'Solved {i}: delta={delta} in {time.time() - start}s           <----------==')\n    return ret","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def ortools():\n    os.nice(10)\n    df = pd.read_csv(TEST_CSV, index_col='id')\n    submission = df.copy()\n    submission.drop(['delta'], inplace=True, axis=1)\n    scores = individual_scores(df, submission)\n    submission['score'] = scores\n    print('Leaderboard score:', sum(scores) / len(scores))\n\n    submission.loc[df.index, 'alive2'] = df.loc[df.index][1:].mean(axis=1) * df.loc[df.index].delta\n\n    to_solve_df = submission.loc[submission.score > 0]\n    to_solve_df = to_solve_df.sort_values('alive2')\n    to_solve = [(idx, int(round(score * N * N))) for idx, score in zip(to_solve_df.index.values, to_solve_df.score.values)]\n\n    submission.drop(['score', 'alive2'], axis=1, inplace=True)\n\n    start_time = time.time()\n    print(f'Starts solving {len(to_solve)} puzzles')\n\n    solutions = Parallel(n_jobs=2)(\n        delayed(load_solve)(df, i, e) for i, e in to_solve\n    )\n\n    print(f'Improved {sum([x is not None for x in solutions])} puzzles in {time.time() - start_time:.2f} seconds 🔥🔥🔥')\n\n    for (i, _), s in zip(to_solve, solutions):\n        if s is not None:\n            submission.loc[i] = s.reshape(N*N)\n\n    submission.to_csv(OR_OUTPUT_CSV)\n    pd_leaderboard_score(df, submission)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Run GA and OR in parallel"},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\n\nd = multiprocessing.Process(name='ortools', target=ortools)\nd.daemon = True\nd.start()\n\ngenetic()\n\nd.join()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Merge solutions"},{"metadata":{"trusted":false},"cell_type":"code","source":"df = pd.read_csv(TEST_CSV, index_col='id')\n\ndf1 = pd.read_csv(OR_OUTPUT_CSV, index_col='id')\ndf1['score'] = individual_scores(df, df1)\n\ndf2 = pd.read_csv(GA_OUTPUT_CSV, index_col='id')\ndf2['score'] = individual_scores(df, df2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"for i in df2.index:\n    if df1.loc[i].score > df2.loc[i].score:\n        df1.loc[i] = df2.loc[i]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"df1.drop('score', axis=1, inplace=True)\ndf1 = df1.astype('int')\npd_leaderboard_score(df, df1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"df1.to_csv(OUTPUT_CSV)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"df1.head()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}