{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Introduction\n\nThis notebook provides background ideas of [the cost minimization notebook](https://www.kaggle.com/saitodevel01/indoor-post-processing-by-cost-minimization).\nHonestly, I was a little wondering about the reaction of the notebook.\nThis method is equivalent to the Kalman smoother under a probabilistic model, neither magic nor trick.\nBut formulation as an optimization problem is more flexible, it is easy to combine other heuristics such that\n[\"Snap to Grid\"](https://www.kaggle.com/robikscube/indoor-navigation-snap-to-grid-post-processing).\nAnd, understanding the probabilistic model behind the cost function allows you to improve post-processing method.","metadata":{}},{"cell_type":"markdown","source":"## Prediction Accuracy: Wi-Fi vs Sensors\n\nIn the early stage of the competition, the accuracy of the wifi model of the public notebooks was about 8m.\nAnd, in the late stage of the competition, [this discussion](https://www.kaggle.com/c/indoor-location-navigation/discussion/232995)\nsuggests that the accuracy of top participants' wifi models (without post-processing) is about 5 to 6m.\nOn the other hand, predicting distances between waypoints using sensor data is much more accurate than pure wifi model.\nThis figure shows a histogram of the prediction error of $\\Delta \\hat{X}$\ndefined in [my notebook](https://www.kaggle.com/saitodevel01/indoor-post-processing-by-cost-minimization).\nThe average error is about 2.7m, and you can get much better results using machine learning.\nThis suggests that, if the position of the previous waypoint is known,\npredicting the position of the next waypoint from sensor data is more accurate than predicting from wifi model.\nA problem is, in test-case, any waypoint is not given, so you must predict at least initial position by other method anyway.\nAnd, using only sensor data will result in significantly worse results because the prediction error is accumulated.\n","metadata":{}},{"cell_type":"code","source":"!pip install kaleido","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import re\nimport json\nimport glob\nfrom dataclasses import dataclass\nimport multiprocessing\nimport numpy as np\nimport pandas as pd\nfrom scipy.interpolate import interp1d\nimport scipy.signal as signal\nimport plotly.graph_objs as go\nfrom plotly.subplots import make_subplots\nfrom PIL import Image\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\n\nINPUT_PATH = '../input/indoor-location-navigation'\n\n@dataclass\nclass ReadData:\n    startTime   : int\n    endTime     : int\n    acce        : np.ndarray\n    acce_uncali : np.ndarray\n    gyro        : np.ndarray\n    gyro_uncali : np.ndarray\n    magn        : np.ndarray\n    magn_uncali : np.ndarray\n    ahrs        : np.ndarray\n    wifi        : np.ndarray\n    ibeacon     : np.ndarray\n    waypoint    : np.ndarray\n\nre_LINE_HEADER = re.compile('\\\\d{13}\\tTYPE_')\ndef to_logical_lines(line):\n    mutch_list = list(re.finditer(re_LINE_HEADER, line))\n    assert(len(mutch_list) > 0)\n    index_list = [m.start() for m in mutch_list] + [len(line)]\n    logical_lines = (line[index_list[i]:index_list[i+1]] for i in range(len(index_list)-1))\n    for logical_line in logical_lines:\n        yield logical_line\n    return\n\ndef read_data_file(data_filename, is_test=False):\n    acce = list()\n    acce_uncali = list()\n    gyro = list()\n    gyro_uncali = list()\n    magn = list()\n    magn_uncali = list()\n    ahrs = list()\n    wifi = list()\n    ibeacon = list()\n    waypoint = list()\n    startTime = None\n    endTime   = None\n\n    with open(data_filename, 'r', encoding='utf-8') as file:\n        lines = file.readlines()\n\n    for line in lines:\n        line = line.strip()\n        if line.startswith('#\\tstartTime'):\n            startTime = int(re.split('[:\\t]', line)[2])\n            continue\n        if line.startswith('#\\tendTime'):\n            endTime = int(re.split('[:\\t]', line)[2])\n            continue\n        if (not line) or line.startswith('#'):\n            continue\n        for line_data in to_logical_lines(line):\n            line_data = line_data.split('\\t')\n\n            if line_data[1] == 'TYPE_ACCELEROMETER':\n                acce.append([int(line_data[0]), float(line_data[2]), float(line_data[3]), float(line_data[4])])\n                continue\n\n            if line_data[1] == 'TYPE_ACCELEROMETER_UNCALIBRATED':\n                if len(line_data) > 5:\n                    acce_uncali.append([int(line_data[0]),\n                                        float(line_data[2]), # x\n                                        float(line_data[3]), # y\n                                        float(line_data[4]), # z\n                                        float(line_data[5]), # x\n                                        float(line_data[6]), # y\n                                        float(line_data[7]), # z\n                                        int(  line_data[8]), # accuracy\n                                        ])\n                else:\n                    acce_uncali.append([int(line_data[0]),\n                                        float(line_data[2]), # x\n                                        float(line_data[3]), # y\n                                        float(line_data[4]), # z\n                                        ])\n                continue\n\n            if line_data[1] == 'TYPE_GYROSCOPE':\n                gyro.append([int(line_data[0]), float(line_data[2]), float(line_data[3]), float(line_data[4])])\n                continue\n\n            if line_data[1] == 'TYPE_GYROSCOPE_UNCALIBRATED':\n                if len(line_data) > 5:\n                    gyro_uncali.append([int(line_data[0]),\n                                        float(line_data[2]), # x\n                                        float(line_data[3]), # y\n                                        float(line_data[4]), # z\n                                        float(line_data[5]), # x\n                                        float(line_data[6]), # y\n                                        float(line_data[7]), # z\n                                        int(  line_data[8]), # accuracy\n                                        ])\n                else:\n                    gyro_uncali.append([int(line_data[0]),\n                                        float(line_data[2]), # x\n                                        float(line_data[3]), # y\n                                        float(line_data[4]), # z\n                                        ])\n                continue\n\n            if line_data[1] == 'TYPE_MAGNETIC_FIELD':\n                magn.append([int(line_data[0]), float(line_data[2]), float(line_data[3]), float(line_data[4])])\n                continue\n\n            if line_data[1] == 'TYPE_MAGNETIC_FIELD_UNCALIBRATED':\n                if len(line_data) > 5:\n                    magn_uncali.append([int(line_data[0]),\n                                        float(line_data[2]), # x\n                                        float(line_data[3]), # y\n                                        float(line_data[4]), # z\n                                        float(line_data[5]), # x\n                                        float(line_data[6]), # y\n                                        float(line_data[7]), # z\n                                        int(  line_data[8]), # accuracy\n                                        ])\n                else:\n                    magn_uncali.append([int(line_data[0]),\n                                        float(line_data[2]), # x\n                                        float(line_data[3]), # y\n                                        float(line_data[4]), # z\n                                        ])\n                continue\n\n            if line_data[1] == 'TYPE_ROTATION_VECTOR':\n                ahrs.append([int(line_data[0]), float(line_data[2]), float(line_data[3]), float(line_data[4])])\n                continue\n\n            if line_data[1] == 'TYPE_WIFI':\n                sys_ts = line_data[0]\n                ssid = line_data[2]\n                bssid = line_data[3]\n                rssi = line_data[4]\n                lastseen_ts = line_data[6]\n                wifi_data = [sys_ts, ssid, bssid, rssi, lastseen_ts]\n                wifi.append(wifi_data)\n                continue\n\n            if line_data[1] == 'TYPE_BEACON':\n                ts = line_data[0]\n                uuid = line_data[2]\n                major = line_data[3]\n                minor = line_data[4]\n                rssi = line_data[6]\n                if is_test:\n                    real_ts = line_data[9]\n                    ibeacon_data = [ts, '_'.join([uuid, major, minor]), rssi, real_ts]\n                else:\n                    ibeacon_data = [ts, '_'.join([uuid, major, minor]), rssi]\n                ibeacon.append(ibeacon_data)\n                continue\n\n            if line_data[1] == 'TYPE_WAYPOINT':\n                waypoint.append([int(line_data[0]), float(line_data[2]), float(line_data[3])])\n                continue\n\n    acce = np.array(acce)\n    acce_uncali = np.array(acce_uncali)\n    gyro = np.array(gyro)\n    gyro_uncali = np.array(gyro_uncali)\n    magn = np.array(magn)\n    magn_uncali = np.array(magn_uncali)\n    ahrs = np.array(ahrs)\n    wifi = np.array(wifi)\n    ibeacon = np.array(ibeacon)\n    waypoint = np.array(waypoint)\n\n    return ReadData(startTime, endTime, acce, acce_uncali, gyro, gyro_uncali, magn, magn_uncali, ahrs, wifi, ibeacon, waypoint)\n\ndef split_ts_seq(ts_seq, sep_ts):\n    \"\"\"\n\n    :param ts_seq:\n    :param sep_ts:\n    :return:\n    \"\"\"\n    tss = ts_seq[:, 0].astype(float)\n    unique_sep_ts = np.unique(sep_ts)\n    ts_seqs = []\n    start_index = 0\n    for i in range(0, unique_sep_ts.shape[0]):\n        end_index = np.searchsorted(tss, unique_sep_ts[i], side='right')\n        if start_index == end_index:\n            continue\n        ts_seqs.append(ts_seq[start_index:end_index, :].copy())\n        start_index = end_index\n\n    # tail data\n    if start_index < ts_seq.shape[0]:\n        ts_seqs.append(ts_seq[start_index:, :].copy())\n\n    return ts_seqs\n\n\ndef correct_trajectory(original_xys, end_xy):\n    \"\"\"\n\n    :param original_xys: numpy ndarray, shape(N, 2)\n    :param end_xy: numpy ndarray, shape(1, 2)\n    :return:\n    \"\"\"\n    corrected_xys = np.zeros((0, 2))\n\n    A = original_xys[0, :]\n    B = end_xy\n    Bp = original_xys[-1, :]\n\n    angle_BAX = np.arctan2(B[1] - A[1], B[0] - A[0])\n    angle_BpAX = np.arctan2(Bp[1] - A[1], Bp[0] - A[0])\n    angle_BpAB = angle_BpAX - angle_BAX\n    AB = np.sqrt(np.sum((B - A) ** 2))\n    ABp = np.sqrt(np.sum((Bp - A) ** 2))\n\n    corrected_xys = np.append(corrected_xys, [A], 0)\n    for i in np.arange(1, np.size(original_xys, 0)):\n        angle_CpAX = np.arctan2(original_xys[i, 1] - A[1], original_xys[i, 0] - A[0])\n\n        angle_CAX = angle_CpAX - angle_BpAB\n\n        ACp = np.sqrt(np.sum((original_xys[i, :] - A) ** 2))\n\n        AC = ACp * AB / ABp\n\n        delta_C = np.array([AC * np.cos(angle_CAX), AC * np.sin(angle_CAX)])\n\n        C = delta_C + A\n\n        corrected_xys = np.append(corrected_xys, [C], 0)\n\n    return corrected_xys\n\n\ndef correct_positions(rel_positions, reference_positions):\n    \"\"\"\n\n    :param rel_positions:\n    :param reference_positions:\n    :return:\n    \"\"\"\n    rel_positions_list = split_ts_seq(rel_positions, reference_positions[:, 0])\n    if len(rel_positions_list) != reference_positions.shape[0] - 1:\n        # print(f'Rel positions list size: {len(rel_positions_list)}, ref positions size: {reference_positions.shape[0]}')\n        del rel_positions_list[-1]\n    assert len(rel_positions_list) == reference_positions.shape[0] - 1\n\n    corrected_positions = np.zeros((0, 3))\n    for i, rel_ps in enumerate(rel_positions_list):\n        start_position = reference_positions[i]\n        end_position = reference_positions[i + 1]\n        abs_ps = np.zeros(rel_ps.shape)\n        abs_ps[:, 0] = rel_ps[:, 0]\n        # abs_ps[:, 1:3] = rel_ps[:, 1:3] + start_position[1:3]\n        abs_ps[0, 1:3] = rel_ps[0, 1:3] + start_position[1:3]\n        for j in range(1, rel_ps.shape[0]):\n            abs_ps[j, 1:3] = abs_ps[j-1, 1:3] + rel_ps[j, 1:3]\n        abs_ps = np.insert(abs_ps, 0, start_position, axis=0)\n        corrected_xys = correct_trajectory(abs_ps[:, 1:3], end_position[1:3])\n        corrected_ps = np.column_stack((abs_ps[:, 0], corrected_xys))\n        if i == 0:\n            corrected_positions = np.append(corrected_positions, corrected_ps, axis=0)\n        else:\n            corrected_positions = np.append(corrected_positions, corrected_ps[1:], axis=0)\n\n    corrected_positions = np.array(corrected_positions)\n\n    return corrected_positions\n\n\ndef init_parameters_filter(sample_freq, warmup_data, cut_off_freq=2):\n    order = 4\n    filter_b, filter_a = signal.butter(order, cut_off_freq / (sample_freq / 2), 'low', False)\n    zf = signal.lfilter_zi(filter_b, filter_a)\n    _, zf = signal.lfilter(filter_b, filter_a, warmup_data, zi=zf)\n    _, filter_zf = signal.lfilter(filter_b, filter_a, warmup_data, zi=zf)\n\n    return filter_b, filter_a, filter_zf\n\n\ndef get_rotation_matrix_from_vector(rotation_vector):\n    q1 = rotation_vector[0]\n    q2 = rotation_vector[1]\n    q3 = rotation_vector[2]\n\n    if rotation_vector.size >= 4:\n        q0 = rotation_vector[3]\n    else:\n        q0 = 1 - q1*q1 - q2*q2 - q3*q3\n        if q0 > 0:\n            q0 = np.sqrt(q0)\n        else:\n            q0 = 0\n\n    sq_q1 = 2 * q1 * q1\n    sq_q2 = 2 * q2 * q2\n    sq_q3 = 2 * q3 * q3\n    q1_q2 = 2 * q1 * q2\n    q3_q0 = 2 * q3 * q0\n    q1_q3 = 2 * q1 * q3\n    q2_q0 = 2 * q2 * q0\n    q2_q3 = 2 * q2 * q3\n    q1_q0 = 2 * q1 * q0\n\n    R = np.zeros((9,))\n    if R.size == 9:\n        R[0] = 1 - sq_q2 - sq_q3\n        R[1] = q1_q2 - q3_q0\n        R[2] = q1_q3 + q2_q0\n\n        R[3] = q1_q2 + q3_q0\n        R[4] = 1 - sq_q1 - sq_q3\n        R[5] = q2_q3 - q1_q0\n\n        R[6] = q1_q3 - q2_q0\n        R[7] = q2_q3 + q1_q0\n        R[8] = 1 - sq_q1 - sq_q2\n\n        R = np.reshape(R, (3, 3))\n    elif R.size == 16:\n        R[0] = 1 - sq_q2 - sq_q3\n        R[1] = q1_q2 - q3_q0\n        R[2] = q1_q3 + q2_q0\n        R[3] = 0.0\n\n        R[4] = q1_q2 + q3_q0\n        R[5] = 1 - sq_q1 - sq_q3\n        R[6] = q2_q3 - q1_q0\n        R[7] = 0.0\n\n        R[8] = q1_q3 - q2_q0\n        R[9] = q2_q3 + q1_q0\n        R[10] = 1 - sq_q1 - sq_q2\n        R[11] = 0.0\n\n        R[12] = R[13] = R[14] = 0.0\n        R[15] = 1.0\n\n        R = np.reshape(R, (4, 4))\n\n    return R\n\n\ndef get_orientation(R):\n    flat_R = R.flatten()\n    values = np.zeros((3,))\n    if np.size(flat_R) == 9:\n        values[0] = np.arctan2(flat_R[1], flat_R[4])\n        values[1] = np.arcsin(-flat_R[7])\n        values[2] = np.arctan2(-flat_R[6], flat_R[8])\n    else:\n        values[0] = np.arctan2(flat_R[1], flat_R[5])\n        values[1] = np.arcsin(-flat_R[9])\n        values[2] = np.arctan2(-flat_R[8], flat_R[10])\n\n    return values\n\n\ndef compute_steps(acce_datas):\n    step_timestamps = np.array([])\n    step_indexs = np.array([], dtype=int)\n    step_acce_max_mins = np.zeros((0, 4))\n    sample_freq = 50\n    window_size = 22\n    low_acce_mag = 0.6\n    step_criterion = 1\n    interval_threshold = 250\n\n    acce_max = np.zeros((2,))\n    acce_min = np.zeros((2,))\n    acce_binarys = np.zeros((window_size,), dtype=int)\n    acce_mag_pre = 0\n    state_flag = 0\n\n    warmup_data = np.ones((window_size,)) * 9.81\n    filter_b, filter_a, filter_zf = init_parameters_filter(sample_freq, warmup_data)\n    acce_mag_window = np.zeros((window_size, 1))\n\n    # detect steps according to acceleration magnitudes\n    for i in np.arange(0, np.size(acce_datas, 0)):\n        acce_data = acce_datas[i, :]\n        acce_mag = np.sqrt(np.sum(acce_data[1:] ** 2))\n\n        acce_mag_filt, filter_zf = signal.lfilter(filter_b, filter_a, [acce_mag], zi=filter_zf)\n        acce_mag_filt = acce_mag_filt[0]\n\n        acce_mag_window = np.append(acce_mag_window, [acce_mag_filt])\n        acce_mag_window = np.delete(acce_mag_window, 0)\n        mean_gravity = np.mean(acce_mag_window)\n        acce_std = np.std(acce_mag_window)\n        mag_threshold = np.max([low_acce_mag, 0.4 * acce_std])\n\n        # detect valid peak or valley of acceleration magnitudes\n        acce_mag_filt_detrend = acce_mag_filt - mean_gravity\n        if acce_mag_filt_detrend > np.max([acce_mag_pre, mag_threshold]):\n            # peak\n            acce_binarys = np.append(acce_binarys, [1])\n            acce_binarys = np.delete(acce_binarys, 0)\n        elif acce_mag_filt_detrend < np.min([acce_mag_pre, -mag_threshold]):\n            # valley\n            acce_binarys = np.append(acce_binarys, [-1])\n            acce_binarys = np.delete(acce_binarys, 0)\n        else:\n            # between peak and valley\n            acce_binarys = np.append(acce_binarys, [0])\n            acce_binarys = np.delete(acce_binarys, 0)\n\n        if (acce_binarys[-1] == 0) and (acce_binarys[-2] == 1):\n            if state_flag == 0:\n                acce_max[:] = acce_data[0], acce_mag_filt\n                state_flag = 1\n            elif (state_flag == 1) and ((acce_data[0] - acce_max[0]) <= interval_threshold) and (\n                    acce_mag_filt > acce_max[1]):\n                acce_max[:] = acce_data[0], acce_mag_filt\n            elif (state_flag == 2) and ((acce_data[0] - acce_max[0]) > interval_threshold):\n                acce_max[:] = acce_data[0], acce_mag_filt\n                state_flag = 1\n\n        # choose reasonable step criterion and check if there is a valid step\n        # save step acceleration data: step_acce_max_mins = [timestamp, max, min, variance]\n        step_flag = False\n        if step_criterion == 2:\n            if (acce_binarys[-1] == -1) and ((acce_binarys[-2] == 1) or (acce_binarys[-2] == 0)):\n                step_flag = True\n        elif step_criterion == 3:\n            if (acce_binarys[-1] == -1) and (acce_binarys[-2] == 0) and (np.sum(acce_binarys[:-2]) > 1):\n                step_flag = True\n        else:\n            if (acce_binarys[-1] == 0) and acce_binarys[-2] == -1:\n                if (state_flag == 1) and ((acce_data[0] - acce_min[0]) > interval_threshold):\n                    acce_min[:] = acce_data[0], acce_mag_filt\n                    state_flag = 2\n                    step_flag = True\n                elif (state_flag == 2) and ((acce_data[0] - acce_min[0]) <= interval_threshold) and (\n                        acce_mag_filt < acce_min[1]):\n                    acce_min[:] = acce_data[0], acce_mag_filt\n        if step_flag:\n            step_timestamps = np.append(step_timestamps, acce_data[0])\n            step_indexs = np.append(step_indexs, [i])\n            step_acce_max_mins = np.append(step_acce_max_mins,\n                                           [[acce_data[0], acce_max[1], acce_min[1], acce_std ** 2]], axis=0)\n        acce_mag_pre = acce_mag_filt_detrend\n\n    return step_timestamps, step_indexs, step_acce_max_mins\n\n\ndef compute_stride_length(step_acce_max_mins):\n    K = 0.4\n    K_max = 0.8\n    K_min = 0.4\n    para_a0 = 0.21468084\n    para_a1 = 0.09154517\n    para_a2 = 0.02301998\n\n    stride_lengths = np.zeros((step_acce_max_mins.shape[0], 2))\n    k_real = np.zeros((step_acce_max_mins.shape[0], 2))\n    step_timeperiod = np.zeros((step_acce_max_mins.shape[0] - 1, ))\n    stride_lengths[:, 0] = step_acce_max_mins[:, 0]\n    window_size = 2\n    step_timeperiod_temp = np.zeros((0, ))\n\n    # calculate every step period - step_timeperiod unit: second\n    for i in range(0, step_timeperiod.shape[0]):\n        step_timeperiod_data = (step_acce_max_mins[i + 1, 0] - step_acce_max_mins[i, 0]) / 1000\n        step_timeperiod_temp = np.append(step_timeperiod_temp, [step_timeperiod_data])\n        if step_timeperiod_temp.shape[0] > window_size:\n            step_timeperiod_temp = np.delete(step_timeperiod_temp, [0])\n        step_timeperiod[i] = np.sum(step_timeperiod_temp) / step_timeperiod_temp.shape[0]\n\n    # calculate parameters by step period and acceleration magnitude variance\n    k_real[:, 0] = step_acce_max_mins[:, 0]\n    k_real[0, 1] = K\n    for i in range(0, step_timeperiod.shape[0]):\n        k_real[i + 1, 1] = np.max([(para_a0 + para_a1 / step_timeperiod[i] + para_a2 * step_acce_max_mins[i, 3]), K_min])\n        k_real[i + 1, 1] = np.min([k_real[i + 1, 1], K_max]) * (K / K_min)\n\n    # calculate every stride length by parameters and max and min data of acceleration magnitude\n    stride_lengths[:, 1] = np.max([(step_acce_max_mins[:, 1] - step_acce_max_mins[:, 2]),\n                                   np.ones((step_acce_max_mins.shape[0], ))], axis=0)**(1 / 4) * k_real[:, 1]\n\n    return stride_lengths\n\n\ndef compute_headings(ahrs_datas):\n    headings = np.zeros((np.size(ahrs_datas, 0), 2))\n    for i in np.arange(0, np.size(ahrs_datas, 0)):\n        ahrs_data = ahrs_datas[i, :]\n        rot_mat = get_rotation_matrix_from_vector(ahrs_data[1:])\n        azimuth, pitch, roll = get_orientation(rot_mat)\n        around_z = (-azimuth) % (2 * np.pi)\n        headings[i, :] = ahrs_data[0], around_z\n    return headings\n\n\ndef compute_step_heading(step_timestamps, headings):\n    step_headings = np.zeros((len(step_timestamps), 2))\n    step_timestamps_index = 0\n    for i in range(0, len(headings)):\n        if step_timestamps_index < len(step_timestamps):\n            if headings[i, 0] == step_timestamps[step_timestamps_index]:\n                step_headings[step_timestamps_index, :] = headings[i, :]\n                step_timestamps_index += 1\n        else:\n            break\n    assert step_timestamps_index == len(step_timestamps)\n\n    return step_headings\n\n\ndef compute_rel_positions(stride_lengths, step_headings):\n    rel_positions = np.zeros((stride_lengths.shape[0], 3))\n    for i in range(0, stride_lengths.shape[0]):\n        rel_positions[i, 0] = stride_lengths[i, 0]\n        rel_positions[i, 1] = -stride_lengths[i, 1] * np.sin(step_headings[i, 1])\n        rel_positions[i, 2] = stride_lengths[i, 1] * np.cos(step_headings[i, 1])\n\n    return rel_positions\n\ndef compute_step_prediction(example):\n    T_ref   = example.waypoint[:, 0]\n    xy_true = example.waypoint[:, 1:]\n    delta_xy_true = np.diff(xy_true, axis=0)\n        \n    acce_datas = example.acce\n    ahrs_datas = example.ahrs\n    step_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\n    headings = compute_headings(ahrs_datas)\n    stride_lengths = compute_stride_length(step_acce_max_mins)\n    step_headings  = compute_step_heading(step_timestamps, headings)\n    rel_positions  = compute_rel_positions(stride_lengths, step_headings)\n    if T_ref[-1] > rel_positions[-1, 0]:\n        rel_positions = [np.array([[0, 0, 0]]), rel_positions, np.array([[T_ref[-1], 0, 0]])]\n    else:\n        rel_positions = [np.array([[0, 0, 0]]), rel_positions]\n    rel_positions = np.concatenate(rel_positions)\n    T_rel = rel_positions[:, 0]\n    delta_xy_pred = np.diff(interp1d(T_rel, np.cumsum(rel_positions[:, 1:], axis=0), axis=0)(T_ref), axis=0)\n    \n    return delta_xy_true, delta_xy_pred\n\ndef compute_step_prediction_error(filename):\n    example = read_data_file(filename)\n    delta_xy_true, delta_xy_pred = compute_step_prediction(example)\n    delta_xy_err = delta_xy_pred - delta_xy_true\n    err_distance = np.sqrt(np.sum(delta_xy_err**2, axis=1))\n    return err_distance\n\ndef visualize_trajectory(waypoint,\n                         trajectries,\n                         floor_plan_filename,\n                         width_meter,\n                         height_meter,\n                         radius=6.0):\n    fig = go.Figure()\n\n    fig.add_trace(\n        go.Scattergl(\n            x=waypoint[:, 0],\n            y=waypoint[:, 1],\n            mode='markers',\n            marker=dict(size=6, color='blue', symbol='circle'),\n            name='true waypoint',\n        ))\n\n    for i in range(waypoint.shape[0]):\n        x0 = waypoint[i, 0] - radius\n        y0 = waypoint[i, 1] - radius\n        x1 = waypoint[i, 0] + radius\n        y1 = waypoint[i, 1] + radius\n        fig.add_shape(\n            type=\"circle\",\n            xref=\"x\", yref=\"y\",\n            x0=x0, y0=y0, x1=x1, y1=y1,\n            line_color=\"LightSeaGreen\",\n        )\n\n    for index, trajectory in enumerate(trajectries):\n        fig.add_trace(\n            go.Scattergl(\n                x=trajectory[:, 0],\n                y=trajectory[:, 1],\n                mode='lines + markers',\n                line=dict(shape='linear', color='green', width=1, dash='solid'),\n                name=f'trajectory_{index}',\n            ))\n\n    # add floor plan\n    floor_plan = Image.open(floor_plan_filename)\n    fig.update_layout(images=[\n        go.layout.Image(\n            source=floor_plan,\n            xref=\"x\",\n            yref=\"y\",\n            x=0,\n            y=height_meter,\n            sizex=width_meter,\n            sizey=height_meter,\n            sizing=\"contain\",\n            opacity=1,\n            layer=\"below\",\n        )\n    ])\n\n    # configure\n    fig.update_xaxes(autorange=False, range=[0, width_meter])\n    fig.update_yaxes(autorange=False, range=[0, height_meter], scaleanchor=\"x\", scaleratio=1)\n    fig.update_layout(\n        autosize=True,\n        width=900,\n        height=50 + 900 * height_meter / width_meter,\n        template=\"plotly_white\",\n        margin=dict(\n            l=0,\n            r=0,\n            b=0,\n            t=0,\n            pad=0\n        ),\n    )\n    return fig","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = pd.read_csv(f'{INPUT_PATH}/sample_submission.csv')\ntmp = sub['site_path_timestamp'].apply(lambda x: pd.Series(x.split('_')))\nsites = np.unique(tmp[0])\nprocesses = multiprocessing.cpu_count()\nERR = []\nfor site in sites:\n    filelist = glob.glob(f'{INPUT_PATH}/train/{site}/*/*.txt')\n    with multiprocessing.Pool(processes=processes) as pool:\n        err = pool.imap_unordered(compute_step_prediction_error, filelist)\n        # err = tqdm(err)\n        err = list(err)\n    ERR.extend(err)\nERR = np.concatenate(ERR, axis=0)\nERR_mean = np.mean(ERR)\nplt.figure(figsize=(6, 4))\nplt.hist(ERR, bins=100)\nplt.grid(True)\nplt.xlim([0, 30])\nplt.ylim([0, 6000])\nplt.plot([ERR_mean, ERR_mean], [0, 6000], color='red')\nplt.xlabel('prediction error [m]')\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"site  = \"5cd56865eb294480de7167b6\"\nfloor = \"F2\"\npath  = \"5cfdd9006fa436000a02de9e\"\n\nfloor_plan_filename = f'{INPUT_PATH}/metadata/{site}/{floor}/floor_image.png'\nfloor_json_filename = f'{INPUT_PATH}/metadata/{site}/{floor}/floor_info.json'\npath_filename       = f'{INPUT_PATH}/train/{site}/{floor}/{path}.txt'\n\nwith open(floor_json_filename) as f:\n    json_data = json.load(f)\nwidth_meter  = json_data['map_info']['width']\nheight_meter = json_data['map_info']['height']\n\nexample = read_data_file(path_filename)\nwaypoint = example.waypoint[:, 1:]\n_, delta_xy_pred = compute_step_prediction(example)\n\ntrajectries_1 = []\nfor i in range(delta_xy_pred.shape[0]):\n    trajectry = waypoint[i, :] + np.concatenate([np.zeros((1, 2)), delta_xy_pred[i:i+1, :]])\n    trajectries_1.append(trajectry)\nfig_1 = visualize_trajectory(\n    waypoint,\n    trajectries_1,\n    floor_plan_filename,\n    width_meter,\n    height_meter,\n)\nfig_1.write_image('trajectries_1.png')\n\ntrajectries_2 = [waypoint[0, :] + np.cumsum(np.concatenate([np.zeros((1, 2)), delta_xy_pred]), axis=0)]\nfig_2 = visualize_trajectory(\n    waypoint,\n    trajectries_2,\n    floor_plan_filename,\n    width_meter,\n    height_meter,\n)\nfig_2.write_image('trajectries_2.png')\n\nfig = plt.figure(figsize=(16.0, 8.0))\nax1 = fig.add_subplot(1, 2, 1)\nax1.imshow(np.array(Image.open('trajectries_1.png')))\nax1.tick_params(labelbottom=False,\n                labelleft=False,\n                labelright=False,\n                labeltop=False)\nax1.spines['right'].set_visible(False)\nax1.spines['left'].set_visible(False)\nax1.spines['top'].set_visible(False)\nax1.spines['bottom'].set_visible(False)\nax1.tick_params('x', length=0, which='major')\nax1.tick_params('y', length=0, which='major')\n\nax2 = fig.add_subplot(1, 2, 2)\nax2.imshow(np.array(Image.open('trajectries_2.png')))\nax2.tick_params(labelbottom=False,\n                labelleft=False,\n                labelright=False,\n                labeltop=False)\nax2.spines['right'].set_visible(False)\nax2.spines['left'].set_visible(False)\nax2.spines['top'].set_visible(False)\nax2.spines['bottom'].set_visible(False)\nax2.tick_params('x', length=0, which='major')\nax2.tick_params('y', length=0, which='major')\n\nfig.tight_layout()\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Caption: Visualizations of predictions based on $\\Delta \\hat{X}$ defined in [my notebook](https://www.kaggle.com/saitodevel01/indoor-post-processing-by-cost-minimization).\n(Left) Single step prediction. (Right) Multiple step prediction.\nThe radius of circles centered on waypoints is 6m.","metadata":{}},{"cell_type":"markdown","source":"## Problem Formulation as Linear Dynamical System\n\nLinear Dynamical System (LDS) or Linear Gaussian State Space Model (LGSSM) is defined by following equation and graphical model.\n$$\n\\begin{array}{ll}\nz_n = A z_{n-1} + w_n & w_n \\sim \\mathcal{N}(w_n \\mid 0, \\Gamma) \\\\\nx_n = C z_n     + v_n & v_n \\sim \\mathcal{N}(v_n \\mid 0, \\Sigma)\n\\end{array}\n$$\nwhere, $z$ is the latent variable and $x$ is the observed variable.\n$A, C, \\Gamma, \\Sigma$ are parameters of LDS.\n\n![graphical_model.png](attachment:12eebeca-d63c-4b8c-96ad-63f9f9296fa7.png)\n\nThe initial distribution of the latent variable is also given as a gaussian distribution.\n$$\np(z_1) = \\mathcal{N}(z_1 \\mid \\mu_0, P_0)\n$$\n$P_0^{-1} = 0$ means that there is no prior information.\n\nThere is another variant of LDS with observable external input $u$.\n$$\n\\begin{array}{ll}\nz_n = A z_{n-1} + B u_n + w_n & w_n \\sim \\mathcal{N}(w_n \\mid 0, \\Gamma) \\\\\nx_n = C z_n     + D u_n + v_n & v_n \\sim \\mathcal{N}(v_n \\mid 0, \\Sigma)\n\\end{array}\n$$\n\nFor this competition, transition process between waypoints can be modeled as the following LDS,\n$$\n\\begin{array}{ll}\nz_n = z_{n-1} + \\Delta z_n + w_n & w_n \\sim \\mathcal{N}(w_n \\mid 0, \\Gamma) \\\\\nx_n = z_n     + v_n              & v_n \\sim \\mathcal{N}(v_n \\mid 0, \\Sigma)\n\\end{array}\n$$\nTherefore, true waypoint position (the latent variable of LDS) can be estimated by inference of LDS called the Kalman filter and smoother.\n\n| symbol | description | |\n|---|---|:---:|\n| $z_n$        | true waypoint position | (unknown) | \n| $\\Delta z_n$ | relative position prediction of sensor model between two waypoints | (known) |\n| $w_n$        | prediction error of sensor model | (unknown) |\n| $x_n$        | prediction of wifi model | (known) |\n| $v_n$        | prediction error of wifi model | (unknown) |","metadata":{},"attachments":{"12eebeca-d63c-4b8c-96ad-63f9f9296fa7.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAjMAAACICAYAAAASwGpdAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydd1gUV9uHf7P0Kr2IBVSaBbsCAgpqsKCCHRUsURPN90ajsfe8MYk9akwvKqgxib03ir0XBAVsFBFEUTrS9vn+MOwrSlmY2Zndde7r4roizJ7nB3fO7NmZM+cwREQQEREREREREVFNwiVCJxARERERERERYYM4mBERERERERFRacTBjIiIiIiIiIhKoylE0fz8fBQUFKCgoACmpqYwMDCAtra2EFFEIPpQNkQfyoXoQ7kQfSgXyuJDoYMZqVSKmJgYREZG4sKFC0hMTERiYiKKioreOdbCwgIuLi5o1aoVvL294evri4YNGyoy3ntHVT4SEhLw6tWrd461tLSEs7MzWrVqBR8fH/j6+sLW1laA1OqLVCrFrVu3EBkZiYsXL9bqw8XFBS1bthR9KIi3fSQkJCAxMbFGH2+er0Qf3CKVSnHz5k1ERUXVyYePjw969Ogh+uCYCh+RkZG4dOlSnXz4+vrCxsZGofkYRTzNFBcXhy1btmDbtm148uQJzM3N4eXlBVdXVzg5OcHa2hoGBgYwMDBAdnY28vPzkZqaivj4eMTExODy5csoKSmBu7s7QkNDERwcDBMTE65jvjfExcVh8+bN2L59O548eQILCwt069atRh8pKSlISEhATEwMLl26hNLSUnh4eCAkJET0wZIKH9u2bUN6enolH87OzrCyspLbR2hoKEaOHCn6YEFsbKysf1T48PLygouLS40+3jxflZaWwtPTEyEhIaIPllTn483zlb6+vlw+KvpHgwYNhP61VJbbt2/LfGRkZMDS0vKd94+3fSQnJ8vOVzz5CAdxyNmzZ6lfv37EMAzZ29vTokWL6ObNm1ReXl6ndvLz8+nw4cM0ZswYMjAwICMjI5o1axalp6dzGVftOXPmjMyHg4MDLV68uN4+Dh06VMnH7NmzKSMjQ0HJ1ZMzZ85Q3759CQAnPkaPHk36+vqij3py+vRp6tOnDwGgZs2a0ZIlS+jWrVsklUrr1E5+fj4dPHhQ5sPY2JjmzJlDT58+VVBy9USRPubOnSv6qCPR0dHk7+/PmY9Ro0Yp0kcYJ4OZ1NRUGjJkCAEgLy8vOnToUJ1/4erIzc2lVatWka2tLRkaGtKqVauopKSEk7bVlZSUlEo+Dh8+zLkPGxsb0YecpKSk0ODBgwkAeXt7c+ojJyeHVq5cSTY2NmRkZESrV68WfdRCcnIyBQUFEQDy8fGhI0eOcOpjxYoVZG1tTUZGRrRmzRoqLS3lpG11JTk5mQIDA2U+jh49qjAfa9euFX3UQlJSksxH9+7dOffxzTffkLW1NRkbG3Ppg/1g5rfffiNDQ0NydHSko0ePchGqSoqKimjZsmWkp6dHbm5uFBcXp7Baqsyvv/7Ki4/CwkJaunQp6erqkpubG925c0dhtVSZX375hQwMDMjJyYmOHTumsDpv+mjbtq3ooxp+/vlnMjAwIGdnZzp+/LjC6hQWFtKSJUtIV1eX2rVrR3fv3lVYLVXmp59+kvk4ceKEwuoUFhbS4sWLZT7i4+MVVkuV+fHHH3nxUVBQQIsXLyYdHR1q3749Fz7qP5gpLCykMWPGEMMwNGfOHHr16hXbMHLx4MED8vDwIAMDAwoLC+OlpipQ4UMikdDcuXN583H//n3RRxUUFBTQ6NGjSSKR0Lx583j14e7uTgYGBhQeHs5LTVWgoKCAgoODSSKR0Pz583nzce/ePeratSsZGhrStm3beKmpCuTn58t8LFiwgIqLi3mpm5iYSF26dCFDQ0Pavn07LzVVgfz8fBo5cqRgPoyMjGjHjh1smqrfYObly5fk5eVF5ubmdOTIETYB6kVJSQl9/vnnxDAMffnll7zXVzZevHhB3bp1I3Nzc4VejamOkpISmjlzJjEMQ8uXL+e9vrLxpg9FXo2pjpKSEpoxYwYxDENfffUV7/WVjaysLPL09CQLCwuFXo2pjpKSEvrss8+IYRj6+uuvea+vbGRlZZGHhwdZWFgo9NN/dRQXF9P06dOJYRj65ptveK+vbGRlZZG7uztZWloK5mPatGnEMAytWLGivs3UfTCTlZVFbm5u1LhxY8EvZX///fekoaFBn3/+uaA5hESZfGzatIkkEgnNmjVL0BxCkpWVRW3atKEmTZoIfmuhwsfs2bMFzSEkz58/p9atW1PTpk0Fv7WwceNG2ZXT95Vnz55Rq1atlMLHhg0bZFdO31cqfNjb21NCQoKgWdavXy+7cloP6jaYKSgoIE9PT2rSpAklJyfXpyDnbNu2jSQSCZsRncpSUFBAHh4e1LRpU0pJSRE6DhERhYeHk0QioZUrVwodhXfy8/OVzkdYWBhJJBJatWqV0FF4Jz8/n7p27UoODg6UmpoqdBwiItqyZQsxDENr1qwROgrv5OXlKZ2PzZs3E8MwtHbtWqGj8E5eXh516dKFmjVrRo8fPxY6DhER/fHHH8QwDK1bt66uL63bYGbYsGFkYWEh+CfOt9mwYQMxDEO7d+8WOgqvDB06lCwtLQX/hPM269evJ4ZhaM+ePUJH4Q2pVEpDhgwhS0tLwT/hvM23335LDMPQ3r17hY7CG1KplIKCgsjKyooSExOFjlOJtWvXEsMwtG/fPqGj8IZUKqXAwECl9LFmzRqSSCS0f/9+oaPwhlQqpUGDBpG1tTXdu3dP6DiVWL16NUkkEjpw4EBdXib/YKbiEqkQ99Tk4aOPPiITExN68OCB0FF4oeIS6cmTJ4WOUiWTJ08mExMTevjwodBReKHiEqmy+pg0adJ75WPdunWkoaFBp06dEjpKlXz44Ydkamr63vhYu3YtaWpq0unTp4WOUiUTJkwgU1NTevTokdBReGHNmjWkqalJZ86cETpKlYwfP76uPuQbzCQmJpKuri4tW7as3uEUTVFREbVt25a8vb05eyZeWYmPjycdHR364osvhI5SLYWFheTm5kY+Pj7vjY///ve/QkeplsLCQmrTpg11795d7X3cvXuXdHR0lHoyekFBAbVu3Zp8fX3V3sedO3dIW1tbqSejFxQUUKtWrd4rH8o8Gb3Ch5+fn7wvkW8w069fP2rdurXSL8Z17do10tDQoM2bNwsdRaH07duX2rRpo/SLP129epU0NDRoy5YtQkdRKD179lQpH1u3bhU6ikLp2bMndejQgcrKyoSOUiNXrlwhiUSi9o/Q+/n5qZQPdX+E3s/Pjzp27Kj0Pi5fvkwSiUTeR+hrH8wcOHCAGIahc+fOsU/HA5988gnZ2NhQQUGB0FEUwr59+4hhGDp//rzQUeRi6tSpZGNjQ4WFhUJHUQh79+4liURCFy5cEDqKXEyZMoVsbW3V1seePXtIIpHQxYsXhY4iFx999BE1bNiQioqKhI6iEHbv3k0SiYQuXbokdBS5mDx5MtnZ2amtj127dpFEIqHLly8LHUUuJk2aRHZ2dvKsC1X7YMbT05MGDhzITTIeePr0Kenr69P69euFjqIQPDw8KDAwUOgYcpORkUH6+vq0YcMGoaMoBA8PDwoKChI6htxkZGSQnp4ebdy4UegoCqFz5840ePBgoWPITXp6Ounp6dGmTZuEjqIQOnXqREOGDBE6htxU+Pj++++FjqIQOnXqREOHDhU6htxU+Pjhhx9qO7TmwUx0dDQBUJlPnRVMmzaNGjdurPSX0epKZGQkAVCZTzkVfPrpp9SkSRO18xEREUEAVOZTTgX/+c9/1NLHqVOnCABduXJF6Ch14pNPPiF7e/s6bziq7Jw8eZIA0NWrV4WOUiemTp1KDg4OaufjxIkTBICuXbsmdJQ6MWXKFGrWrFltPmoezIwdO5Y6d+7MbTIeuHfvHjEMI8hquIokJCSEunbtKnSMOlPhQ4jVcBXJmDFjyN3dXegYdSYxMZEYhhFkNVxFMnr0aPLw8BA6Rp1JSEggAEr7JFx9GTVqFHl6egodo87cuXOHACjtk3D1JTg4mLp16yZ0jDpT4SMiIqKmw8IkqIaCggLs3r0boaGh1R3CORkZGYiKimLdTosWLeDu7o6wsDD2oZSE/Px83nzk5eXhp59+wty5c/Hrr7+isLCQVXstWrRA165d1c7Hnj17ePGRnZ2NNWvWYNq0aTh+/DjKy8tZtefo6IguXbqolY+8vDzs3buX1/MVAGRlZeHrr79m1YaTkxM6d+4s+mDBoUOHsGPHDtnXypUr633ecnV1RadOndTKR25uLvbt28d7/7h16xY2btyIn376CY8fP65XG66urujYsWPtPqob5hw+fJgYhqH09HTuh1pvkZmZSTNnziQ9PT369NNPOWlz3bp1ZGZmpjaXCg8ePEgMw9DTp08VWic+Pp5sbGzI0dGRtLW1CQA1b96c9f8Ha9euJXNzc7XxUTExPjMzU6F1srKyqHnz5hQSEkJ+fn4kkUioS5curNtds2aNWvnYv38/SSQSevbsGa91AwMDydramnU7q1atIktLS7V5LHjfvn28+bh79y4xDEMAZF8jR45k1ebKlSvJyspKbXxUPKjw/PlzXuo9e/aMPvzwQ+rbty8nuwWsWLGiNh/VX5mJiIhAy5YtYWNjU6/RVF1ISkpCaGgoioqKOGvTz88PL168wM2bNzlrky1lZWX1/lQdERGB1q1bw8rKiuNUlfnss89w7NgxJCYm4vHjx5g4cSIePHiABQsWsGrXz88PWVlZuHXrFkdJ2cPWR5s2bWBpaclxqsr89ddfuHz5MrZu3YpTp05h6dKluHz5Ms6dO8eq3QofMTExHCVlDxc+LCwsOE5VPb/88gvi4uI4acvPzw/Pnj3D7du3OWmPC9j6cHNz48XH2rVrERERgZSUFNnXH3/8wapNPz8/ZGZmIjY2lqOU7GHro23btjA3N+c41bskJSXB1dUVxcXFOHz4MJo0acK6zQofNfW3agcz58+fR/fu3VmHkIfOnTvDxcWF0zbbtGkDc3Nz1id9LklNTYWdnR2mTZuGixcvgojkfi0fPq5du4bRo0fDzc0NAGBpaYkvvvgCEokE58+fZ9W2m5sbzMzMlMpHSkqKUvsoKSmBv78/zMzMZN+ruExsbGzMqm03NzeYmpoqlY/k5ORKPuoCn+crAEhMTMSNGzcQEBDASXvt2rWDiYmJUvlISkqCnZ0dpk+fjkuXLtXptXz5yMjIQExMDFq0aIHGjRvLvnR1dVm12759ezRo0ECpfDx69AiNGjVSah8lJSUYPnw4zMzM8OOPP3LWbvv27WFsbFyjD83qfhAfH4/Ro0ezDhEWFlblaLJNmzbo2LEj6/arg2EYuLq6Ij4+XmE16sPTp0+xadMmbNiwAXZ2dhg7diyCg4PRunXrGl8XHx+PsWPHsq5fkw97e3t06NCh0vdtbW3RsWNHaGpW+7+KXKiKj3HjxiE4OBitWrWq8XUJCQkYP3486/q19Q8HB4dK34+JiUFAQADatGnDqq5EIlErH/Hx8fjwww9Z15fnfFVaWoqFCxfit99+w5IlS1jXBF77cHFxUUof3333HdavXw87OzuMHz8ewcHBaNmyZY2vS0hIwKRJk1jXr83Hxo0bcenSJTRu3BgODg5YvHgxxo4dC4ZhWNVVVh8ZGRkyH40aNZL1j5p8EBESEhLw0Ucfsa5fm48FCxbgypUr+PXXX2FgYMC6XgUaGhq1+6jq5lNmZiZns+vbt29Pe/bsoZiYGLp58ya5urqSnp7eO5sjFhcXEwDO5swQEU2cOLEuyyErnIcPH1a6rwtANi+lRYsWtGTJkio3/crIyJBnNrdcyOvjTWxsbDjZOuHDDz+knj17sm6HKx48eFAvH+np6QSAIiMjWWeQ14dUKqWdO3dSy5YtOdtxeMKECdSrVy9O2uKC+/fv18vHkydPCABFRUWxziCPj4ULF8oWEf3ss884mTNDRDRu3Dj64IMPOGmLC+7du1erj/v377/zurS0NAJA0dHRrDPU5uPYsWM0a9Ys8vLyIi0tLQJAvXr14mTZgbFjx5K/vz/rdrgiMTGxXj4eP35MADjZF6s2H3Z2dqSpqUnTpk0jX19fMjAwIG9vb04eBw8NDaU+ffpU9+OqH82OjY0lABQXF8c6wJtL2f/4448EgFavXv3OcYoYzCxatIhat27NWXtsqWow8+ZXRWd0c3Ojb7/9VjbpNiYmhgBwslu5vD4qiI6OpkaNGlFeXh7r2gsXLqQ2bdqwbocrqhrMyOPj1q1bBICT3crl8ZGfn0+TJk0ifX19AkAmJiacrG2zYMECcnNzY90OV1Q1mJHHx82bNwkAJ7uV1+YjKiqKli5dKvs3l4OZefPmUbt27ThpiwuqGsy8+aWpqVnJR0ZGBhER3bhxgwBwsjt2Xc5XN2/eJBcXFwLAyb5Dc+fOpfbt27NuhyuqGszI4+P69esEgJPdsWvyUTFoateuHWVlZRHR62UHbG1tydDQkB4/fsyq9pw5c6hDhw7V/TisynsHeXl5AAAjI6OqflwnKu7xp6amYtasWfD09MRnn33Gul15MDIyQk5ODq5du8ZLvdp48uRJjT8vLS0FANy+fRszZ87EjBkz4OvrCw8PDwD8+ygvL8fixYuxf/9+GBoasq6tbD7S0tJq/HlVPvz8/ODu7g6APx8GBgb4+eef8eOPP2LDhg34/PPPMWXKFFy9epVVbWXzUdujm9X56Nq1KwDF+8jOzsZ3332HHTt2sK5TFcrmIzU1tcafl5WVAajso2fPnujSpQsA/s9Xbdu2xbVr1+Ds7IwdO3Zg7ty5rGorm4+UlJQafy60j+vXrwMAAgMDZfP8nJycsHbtWgQHB+P777/H8uXL613byMhINjapkqqGOBUrm3L52Gnfvn1JT0+v2k9Pirgy891331GDBg1qHM2q0tfff//N2d+mNh9ERNOnT6d9+/ZxVnPjxo1kYmIi+N+Rq69du3Zx9reRx0cFw4YNI4lEIs9+JTWyYcMGtfKxZ88eVn+PN6nKx7hx42jFihW0a9cu2VdAQAA1aNCAdu3axXqRtW+//ZbMzMwE/zuqio/q+OSTT0hPT491zYrlPYT+O3L1xeW5vCofFTsGfPfdd5WOTUpKIgCst7VYu3Yt2dnZVffjqq/M6OvrAwBnj0pv3boVR44cwerVq+Hk5MRJm/JQUFAAExMT2YhRaFJTU9GjR49aj9PQ0IBUKoW2tjYCAwPRsWNHzJ49W3ZFgC3y+Pj555/Rvn17DBw4kJOawP98KMsnnfr4CAoKQocOHTB79mzZFQG21LV/9O7dG5GRkdDR0WFVt6CgAKampkrjIyUlBb6+vrUeV52Pik+gbKnOx7Nnz3DixIlKx+bk5KCwsBCffvopWrVqBT8/v3rXrfBx5cqVerfBJcnJyXL9PhU+dHR0EBgYiPbt22POnDmC9Q8XFxdO3mcKCgpgZmamND6SkpLQs2fPWo9700dQUBDatWuHOXPmKLx/VPz32+eTJk2aQEtLi/WVoYKCgponFVc1xImLiyMAFBsby2okRfR6sqSpqSl5enpWWqDr7f2eFHFlRtnmaNQ0Z4ZhGNLQ0CCJREK+vr60ZcsWys3NJSKi27dvEwC6c+cO6wzy+Ni9ezf9+OOP77yW7QRLZZujUdOcmZp8cDmHSd7+8SaffvopjR8/nnXt+fPnU9u2bVm3wxU1zZmpyQeXc5jq6mPWrFmczZlRtjkaNc2Z0dDQIE1NzUo+KubVcTmHqT79w8/PjxYvXsy6di1zNHinpjkzNfngcg5TbT78/f3J1dW10msqtuv45ZdfWNWePXs2dezYsbofV71oXsVCefVdfvhNpk6dilevXuGPP/6ARPK6XElJCbZt21bpuJcvXwIAXr16xbpmBWlpabws+ldfGIaBlpYWGIZB586dsWbNGqSnpyMiIgKhoaGykSyfPk6ePIkVK1agtLQU3333newxwI8++oj1Amuij3epyUdRURGWL19eaeGurKws3LhxA+vWrWNdW/TxLvKerxSBKvno2LEjVq9ejYyMDJmPinl1fPlITEzE9OnTcePGDdnxcXFxKCgowMKFC1nXFn28S239Y82aNUhNTa20LllkZCRcXV0xbtw4VrVr81HlbSYzMzNYWFggMTER/v7+9S6+e/du7NmzB87Ozti4cSOA14OVq1evyia1AsCRI0ewZcsWAMDevXvRuXNnBAQEsP4fKT4+Hp06dWLVhiKQSCSQSqVo1aoVxo0bhxEjRqBRo0bVHm9hYQFzc3MkJiaid+/e9a5bm4/r168jMDAQBQUF7yzKpKurW+uE2dqIj4/n7FInl7zpY/z48RgxYgTs7OyqPd7S0hJmZmZITExEr1696l23Nh9SqRS7du3CokWL0KlTJ/Tp0wcWFhY4fPgwJxOy4+PjObt1ySUVPlq3bi3rHzX5sLKygqmpKRITE+W6DF8d8p6vFEV8fDy8vLwUXqeu1NWHtbU1TExMkJiYyOq2W20+8vPzsXnzZqxfvx6+vr7o0qULzMzMEBkZCS0trXrXrSA+Ph4+Pj6s2+GaCh9t27ZFaGgoRowYgYYNG1Z7vI2NjcyHPLdxq0Oe/tGqVSucO3cOM2bMQLdu3aCjo4MLFy7g1KlTrNcqi4+Pr3laQHXXbHx8fGjSpEmsLgsJSVlZGZmYmND3338vdBQZDx8+JEdHR1qyZEmdL4l7eXnRRx99pKBkiqesrIwaNGhAP/zwg9BRZDx48ICcnJxo6dKldb4k3q1bN/r4448VlKwyL1++pIKCAk7bLCsrI2Nj4ypvJwrF/fv36+3D09OTpkyZoqBkiqe0tJSMjIzo559/FjqKjHv37sl81PUWhYeHB02dOlVByf7Hq1evKDExkfVjv29TWlpKhoaGrG+NcEliYiI5OzvTsmXL6uzD3d2dPvnkEwUlq5q0tDR68eIFJ21V+Pj111+rO6TqCcAA4OXlhX/++YfVSEpIbty4gezsbHh7ewsdRUbTpk2RmJhYr9d6eXlh7969HCfij+vXryMnJ0epfNjb2yMhIaFer/Xy8sL+/fs5TlQ1JiYmnLd57do15ObmKpUPBwcHVj4OHjzIcSL+uHr1KvLy8pTKR7NmzVj5OHz4MMeJ3kVHRweOjo6ct3vlyhXk5+crlY/mzZvXe0ViLy8vHD16lONENVPT1aK6cvny5Vp9VLs3k6+vLxITE2tda0BZOXXqFKysrGpdBp1PKu4x1gdfX1/Ex8ezvtUjFKdOnYK1tXWty6DzCVsfd+/eVWkfNjY2cHV1FTqKDC581LaWk7Jy6tQpNGzYkPM96tjA1sedO3eQnp7OYSL+qPDh7OwsdBQZbH3ExcUhIyODw0T8cerUKdjZ2dX4lFq1fx1vb2+Ymprizz//VEg4RfPXX39hwIABrPfoUBZ8fHxgYmIi+lASfHx80KBBA+zcuVPoKPVC3Xx0794dxsbGKu9DXejRoweMjIxU2geXy1IITY8ePWBoaKjWPqodzOjo6GD48OEICwvjPJiiiYuLw/Xr1xESEiJ0FM7Q1dXFsGHDVNJHbGwsbty4oVY+9PT0VNbH7du3cfPmTbXzMXToUJX0cevWLcTExIg+lISbN2/i9u3bauVDX19fZX3cuHEDsbGxtfuoadLN5cuXCQDrlS35ZuLEidSiRQuSSqVCR+GUixcvEsDNBod88uGHH5Kjo6Pa+bhw4QIB3GxwyCcTJkwgJycntfNx/vx5ArjZ4JBPxo0bRy4uLmrn49y5cwRws8Ehn4wdO5ZcXV3VzsfZs2cJAJ05c0boKHUiNDRUHh9VbzT5Jn5+fkq183RtpKSkkLa2dk2znlWaHj16KNVOx7VR4eO3334TOopC6N69O/Xu3VvoGHJT4eP3338XOopC8PHxUaqdp2sjOTmZtLW1afPmzUJHUQje3t5KtfN0bSQnJ5OWllalDRXVCW9v75p2nlY6Hjx4QJqamrR169baDq19MBMZGUkA6MiRI9ykUzAhISFkb29PJSUlQkdRCKdOnSIAdOzYMaGjyMWYMWPIwcFBbX2cPHmSANDx48eFjiIXo0ePpmbNmqmtjxMnThAAOnHihNBR5CI4OJiaN29OpaWlQkdRCMePHycAdPLkSaGjyMXIkSOpRYsWauvj2LFjKnW3ZcSIEeTo6CiPj9oHM0REQ4cOpRYtWlBRURH7dArkzJkzxDAMpxucKSODBw9WKR979+4VOopCCQoKIkdHR9YbPyqa06dPE8MwnG44p4wEBgaqhI/o6GhiGIb2798vdBSFMmjQIJXyceDAAaGjKJSBAweSk5OT0vuo+KAopw/5BjOpqalkZGRE06dPZ5dOgeTk5FCLFi2of//+QkdROCkpKWRoaEifffaZ0FGqJTs7m5o3b04BAQFCR1E4ycnJZGhoSDNmzBA6SrVU+BgwYIDQURROUlISGRoa0ueffy50lGp5+fIlNWvWjAYNGiR0FIWTlJREBgYGNGvWLKGjVMvLly/JwcGBAgMDhY6icB49ekQGBgY0e/ZsoaNUy4sXL8je3p6CgoLkfYl8gxkiom3btin1VY/hw4eTjY0NpaenCx2FF8LDw5X2qodUKqVhw4aRjY0NZWRkCB2HF8LCwpT2qodUKqWhQ4eSra3te+Njy5YtSvspWyqV0uDBg6lhw4b09OlToePwwubNm4lhGDp48KDQUd5BKpVSUFAQ2dnZUWZmptBxeOGPP/4ghmHo0KFDQkd5B6lUSoGBgXX1If9ghoho8uTJZGRkRFeuXKl7QgWyePFi0tTUVJn7gFwxadIkpfSxaNEi0tTUpIiICKGj8MrEiRPJyMiIrl69KnSUSixcuPC99DFhwgQyNjama9euCR2lEvPnzyctLS2VewqOLePHjydjY2O6fv260FEqMW/ePNLS0lK5p+DYMm7cOGrQoAHduHFD6CiVmDt3LmlpadX1Kbi6DWZKSkqoT58+ZGlpSXfu3KlbQgWxceNGYhhGbZ+WqYkKH1ZWVkrjY8OGDe+1D39/f7KysqK7d+8KHYeIiNavX08Mwx3TTtsAACAASURBVKjt00s1UVxcTL179yZra+s674WmKL799ltiGEZtn16qieLiYurVq5dS+Vi3bh0xDKO2Ty/VRIUPGxubOu+FpijWrl1LDMPI8/TS29RtMENElJ+fT926dSNzc3O6cOFCXV/OKcuWLSOGYWjFihWC5hCSvLw8mY+LFy8KmmXp0qXEMAytXLlS0BxCkpeXR56enmRhYUGXLl0SNMuSJUuIYRhatWqVoDmEJDc3lzw8PMjS0pIuX74saJbFixcTwzC0evVqQXMISW5uLrm7uyuFj0WLFhHDMLRmzRpBcwhJbm4ude3alSwtLQW9wi+VSmnhwoXEMAytXbu2Pk3UfTBDRFRQUEABAQFkYGBQnxEUa/Lz8ykkJIQ0NDTop59+4r2+svGmj7CwMN7rv+lDmXb9FYqCggLq37+/YD7y8vJozJgxpKmpqVS7/gpFfn4+9evXjwwNDWnbtm2818/Ly6PRo0eTpqbme3nF8m3y8/Opb9++ZGhoSNu3b+e9fl5eHo0aNYo0NTXfyyuWb5Ofn099+vQRzEdubi4FBwez9VG/wQzR6y25Z86cSQzD0Pjx4+nly5f1bapOXL16lVxdXcnCwkJl1r7hg9LSUpoxYwYxDEMTJkzgzceVK1dEH1VQWlpKn332GTEMQx9++CFlZ2fzUvfKlSvk4uJClpaWdPToUV5qqgIlJSU0ffp0YhiGJk6cyJuPy5cvk7OzM1laWqrM2lB8UFJSQtOmTSOGYWjSpEm8+bh06RI5OzuTlZWVyqwNxQclJSX06aefEsMwNHnyZMrJyeGl7qVLl8jJyYmsrKzYrg1V/8FMBfv37ycrKyuytramLVu2UHl5Odsmq+T58+c0depU0tDQoB49elBqaqpC6qg6+/btIysrK7KxsVG4jylTpog+amHv3r0yH1u3blWYj2fPntGUKVNIIpGQr68vPX78WCF1VJ09e/aQpaUl2draUlhYmEJ9fPzxxySRSMjPz0/0UQ27d++W+QgPD1eoj48++ogkEgn17NmT0tLSFFJH1dm9ezdZWFgo3EdmZiZNnjyZSx/sBzNEr58JrxhouLi40JYtWzhbkCctLY1mz55NhoaGsgGTuu2ZwTUvXryQvbG5urpy7mPWrFmijzrwto+tW7dy5uPx48f0+eefy3xs3bpV9FELWVlZsoFGy5YtKSwsjIqLizlp+00fNjY2FBYWJvqohaysLNlAQxE+Zs6cSQYGBmRjY0Ph4eGij1rIysqSDTRatWpF4eHhnPlITU2V+agYMHHkg5vBTAXx8fE0duxY0tLSIlNTU/r444/p+PHjVFBQUKd20tLSaPPmzeTv708aGhpkbW1Nq1atory8PC7jqj13796l0NBQ0tTUJDMzM8582NjYiD7qwds+pkyZQidOnKizj8ePH9Mff/xBH3zwgczH6tWrRR915M6dOxQSEvKOj8LCwjq187YPW1tbWr16NeXn5ysouXryto+pU6fW28fvv/9OvXv3lvlYs2aN6KOOxMXFyebemZub09SpU+nkyZPK6iOMISLiesvutLQ0bN++HeHh4YiJiYGOjg46deoEV1dXODs7w8rKCoaGhjAwMEBubi5ycnKQkpKChIQExMTEID4+Htra2vD390dISAgGDBgAXV1drmO+N1T4CAsLw+3bt2U+WrZsCScnpzr5CA0NRUBAgOiDBWlpadi2bRvCw8Pl9pGcnIzExETcunULCQkJ0NHRqeRDR0dH6F9LZXn8+LHMR2xsLHR0dNC5c2e4urqKPgSAKx99+vRBSEiI6IMlbHzcvHkTiYmJMh+hoaHo37+/InyEK2Qw8yZpaWmIiIjAhQsXkJCQgMTERDx//hyvXr2SHWNsbIyGDRvC1dUVLVu2hI+PD7y8vKCvr6/IaO8log/lQl4fdnZ2cHFxQcuWLdG9e3d069ZN9KEAHj9+jIiICFy8eLFaHw0aNEDDhg1FHzxQ4ePChQtITExEQkICsrKyqvXRqlUr+Pj4iD4URH18VPQPPT09RUZT/GCmOnbs2IFRo0ZBoPIib7F9+3aMGTMGUqlU6CgiAMLDwzFu3DiUlZUJHUUEQFhYGCZMmIDS0lKho4gA2Lp1KyZOnIiSkhKho4gA2LJlCyZPnozi4mKhIoRLhKqsoaEhVGmRKhB9KBeiD+VC9KFciD6UC2XwIdhgRkRERERERESEC8TBjIiIiIiIiIhKIw5mRERERERERFQacTAjIiIiIiIiotKIgxkRERERERERlUYczIiIiIiIiIioNJp8FCkqKkJ6enql7z19+hQA8PDhw0rf19DQQNOmTfmI9d5SlY/MzEwAog8hEH0oF9X5ICLRhwCIPpSLuvjQ1NREkyZNeMnFy6J5WVlZsLGxkWsBsD59+uDIkSOKjvRe8/z5c9ja2srlo2/fvjh8+DAPqd5fnj17BltbW5SXl9d6bP/+/XHw4EEeUr2/1MVHQEAADhw4wEOq95enT5/Czs5O9KEk1MXHgAEDsH//fh5S8bRonrm5OXr16gWJpPZywcHBPCR6v7GwsEDPnj1FH0qCpaWl3D5GjhzJQ6L3G0tLS/j6+oo+lARra2t0795dPF8pCdbW1vDx8VE6H7zNmRkzZkytx2hrayMwMJCHNCLy+NDR0UFQUBAPaUTk9SH2D34ICQmp9RgdHR0MGjSIhzQi8vjQ1dXFwIEDeUgjoow+eBvMBAYGQltbu9qfa2pqYuDAgTA2NuYr0ntNUFCQXD4MDQ15TPX+EhgYCC0trWp/rqmpiUGDBok+eCIoKAiamtVPKdTU1ERgYKDogycGDx4s+lAihgwZUqMPLS0tBAUFwcDAgLdMvA1mDAwMMHDgwGpP2OXl5Rg9ejRfcd57DAwMMGDAANGHkmBkZISAgADRh5Ig+lAujI2N0a9fv2rfQEUf/GJsbIy+fftW66O0tJR3H7w+mj169Ohqd53V09NDnz59+Izz3jN69OhqJwEbGBiIPnimJh+Ghobw9/fnOdH7jehDuRgzZky1k04NDQ3xwQcf8Jzo/aYmH8bGxujduzeveXgdzPTp06fK20haWloYMWIEdHV1+Yzz3tO3b98qL8tqaWlh+PDh0NHRESDV+0v//v2r9TFixAjRB8/079+/ysvkWlpaGDlyZI23aUW4p3///tDX13/n+1paWggODhZ98ExAQIBS+eB1MKOtrY1hw4a9c+m2tLQUo0aN4jOKCEQfyoa2tjaGDh36zklA9CEMOjo6og8lQldXF0OGDBHPV0qCrq4uBg8erDQ+eF8BeNSoUe/cajI3N4evry/fUURQvY8ePXoIE+g9Z9SoUSgpKan0PQsLC/j4+AiU6P2mKh+Wlpbw9vYWKNH7TVXnK0tLS3h5eQmU6P1GmXzwPpjp0aMHLC0tZf/W1tbGmDFjoKGhwXcUEQC+vr7v+AgJCRF9CISfnx8sLCxk/9bS0kJoaKjoQyD8/Pxgbm4u+7foQ1h69eoFMzMz2b+1tbUxduxY0YdAVOVj3Lhxcq1BwzW8V5RIJBgzZozs0lRJSYm40JGASCQSjB49WvShJFT4qLi1UVpaKvoQEA0NDdGHEqGhoYFRo0bJfIjnK2HR1NREcHCwUrx/CLLRZHBwsOzSVKNGjdClSxchYoj8y5s+GjdujM6dOwuc6P1m5MiRslsbTZo0QadOnQRO9H7zpo+mTZuiY8eOAid6vwkODpb5sLe3R4cOHQRO9H7z5vuHg4MD2rdvL0gOQQYznTt3hoODAwBg3LhxYBhGiBgi/9KlSxfZ5myiD+Fxd3ev5ENEWDw8PGSb5Yk+hMfDwwONGjUCIPpQBjw9PWFnZwdAWB+87JoNvN5sMiEhAffu3UNWVhbs7e3x6NEjFBQU4LvvvoONjQ2cnJzg7OwsPoLKA2/7cHBwQHJyMvLz82U+nJ2d4eTkJPrggaysLMTHx+P+/fuy/iH6EI4KH/fu3cOLFy9gb2+PlJQU5OfnY9OmTbC2thZ98Mjz589l56sXL16gWbNmePz4MfLy8rBp0ybZ+4fogx+q8pGWllbJh7OzMxwdHXnzobBds/Py8nDgwAGcPHkSJ05G4nFqEgBAS1sP+kamYCSayM95BiNTa5S+KkBh3gtIpeWQaGigQ4dO6NXTF3379oW3t7d4pYAD3vRx/GQE0lKTAfzPh0RDC/nZmTB8y4eGhibad+iIXj190a9fP3h5eYk+OCAvLw/79+9/3T9ORdbZR+9efujbt6/ogyPq6qMgNwtEUmhoaKJDx/+dr0Qf3PCmj+MnI/DkcQqAt3zkZMLQpHof/fr1Q7du3UQfHJCbm1upf7ztg5FooiD3OQxNrKr0UXG+UqCPcM4HMxcuXMCmTd9j9549KC0phXWzTrB29IZ1sy5oYNUChqZ2wL+/zOM7p9CoZU8AQHlZMXIzH+JlejzS751H5v0zyEq/h8ZN7DF+XCimTp0Ka2trLqO+F1T42LV7N8pKy0QfAnP+/HlZ/ygrLYNN886wauEF62Zd0cCqeR18nEZW+n00aeog82FlZSXkr6aSVOXDuoU3rJp1kcvHi/S7yJD1D9EHW86dO4fvv/9B5sO6eWfYyOkjJ/MBXqbHV/LR1L4Zxo8LxZQpU0Qf9eDs2bMyH+XlUtg06wzrFl519vH03mm8yHgg8zF16tRKT9FyAHeDmdOnT2PR4qU4HR0Jm2Yd0bzLSDTrEAQdA9N6t/nyyV3cu/wXHl7+E6XFeZg08UPMnz8ftra2XERWa6Kjo7F4ybLXPhw6onlXjnxc2omHV3ei9FUeJk+aiPnz58PGxobD5OpJdHQ0Fi1eijOno2DTrNPr/tExCDr6JvVu88WTO7h/6S/RRz2IiorC4iXLRB9KQlRUFBYvXoozZ6Jh27wTmncJhkOHQHY+0uJw7/JfeHRlJ8pKCmQ+xA9htRMZGYlFi5fi3NnT3Pq4tBOPrv6FspICfDR5EubNm8eVD/aDmYyMDMyc+Tl27NiORi4+cPOfCVvHblyEk1FW+gqJ57ch9tR6UEk+li//L6ZOnSquLVAFFT62b9+Gxq7dFegjHLGnNog+aiE9PV3WPxq37AG3D2YozsfJ9aDSAnz11ZeYMmWK6KMK0tPTMWPGTPz55w40btkDbf1nwqaFJ6c1ykpfIeFcGOJObQDKCrF8+X9FH9Xwpo8mLX3h5j9DQT62/uujCF9/vRwfffSR6KMKnjx5ghkzZmLnzj/RpJUf3D6YAZsWHpzWKCspQsL5sEo+Pv74Y7Zr07AbzBw4cAChY8ejnNFGh4BFcOw6gk2YWikrfYWY4+sRc2I92rVri792/il7KkrktY+Q0HGQSnTQJWg5HDoMUmi9t338/ddO2NvbK7SmKrF//36Ejh0Pkuii8+DlcGg/UKH1RB81w7uPkiLEnNiAmBPr0b5dO/z115+ijzcQ0od7167488/tsqdwRIB9+/Zh7LgJYLSN0XXI12jcWrEbd3Lso36DGalUijlz5mDNmjVw9hgF92HfQFP73Q2nFMXLJ3cR9ccElBdl4a+dO3jfnVPZqOTDczTch34DTW093uq/eHIH0X98COmrF/hr5w706tWLt9rKiFQqxezZs7F27VrBffz915/o2bMnb7WVEalUilmzZmHdunVw9hwD96Ff8+sjLe61j+KXog8A5eXlmDVrFr799lu4dAtB1yFf8ewjFlF/TARKsvH3X3/Cz8+Pt9rKSHl5OT7//HOsX78erl5j0WXIcmhq8bfpc9bj24jePAkoycY/f++s79ZGdR/MlJSUIDR0LHbt3gOvUd+iRZfh9SnMmrKSQpzd/hmSbuzDli2b39uNxkpKShASEoo9e/ehW/A6wX0k39yPLVs2v7ercpaUlGDMmBDs3bcf3UZ9ixadhwmS47WP6Ui+eQBbt27ByJEjBckhNG/68Bq1Hs07DxUkR2lxAc5tn47kWwcRFrYVI0Yo9iq2slJcXIwxY0Kwb/8BwX2c3fYpUmIOIzw8DMOHC3PeFJri4mIEjxqNgwcPw2vMBjTvOFiQHKWv8nF226dIjT2KsLCt9fERrrF06dKlchcsLUVg0GAcPX4SvT/egaZufetakDMkGlqwbxuAkqI8/Lh2Pho1avTerQRZWlqKQYFBOHbiFHp//KdS+CguysGPa+ejSZMmgq0EKRQVPo6fiHjdP9r0ESyLzEdhNn5YOx9NmzZFu3btBMsjBCUlJRg4KBDHT0ai98c70ERAHxqa2rBvNwCvCl/ihzXvt48Tp6LwwZQ/0aSNv2BZZD4KXvuwt7d/L30EDBiIyKgz6D3lTzRR8G2lmtDQ1IZ9+wF4lZ+FH9YsgIODA9q2bVuXJmLkXjSPiDB+/AREREajz//tgaW9EgwcGAZdBn8BDW09TJ78ESwsLDBokGLniSgLUqkU48aNR2TUafT5z15YNlWCgQPDoOvg/0JTSxeTJk2Gubk5Bg5U7H1wZUEqlWLs2HGIjDqNvp/ug0UTJTgxMgy6DvkSGlq6mDhxEszNzTFgwAChU/FChY/o02fR5z97lcaH+7+X8CdOnAQLCwsEBAQInYoXpFIpQkPH4vSZc+jzn32waFKnNyqFwDASuA/9Cppaevjww4mwsLBA//79hY7FC1KpFCEhoTh77iL6fLoPFo3dhI70r4+vofGvD3Nz8zr5kHv68LJly7Bz51/wm7hVOQYyb9AxYB6cu4VgxMhg3Lp1S+g4vLBs2TL89fc/6DkpTDkGMm/QMWA+HN1HYfiIkYiJiRE6Di8sXboUf/+zCz0nhSnHG+cbdBqwAI7uwRg+YiRu374tdBxeWLJkCf7ZtRt+E5XRx0I4uo/EsOEjEBsbK3QcXli8eDF27d6DnpPClWIg8yadBi5Ei64jMHTY8PfGx8KFC7Fn7z70nByuFAOZN+k8cBGadx6GYcNHIC4uTu7XyTVnJioqCj179oLH8BVw9R7PKqiiIJLi+Kah0CxOw80b12BsbCx0JIVR4cNzxEq4eI0TOk6VkLQcxzYNhXZpOm5cv6rWPiIjI9GrV294jlgFF6+xQsepEqm0DEc3BEJX+hw3b1yDkZGR0JEUhszHyNVw6RYqdJwqee1jEPToBW5cv/pe+OgWvAbOniFCx6mS98lHREQEevXuDe9R38LJY7TQcapEKi3DkfUDoY+X8voIr/XKTG5uLkYGj4ZD+wClHcgAry9ReYd8j8ysHMyc+bnQcRRGTk7Ovz4GKO1ABgAYiQa8Q77H02cvMWv2bKHjKAyZjw4DlXYgAwASiSa6j/sZTzNfqLWP7OxsjBg5Cs06BirtQAao8PELMjKzMHvOHKHjKIwKH807BintQAb418fYn5H+9DnmzJ0rdByF8fLlS4wYOQotOg1R2oEM8NpHj3G/IP3pc8ybN0++19R2wOLFi5FX8Aoew1exDqho9BvYoMuQr/Db77/hwoULQsdRCP/zsVLoKLViYGKLLoOX45dffsHFixeFjqMQFi1ahPzCYngMWyF0lFoxMGmIzkOW4+eff1ZbHwsXLkRRcbmK9I+G6Bz0JX766SdcunRJ6DgKYcGCBSgqkcJ9uAr0D1M7dBm8HD/++KNa+3hVSvBQER+dg77EDz/8gMuXL9d6fI23me7cuYM2bm7oNnItnD3HcBpUkRz7bjAs9Apw7epltdpkLC4uDm5t28IreJ1Sj6orQYSj3wXB2rAYV6+o1wkiNjYWbdu1g1fwt3DyUJGlAYhwdGMgbIxLceWyeg1oYmJi0L59B3iP2ajwBTw5gwhHNg5CwwbluHxJvT6AvfbRHj4hmwRbMqLOEOHIhoGwMyVcunhe6DSccuvWLbTv0AHdVczH4fUD0MgMtfmo+TbT8uVfwczGUXVO1P/SadBS3LxxDUeOHBE6CqcsX/4VzG2d4OSuQj4YBp0GLcP1a1dw9OhRodNwymsfznByV6E1dRgGnQKX4eqVSzh27JjQaThl+fKvYNHIFY6qcqIGZP3jyuWLOH78uNBpOOXLL5fDolErwdZaqhcMg06BX+DypQs4ceKE0Gk45csvl8OycWuV89H5Xx8nT56s+dDqrsw8evQIjo5O8A75TrV++X85+dMo2Bjk48L5s0JH4YSHDx/CyckZPiGbBFtoig0nfhwJO+NXOHf2tNBROOHBgwdwcnZG99Af0LzTEKHj1JkTP4xAI5MSnD0TLXQUTrh//z6cXVzQY+xPaNYxSOg4deb498PRxKwMZ05HCR2FE+7duwcXF1f0GP8zmnUIFDpOnTn+/TA0NZfidHSk0FE4ITExEa6urvAd/6vCt7lRBMe/Hwp7CyA6KqK6Q6q/MvP777/DoIGlSp4YAKCV3ye4eOFcnR7tUmZ+//13GJhYqa4P309w/twZ3L17V+gonPD777/DyMRGJU/UwOv+ce7sacTHxwsdhRN+++03GJnaKnx/H0XRyu8TnD0TjYSEBKGjcMJvv/0GI7OGKu3jzOkoJCYmCh2FE3777TcYmzeCfXvVXGeqle8nOB0dWaOPKgczRIStYdvg0Gk4JBK519VTKmxbeMLEqinCwsKEjsIaIkJY+HY4dBwORqKaO702dPKCiWUTtfJh32mYyvqwdfJCA3M7tfAhlUpl5ytV9dHQ2VutfISFb0ezTsPBMKx2QhaMhs4+MDZvqFY+HFTZh0t3GJs3RHh4eLXHVPmbXblyBSnJj3i/fF6Um4n0e+e4aYxh0LR9EHbs/Jub9gTk8uXLr3105sdH6at8xJ/djCv7vkDC+TCUlRSxb1SNfFy6dAmpKUlo0Ymf230lRTm4fWoTLvw9D2l3I0HSctZtMowETTsMUQsfFy9exJO0VN58vElxwQvcOvYt63YYRoKm7Qdjx5+q7+PChQt4kpbK6+3w1NjjeHB1l+wr5sQGVuctmQ816B/nz59H+pPHgkxPeJEWi7ioXxB/djMKsp/Uux2GkaBpu6Aa+0eVg5lTp07ByNQGZg1b1rt4XXiV/xyXdi/GzsXtkXTzIGftNnL1Q0rSQyQlJXHWphDw6SPn6X38tawzbp/ahNsRP+DMtunYvdwbRbmZrNu2c/FF0sP7SE5O5iCpcLz2YQvThq4Kr1Vc8BJ7v+mJrMexeJl+F0c3Dcf+NdzsMdTI1RePHtxTCx8NzO1gYuvMe+3T26YhNvJHTtqyc/XFwweJePz4MSftCYXMh40TL/Wyn97DsR9HIfKPybKvrNTbrHfibuTqiwf3EtTDh0UjmFg78lbzVX4Wzmybhiv7/oumbn3h4jUOBiYNWbXZyNUX9+/FIy0trcqfVzmYiYiIgnULL4Cnx5rzslLh2HUEykpfcdqulUMnaOvoITJStSdxRUREwdrRm5daF3ctQN//+wfDllzBqOW34ewZgtznj3Bl/5es27Zu1gVa2rqIiopiH1RAIiKiYMOTj4fX92LQ7BPoMfYH9Pt0Lzr0n4NnSdfx9CH7x9ytm3WFlrYuoqNVexLwqYrzFc/En9uK7HTu5hxZN3/tQ9XPV6d4PF8BQOyp79F/2j4Efxkj+/IJ+Y51u9bN3aGppaPy56vX/YM/H3lZKfjnvx4oLyuB/9SdMDRrxEm7tfmocjBz89YtWDTlb/8ly6btFTJq1NDUgXmjViq/P9DNW7dgyYOP5ym30LzzUJjZtQIA6BpaoGPAPDCMBJkPa1+0qDY0tHRgYddSLXzw0T+kZSVo1NIPOgamsu9VrJ+ipct+uXUNLR2Y27mqvI9bt27BomlHXmvmZD5AVuptNG7N3c7Pmlq6MG/oovI+YmJieHv/KMrNxIu0OBhbOsDA1E72paGlw7ptTS1dWNiphw++9lOUlpUg4rcJ0NE3gVfwGk7b1tTWg3lD52p9vDO7Nzs7G8+fPUUn6xasCr94cgfPU15v+iiRaMDO1RfPU26hKO8ZJBqaaNYhEBINLVY15MHQojnuxqvuEwLZ2dnIep6JBjz4MDJv/M6mY/oNrGHRpC0YjiaCG1i2UGkfL168wIusZzDho39oasPIvGnl16XFoUlrf85uORqquI+srCxkv8zipX9UnK+k5aW4emA5fEZvwLVD37D+Hd5E1X08e/YM2S+zWH84lddHXNTPyEy6hh0L2sDIvCna95sFp64jOburYGCh2j4yMzORk/0CDaz48XHlwHI8S74B79Hroamtz8WvUIma+sc7V2YePnwIADCycGBV1KxhSzBgcDrs//D4biT0jCzBSCS4d3EHGrfsyctABgCMLZvh3r0HvNRSBA8evM5uzIMPHQOzKk8C+S/T0LhVL1b1KzC2dFBpH7L+YWnPqp069w8iPLy+F5f3fYFuwatZ1X4TYwvV9iHrHzz6uHF4FVr7fgwtXUNWNavCSF36h4U9q3bk9WHr2A1uvf4Dm+buKMh+gtNh/4cj3w3hZJI8UHG+us9JW0JQ4YOv/vHw6i5IJJp4+eQODq0fhM2fNcbBdQF4nsrN1a3X56uqfbwzmMnJyQEA6Og3YF3Y0X0kWnQZjkc39iEn8yHuRP0Kvw9/e/2myRM6+ibIzc3hrR7XCO0j4/55SCSaaO03hXV94LWPHHXwoWfCui15fZSVFOLsjhk4HfZ/yE5PwK4vvfAs+Qbr+gCgrW8i+51Ukf/1D358pN87B0ZDE9bNurCuVxU6og8Z8viwc/VFl6ClCJhxCIFzTsLE2hFp8dGIObmRdX2gwkcuJ20JAZ8+CrLTUZCdDtOGrmjfdxb6T9uHoHlRyH32EIfWBaAgO511Bu0afLwzmMnPzwcAaOlw86nDY+hX0NZrgP2r/eHkORp6RpactCsvWjqGyMtT3f8ZK3xo6hhw0l5dfJC0HNcOfoMPPt4GLY7qa+kYIl8NfHD195DHh6a2PrxGrcPYtSlwH7ocpcX5OPfnTE7qa+saIj9f9X1wdUm7Jh8lRTm4E/0r2vWZwUmtqlCX/qGpo3gfb2Nm1xqBcyNhYNIQD67u5qS+lq7q+2AYhpf+kZX6+jZU07b9ZPP8Glg1R9chX6K0uAB3T//GJahZ1gAAFCVJREFUun5NPt4ZzJSXv748x9UGjToGpug0YAGKC16grLiAkzbrAiORQFrOzSVHIfifD24WO6qLj0u7F6O131SYvzWPhg2MREMtfEDCvw+GkaC178ewbxeArNTbKC8rZl2fYTT+9zupIOXl5WAYhrPF8mrycfGfBbBs2h4pMUeRdPMgkm4eRG7mQ5SXFSPp5kE8STjDuj4j0UA5R7dIhEDmQ4DzFfB6kmjTtn2Rm8nNrTqGUX0fAD8+tPWMAQC6BuaVvm/t0BnA62U/2CKpoX+88xsaGr6+IlNaUsi6MAAQSZEaexxWDp1w4e95nKxXUhdKi/NhYMj+yQ+hqPBRxrOP+LNbYN7YDU3duFnTpILSV2rio1i4/mHn0gM6BqbQ0GT/xEZpcT4MVdwHEfHioyg/C3FRP+PC33NlX+n3zqL0VT4u/D0Xt06sZ12/9JWa+OBioU3Ur3+YWDvB2IrdhPAK1KN/SFFWqngfDf79mz9PvVk5g2kjSDS0OJljVlLD+8c7gxkjo9cHlr7KY10YAGIjfkBTt37wHf8zystKcJajy+PyUlKUJ/udVJGK7CVF/PlIunUIAMkeA66Ai9WZS1+phw8h+8fLJ3fRpA03jwSXvMqDkZExJ20Jgax/8ODDf8oOBC+PrfTl6j0BuobmCF4ei77/9w/r+qVq4kPI/pF06yCatu3LSX216R88vH/oGVuhkasfMh9drfSanGcPIS0vhXXzrqzrl77Kg3E1Pt4ZzDRu3BgAkP8ilXXhl0/uIj3xHBzdR75+bK7v50i+dRj3L//1zrHFRa8nKnFx6fxN8l+komnTprUfqKTw7SMtPhoxx9dDWl6GO9G/4k70r4iN/Alnd8zAi7Q7rDPkvUhRCx95PPgoK32Fm0fX4uWT/23OWVzwAlmPb8N9yHLW9QEgPysFTZs24aQtIRDqfKUo8kQfMmrzkZP5ABf/mY+sN56UeZkej7LiQrTvw82HZrF//A95+kfXIf9Fwcs0PH1jXbL0xDMwsXGCo3sw6ww19Y93Fg9p2LAhDA2Nkf30PmxaeNa7aHriWZwO/w/s2w8EiACGgaGpHQDgzPbpKC8rhrNnCAAgNe4k7l36EwCQfOsQLJu2R5PW/tAztqp3/Qrynt2Ht4/il51XFHZ2djA0NEZO5n3YtPCodzvy+DBv3BYnfhqDspJCZCZdq/R6DS0djFrOfgfy/Gf30cO3Det2hMLOzg76BobIeXofNs3d692OPD6adxqCRzcP4OrBr2DZpB0atewJXUNz+E/dydkE5Lxn9+Hbsy0nbQlB48aNZT7YPGFUl/OVIsl/fh8te7VXeB1F0aRJE+jq6SP76X1Y/TtXoj7Ie75KvLgDsZE/wdbJC1b2HaGjb4J+0/dztvRH/rP7aPkBvwsyckmFj5yn92Hl0Kne7cjbP0xtXTBg5hFc3LUQ1s27QkNTB5kPr6Dfp3s52bS64Pl9tPSv+v8rhojo7W927uqBQt1W8ByxknVxQSHCjvnOWP7FIkybNk3oNPWmcxd3FOq3hudw1fZBJMWOuU745qtl+M9//iN0nHrTqXNXFBm4wXP4Cl7qlRTlQKKhzXqvmbep8LHi6y/wf//3f5y2zSftO3ZGqXEHeAz7WugorCCSYvtcR6z65kt88sknQsepN+06dEK5SSe4D/1K4bXKy4qR/yINmtp6MDCx5bRtkpZj+1xHrF75FaZOncpp23zStl1HSM26wH0oN1dz5aUwJwMaWrqcPBYOvPaxbU4LrF39DaZMeWepkPAqpzj79fBB5gOOdq8WkKy0WBTkZqF79+5CR2GFbw8fPFMHH49jUZj/UvRRR7T1GnA+kAGArNTbauHDT136R2oMivKz1cIHX+8fGpo6aGDVjPOBDAA8T41BUUGO6vvwFeb9XL+BDWcDGeD1djuvCnOr9VH1YMbPD1lPEjhZ5EZIniSchompOdzcuHu0WAj8/PzwPC0BhTkZQkdhxZOEaJiZW6JNG9W9zQRU+IhHYc5ToaOwosJH69athY7CCj8/Pzx7fAdFec+EjsKKtPhomFtYoVWrVkJHYYWfnx+epcapvI8nCdGwsLRGy5bcbB0iFP/z8VzoKKx4knAallY2cHWtetpIlYOZ7t27w9DQGI+u71VoOEWTfGMvBg4IgISjNUGEosLHQ3XwMTCAszWMhKJHjx4wMDBSCx+DBg1QeR++vr7Q1zfEw2sq7uOmevl4dH2f0FFYoS79w8/PD3r6+irvI+nGnhp9VPkur6urixEjhuHRVf5m8XNNTuYDZDy6jrFjQ4WOwho9PT0MHz5UpX1kZyTiadJNjA0VfSgD2ekJeJp8Sy186OvrY9iwISrvIzM5Ri18GBgYYMjQwXiowj5epscjM+W2+vgYPFil+8fL9Hg8S42t0Ue1lyzGjx+Pp8kxyHhwUSHhFE1c1E9o3MQePXr0EDoKJ4wfPx5Pk27h6cNLQkepF3eif0FT+2bw8fEROgonvPZxs9IjiKpE3L8+vL29hY7CCePHj0fGo+vvrHGhKsRF/wx7h+Zq42PC+PHIeHjtnaciVYU7UT/DoVkLeHl5CR2FEyZMGI/0h1fxLOm60FHqRVzkT2jW3BHdunWr9phqBzOenp7w7OaN28fXKSScIinKzcS9C9sxb+5slb/FVEG3bt3g4emFGBX0UZjzFPcuqpcPLy8vlfZx/9IOzJ83R218+Pj4wN2jm0r6KMh+gnsXX/tQ9VsaFXTv3h1d3T1V8v1DHX306NEDXd09EXNCBX28TMP9S3/W6qPGM9niRQuQHHsSGfcvcB5Qkdw4vBLm5mYYP3680FE4ZfGiBUiJPalyV8tuHF4BCwsLjBs3TugonLJo4XykxJ5Quatl1w99A0tLS4wdO1boKJyycME8JN8+hsxHV4SOUiduHF4JKysrtfSRFHNU5a6W3Ti8AtY2NghVg1tMb7Jg/lwk3Tqicj6uH14BG1tbhITUvM5TjYMZf39/9O3XHxf/+hzS8lJOAyqK5ym3EH9uK1at/Aa6urpCx+GUPn36wL9PX5XzkXA+HKtXrYCODvu9hJSJvn37wr9PX1zYOVOFfNxEwgX19NG/f3984N8H5/+cAam0TOg4cvHaxzasWb0S2traQsfhlICAAPT+wB/n//xMZXw8S76BhAvb1dLHgAED0KOHHy7sVJ3+8Sz5BhIv7pDLR5WL5r3JgwcP0LJVa7T5YAZnS0QrivKyYhxa2wfOTUwRHR2pNpcI3+T+/fto1boN3D6YiXZ9Zggdp0bKS4txaK0/nO3NER0VoZY+7t27h9Zt3ODm/zna+X8mdJwaKS8txsE1H6BlM0tERUUIHUchJCYmok0bN7TtOxttP5gudJwaee2jN1o1t0Zk5Cmh4yiEhIQEuLm1Rdu+c9D2A+VeuLS8tBgHVvdCG0dbREScFDqOQoiPj4ebW1u07z8Pbr0/FTpOjZSVvsLB1b3h5tQQp06dqO3wcI2lS5curekIMzMzGBsZYeumJbB19IKhWWPOwnLNxX/m4/mD8zhy5BDMzc1rf4EKYmZmBiNDQxXxMQ9Zjy7iyJFDMDMzEzqOQjA3N4ehgcFrH07eMDRrJHSkarnw91y8SLqEI0cOwdTUVOg4CsHc3Pz/27vzoKbPPI7jnxRlISigyyHBKqASvLq7al0Qj50phxzOKFJnXaOcHmhXraK4jii6dqYcauUIIAgB7AIqulNLCKjchO66SAUREw9ESAIqRwJIEI37B0fFQiWQ5Bfh9/o/eb4z73mSZ4YfT0ClaiOZGQQafSUmTVHnHgForr2FLPaPY7aHgYEBtLW1kBJ9vLeHKdEjDYl78SBanv5vTO8PAwMDaGn19DC1WgUdde6RfgCt9WXIysqEvv4HL9+rGNbTf7t374aLqysKWL4K+cEqZeD/lIrqwgScPx+H2bMV8/Pv6mrPnj1wdnFR7x6l/0J1USISEuIxa9YsosdRqr1798LJ2RmFLF+0N9cTPc6g+KXf434xC4mJ52FhYUH0OEq1b98+ODk5oSDRR2178LgXcL84Caxx0GP//v1wdHREQaIPOloERI8zKB43BbySZCSxEmBubk70OErl7+8PBwcH5Cf6oKNVSPQ4g7pfkgweNwVJrASYmZkN6zXDOsxQKBQkJ7Ewc7oxcphfQtquXjcJ1lVdR0nq1zh8+DA2bNhA9DhKR6FQkJKchBmmRr09mogeaYCnd7NRnPo1jhw5And3d6LHUbq+Hp/SDHGd6a5+PSo5KE7dh8DAQKxfv57ocZSOQqHgQkoyppsY4Hq0Gu6PSg64aftx9OhRuLm5ET2O0lEoFHx/IQXTTX6PHKY7ujqaiR5pgNoKDrhp/jh27BjWrVtH9DhK19eDZqSP60x3dHW0ED3SALUVWeCm+SMoKAhr164d9us++MzMu0QiEWyWLUfHay047LwEqt60EQ2rSLUVWchP3IrNjL8hPj5uTD6XMRSRSARrG1u8fKOtPj3usJHP2gaPLQycOxc7rnoIhULYLFuOThkV9n4X1aLHkzuZKGBtg6fnFpyLjSV6HJUSCoWwtrGF9K0O7P0ugapnTPRI/T28vDwQGxND9Dgq9UuPSb37Qw16/PwjCpK2w9vbEzHR0USPo1ICgQDWNrZ4RdGFvd9FaOsaET0SasqvoTB5x0h6DP5Dk0MxMTFBYUEepuq8BfuME1pFPPkmVbDqokTcjPOAt5fHuPviBHp6FBXm/9KjgU/oPNVFCbgZ7wkfb0/ExESPux40Gg1FhfnQ134D9hlnteiRG+8FXx/vcfdBDbzfwwmtjQ8Inede4Xnkxnthq68PoplMQmchQl8PPa1uZH3nTHyPgnjknvfG9m2+YEZFEToLEUxNTVFUmI/Jml3IPOMEceNDQuepyo9DXoIPdmzfOqIect+YNWPGDHBLimA1+1P8EGaHB/9Jl3vR0Xr96iUKU3aCm34AQUHHEBMTDQ0NDZXPoQ7e7XEt1A4P/6v6K6u7uzpQmOwHbvpBHD8ehOho5rjuUcotBn2WKbE9kvzAvRiAEyeOg8mMGjOX48lr5syZ4JYUwdKChmshXxDWoyBpB0ovHcKJE8cRFRU5rnuUcosxx9ykp8etSyqfoVvajgLWdpRe/gdOnvwnIiIixm0PMzOznh5m0/BD6Bd4dOuyymfolrYjn7UNP2UcxjffnER4ePiIenzwv5kGQ6VS4bFlM9rbJEiLOQKxqBpGFkuhqTVZ7gHkVV+di9w4BjqeVSMj4zJ8fHyUvqa66+vRJhH39Gi4DyNzFfY4x8DL5/fJHr2oVCo8PLZAIm5FWmyganvcu4ncOAY6m3i4knEZ3t7eSl9T3fXtD7G4BemxgRA38FTWo67qRm8PPq5eyRhzF3mOxK96NPJhbLEUE7UmKX3tuqrryI3fDGnzQ1y9kjHmLvIcCR0dHXh4bEFLcxPSzwVC0vgARirscTOOga6WR/j31SujuTiyQq5nZgaTk5ODHX67IBSKsMBuN+b/ZSs0tfVG85aDahZUoZwdjJqfM7F2nRsiI8Jhaqq+/1ZGlOzsbPjt/OqdHtugqa2r8HWaBVUoz/wWNXfYWOe2HhHhZ8keg8jOzsYOv10QNTT29Fi1VXk9eveH23p3RISfBY1GU/g6HzsOhwO/nV+poMddlGcGo+YOm+zxG97tsdBuD+at8lVKj6b6SpSzg/HkThbcv9yA8LPfwcTEROHrfOzYbDZ27vo7Gp89xwK7PZi/yhcTlXDob6qvRHnmt3hSwVFUjwujPswAgFQqRVhYGEJDT6H7jQxzlnnC0noj9KdZjup938reoL46D7wSFmorOJg3fyHCQoOxevXq0Y48pr3fw9LWC3OsN0LfeM6o3lcmew1BdR54JUmoreBg/oLPEBYaDEdHRwVNPjZJpVKEhoYiLOy0cnoUs1BbmU32GKaBPd729virwnssWPgHhIUGw8HBQUGTj01SqRQhISEIO3Uar98AlrZesLTeCD3j0V2xIZO9huBebs/nVW+PU2EhsLe3V9DkY1NnZydCQkJw6vSZnh7LvWD5Z0X2YKG2MgcLP/sjToWFwM7OThFjK+Yw00cikSAqKgqRUdEQCupgbPYnmM6zB42+AoYzF0Nj4oevT+9sewHRg2IIeUUQ3M1Cu/gZbJYtR8BBf6xZs2bcPVQ6GhKJBJGRkYhixii8x6GAA3B1dSV7yEEsFvfvD5GwfmAPs8XQmDCcHs8helACIa8I9XfZ6BA/xzLbFQg46E/2kJNYLO7fHyJhPaaZLwJtrp38PfjFEPKLUH83q7/HoYADcHFxIXvIYegeK2Fotkj+HpVsdEhewHb5SgQc9Cd7yKm1tbW/R4NIgGkWi3t6WK6QvwevsGd/9PY4FHAAzs7Oiuyh2MNMH5lMhvz8fKSlpYGTfQN1T2vwySca0DecgUmGs6BJnYqJv9PBBE0qurva8apTgi6JCK2ND9EheQENjQlYtHgJXF2cwGAwxvylUsrW1yM1NRWc7Buor3vy2z1eitHV1vCrHmtcnbFp0yayxyjJZDLk5eX174/h9JC2iSBufIgOSRM0NCZg8ZLP+/fHWL/kS9mG12MSJmhqkz1UQCaTITc3F2lpacjOuTmgx2TD2ZhInTKsHmtcncFgMIZ96RppcO/24OTcgKCu9oM9utr6vs97eiz5fGn//lBSD+UcZt5XU1ODsrIy8Hg88Pl8NDW3oK2tHe3tHdDT08XUKXowNjYGnU6HlZUVrK2toaur+L+bknoMt4eVlRXodDpsbGwwebLyH5Ycrx4/foyysjLw+XzweDw0t7QO2aNvf5A9lIfsoV7IHurl0aNHuH37dv/3x1A95s6dCzqdrqoeqjnMkEgkEolEIimJfJfmkUgkEolEIqkb8jBDIpFIJBLpo0YeZkgkEolEIn3U/g90/GC2clBHoQAAAABJRU5ErkJggg=="}}},{"cell_type":"markdown","source":"## Inference of LDS\n\nThe recurrence formula of the Kalman filter and smoother is given by following equations.\n(See [the Bishop's PRML book](https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/) for derivation).\n\n#### Kalman filter\n$$\np(z_n \\mid x_{1:n}) = \\mathcal{N}(z_n \\mid \\mu_n, V_n) \\\\\np(z_{n+1} \\mid x_{1:n}) = \\mathcal{N}(z_{n+1} \\mid A \\mu_n, P_n)\n$$\n\n$$\nK_n   = P_{n-1} C^T (C P_{n-1} C^T + \\Sigma)^{-1} \\\\\n\\mu_n = A \\mu_{n-1} + K_n (x_n - C A \\mu_{n-1}) \\quad (n \\neq 1) \\\\\nV_n   = (I - K_n C) P_{n-1} \\\\\nP_n   = A V_n A^T + \\Gamma  \\\\\n\\mu_1 = \\mu_0 + K_1 (x_1 - C \\mu_0)\n$$\n\n#### Kalman smoother (Rauch-Tung-Striebel smoother)\n$$\np(z_n \\mid x_{1:N}) = \\mathcal{N}(z_n \\mid \\hat{\\mu}_n, \\hat{V}_n) \n$$\n\n$$\nJ_n = V_n A^T (P_n)^{-1} \\\\\n\\hat{\\mu}_n = \\mu_n + J_n (\\hat{\\mu}_{n+1} - A \\mu_n) \\quad (n \\neq N) \\\\\n\\hat{V}_n   = V_n + J_n (\\hat{V}_{n+1} - P_n) J_n^T   \\quad (n \\neq N) \\\\\n\\hat{\\mu}_N = \\mu_N \\\\\n\\hat{V}_N   = V_N\n$$\n\nThe distribution of the latent variable obtained by the Kalman smoother is the posterior distribution with all observation (including future information).\n\nAnd, the joint distribution of the latent variables can also be expressed by Bayes' theorem as follows,\n$$\np(z_1, \\cdots, z_N \\mid x_1, \\cdots, x_N) \\propto p(z_1, \\cdots, z_N, x_1, \\cdots, x_N) \\\\\n= p(z_1) p(x_1 \\mid z_1) \\prod_{n=2}^{N} p(z_n \\mid z_{n-1}) p(x_n \\mid z_n)\n$$\n\n$$\n\\ln p(z_1, \\cdots, z_N \\mid x_1, \\cdots, x_N)\n= - \\frac{1}{2} \\left\\{\n(z_1 - \\mu_0)^T P_0^{-1} (z_1 - \\mu_0) + (x_1 - C z_1)^T \\Sigma^{-1} (x_1 - C z_1) + \n\\sum_{n=2}^{N} \\left[ (z_n - A z_{n-1})^T \\Gamma^{-1} (z_n - A z_{n-1}) + (x_n - C z_n)^T \\Sigma^{-1} (x_n - C z_n) \\right]\n\\right\\} + \\text{const.}\n$$\n\nSince the joint distribution of the latent variables is a multivariate gaussian distribution,\nthe logarithm of the probability density function is quadratic with respect to the latent variables.\nTherefore, the logarithm of the probability density function can be written as\n$$\n\\ln p(z_1, \\cdots, z_N \\mid x_1, \\cdots, x_N) = - \\frac{1}{2} (Z - Z^{*})^T Q (Z - Z^{*}) + \\text{const.}\n$$\n\nHere, $Z^{*}$ is the posterior mean of the joint distribution, and $Z^{*}$ is obtained by minimizing the following quadratic function.\n$$\nL(z_1, \\cdots, z_N) = (z_1 - \\mu_0)^T P_0^{-1} (z_1 - \\mu_0) + (x_1 - C z_1)^T \\Sigma^{-1} (x_1 - C z_1) + \n\\sum_{n=2}^{N} \\left[ (z_n - A z_{n-1})^T \\Gamma^{-1} (z_n - A z_{n-1}) + (x_n - C z_n)^T \\Sigma^{-1} (x_n - C z_n) \\right]\n$$\n\nFor this competition, by assuming $P_0^{-1} = 0$, the logarithm of the joint distribution is expressed as follows,\n$$\n\\ln p(z_1, \\cdots, z_N \\mid x_1, \\cdots, x_N, \\Delta z_2, \\cdots, \\Delta z_N) = - \\dfrac{1}{2} \\left\\{\n(x_1 - z_1)^T \\Sigma^{-1} (x_1 - z_1) + \\sum_{n=2}^{N} \\left[ (z_n - (z_{n-1} + \\Delta z_n))^T \\Gamma^{-1} (z_n - (z_{n-1} + \\Delta z_n)) + (x_n - z_n)^T \\Sigma^{-1} (x_n - z_n) \\right]\n\\right\\} + \\text{const.}\n$$\n\nThe advantage of using Kalman filter and smoother is that the computational complexity of inference is reduced from cubic to linear with respect to length of the sequence ($N$).\nHowever, by using a sparse matrix solver, it is possible to obtain linear order complexity even with the direct method.\n[The cost minimization notebook](https://www.kaggle.com/saitodevel01/indoor-post-processing-by-cost-minimization) solves the posterior mean directly.","metadata":{}},{"cell_type":"markdown","source":"## Parameters of LDS\n\n$\\Gamma$ and $\\Sigma$ are the parameters of the above LDS model.\nThe definition of these parameters is the variance of the prediction error.\n$$\n\\Gamma = \\mathbb{E}[w_n w_n^T] \\\\\n\\Sigma = \\mathbb{E}[v_n v_n^T]\n$$\nThese parameters determine which of two models, wifi and sensor, to prioritize.\n\nAnd treating these parameters as constant values doesn't lead to the best results\nbecause the variance of wifi models is greatly depending on paths and timestamps.\nI think proper error estimation of machine learning models is another key technique to capture this competition.","metadata":{}}]}