{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom numba import njit, prange\nfrom scipy.signal import convolve2d\nfrom tqdm.auto import tqdm\nimport matplotlib.pyplot as plt\nimport pycosat\nimport random\nimport time\nimport tensorflow as tf\nnew_model = tf.keras.models.load_model(\"../input/testtest2/model_iterative_1step_40epochs_new.h5\",compile=False)\nnew_model.summary()\n\nimport multiprocessing\nfrom multiprocessing import Process\nfrom time import sleep\nimport copy","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"window = np.ones((3, 3))\n\ndef plot_3d(solution_3d: np.ndarray, title, size=3, max_cols=10, has_target=False):\n    N = len(solution_3d)\n    if N <= 0:\n        return\n    cols = min(N, max_cols)\n    rows = (N - 1) // max_cols + 1\n    plt.figure(figsize=(cols*size, rows*size))\n    plt.suptitle(title)\n    for t in range(N):\n        board = solution_3d[t]\n        plt.subplot(rows, cols, t + 1)\n        plt.imshow(board, cmap='binary')\n        plt.title('target' if t == 0 and has_target else f'state {t}')\n    plt.show()\n\ndef csv_to_numpy_list(df) -> np.ndarray:\n    return df[[ f'stop_{n}' for n in range(25**2) ]].values.reshape(-1,25,25)\n\ndef life_step(X: np.ndarray):\n    C = convolve2d(X, window, mode='same', boundary='wrap')\n    return (C == 3) | (X & (C == 4))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For SAT solvers we only need Test data: there is nothing to train and we always know if solution is correct."},{"metadata":{},"cell_type":"markdown","source":"* Ok, let's try out some of the hepler functions above and look at some sample test board."},{"metadata":{},"cell_type":"markdown","source":"Ok, finally we are getting to the solver itself.\n\nLooks like pycosat requires input in [Conjunctive normal form](https://en.wikipedia.org/wiki/Conjunctive_normal_form) and does not support weights, optional clauses or native cardinality constraints. It can present some problems later when we would want to limit the search space or prefer some solutions to another. But for now the simpler the better.\n\nWhile trying to build CNF expression for Game of Life, I found this article, that helped me to reduce number of clauses and variables a lot:\nhttps://sat-smt.codes/SAT_SMT_by_example.pdf (Pages 501-504)"},{"metadata":{"trusted":true},"cell_type":"code","source":"SIZE = 25\n\nempty_board = np.zeros((SIZE,SIZE), dtype=bool)\n\n# Computes number of variables.\n# I negated variable states after some testing, because solver\n# seems to prefer True as default state for solution variables,\n# while I want cells to be empty whenever possible\ndef v(c):\n    return -(SIZE * c[0] + c[1] + 1)\n\ndef dead_clauses(res, c, x):\n    # if cell is dead, there was not exactly 3 alive neighbours (56 clauses)\n    for i1 in range(0, 6):\n        for i2 in range(i1+1, 7):\n            for i3 in range(i2+1, 8):\n                a = [v(x[i]) for i in range(8)]\n                a[i1], a[i2], a[i3] = -a[i1], -a[i2], -a[i3]\n                #print(f\"1/56={a}\")\n                res.append(a)\n    # if cell is dead and was alive, was not 2-3 alive neighbours (28 clauses)\n    for i1 in range(0, 7):\n        for i2 in range(i1+1, 8):\n            a = [v(x[i]) if i < 8 else -v(c) for i in range(9)]\n            a[i1], a[i2] = -a[i1], -a[i2]\n            res.append(a)\n\ndef live_clauses(res, c, x):\n    # if cell is alive, there was less then 4 alive neighbours (70 clauses)\n    \n    #cell live      prev   dead      prev neighbours 3\n    #                       live         2 or 3\n    for i1 in range(0, 5):\n        for i2 in range(i1+1, 6):\n            for i3 in range(i2+1, 7):\n                for i4 in range(i3+1, 8):\n                    #from each 4 at least 1 was dead\n                    res.append([-v(x[i1]), -v(x[i2]), -v(x[i3]), -v(x[i4])])\n    # if cell is alive and was dead, there was more than 2 alive (less than 6 dead) neighbours (28 clauses)\n    \n    for i1 in range(0, 7):\n        for i2 in range(i1+1, 8):\n            a = [v(x[i]) if i < 8 else v(c) for i in range(9) if i != i1 and i != i2]    # v(c)\n            res.append(a)\n    \n    # if cell is alive, there was more than 1 alive (less than 7 dead) neighbours (8 clauses)\n    for i1 in range(0, 8):\n        a = [v(x[i]) for i in range(8) if i != i1]\n        res.append(a)\n\ndef board_clauses(board, use_opt = True):\n    res, opt1, opt2 = [], [], []\n    for i in range(SIZE):\n        for j in range(SIZE):\n            x = [((i + k % 3 - 1) % SIZE, (j + k // 3 - 1) % SIZE) for k in range(9) if k != 4]\n            #print(x)\n            if board[i,j]:\n                live_clauses(res, (i, j), x) # 106 clauses  x means coord of neighbours\n            else:\n                dead_clauses(res, (i, j), x) # 84 clauses\n                if use_opt:\n                    y = [((i + k % 5 - 2) % SIZE, (j + k // 5 - 2) % SIZE) for k in range(25)]\n                    if sum(board[ii,jj] for ii,jj in y) < 1: # No alive neighbours\n                        res.append([-v((i, j))]) # Very dead space MUST stay dead! (1 clause)\n                    elif sum(board[ii,jj] for ii,jj in x) < 1: # No alive neighbours\n                        opt1.append([-v((i, j))]) # Dead space should stay dead! (1 clause)\n                    elif sum(board[ii,jj] for ii,jj in x) < 2: # Too few alive neighbours\n                        opt2.append([-v((i, j))]) # Dead space should stay dead! (1 clause)\n           \n    return res, opt1, opt2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Depth(delta)-1 tasks can be solved relatively easy:"},{"metadata":{},"cell_type":"markdown","source":"It takes just a bit more than half a second to solve delta-1 task in average.\n\nNext, I've built step-by-step solver for individual tasks, with ability to gradualy exclude optional clauses and fallback to previous step when there is no solution.\n\nLet's try it out for first tasks."},{"metadata":{"trusted":true},"cell_type":"code","source":"def satsolver(n, ret):\n    T = min(deltas[n], 5)\n    board = np.tile(empty_board, (T+1, 1, 1))\n    board[0] = boards[n]\n    solvers = [None for _ in range(T)]\n    opt = [None for _ in range(T)]\n    os = [0 for _ in range(T)]\n    oe = [0 for _ in range(T)]\n    t = 0\n    while 0 <= t and t < T:\n        if solvers[t] is None:\n            clauses, opt1, opt2 = board_clauses(board[t])\n            solution = pycosat.solve(clauses)\n            if not isinstance(solution, str):\n                if t == T - 1:\n                    #print(t, '!!', end=\" \")\n                    t += 1\n                    board[t] = np.array(solution[:SIZE**2]).reshape(SIZE,SIZE) < 0\n                    continue\n                else:\n                    #print(t, '??', end=\" \")\n                    random.shuffle(opt1)\n                    random.shuffle(opt2)\n                    opt[t] = opt1 + opt2\n                    os[t] = len(opt[t])+1     \n                    oe[t] = len(opt[t])+1\n                    solvers[t] = pycosat.itersolve(clauses + opt[t])\n                    \n        try:\n            #print(\"trytrytry\")\n            solution = next(solvers[t])\n            if oe[t] - os[t] > 1:\n                os[t] = (os[t]+oe[t])//2\n                solvers[t] = pycosat.itersolve(clauses + opt[t][:(os[t]+oe[t])//2])\n            else:\n                t += 1\n                board[t] = np.array(solution[:SIZE**2]).reshape(SIZE,SIZE) < 0\n        except Exception as err:\n            #print(\"except\")\n            if solvers[t] is not None and (os[t]+oe[t])//2 > 0:\n                if os[t] == oe[t]:\n                    os[t] = 0\n                    oe[t] = len(opt[t])+1\n                elif oe[t] - os[t] > 1:\n                    oe[t] = (os[t]+oe[t])//2\n                else:\n                    os[t] -= 1\n                    oe[t] -= 1\n                solvers[t] = pycosat.itersolve(clauses + opt[t][:(os[t]+oe[t])//2])\n            else:\n                solvers[t] = None\n                opt[t] = None\n                t -= 1   \n\n    if t == T:\n        pred_start = 1 * board[T].ravel()\n        for k in range(625):\n            ret[k] = pred_start[k]\n#         ret = copy.deepcopy(pred_start)\n#         ret = pred_start\n    \ndef unet(n):\n    for i in range(deltas[n]):\n        if i == 0:\n            preds = new_model.predict(boards[n].flatten().reshape(1,25,25,1))\n        else:\n            preds = tf.where(preds>0.5,1,0)\n            preds = new_model.predict(preds)\n\n    preds = tf.cast(tf.where(preds>0.5,1,0), tf.uint32)\n    preds=tf.reshape(preds,[-1])\n    pred_start = np.array(preds)\n    #print(pred_start.shape)\n    test_result_df.loc[test_df.index[n]] = pred_start\n    err = get_scores_on_1_game(pred_start, n)\n    print(f\"N={n}, delta={deltas[n]}, Method: CNN, err={err}\", end='')\n    errs.append(err)\n\ndef apply_all_dead(n):\n    \n    pred_start =np.zeros((625,1)).astype(\"int64\")\n    pred_start=pred_start.reshape(-1)\n    #print(pred_start.shape)\n    test_result_df.loc[test_df.index[n]] = pred_start\n    err = get_scores_on_1_game(pred_start, n)\n    print(f\"N={n}, delta={deltas[n]}, Method: Apply All Dead, err={err}\", end='')\n    errs.append(err)\n    \n    \ndef run_with_limited_time(func, n, kwargs, time):\n    \"\"\"Runs a function with time limit\n\n    :param func: The function to run\n    :param args: The functions args, given as tuple\n    :param kwargs: The functions keywords, given as dict\n    :param time: The time limit in seconds\n    :return: True if the function ended successfully. False if it was terminated.\n    \"\"\"\n    \n    ret = multiprocessing.Array(\"i\", 625)\n    p = Process(target=func, args=[n, ret], kwargs=kwargs)\n#     p = Process(target=func, args=(n,), kwargs=kwargs)\n    p.start()\n    p.join(time)\n    if p.is_alive():\n        p.terminate()\n        return False\n#     print(\"\\n haha\\n\",np.array(ret))\n    pred_start = np.array(ret)\n    err = get_scores_on_1_game(pred_start, n)\n    print(f\"N={n}, delta={deltas[n]}, Method: SAT, err={err}\", end='')\n    test_result_df.loc[test_df.index[n]] = pred_start\n    errs.append(get_scores_on_1_game(pred_start, n))\n    return True\n\ndef get_scores_on_1_game(board, n):\n    board=np.array(board)\n    for j in range(deltas[n]):\n        board = life_step(board.reshape(SIZE,SIZE))\n    d = np.sum(board ^ boards[n])\n    return d/625","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission_df = pd.read_csv('../input/conways-reverse-game-of-life-2020/sample_submission.csv', index_col='id')\ntest_df = pd.read_csv('../input/conways-reverse-game-of-life-2020/test.csv', index_col='id')\ndeltas = test_df['delta'].values\nboards = csv_to_numpy_list(test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"delta5=test_df.iloc[[0,100,200]]\nprint(delta5)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_size = 1000\ntimes=[]\nerrs=[]\ntest_result_df = sample_submission_df.copy().iloc[9000:10000]\nfor n in tqdm(range(9000,10000), total=test_size):\n\n    start=time.time()    \n    if deltas[n]==5:\n        unet(n)\n        processtime = time.time()-start\n\n        times.append(processtime)\n        print(f\" time={processtime}\")\n        continue\n        \n    success=run_with_limited_time(satsolver, n, {}, 30)\n    if success==False:\n        #apply_all_dead(n)\n        unet(n)\n\n    end=time.time()\n    processtime = end-start\n    times.append(processtime)\n    print(f\" time={processtime}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(times))\nprint(len(errs))\nprint(test_result_df.shape)\n#print(test_result_df.iloc[997:1001])\n\n#test_result_df=test_result_df.iloc[1000:2000]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_result_df[\"time\"] = np.array(times)\ntest_result_df[\"err\"] = np.array(errs)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_result_df.head()\ntest_result_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_result_df.to_csv(\"result_combined_9000-10000.csv\", index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.mean(np.array(times)))\nprint(np.mean(np.array(errs)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_result_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_time=test_result_df['time'].mean()\nmean_err=test_result_df['err'].mean()\nprint(mean_time,mean_err)\nm=test_result_df['time'].max()\nprint(m)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import f1_score,precision_recall_fscore_support,accuracy_score\n\nN=10000\nfinal=[]\nf1score=[]\nfor i in tqdm(range(N),total=N):\n    \n#     y_pred=np.array(test_result_df.iloc[i][0:625].astype(\"int64\"))\n#     for j in range(deltas[i]):\n#         y_pred = life_step(y_pred.reshape(SIZE,SIZE))\n    \n#     y_pred=y_pred.reshape(-1)\n\n    y_pred=np.zeros((625,1))\n    y_true=np.array(boards[i]).reshape(-1)\n#     print(y_pred.shape)\n#     print(y_true.shape)\n    re=accuracy_score(y_true,y_pred)\n    #f1=f1_score(y_true,y_pred,average='binary')\n    #print(re)\n    final.append(re)\n    \n#print(final)\nprint(np.mean(final))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def confirm_on_test_result(test_index):\n    result = np.array(test_result_df.iloc[test_index])\n    pred = result[:-2].astype('int64')\n    true_err = get_scores_on_1_game(pred, test_index)\n    print(true_err)\n    return true_err==result[-1]\n\nprint(confirm_on_test_result(50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_result_df.iloc[0:5]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aaa_df=pd.read_csv(\"../input/testtest/result_cnn_0-10000.csv\")\naaa_df.head()\naaa_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.iloc[idx, 1:626]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# visulization\n\nidx = 15\n\nfig, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize=(10,10))\n#fig.suptitle(f\"idx={idx}\")\nax1.imshow((1-test_df.iloc[idx, 1:626].values).reshape(25, 25), cmap=\"gray\")\nax1.set_title(\"End Board\")\n\nprint(aaa_df.iloc[idx, 1:626])\n\nax2.imshow((1-aaa_df.iloc[idx, 1:626].values).reshape(25, 25), cmap=\"gray\")\nax2.set_title(\"Predicted Start Board\")\n\n\ny_pred=np.array(aaa_df.iloc[idx][1:626].astype(\"int64\"))\nfor j in range(deltas[idx]):\n    y_pred = life_step(y_pred.reshape(SIZE,SIZE))\n    \ny_pred=y_pred.reshape(-1)\ny_pred=np.where(y_pred>0.5,0,1)\nprint(y_pred.reshape(25,25))\n\nax3.imshow(y_pred.reshape(25, 25), cmap=\"gray\")\nax3.set_title(f\"Predicted Start Board After {deltas[idx]} Steps \")\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"time_df=aaa_df[aaa_df[\"err\"]==0.0]\nprint(time_df)\n# m=time_df[\"time\"]\n# print(m)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"unet_df=aaa_df[aaa_df[\"err\"]!=0.0]\nprint(unet_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import f1_score,precision_recall_fscore_support,accuracy_score\n\nN=372\nfinal=[]\nf1score=[]\nfor i in range(N):\n    \n#     index=np.array(aaa_df.iloc[i].index.astype(\"int64\"))\n    index=int(unet_df.iloc[i]['id'])-50000\n\n#     for j in range(deltas[5000+i]):\n#         y_pred = life_step(y_pred.reshape(SIZE,SIZE))\n    #print(index)\n    y_pred=np.zeros((625,1))\n    y_true=np.array(boards[index]).reshape(-1)\n    #print(y_pred)\n    #print(y_true)\n    re=accuracy_score(y_true,y_pred)\n    #f1=f1_score(y_true,y_pred,average='binary')\n    #print(re)\n    final.append(re)\n    \nprint(final)\nprint(np.mean(final))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nscores=[]\n\n\nfor i in range (0,100):\n    board=aaa_df[i][0:625].astype(\"int64\")\n    board=np.array(board)\n    for j in range(deltas[i]):\n        board = life_step(board.reshape(SIZE,SIZE))\n    d = np.sum(board ^ boards[i])\n    scores.append(d/625)\n\nprint(scores)\nprint(np.mean(scores))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"limit=10  error=0.115696   0.114 with 1000 data\n\nlimit=1000  error=0.1141\n\nlimit=10000  error=0.1083\n\nlimit=20000  error=0.1033\n\nwithout limit error=0.091 \n\nT=1   average 0.5\n\nT= 1 2  average 5.\ntime=60 error=0.15"},{"metadata":{},"cell_type":"markdown","source":"Finally, let's save resulting submission"},{"metadata":{},"cell_type":"markdown","source":"![](http://)So it does solve some tasks already!\n\nObviously performance is not good enough and a lot of optimization needs to be done for deep backwards search (3+ steps).\n\nBut depth-1 tasks can be solved relatively easy:"},{"metadata":{"trusted":true},"cell_type":"code","source":"start,end=300,400\n\nresult=sample_submission_df\nresult[\"time\"]=np.zeros([50000,1])\nresult[\"time\"][start:end]=np.array(times)\n\nresult[start:end].to_csv(\"result300-400.csv\", index=True) #!!!!!!!!!!!!!\n\nsubmission_df = pd.read_csv('../input/testtest/result0-100.csv', index_col='id')\nsubmission_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nN = 10\nscore = 0\nn=1\nstart=time.time()\nfor n in tqdm(range(N), total=N):\n    print(\"xXXXXXXXXXxxx\")\n    clauses, opt1, opt2 = board_clauses(boards[n], use_opt = False)\n    print(\"eeeeeeeeeeeeeeee\")\n    solution = pycosat.solve(clauses)\n    print(solution)\n    if isinstance(solution, str):\n        print(f'{n} not solved!')\n        continue\n    board = np.array(solution[:SIZE**2]) < 0\n    sample_submission_df.loc[test_df.index[n]] = 1 * board\n    \n    board = life_step(board.reshape(SIZE,SIZE))\n    #board = life_step(board.reshape(SIZE,SIZE))\n\n    d = np.sum(board ^ boards[n])\n    score += d / 625\n#print(len(clauses))\n#print(f\"clauses={clauses}\")\n#print(solution)\n#print(board)\nend=time.time()\nptime=end-start\nprint(ptime,ptime/N)\nprint(score/N)","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}