{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 3rd Place solution (part): Simulated Annealing\n\nFirst of all, I appreciate to staff planned and organized this competition and also thanks to competitors fought together on this problem. Thanks to your effort, this competition have been a lot of fun.\n\nI mainly used simulated annealing for this problem. This algorithm is well-known and easy. The implementation of my (early) solution is below."},{"metadata":{},"cell_type":"markdown","source":"### prepare for execution\n\n- import\n- setting constants\n- forwarding function\n- score function"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom scipy.signal import convolve2d\nfrom time import perf_counter\nfrom tqdm.auto import tqdm\nimport os\n\nSIZE = 25\nNUM_PROBLEM = 10 # for debug. 50000 for submission\nSOLVE_TIME = 10 # must be short enough to submit\n\nconway_kernel = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], dtype=np.uint8)\nconway_transition = np.array([[0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0, 0, 0, 0]], dtype=np.uint8)\n\ndef game_forward(x: np.ndarray):\n    adjacent_sum = convolve2d(x, conway_kernel, mode=\"same\", boundary=\"wrap\")\n    return np.where(x == 0, conway_transition[0][adjacent_sum], conway_transition[1][adjacent_sum])\n\ndef calc_score(gt, pred):\n    assert gt[\"id\"] == pred[\"id\"], (gt, pred)\n    table = np.copy(pred[\"start\"])\n    for _ in range(gt[\"delta\"]):\n        table = game_forward(table)\n    return np.sum(table == gt[\"stop\"])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"execution\n\n- load dataset\n- solve using simulated annealing\n- write submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_dataset(csv_file_name):\n    with open(csv_file_name, \"r\") as f:\n        df = pd.read_csv(f)\n    data = np.zeros(NUM_PROBLEM, dtype=[\n        (\"id\", np.uint64),\n        (\"delta\", np.int64),\n        (\"stop\", np.uint8, (SIZE, SIZE))])\n    for i in tqdm(range(NUM_PROBLEM), desc=os.path.basename(csv_file_name)):\n        data[i][\"id\"] = df.iloc[i, 0]\n        data[i][\"delta\"] = df.iloc[i, 1]\n        data[i][\"stop\"] = np.array(df.iloc[i, 2:], dtype=np.uint8).reshape(SIZE, SIZE)\n    return data\n\ndataset = load_dataset(\"../input/conways-reverse-game-of-life-2020/test.csv\")\n\ndef solve(dataset):\n    submission = np.zeros(NUM_PROBLEM, dtype=[\n        (\"id\", np.uint64),\n        (\"score\", np.uint64),\n        (\"start\", np.uint8, (SIZE, SIZE))])\n\n    begin_temperature = 20\n    end_temperature = 0.01\n    temperature_coeff = np.log(end_temperature/begin_temperature)/SOLVE_TIME\n    progress_bar = tqdm(dataset)\n    score_sum = 0\n    for data_idx, data in enumerate(progress_bar):\n        calc_start = perf_counter()\n        submission[data_idx][\"id\"] = data[\"id\"]\n        submission[data_idx][\"start\"] = data[\"stop\"]\n        submission[data_idx][\"score\"] = calc_score(data, submission[data_idx])\n\n        duration = 0\n        while duration < SOLVE_TIME and submission[data_idx][\"score\"] < SIZE * SIZE:\n            temperature = begin_temperature * np.exp(temperature_coeff * duration)\n            i, j = np.random.randint(SIZE, size=2)\n            submission[data_idx][\"start\"][i][j] = 1 - submission[data_idx][\"start\"][i][j]\n            tmp_score = calc_score(data, submission[data_idx])\n            if tmp_score > submission[data_idx][\"score\"] or np.random.random() < np.exp((tmp_score - submission[data_idx][\"score\"]) / temperature):\n                submission[data_idx][\"score\"] = tmp_score\n            else:\n                submission[data_idx][\"start\"][i][j] = 1 - submission[data_idx][\"start\"][i][j]\n            duration = perf_counter() - calc_start\n        score_sum += submission[data_idx][\"score\"]\n        progress_bar.set_description(\"score %f\" % (1 - score_sum/(SIZE*SIZE*(data_idx + 1))))\n    return submission\n\nsubmission = solve(dataset)\n\nwritefile = \"submission.csv\"\nwith open(writefile, \"w\") as f:\n    print(f\"id,{','.join(['start_{0}'.format(i) for i in range(SIZE*SIZE)])}\", file=f)\n    for v in tqdm(submission, desc=os.path.basename(writefile)):\n        print(f\"{v['id']},{','.join(np.reshape(v['start'], -1).astype(str))}\", file=f)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## optimization\n\nThis solver is too slow (about 2.5k ~ 7k iterations/sec) to achive good score. Then, I tried to execute as fast as possible.\n\n- **Use C++**. Programs in python can run faster, but C++ is easier to optimize.\n- **Partial update and compare**. When one cell of a table is swapped, the range that changes is not more than delta. So memorizing all result of updated table, this computation can be faster.\n- **Multithreading**. I use 2 PCs (Ryzen 5 3600 (6C12T) and Core i7-8500U (4C8T)). Parallelize task-wise is easy.\n- **Bitwise update and compare**. Theoretically, each cell can be compute using only 1 bit, and I implemented in practice. The core function is `pointwise_update()`."},{"metadata":{"trusted":true},"cell_type":"code","source":"def sum_of_two_is_0(a: int, b: int):\n    return 1 - (a | b) # in C++, ~(a|b)\n\ndef sum_of_two_is_1(a: int, b: int):\n    return a ^ b\n\ndef sum_of_two_is_2(a: int, b: int):\n    return a & b\n\ndef sum_of_three_geq_2(a: int, b: int, c: int):\n    return (a & b) | ((a ^ b) & c)\n\ndef sum_of_three_mod_2(a: int, b: int, c: int):\n    return a ^ b ^ c\n\ndef sum_of_four_is_1(a: int, b: int, c: int, d: int):\n    return (sum_of_two_is_0(a, b) & sum_of_two_is_1(c, d)) | (sum_of_two_is_1(a, b) & sum_of_two_is_0(c, d))\n\nfor a in [0, 1]:\n    for b in [0, 1]:\n        assert sum_of_two_is_0(a, b) == (1 if a+b==0 else 0)\n        assert sum_of_two_is_1(a, b) == (1 if a+b==1 else 0)\n        assert sum_of_two_is_2(a, b) == (1 if a+b==2 else 0)\n        for c in [0, 1]:\n            assert sum_of_three_geq_2(a, b, c) == (1 if a+b+c>=2 else 0)\n            assert sum_of_three_mod_2(a, b, c) == (a+b+c)%2\n            for d in [0, 1]:\n                assert sum_of_four_is_1(a, b, c, d) == (1 if a+b+c+d==1 else 0)\n\ndef pointwise_update(x: np.ndarray):\n    assert x.shape == (3, 3)\n    ul, u, ur = x[0]\n    l,  a, r  = x[1]\n    dl, d, dr = x[2]\n    sum_of_remaining_term_geq_2 = sum_of_three_geq_2(\n        sum_of_two_is_1(l, r),\n        sum_of_three_mod_2(ul, u, ur),\n        sum_of_three_mod_2(dl, d, dr))\n    sum_of_neighbor_is_2_or_3 = sum_of_four_is_1(\n        sum_of_two_is_2(l, r),\n        sum_of_three_geq_2(ul, u, ur),\n        sum_of_three_geq_2(dl, d, dr),\n        sum_of_remaining_term_geq_2)\n    sum_of_neighbor_mod_2 = sum_of_three_mod_2(\n        sum_of_two_is_1(l, r),\n        sum_of_three_mod_2(ul, u, ur),\n        sum_of_three_mod_2(dl, d, dr))\n    return sum_of_neighbor_is_2_or_3 & (a | sum_of_neighbor_mod_2)\n\nfrom itertools import product\nfor v in product([0, 1], repeat=9):\n    x = np.array(v).reshape(3, 3)\n    assert pointwise_update(x) == ((sum(v) - x[1][1] == 2 and x[1][1] == 1) or sum(v) - x[1][1] == 3), x\n\ndef get_neighbor(x: np.ndarray, i: int, j: int):\n    return np.array([[x[(i+di)%SIZE][(j+dj)%SIZE] for dj in [-1, 0, 1]] for di in [-1, 0, 1]])\n\ndef game_forward_bitwise(x: np.ndarray):\n    return np.array([[pointwise_update(get_neighbor(x, i, j)) for j in range(SIZE)] for i in range(SIZE)])\n    \ntable = np.random.randint(2, size=(SIZE, SIZE))\nassert (game_forward(table) == game_forward_bitwise(table)).all()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- parallelize using 256bit **SIMD** (avx2) instead of 32bit representation.\n\nNow that we have been able to update to bitwise, we can update 8 tables at the same time by just transforming it appropriately. The implementation is below.\n\n```cpp\nconstexpr size_t MAX_DELTA = 5;\nconstexpr size_t SIZE = 25;\nconstexpr size_t NUM_PROC = 32;\n\nstruct aligned_table{\n    alignas(32) uint32_t data[MAX_DELTA+1][SIZE][NUM_PROC];\n};\n\nconstexpr uint32_t VALID_BITS = ((1 << SIZE) - 1) << 1; // 00000011 11111111 11111111 11111110\nusing bitarray_x8_t = __m256i;\n\n// after swapped tables.data[0][i0][proc..proc+7] in other functions, this function can update tables.data[1..delta][i0-delta..i0+delta][proc..proc+7].\nvoid update_partial_x8(aligned_table &tables, const int i0, const int proc, const int delta){\n    const bitarray_x8_t valid_bits_x8 = _mm256_set1_epi32(VALID_BITS);\n    const size_t buffer_offset = 6;\n    const size_t buffer_size = 2 * buffer_offset + 1;\n    bitarray_x8_t last_table[buffer_size];\n    bitarray_x8_t l_and_r[buffer_size];\n    bitarray_x8_t l_xor_r[buffer_size];\n    bitarray_x8_t xor_a_l_r[buffer_size];\n    bitarray_x8_t sum_a_l_r_geq_2[buffer_size];\n\n    for(int d=1; d<=delta; ++d){\n        for(int di=-d-1; di<=d+1; ++di){\n            const auto i_data = (i0 + SIZE + di) % SIZE;\n            const auto i = buffer_offset + di;\n            last_table[i] = _mm256_load_si256((bitarray_x8_t *)(tables.data[d-1][i_data] + proc));\n            // left rotate\n            const auto ll = _mm256_slli_epi32(last_table[i], 1), lr = _mm256_srli_epi32(last_table[i], SIZE-1);\n            const auto l = _mm256_and_si256(_mm256_or_si256(ll, lr), valid_bits_x8);\n            // right rotate\n            const auto rl = _mm256_slli_epi32(last_table[i], SIZE-1), rr = _mm256_srli_epi32(last_table[i], 1);\n            const auto r = _mm256_and_si256(_mm256_or_si256(rl, rr), valid_bits_x8);\n            l_and_r[i] = _mm256_and_si256(l, r);\n            l_xor_r[i] = _mm256_xor_si256(l, r);\n            xor_a_l_r[i] = _mm256_xor_si256(l_xor_r[i], last_table[i]);\n            sum_a_l_r_geq_2[i] = _mm256_or_si256(l_and_r[i], _mm256_and_si256(l_xor_r[i], last_table[i]));\n        }\n        for(int di=-d; di<=d; ++di){\n            const auto i_data = (i0 + SIZE + di) % SIZE;\n            const auto i = buffer_offset + di;\n            const auto u_xor_d = _mm256_xor_si256(xor_a_l_r[i+1], xor_a_l_r[i-1]);\n            const auto u_and_d = _mm256_and_si256(xor_a_l_r[i+1], xor_a_l_r[i-1]);\n            const auto next_sum_mod_2 = _mm256_xor_si256(u_xor_d, l_xor_r[i]);\n            const auto sum_u_d_lr_geq_2 = _mm256_or_si256(_mm256_and_si256(u_xor_d, l_xor_r[i]), u_and_d);\n            const auto next_sum_is_2_or_3 = _mm256_or_si256(\n                _mm256_andnot_si256(\n                    _mm256_or_si256(sum_a_l_r_geq_2[i+1], sum_a_l_r_geq_2[i-1]),\n                    _mm256_xor_si256(sum_u_d_lr_geq_2, l_and_r[i])\n                ),\n                _mm256_andnot_si256(\n                    _mm256_or_si256(sum_u_d_lr_geq_2, l_and_r[i]),\n                    _mm256_xor_si256(sum_a_l_r_geq_2[i+1], sum_a_l_r_geq_2[i-1])\n                )\n            );\n            const auto new_array = _mm256_and_si256(next_sum_is_2_or_3, _mm256_or_si256(last_table[i], next_sum_mod_2));\n            _mm256_store_si256((bitarray_x8_t *)(tables.data[d][i_data] + proc), new_array);\n        }\n    }\n}\n```\n\nAt this point, updating the table was no longer a bottleneck, and we needed to speed up other parts such as random sampling.\n\nApply these accelerations, my solver can compute about **40M iterations/sec.** (sum to all processors. i.e. 2M iter/proc/sec.)\n\n(I wonder how fast it runs if this algorithm is implemented in GPU, but I don't have enough time to do)"},{"metadata":{},"cell_type":"markdown","source":"## time-homogeneous parallel simulated annealing\n\nTo reduce the difficulty of temperature control, I used the idea of **time-homogeneous parallel** simulated annealing, describing in https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.10.7865&rep=rep1&type=pdf .\n\nI set the temperatures `np.logspace(np.log10(20), np.log10(0.1), 32)`, and ran this solver about 25 days, then achive a score about 0.015.\n\n- `25 * 86400(sec/day) * 20(processors) / 50000 (problems) = 864(cpu-sec/problem)`\n- `864 * 2M = 1.728 B(iter/problem)`\n    - some of problems solve easily, then I'd executed about 3B iterations for each problems\n\n## combine with other solvers\n\nIn addition to this, I used z3solver https://www.kaggle.com/jamesmcguigan/game-of-life-z3-constraint-satisfaction for unsolved delta=2 problems, because discussion https://www.kaggle.com/c/conways-reverse-game-of-life-2020/discussion/190796 says that these solver can easily solve low-delta problems. Using this, about 800 problem solved in 4 days, then I can achive a score 0.01315, what is my best score.\n\nAfter all, my submission is combine with submission of my team mate @markuskarmann, and we could achive the great (and so amusing) score of 0.01234.\n\nhis solution is written in https://www.kaggle.com/markuskarmann/3rd-place-solution-part , it is so interesting."},{"metadata":{},"cell_type":"markdown","source":"## statistics\n\nAs shown the following graph, many problems remain unsolved. I'm wondering if it will be better when combined with many other solvers."},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import Image\nImage(\"../input/conway-final-submission-and-stats/score_nonsolve_20201130.png\")","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}