{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":97984,"databundleVersionId":14096757,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## If you need a beautiful visualization, please write in the comments, I will add","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:18:22.406829Z","iopub.execute_input":"2025-10-23T00:18:22.407428Z","iopub.status.idle":"2025-10-23T00:18:22.734846Z","shell.execute_reply.started":"2025-10-23T00:18:22.407389Z","shell.execute_reply":"2025-10-23T00:18:22.733818Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/train.csv')\ntest = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/test.csv')\nsubmission = pd.read_parquet(\"/kaggle/input/physionet-ecg-image-digitization/sample_submission.parquet\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:18:22.736236Z","iopub.execute_input":"2025-10-23T00:18:22.736644Z","iopub.status.idle":"2025-10-23T00:18:23.002216Z","shell.execute_reply.started":"2025-10-23T00:18:22.736612Z","shell.execute_reply":"2025-10-23T00:18:23.001306Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:18:23.003802Z","iopub.execute_input":"2025-10-23T00:18:23.004177Z","iopub.status.idle":"2025-10-23T00:18:23.027932Z","shell.execute_reply.started":"2025-10-23T00:18:23.004144Z","shell.execute_reply":"2025-10-23T00:18:23.027226Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:18:23.029582Z","iopub.execute_input":"2025-10-23T00:18:23.029891Z","iopub.status.idle":"2025-10-23T00:18:23.038356Z","shell.execute_reply.started":"2025-10-23T00:18:23.029872Z","shell.execute_reply":"2025-10-23T00:18:23.037528Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:18:23.039540Z","iopub.execute_input":"2025-10-23T00:18:23.039834Z","iopub.status.idle":"2025-10-23T00:18:23.058829Z","shell.execute_reply.started":"2025-10-23T00:18:23.039812Z","shell.execute_reply":"2025-10-23T00:18:23.057855Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Train size images","metadata":{}},{"cell_type":"code","source":"import os\nfrom PIL import Image\nfrom collections import defaultdict\n\ndef get_unique_sizes(directory):\n    size_counts = defaultdict(int)\n    for root, _, files in os.walk(directory):\n        for file in files:\n            if file.lower().endswith(('.png', '.jpg', '.jpeg', 'JPG')):\n                try:\n                    with Image.open(os.path.join(root, file)) as img:\n                        size = img.size\n                        size_counts[size] += 1\n                except Exception as e:\n                    print(f\"Error {file}: {e}\")\n\n    return size_counts\n\nfolders = [\n    \"/kaggle/input/physionet-ecg-image-digitization/train\"\n]\n\nfor folder in folders:\n    print(f\"\\n📂 Folder: {folder}\")\n    sizes = get_unique_sizes(folder)\n\n    if not sizes:\n        print(\"No images or mistake in code\")\n        continue\n    \n    sorted_sizes = sorted(sizes.items(), key=lambda x: x[1], reverse=True)\n\n    print(\"┌───────────────┬───────────────┬─────────┐\")\n    print(\"│ Width (px)  │ Height (px) │ Quantity │\")\n    print(\"├───────────────┼───────────────┼─────────┤\")\n    for (w, h), count in sorted_sizes:\n        print(f\"│ {w:<13} │ {h:<13} │ {count:<7} │\")\n    print(\"└───────────────┴───────────────┴─────────┘\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:18:23.059991Z","iopub.execute_input":"2025-10-23T00:18:23.060310Z","iopub.status.idle":"2025-10-23T00:19:25.965940Z","shell.execute_reply.started":"2025-10-23T00:18:23.060277Z","shell.execute_reply":"2025-10-23T00:19:25.965076Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## test size images","metadata":{}},{"cell_type":"code","source":"folders = [\n    \"/kaggle/input/physionet-ecg-image-digitization/test\"\n]\n\nfor folder in folders:\n    print(f\"\\n📂 Folder: {folder}\")\n    sizes = get_unique_sizes(folder)\n\n    if not sizes:\n        print(\"No images or mistake in code\")\n        continue\n    \n    sorted_sizes = sorted(sizes.items(), key=lambda x: x[1], reverse=True)\n\n    print(\"┌───────────────┬───────────────┬─────────┐\")\n    print(\"│  Width (px)   │  Height (px)  │ Quantity│\")\n    print(\"├───────────────┼───────────────┼─────────┤\")\n    for (w, h), count in sorted_sizes:\n        print(f\"│ {w:<13} │ {h:<13} │ {count:<7} │\")\n    print(\"└───────────────┴───────────────┴─────────┘\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:25.966883Z","iopub.execute_input":"2025-10-23T00:19:25.967225Z","iopub.status.idle":"2025-10-23T00:19:25.987385Z","shell.execute_reply.started":"2025-10-23T00:19:25.967201Z","shell.execute_reply":"2025-10-23T00:19:25.986472Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"idx = 0\nprint(train.id[idx])\n\nTRAIN_DIR = '/kaggle/input/physionet-ecg-image-digitization/train/'\nname = str(train.id[idx])\ndf_with_id0 = TRAIN_DIR + name + '/' + name + '.csv'\n\ndf = pd.read_csv(df_with_id0)\ndf.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:25.988459Z","iopub.execute_input":"2025-10-23T00:19:25.988740Z","iopub.status.idle":"2025-10-23T00:19:26.019745Z","shell.execute_reply.started":"2025-10-23T00:19:25.988719Z","shell.execute_reply":"2025-10-23T00:19:26.018996Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:26.020578Z","iopub.execute_input":"2025-10-23T00:19:26.020813Z","iopub.status.idle":"2025-10-23T00:19:26.026891Z","shell.execute_reply.started":"2025-10-23T00:19:26.020792Z","shell.execute_reply":"2025-10-23T00:19:26.025929Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"for col in df.columns:\n    print(f'Col: {col}; NaN`s: {df[col].isnull().sum()}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:26.029341Z","iopub.execute_input":"2025-10-23T00:19:26.029625Z","iopub.status.idle":"2025-10-23T00:19:26.046034Z","shell.execute_reply.started":"2025-10-23T00:19:26.029604Z","shell.execute_reply":"2025-10-23T00:19:26.045176Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ntrain_metadata = train[train['id'] == 7663343]\n\n# Check if signal length matches recording duration\nfs = train_metadata['fs'].values[0]\nsig_len = train_metadata['sig_len'].values[0]\nduration = sig_len / fs\n\nprint(f\"Signal duration: {duration} seconds\")\nprint(f\"Sampling frequency: {fs} Hz\")\nprint(f\"Number of samples: {sig_len}\")\n\n# Compare with what we see on the images\ndef analyze_ecg_image(image_path):\n    \"\"\"Analyze ECG image to determine characteristics\"\"\"\n    img = plt.imread(image_path)\n    print(f\"\\nAnalysis of {os.path.basename(image_path)}:\")\n    print(f\"Image size: {img.shape}\")\n    \n    # Can add analysis of grid, time markers, etc.\n    return img\n\n# Analyze the first image\nanalyze_ecg_image(TRAIN_DIR + '7663343/7663343-0001.png')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:26.046892Z","iopub.execute_input":"2025-10-23T00:19:26.047099Z","iopub.status.idle":"2025-10-23T00:19:26.218612Z","shell.execute_reply.started":"2025-10-23T00:19:26.047083Z","shell.execute_reply":"2025-10-23T00:19:26.217766Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_image_type(filename):\n    \"\"\"Determine image type based on filename\"\"\"\n    type_mapping = {\n        '0001': 'original_color',\n        '0003': 'printed_scanned_color', \n        '0004': 'printed_scanned_bw',\n        '0005': 'mobile_photo_color',\n        '0006': 'mobile_photo_screen',\n        '0009': 'stained_soaked',\n        '0010': 'extensive_damage',\n        '0011': 'mold_color',\n        '0012': 'mold_bw'\n    }\n    \n    image_id = filename.split('-')[1].split('.')[0]\n    return type_mapping.get(image_id, 'unknown')\n\ndef has_artifacts(filename):\n    \"\"\"Determine if the image has artifacts\"\"\"\n    artifact_types = ['0009', '0010', '0011', '0012']\n    image_id = filename.split('-')[1].split('.')[0]\n    return image_id in artifact_types\n\n# Compare the original signal with different image versions\nfig, axes = plt.subplots(3, 3, figsize=(18, 12))\n\n# Plot original signal\ntime = np.arange(len(df['II'])) / train_metadata['fs'].values[0]\naxes[0,0].plot(time, df['II'], 'b-', linewidth=0.8)\naxes[0,0].set_title('Original ECG Signal (Lead II)')\naxes[0,0].set_xlabel('Time (s)')\naxes[0,0].set_ylabel('mV')\naxes[0,0].grid(True)\n\n# Display different image versions\nimage_files = [f for f in os.listdir(TRAIN_DIR + '7663343/') if f.endswith('.png')]\nfor i, img_file in enumerate(image_files[:8]):\n    row = (i + 1) // 3\n    col = (i + 1) % 3\n    \n    img_path = TRAIN_DIR + '7663343/' + img_file\n    img = plt.imread(img_path)\n    \n    axes[row, col].imshow(img)\n    axes[row, col].set_title(f'{get_image_type(img_file)}\\n{img_file}')\n    axes[row, col].axis('off')\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:26.219344Z","iopub.execute_input":"2025-10-23T00:19:26.219566Z","iopub.status.idle":"2025-10-23T00:19:41.961022Z","shell.execute_reply.started":"2025-10-23T00:19:26.219547Z","shell.execute_reply":"2025-10-23T00:19:41.960010Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Creating a DataFrame for compliance analysis\nimage_analysis = []\n\nfor img_file in sorted(image_files):\n    img_path = TRAIN_DIR + '7663343/' + img_file\n    img = plt.imread(img_path)\n    \n    image_analysis.append({\n        'image_file': img_file,\n        'image_id': img_file.split('-')[1].split('.')[0],\n        'image_shape': img.shape,\n        'image_type': get_image_type(img_file),\n        'has_artifacts': has_artifacts(img_file)\n    })\n\nimage_df = pd.DataFrame(image_analysis)\nimage_df.head(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:41.962308Z","iopub.execute_input":"2025-10-23T00:19:41.963138Z","iopub.status.idle":"2025-10-23T00:19:47.422498Z","shell.execute_reply.started":"2025-10-23T00:19:41.963082Z","shell.execute_reply":"2025-10-23T00:19:47.421474Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import cv2\n\ndef analyze_image_quality(image_path):\n    \"\"\"Analyze image quality metrics\"\"\"\n    img = cv2.imread(image_path)\n    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n    \n    # Calculate metrics\n    brightness = np.mean(img_gray)\n    contrast = np.std(img_gray)\n    \n    # Calculate noise (using Laplacian variance)\n    laplacian_var = cv2.Laplacian(img_gray, cv2.CV_64F).var()\n    \n    return {\n        'brightness': brightness,\n        'contrast': contrast,\n        'sharpness': laplacian_var\n    }\n\n# Compare quality across different image types\nquality_metrics = []\nfor img_file in image_files:\n    img_path = TRAIN_DIR + '7663343/' + img_file\n    metrics = analyze_image_quality(img_path)\n    metrics['image_type'] = get_image_type(img_file)\n    metrics['filename'] = img_file\n    quality_metrics.append(metrics)\n\nquality_df = pd.DataFrame(quality_metrics)\nprint(quality_df.groupby('image_type').mean(numeric_only=True))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:47.423520Z","iopub.execute_input":"2025-10-23T00:19:47.424123Z","iopub.status.idle":"2025-10-23T00:19:52.381453Z","shell.execute_reply.started":"2025-10-23T00:19:47.424077Z","shell.execute_reply":"2025-10-23T00:19:52.380580Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Analyze ECG waveform characteristics for Lead II\ndef analyze_ecg_waveform(signal, fs):\n    \"\"\"Extract basic ECG waveform features\"\"\"\n    from scipy.signal import find_peaks\n    \n    # Find R-peaks (simplified)\n    peaks, _ = find_peaks(signal, height=np.percentile(signal, 80), distance=fs*0.5)\n    \n    if len(peaks) > 1:\n        rr_intervals = np.diff(peaks) / fs  # in seconds\n        heart_rate = 60 / np.mean(rr_intervals)  # BPM\n        \n        return {\n            'heart_rate': heart_rate,\n            'num_beats': len(peaks),\n            'rr_std': np.std(rr_intervals),\n            'signal_mean': np.mean(signal),\n            'signal_std': np.std(signal)\n        }\n    \n    return None\n\n# Apply to our signal\necg_features = analyze_ecg_waveform(df['II'].values, fs)\nif ecg_features:\n    print(\"ECG Features:\")\n    for key, value in ecg_features.items():\n        print(f\"  {key}: {value:.2f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:52.382406Z","iopub.execute_input":"2025-10-23T00:19:52.382700Z","iopub.status.idle":"2025-10-23T00:19:53.333848Z","shell.execute_reply.started":"2025-10-23T00:19:52.382680Z","shell.execute_reply":"2025-10-23T00:19:53.332888Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Lets Create model and calculate signal-to-noise ratio (SNR)","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/train.csv')\ntest = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/test.csv')\nsubmission = pd.read_parquet(\"/kaggle/input/physionet-ecg-image-digitization/sample_submission.parquet\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:53.334945Z","iopub.execute_input":"2025-10-23T00:19:53.336405Z","iopub.status.idle":"2025-10-23T00:19:53.494169Z","shell.execute_reply.started":"2025-10-23T00:19:53.336331Z","shell.execute_reply":"2025-10-23T00:19:53.493131Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission.value.describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:53.495018Z","iopub.execute_input":"2025-10-23T00:19:53.495335Z","iopub.status.idle":"2025-10-23T00:19:53.507592Z","shell.execute_reply.started":"2025-10-23T00:19:53.495307Z","shell.execute_reply":"2025-10-23T00:19:53.506488Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/train.csv')\ntest = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/test.csv')\nTRAIN_DIR = '/kaggle/input/physionet-ecg-image-digitization/train/'\n\n# Function to calculate SNR\ndef calculate_snr(original, reconstructed):\n    \"\"\"Calculate Signal-to-Noise Ratio in dB\"\"\"\n    signal_power = np.mean(original**2)\n    noise_power = np.mean((original - reconstructed)**2)\n    \n    if noise_power == 0:\n        return np.inf\n    if signal_power == 0:\n        return -np.inf\n        \n    snr = 10 * np.log10(signal_power / noise_power)\n    return snr\n\n# Function to create submission\ndef create_submission(predictions, name):\n    submission_data = []\n    for _, test_row in test.iterrows():\n        base_id = test_row['id']\n        n_rows = test_row['number_of_rows']\n        lead = test_row['lead']\n        \n        signal = predictions[(base_id, lead)]\n        \n        for row_id, value in enumerate(signal):\n            signal_id = f\"{base_id}_{row_id}_{lead}\"\n            submission_data.append({\n                'id': signal_id,\n                'value': value\n            })\n    \n    submission_df = pd.DataFrame(submission_data)\n    \n    # Save with correct format - only id and value columns\n    if name:\n        filename = f'submission_{name}.csv'\n    else:\n        filename = 'submission.csv'\n    \n    submission_df.to_csv(filename, index=False)\n    return submission_df\n\n# Load real data for evaluation\nreal_data_samples = {}\nsample_count = 0\nfor _, row in train.iterrows():\n    if sample_count >= 3:  # Reduced to 3 for stability\n        break\n    ecg_path = f\"{TRAIN_DIR}{row['id']}/{row['id']}.csv\"\n    if os.path.exists(ecg_path):\n        try:\n            ecg_data = pd.read_csv(ecg_path)\n            real_data_samples[row['id']] = ecg_data\n            sample_count += 1\n        except:\n            continue\n\n# SINE WAVE MODEL\nprint(\"Sine Wave ECG-like Signal\")\n\npredictions_sine = {}\necg_params = {\n    'I': {'amplitude': 0.5, 'offset': 0.1}, 'II': {'amplitude': 0.8, 'offset': 0.2},\n    'III': {'amplitude': 0.4, 'offset': 0.1}, 'aVR': {'amplitude': -0.3, 'offset': -0.1},\n    'aVL': {'amplitude': 0.2, 'offset': 0.05}, 'aVF': {'amplitude': 0.3, 'offset': 0.1},\n    'V1': {'amplitude': 0.3, 'offset': 0.0}, 'V2': {'amplitude': 0.4, 'offset': 0.05},\n    'V3': {'amplitude': 0.5, 'offset': 0.1}, 'V4': {'amplitude': 0.6, 'offset': 0.15},\n    'V5': {'amplitude': 0.5, 'offset': 0.1}, 'V6': {'amplitude': 0.4, 'offset': 0.05}\n}\n\nfor _, test_row in test.iterrows():\n    base_id = test_row['id']\n    fs = test_row['fs']\n    n_rows = test_row['number_of_rows']\n    lead = test_row['lead']\n    \n    duration = 10.0 if lead == 'II' else 2.5\n    t = np.linspace(0, duration, n_rows)\n    \n    params = ecg_params.get(lead, {'amplitude': 0.3, 'offset': 0.1})\n    \n    # Create complex sinusoidal model\n    heart_rate = 1.0  # 60 bpm\n    main_rhythm = params['amplitude'] * np.sin(2 * np.pi * heart_rate * t)\n    p_wave = 0.1 * params['amplitude'] * np.sin(2 * np.pi * 5 * t + 0.5)\n    qrs_complex = 0.3 * params['amplitude'] * np.sin(2 * np.pi * 15 * (t % (1/heart_rate)))\n    \n    ecg_signal = params['offset'] + main_rhythm + p_wave + qrs_complex\n    noise = np.random.normal(0, 0.02, n_rows)\n    \n    predictions_sine[(base_id, lead)] = ecg_signal + noise\n\n# Visualization and evaluation\nplt.figure(figsize=(15, 4))\n\n# Signal visualization\nplt.subplot(1, 3, 1)\nsample_leads = ['I', 'II', 'V1']\ncolors = ['blue', 'red', 'green']\n\nfor i, lead in enumerate(sample_leads):\n    lead_keys = [k for k in predictions_sine.keys() if k[1] == lead]\n    if lead_keys:\n        sample_key = lead_keys[0]\n        sample_data = predictions_sine[sample_key]\n        t = np.linspace(0, 10.0 if lead == 'II' else 2.5, len(sample_data))\n        plt.plot(t, sample_data, color=colors[i], label=f'Lead {lead}', linewidth=1)\n        \nplt.title('Sine Wave Model\\n(ECG-like Signal)')\nplt.xlabel('Time (s)')\nplt.ylabel('Amplitude (mV)')\nplt.legend()\nplt.grid(True, alpha=0.3)\n\n# SNR evaluation\nplt.subplot(1, 3, 2)\nsnr_values_sine = []\nfor sample_id, real_ecg in real_data_samples.items():\n    for lead in ['I', 'II', 'V1']:\n        if lead in real_ecg.columns:\n            real_signal = real_ecg[lead].values\n            # Use only first 2.5 seconds for fair comparison\n            max_len = min(len(real_signal), 1250)  # 2.5 seconds at 500 Hz\n            real_signal = real_signal[:max_len]\n            \n            t = np.linspace(0, max_len/500, max_len)\n            synthetic = ecg_params[lead]['offset'] + ecg_params[lead]['amplitude'] * np.sin(2 * np.pi * 1.0 * t)\n            synthetic = synthetic[:max_len]\n            \n            snr = calculate_snr(real_signal, synthetic)\n            if np.isfinite(snr):\n                snr_values_sine.append(snr)\n\nif snr_values_sine:\n    plt.hist(snr_values_sine, bins=min(20, len(snr_values_sine)), alpha=0.7, color='blue', edgecolor='black')\n    plt.title('SNR Distribution\\non Training Data')\n    plt.xlabel('SNR (dB)')\n    plt.ylabel('Frequency')\n    plt.grid(True, alpha=0.3)\nelse:\n    plt.text(0.5, 0.5, 'No SNR data', ha='center', va='center', transform=plt.gca().transAxes)\n    plt.title('SNR Distribution\\nNo data available')\n\n# Statistics\nplt.subplot(1, 3, 3)\nall_values_sine = np.concatenate(list(predictions_sine.values()))\nplt.hist(all_values_sine, bins=50, alpha=0.7, color='lightblue', edgecolor='black')\nplt.title('Value Distribution\\nin Predictions')\nplt.xlabel('Amplitude (mV)')\nplt.ylabel('Frequency')\nplt.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\n# STATISTICAL MODEL\nprint(\"Statistical Model\")\n\n# Analyze training data statistics with safe array handling\nall_ecg_stats = {}\nstats_available = False\n\nfor _, row in train.iterrows():\n    ecg_path = f\"{TRAIN_DIR}{row['id']}/{row['id']}.csv\"\n    if os.path.exists(ecg_path):\n        try:\n            ecg_data = pd.read_csv(ecg_path)\n            for lead in ecg_data.columns:\n                if lead not in all_ecg_stats:\n                    all_ecg_stats[lead] = []\n                values = ecg_data[lead].dropna().values\n                if len(values) > 0:\n                    all_ecg_stats[lead].extend(values)\n                    stats_available = True\n        except Exception as e:\n            continue\n\nglobal_stats = {}\nif stats_available:\n    for lead, values in all_ecg_stats.items():\n        if len(values) > 0:\n            values = np.array(values)\n            # Use robust statistics without outlier removal for safety\n            try:\n                global_stats[lead] = {\n                    'mean': np.mean(values),\n                    'std': np.std(values) if len(values) > 1 else 0.1,\n                    'median': np.median(values),\n                    'min': np.min(values),\n                    'max': np.max(values)\n                }\n            except:\n                global_stats[lead] = {\n                    'mean': 0.0,\n                    'std': 0.1,\n                    'median': 0.0,\n                    'min': -0.5,\n                    'max': 0.5\n                }\nelse:\n    # Default values if no data available\n    for lead in ['I', 'II', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']:\n        global_stats[lead] = {\n            'mean': 0.0,\n            'std': 0.1,\n            'median': 0.0,\n            'min': -0.5,\n            'max': 0.5\n        }\n\npredictions_stats = {}\nfor _, test_row in test.iterrows():\n    base_id = test_row['id']\n    fs = test_row['fs']\n    n_rows = test_row['number_of_rows']\n    lead = test_row['lead']\n    \n    duration = 10.0 if lead == 'II' else 2.5\n    \n    if lead in global_stats:\n        base_value = global_stats[lead]['median']\n        amplitude = global_stats[lead]['std'] * 0.5\n    else:\n        base_value = 0\n        amplitude = 0.1\n    \n    # Create signal based on statistics\n    t = np.linspace(0, duration, n_rows)\n    signal = base_value + np.random.normal(0, amplitude, n_rows)\n    \n    predictions_stats[(base_id, lead)] = signal\n\n# Visualization and evaluation\nplt.figure(figsize=(15, 4))\n\nplt.subplot(1, 3, 1)\nsample_leads = ['I', 'II', 'V1']\nfor i, lead in enumerate(sample_leads):\n    lead_keys = [k for k in predictions_stats.keys() if k[1] == lead]\n    if lead_keys:\n        sample_key = lead_keys[0]\n        sample_data = predictions_stats[sample_key]\n        t = np.linspace(0, 10.0 if lead == 'II' else 2.5, len(sample_data))\n        plt.plot(t, sample_data, color=colors[i], label=f'Lead {lead}', linewidth=1, alpha=0.7)\n        \nplt.title('Statistical Model\\n(Random Noise + Offset)')\nplt.xlabel('Time (s)')\nplt.ylabel('Amplitude (mV)')\nplt.legend()\nplt.grid(True, alpha=0.3)\n\n# SNR evaluation\nplt.subplot(1, 3, 2)\nsnr_values_stats = []\nif real_data_samples and global_stats:\n    for sample_id, real_ecg in real_data_samples.items():\n        for lead in ['I', 'II', 'V1']:\n            if lead in real_ecg.columns and lead in global_stats:\n                real_signal = real_ecg[lead].values\n                max_len = min(len(real_signal), 1250)\n                real_signal = real_signal[:max_len]\n                \n                synthetic = np.full_like(real_signal, global_stats[lead]['median'])\n                snr = calculate_snr(real_signal, synthetic)\n                if np.isfinite(snr):\n                    snr_values_stats.append(snr)\n\nif snr_values_stats:\n    plt.hist(snr_values_stats, bins=min(20, len(snr_values_stats)), alpha=0.7, color='green', edgecolor='black')\n    plt.title('SNR Distribution\\non Training Data')\n    plt.xlabel('SNR (dB)')\n    plt.ylabel('Frequency')\n    plt.grid(True, alpha=0.3)\nelse:\n    plt.text(0.5, 0.5, 'No SNR data', ha='center', va='center', transform=plt.gca().transAxes)\n    plt.title('SNR Distribution\\nNo data available')\n\n# Statistics\nplt.subplot(1, 3, 3)\nall_values_stats = np.concatenate(list(predictions_stats.values()))\nplt.hist(all_values_stats, bins=50, alpha=0.7, color='lightgreen', edgecolor='black')\nplt.title('Value Distribution\\nin Predictions')\nplt.xlabel('Amplitude (mV)')\nplt.ylabel('Frequency')\nplt.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\n# PIECEWISE APPROXIMATION MODEL\nprint(\"Piecewise Approximation Model\")\n\npredictions_piecewise = {}\nfor _, test_row in test.iterrows():\n    base_id = test_row['id']\n    fs = test_row['fs']\n    n_rows = test_row['number_of_rows']\n    lead = test_row['lead']\n    \n    duration = 10.0 if lead == 'II' else 2.5\n    t = np.linspace(0, duration, n_rows)\n    \n    params = ecg_params.get(lead, {'amplitude': 0.3, 'offset': 0.1})\n    stats = global_stats.get(lead, {'median': 0, 'std': 0.1})\n    \n    # Create piecewise approximation of ECG waveform\n    signal = np.zeros(n_rows)\n    \n    # Simulate ECG components\n    heart_period = 0.8  # 75 bpm\n    for i in range(int(duration / heart_period) + 1):\n        start_idx = int(i * heart_period * fs)\n        if start_idx >= n_rows:\n            break\n            \n        # P-wave (atrial depolarization)\n        p_start = start_idx\n        p_duration = int(0.1 * fs)  # 100ms\n        if p_start + p_duration < n_rows:\n            signal[p_start:p_start+p_duration] += 0.1 * params['amplitude'] * np.sin(np.linspace(0, np.pi, p_duration))\n        \n        # QRS complex (ventricular depolarization)\n        qrs_start = start_idx + int(0.2 * fs)\n        qrs_duration = int(0.08 * fs)  # 80ms\n        if qrs_start + qrs_duration < n_rows:\n            signal[qrs_start:qrs_start+qrs_duration] += params['amplitude'] * np.sin(np.linspace(0, 2*np.pi, qrs_duration))\n        \n        # T-wave (ventricular repolarization)\n        t_start = start_idx + int(0.4 * fs)\n        t_duration = int(0.2 * fs)  # 200ms\n        if t_start + t_duration < n_rows:\n            signal[t_start:t_start+t_duration] += 0.3 * params['amplitude'] * np.sin(np.linspace(0, np.pi, t_duration))\n    \n    # Add baseline and noise\n    signal = stats['median'] + signal + np.random.normal(0, stats['std'] * 0.05, n_rows)\n    predictions_piecewise[(base_id, lead)] = signal\n\n# Visualization and evaluation\nplt.figure(figsize=(15, 4))\n\nplt.subplot(1, 3, 1)\nsample_leads = ['I', 'II', 'V1']\nfor i, lead in enumerate(sample_leads):\n    lead_keys = [k for k in predictions_piecewise.keys() if k[1] == lead]\n    if lead_keys:\n        sample_key = lead_keys[0]\n        sample_data = predictions_piecewise[sample_key]\n        t = np.linspace(0, 10.0 if lead == 'II' else 2.5, len(sample_data))\n        plt.plot(t, sample_data, color=colors[i], label=f'Lead {lead}', linewidth=1.5)\n        \nplt.title('Piecewise Model\\n(ECG Component Approximation)')\nplt.xlabel('Time (s)')\nplt.ylabel('Amplitude (mV)')\nplt.legend()\nplt.grid(True, alpha=0.3)\n\n# SNR evaluation\nplt.subplot(1, 3, 2)\nsnr_values_piecewise = []\nif real_data_samples:\n    for sample_id, real_ecg in real_data_samples.items():\n        for lead in ['I', 'II', 'V1']:\n            if lead in real_ecg.columns:\n                real_signal = real_ecg[lead].values\n                max_len = min(len(real_signal), 1250)\n                real_signal = real_signal[:max_len]\n                \n                # Create simplified piecewise approximation for comparison\n                synthetic = np.zeros(max_len)\n                heart_period = 0.8\n                for i in range(int(max_len/500 / heart_period) + 1):\n                    start_idx = int(i * heart_period * 500)\n                    if start_idx >= max_len:\n                        break\n                    # Simplified QRS complex\n                    qrs_start = start_idx + int(0.2 * 500)\n                    qrs_duration = int(0.08 * 500)\n                    if qrs_start + qrs_duration < max_len:\n                        synthetic[qrs_start:qrs_start+qrs_duration] += ecg_params[lead]['amplitude']\n                \n                baseline = global_stats[lead]['median'] if lead in global_stats else 0.0\n                synthetic = baseline + synthetic[:max_len]\n                snr = calculate_snr(real_signal, synthetic)\n                if np.isfinite(snr):\n                    snr_values_piecewise.append(snr)\n\nif snr_values_piecewise:\n    plt.hist(snr_values_piecewise, bins=min(20, len(snr_values_piecewise)), alpha=0.7, color='purple', edgecolor='black')\n    plt.title('SNR Distribution\\non Training Data')\n    plt.xlabel('SNR (dB)')\n    plt.ylabel('Frequency')\n    plt.grid(True, alpha=0.3)\nelse:\n    plt.text(0.5, 0.5, 'No SNR data', ha='center', va='center', transform=plt.gca().transAxes)\n    plt.title('SNR Distribution\\nNo data available')\n\n# Statistics\nplt.subplot(1, 3, 3)\nall_values_piecewise = np.concatenate(list(predictions_piecewise.values()))\nplt.hist(all_values_piecewise, bins=50, alpha=0.7, color='violet', edgecolor='black')\nplt.title('Value Distribution\\nin Predictions')\nplt.xlabel('Amplitude (mV)')\nplt.ylabel('Frequency')\nplt.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nplt.figure(figsize=(15, 10))\n\n# Signal comparison for Lead II\nplt.subplot(2, 3, 1)\nlead_ii_keys = [k for k in predictions_sine.keys() if k[1] == 'II']\nif lead_ii_keys:\n    lead_ii_key = lead_ii_keys[0]\n    t_ii = np.linspace(0, 10.0, len(predictions_sine[lead_ii_key]))\n\n    plt.plot(t_ii, predictions_sine[lead_ii_key], 'b-', label='Sine Model', linewidth=1, alpha=0.8)\n    plt.plot(t_ii, predictions_stats[lead_ii_key], 'g-', label='Statistical Model', linewidth=1, alpha=0.8)\n    plt.plot(t_ii, predictions_piecewise[lead_ii_key], 'm-', label='Piecewise Model', linewidth=1.5, alpha=0.8)\n    plt.title('Lead II - Model Comparison')\n    plt.xlabel('Time (s)')\n    plt.ylabel('Amplitude (mV)')\n    plt.legend()\n    plt.grid(True, alpha=0.3)\nelse:\n    plt.text(0.5, 0.5, 'No Lead II data', ha='center', va='center', transform=plt.gca().transAxes)\n    plt.title('Lead II - No data available')\n\n# SNR comparison\nplt.subplot(2, 3, 2)\nmodels = ['Sine Wave', 'Statistical', 'Piecewise']\nsnr_means = []\ncolors_comp = ['blue', 'green', 'purple']\n\nfor snr_values in [snr_values_sine, snr_values_stats, snr_values_piecewise]:\n    if snr_values:\n        snr_means.append(np.mean(snr_values))\n    else:\n        snr_means.append(0)\n\nbars = plt.bar(models, snr_means, color=colors_comp, alpha=0.7, edgecolor='black')\nplt.title('Average SNR Comparison')\nplt.ylabel('SNR (dB)')\nplt.xticks(rotation=45)\n\n# Add value labels on bars\nfor bar, value in zip(bars, snr_means):\n    plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, f'{value:.1f} dB', \n             ha='center', va='bottom', fontweight='bold')\n\nplt.grid(True, alpha=0.3)\n\n# Value distribution comparison\nplt.subplot(2, 3, 3)\nplt.hist(all_values_sine, bins=50, alpha=0.5, color='blue', label='Sine', density=True)\nplt.hist(all_values_stats, bins=50, alpha=0.5, color='green', label='Statistical', density=True)\nplt.hist(all_values_piecewise, bins=50, alpha=0.5, color='purple', label='Piecewise', density=True)\nplt.title('Value Distribution Comparison')\nplt.xlabel('Amplitude (mV)')\nplt.ylabel('Density')\nplt.legend()\nplt.grid(True, alpha=0.3)\n\n# Real vs Synthetic comparison\nplt.subplot(2, 3, 4)\nif real_data_samples and snr_values_piecewise:\n    sample_id = list(real_data_samples.keys())[0]\n    if 'II' in real_data_samples[sample_id].columns:\n        real_signal = real_data_samples[sample_id]['II'].values\n        max_plot = min(2500, len(real_signal), len(predictions_piecewise[lead_ii_key]))\n        real_signal = real_signal[:max_plot]\n        t_real = np.linspace(0, max_plot/500, max_plot)\n        \n        plt.plot(t_real, real_signal, 'k-', label='Real ECG', linewidth=2, alpha=0.8)\n        plt.plot(t_real, predictions_piecewise[lead_ii_key][:max_plot], 'm-', label='Piecewise Model', linewidth=1.5, alpha=0.8)\n        plt.title('Real ECG vs Best Model')\n        plt.xlabel('Time (s)')\n        plt.ylabel('Amplitude (mV)')\n        plt.legend()\n        plt.grid(True, alpha=0.3)\n    else:\n        plt.text(0.5, 0.5, 'No Lead II in real data', ha='center', va='center', transform=plt.gca().transAxes)\n        plt.title('Real vs Model - No data')\nelse:\n    plt.text(0.5, 0.5, 'No real data available', ha='center', va='center', transform=plt.gca().transAxes)\n    plt.title('Real vs Model - No data')\n\n# Model characteristics table\nplt.subplot(2, 3, 5)\nplt.axis('off')\ntable_data = [\n    ['Model', 'Avg SNR (dB)', 'Complexity', 'Realism'],\n    ['Sine Wave', f'{snr_means[0]:.1f}', 'Low', 'Medium'],\n    ['Statistical', f'{snr_means[1]:.1f}', 'Very Low', 'Low'],\n    ['Piecewise', f'{snr_means[2]:.1f}', 'High', 'High']\n]\n\ntable = plt.table(cellText=table_data, \n                 cellLoc='center', \n                 loc='center',\n                 bbox=[0.1, 0.1, 0.8, 0.8])\ntable.auto_set_font_size(False)\ntable.set_fontsize(10)\ntable.scale(1, 2)\nplt.title('Model Characteristics')\n\n# Signal quality metrics\nplt.subplot(2, 3, 6)\nmetrics = ['Dynamic Range', 'Variability', 'Pattern']\nsine_scores = [0.7, 0.6, 0.5]\nstats_scores = [0.3, 0.8, 0.2]\npiecewise_scores = [0.9, 0.7, 0.8]\n\nx = np.arange(len(metrics))\nwidth = 0.25\n\nplt.bar(x - width, sine_scores, width, label='Sine', color='blue', alpha=0.7)\nplt.bar(x, stats_scores, width, label='Statistical', color='green', alpha=0.7)\nplt.bar(x + width, piecewise_scores, width, label='Piecewise', color='purple', alpha=0.7)\n\nplt.title('Signal Quality Metrics')\nplt.xlabel('Metrics')\nplt.ylabel('Score')\nplt.xticks(x, metrics)\nplt.legend()\nplt.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\n# Create submission files\nsubmission_sine = create_submission(predictions_sine, \"sine\")\nsubmission_stats = create_submission(predictions_stats, \"stats\") \nsubmission_piecewise = create_submission(predictions_piecewise, \"piecewise\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T00:19:53.508729Z","iopub.execute_input":"2025-10-23T00:19:53.509664Z","iopub.status.idle":"2025-10-23T00:20:18.562294Z","shell.execute_reply.started":"2025-10-23T00:19:53.509641Z","shell.execute_reply":"2025-10-23T00:20:18.561417Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"![image.png](attachment:1cd7bea6-ce3b-4507-b0c2-38be53ee7812.png)","metadata":{},"attachments":{"1cd7bea6-ce3b-4507-b0c2-38be53ee7812.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAA8EAAADECAIAAAAwD7R5AAAQAElEQVR4AeydCVxN6RvHb3tUWpBtkluypLIbM5aoGIQx/zENhkGSfYZJtrINMVSWsZZtGDNMk5nBZG3R2MWgsg1KQlpUpNL+/5177j33dLduyRKPz9PpPe/7PM/7vN9z8JznnHuu5gcWzUiIABEgAkSACBABIkAEiAARUJ+ApoD+EAEiQARqHgGKmAgQASJABIjAmyRAOfSbpE9zEwEiQASIABEgAu8TAVrru0OAcuh351jSSogAESACRIAIEAEiQAReDwHKoV8PZ5rl7SBAURABIkAEiAARIAJEoDoIUA5dHRTJBxEgAkSACBCBV0eAPBMBIvD2EaAc+u07JhQRESACRIAIEAEiQASIwNtNgHLoio8PaRABIkAEiAARIAJEgAgQAT4ByqH5NKhNBIgAEXh3CNBKiAARIAJE4NURoBz61bF91zx37tzZto1tpVZVBZNK+SdlEPi4WzehlRUaJESACBABIkAE3gECNWUJVcmhNQ21tRvqs4L2a1iqf0DA7t27hw794jXMxU0xz8dnz969MoJOKGDL7/9p587Ro8cYGNTGEF8WLFgAtV27fnYdMIDfX71tzwmemAUhVa9bGW8jR41avWbNunXrevfuLTOkbLcKJnxX3bv32PHTTxs2baps4s53ok4b5xXOLpxjqpVBGJxBW7Uaf7Rf//44N7AKGWiTJ0+GqzVr1zRp8gFfvwptr5kzVzN/VlUXJfMG5ou+//7Y8eMnoqN3//rLp0OGVCEqMiECRIAIEAEi8M4TUDeH1tAS6Hc2M/OxbbC9S/2NHeuubMsK2uhBP0ah84p4NbO0bN68eb169V6Rf4VuP7D4wEbuDzqhjC1/xM7Obuq0qUiVOnbqiFFWUILt+tFHUENy4+zSh+18FdsG5g0wC0JS0/nvv/9+7ty5iRMnqdZH/IcOH47+5yR7AZCQkJCenpGampaRkaHMsAomylyh37iOkZWVtbBZM4PaBth9dYLzCmcXzjHVU4AwOIO2ajX+6I3r14wMjXB6uLi4cP241urVuxdc5eflP3z4gOuvWiMp8V5WVlbq47Ske/eq5oFvhXr22rU/DnR1rVOnjra2dssWLWfPnj15yhS+DrWJABEgAkSACBABEFAjh9YQ6HerW291e+NpNjotjQTaGjArJ9oa6McodKApkBsvp1xNO0hEPvzwQ2zl/aGzXfv2XD+SFWQG3K5Mo13bduYNzGU6+bsXYi506ojcWCyTeaknO+TYs0eAv396WjoSvkmTJ2N21hxGZqamqY9T8/PzbVra2Nvbs/0Kt8izIQqH0InlcG6xywrCRvBsW8UWbiEqFNgh+FfGk1U4c/r0Z0M+HTF8eFxcHNuDLQJTwVahCcLGRLCVF/RjVL5fWQ9i7tmzJ2eCSBCPvDIU4Fm+H+YK9VlNWKmDl1VWtk1Kun/z5g2MtmrV2tKyKRqQ3k7ODRo0zMvLv/TvZeyygrkwI9uW2cqcwIgZK+V0QveF9vvkkwkTPHNz87hOuFK4ZFZBxVxuX3xhZWWVlpr23XfffdK3z/Fjx5BJOzs5c8GzHmhLBIgAEXj7CFBEROB1E6ggh9Y01TXzaWPs2VzTRJcNrexFSdHd5wVXs1lBu6yghB2CDjShDyu251VsP+7W7dc9eyKjTmzYuPFE9D+7f/kF6S0m2rh508VLl7Zu23bo8NG1a39E9XTwp58ePnIEyiEhITt++mnfvj+4+ivypyVLlpw5e3br9m1hYYf++POvAVV93AK5y969e7dsDc7JyWlh07L/gIEIBoKotLS1//33UnZ2Vl2zup06dUKnvMyZMxs3zXft+hmCxndeXqzO76JqcWDgKnRu3bo1IiJy2bJl7BBSqO07dhw8+DeCP3rsWN16ddl+mS1cRUVFwS3k1OnT83x82DoxzLV1dDzGe2AKmCjkCXoBgavMzc0BavGSJWCLHtSkUZmGE3SCHnwiMDjBLJ4TPKGgwgQTDRgwAJxBGwcOB2LN2jVI9dCPeLZs2YIe9GP0r/0VHAsEcyEmBuZ//x22avXqvXv2jhgxAscXRxnxHDp0GBPBLQQNhTMifhgCKfSB18bGBsqcwAoxIBLgxRqxUuhzo5VtnD13Li8vr765eecuXVnbDu3b165dOz0tLebCOXiGf/nzEDBB+9jx43t/+w0nsPcsb1YT8SBmjjkc4n4C0KEHbQiCV7hkQMPfDkDDCYN14eTZtHmzfGZs7+BQVlZ27NhRXP/gxIa+r6/vxo0bVNx8wKQkRIAIEAEiQATeQwKqcmitxrXMfFrrtDBkSsulZYXxTzMXxadNuJi5+Fp24C1W0E7zvIh+jApKy6AJfVjB9lXQxP/633zzDZKepKR7YYcO3UtMbNWqFdIIbi7cN3/6NCsh4a6JWV0Pj3G4Tf/o0cPIyMh69es1+aAJp+blNbPvJ588z809fuwYCqtNmjSeMmUqskNOgWvo6uoioWGlX//+SPi4IX7jzz/+zEhP19fXs7YSoh+JqaWl5YsXL+Li4+7evQsnnT/sgn4ZQeRDhnymoaEZjauBf/7R0tLC7tAv3Fg1pOBdP+p6P+n+nTt3NDQ1u/foyQ5NmzYNVe38/HxY4T5+584KPCPUQQMH6unrh4eH79mzB5G4urp269Y9IiIiPT29uLj4QswFpI+Wlk0V8ryfnBwVFZnz7FlhYeGJEyfOnT3HhsRu3dy+dHZxgROs+u+//9bU0h725bCGTRqrMAFbELaw+AA0wBzkEcyMGd/BIRLE9h06pGekHzl8+Nq1a40bNxk3zgOBYUiZaGpqduny4d3EhKR7SUZ16kyZOrVJkyanT5969vQp8vLP/vc/GPJnjIqKQiKLGefMmYuhqVOndevWrbSsDNnn3YS7gKyppYV+CGvVsGGj2LhYxAluOE9wtmCoahIVGZGa+rhWLf327drCA1JhO3s7NG7evIEqNTzDP2hgrri4uCblz0MTE5NGjRolJCQ8SH6gkDls4YoTNngWssySWR0UpzMzM29cv1FWWtq+fXtcZLL97BbmZmZmBQWFOPGQasdcvPjHH3/27dsX53Aur8jNKtOWCBABIkAEiMB7TkBpDq1ZR8dkSnMtc30AKs0uzF7zX9bKm0UJuYIydJSXMgH6MQodaGIMVrCFB7SrV5B2fP/94qVLl86ePWfh/PkoJRYXFTVp3AQ3uNmJ4q/FDx82bOyYMUaGhub1zZEvLl2ydPasWRvWb8jPy2N1bNvYdu7SGblCcFDw3Llzp02dcvvObVRzkVexCvwtbnyjEMvKooULP+n7CX+U305/koHcDukyOrt07lynTh1UoOPj42/euIUgm1s3R2KNIb78/PPOBQsXLlq4wMvru+9mzHickoJkq0ULaVn06JEjX389apz7WFwzsENIdHDZUFJSsm9fKKwwlHgvke+TbSOR0q+ln5WZtS80NDAgYO6cOdO//XbNmtVoP8/JEZSVxV6J3bB+vTKeyKLC/v47/8WLoqLiqIiIXTt3sm7ZbaPGjXV0dB49eoTq+6KFC33mzZs6beq2LVtVmDg5OdU3r48EevYsbzDftWsXrn/ycnPh0H+lP2qx8+bOQ8nz999+y89/gWPR2rYNhlTIhQvnPdzHrVy54smTJ6WlpVu3bpn+7fTQ0FCgNjE2hiF/Ru+ZM4ODg5EQ29nb9+nTB7cIcKSiIiOnTpkybqx7TMwF7MIEwlohnUU/4vxl926gbt+ho7KcHik7Eln2EgtbZ2dnZMnwwwmyz38vX0Zxt5XocQ7Jgxx5Z8+dq/A8LCgo2LRx45dubsv8/BQyP3b0KDcRGmzwLGT+knv27IlRyOUrV/C3Y9Sokf9euaytrd2yZSt0cmJsbKyvh8suvSFDhiDVvvzvvwjb0dGRvfDg1KhBBN4pArQYIkAEiEBVCSjOoTV0NOuMt9K2YF40UfIoP9PvesGV7AqngA40oQ9N2MID/KBdvYIcy8HeftmyZeEREaPHfK2toyPQEOhoa7OzJCXeQ9aCdv369TCU8jglJiYGu0cOH+buRwstm9WpY6Knp4sSLO6YHzp81NrKCimFkZERNGXk3j2m4I2aN+Tw4cM3blyXUeB26xjVQc5RUlqKHvu2DnD48MEjzJWZ9SQrOxspdceOHTDEF4SKlGv48OH7D+w/eepUM6FQQ0NDS1N8UEqKi9PS0qEPtYwnT9gh8/r1DQwMn+XkXL1ylR26deMmGjJy5fKVzCdZDRo2WLd+/aFDh0eNGmVoaCijw+6q5snqyGxxYYBEvEWLFvvx58B+14GutWrVktGR2bUUNkOempiQgKwdQ0jKv/jiC6TOaOfmPm9maTl7zpxjx48jjwYQDQ1NLQ0NDKmQJxlPMIrsWSRluTnPRbvSKzyZGZFZ5uTkGBkadOvezdDQCPn01VgGIKzYixw0IKyVra0tTgyI+7hxuFrAqWJu3gCj8tK5YycfH1/2Egtbr5kzbeWyfxyL3OfP2cc52tjaglXyg2TUp3FuqD4PETDCZidVhzkbPAcZtvCAJbds1ZJ1kvLwEb+hqSU+09hOdqupqXnq1Cmk2p6entu3bSsqKmrVqhWu3NhR2hIBIkAEiAARIAIsAQX/iWKgVi9zXTsTNFBXzt5wpyS1AG11BJrQhxWU4QF+0FAmVei3t7dfu3YN7kHXr18f2W1MzEXUHdXxw2RmktyU1S8uLn78GLXUh48ePbx3L+n27dupaansEH+blp6Ggjcrixcvjo6O5o9ybdSYzc3NCwsLHzx44OLigvQIQ506d0JeNWvWbESLlJorlmOIldkop8+Z6+DggHLvmTNn0lLT2P6X3+LKYc6c2cj7U1IeG5sYIzy/Zcs8PSfIeK4azwP798+bNy8qKiorMxurRmU3wD8AB0XGufxuQWGhTCeOS2DAqlFff/3BBx88Tnl8+vRpZLcyOi+zKz+jOt6ys7JxVkBwNO/cuZN4715uHlMyl7e9n5x89OgRcGYlIiIiTe4sOnL4ME6yWrX0bVu3aiFKZ5EQ47qI9abmeag+86otGcHgOu3FiwJcuSXdS8Iu5Pr168+ePcM1G67csEtCBIgAESACRIAIcAQ0uRbX0KyjU6tPAw0tQVlRac7upOJk8SMQnILqBvRhBVt4gB94U61fqdGuH3VFPvrgwcNx7mM9xo27eeOGMvP09Ayk1xYWTdlX8/bv71qvnvjVeGnp6ah9lpSUhoSEot4GQb1tw/oNwUHByryp7scNfRR6TU1NU1PTTp482bFjR0Mjo8ysLCRPbGp1+tQpBIPEGuk131X7Dh1Qut67d+/Qzz//fvGi/Bf5/FGFbTb4OoaGbL0TOWjL1uXuyLNWtm1sO3XqFHfl6mdDPu3bxyUuLk5XV9ehrQM7ym3V58mZoIGqpIND26NHjgwY0H/EiK8ePEjGeu3s7DCkTB4/SsFKm9s0R8DQcfvSbe9vv+ECw3XgQIumFmC1YP78r78ehQpoaam0lgzNKovMjJbNhCg/5+flx1yIfDAckwAAEABJREFUef48R09PjwUI/61at8QtCzQgrNWjlEc4KyC4VNizZ0+Av//1a4rvPwDsMj8/9hIL28CAALbQDld8uRgTU1JSgiuWRg0boSaNyjRG2UOp5nmoDnM2eA4yt+TEBAVP+yAAGbly+fLDRw+1tLUtm1myQ7a2trh/gr8sCJXtoS1LgLZEgAgQASJABBTk0HrtTLTq6wNNYdzTFzGZaKgjzMMbo5uxD2/ACrawgh94Q6O6JCXlcWFhkampydAv3MZ7jh84cBD+y1fo/OixoygfmpmaLly06Pfff58ydYq+nh6riRrt1atXURf0nOC5YMGCFStX+s5fsPyH5cOGDWMV+Fub5jZ7eN+0Ms/Hhxtlh/b98edff+1HivP8+fNfd+9OS31s39YBN8TjY+N8fX2RV0G2bt2a/iQDiWanLp05czQyszK1tLRQn/5y2JeoEzdp0hidqgXB37x5E6seNnxYYOCqoOAtza2by5s4OztPmDhx6jfTfHx8UX62sLBADpeYyORSJaWlMO/eo7u39ywVPIuKiwVlAn19vcGfDpk0aTJ/iiFDhniM9/CeNeubb74ZPXq0qWndgoKC+/fvqzDBpUVWdjaIrVu/Yfny5R7jPa2trY0MDZHe5eXlGRoYoJj91Vdf4VIEM/LnqnJbZkYcaxzx6zduhIWFXbp4CQl9375912/YsGXLli5dPsQuOxFrhWR367ZtuIGwYsUP8+fP9/H1YVN/VqcK2wsxMSjoWjRtigs51KRxcQUnOJTqn4fKmMMPJ2zwHGRuyeHh4ZyO6kZEePiLFy9w2bltx/Yff/xxzFh3HR0dBIlQVRvSKBEgAkSACLz9BCjC6iUgl0NrCPS71kUJuayoNP9kOrIodeZDAm3yXctaTg2Mp9kwaXSZALbwAD/wJqjg0VZ1ZhDrREVGnD17Vl9ff/jw4ePGeTxOfZyf/0I8Vv5XYkJCQIA/yoToRuJy6uTJ5AfS77NYvXrV6dOnUM0d/OmnSDfLykr//PNP1IOhLCOmpqY2vD8fWHzAKbBDlpZNS0tLMdG8uXND94X26OnYpFFj5JRXrl7hNDF6984dTU1NFGv52djfBw6mp6WjE+ls69a2CerVC9etWxcfH29gYODYy7Fxo0YXzp/nJuIaqKwfP3YMtefP/vfZVyNHIkk9ffr05k0boRAZGZGbm9uqVSukShkZGcp4oip56vQpZN4dO3XERDDkBPklkkJjY+OvR48eOHCglpbmoUOHdv/8swqTqKio4M1BmZmZSE/79O1bx8gIHrAQJGdRkVFlZWX9BwyYPmNGUVFRdnY2N9HLNBTOGBgYCJ/r1687c+YMDkfXrl1t27Q5f/58CS4YMCAQsFZZWVlt27ZFsVwotLp169amjRu5Ry9EWpXenDl9OikpSVv01H5cbCxnr/55qIw55woNNngZyOySMaqO4K/Azz/vys/La+vQ9uNu3fR0dU6fPoUg1bElHSJABN5zAg5uQ91E4qDq6xAqBcmOdejmZid7F7VSbkiZCLwaArI5tJaZnlYj5vNhpVmFRXcVPwMqEwmbQGvV1UWurG2up2mmAwXYwgMa8KZlJi4AY7dq8sUXXyDd2bx5E1KZObNnI4FGidfdfaz72LE9e3Qf0L8/UrHJEyd16tiR/ZgaO8vNGzfWrl3by9ERyeIPPyzX0JDm8mmpadO/nT7408G4WY8GdFaJsivWkN2yDuGTL+jEKLb8zu7duo0dMwbJKIZQYsR03T7+eNfOcu+ywCwwGTniKywBaqwg9RwwoD+GZs+a9b/Phnw1YgR02CVwS2Y12RnZIVweYLrRo7+G1aBBA1EMhhUUWE12i1lQQO3j4gwdiLOz03czZqATo8FBwVjvlMmTx4wdfe7sWWU8obnMzw8eoDnOfWzYoUOOPXuwqBHAlEmTBg8exNLD8qEJfQgaCk0wFBr6e/9+/SZPmoRjhyMID/CDfn//lZ99NgSuMDRi+PC+ffpgIkwHQYOdEWqcYKVYL4sCxx0KUIMyFHCG4DwBOrQhymYEhxnTpyMGrB1To63Mys3NbeRXX6FuDW/8ebFbWRnv4YGwu3TuvGLFCs5W2XmI5WBRWBoWyCqDFYjJM1dzyTLBgx6CQSfrnL/F6YETeNLEiexpg5MTQfIVqE0EiAARUETAbvAU/8AARma0VjRehT77QVNEDgMDpit4ZrEKDsmkJhMwMKg947vvfvn11z179+JGq/xS/AMCMMSXDZs22baxZTXNG5gv+v77Y8ePR/9zct++P0aPHsP2v8xWNofWrKujWVsLHksyCkqzC9FQLdIEWiAoeZSftfpWSSrzAUTYwgNs4Q0+0ahGQT6BVFXZI6rsRN2799j3x5+bN20KCAz8ctiXa9b+aGFh8SwnJzZOWgVEcnDs6NFTp06yJm9qiwAiIlAbrtxz51h+hVZIFqEDQUNmdSi+YvlspwqeMIQmtqwmfwtzhfSgrMwE5kgKcewwI9qcsK4wxPVUYwNu5WeEf8QQHh6OqdGWF2VW8pov34MYFJKU96ym5ssHDw8KTxv5kKiHCBCBd50ArY8IvHkCuCUeFLxlxIgRLVu2tLGxaaDoZVnNLC0xxBdhs2YGtQ0QvdDKau3aHwe6upqZmSEXt2xmOWnypEWLF2PoZUQuh66lLdBi6rUlKdJnJDTr6JhMb6FrZywzk7IEmlUTe9DS0IRPtus1bpGY/vHHvvz8/F69enl7z2rXrh3ucaPGhrvqrzEKmooIEAEiQASIABEgAkTgpQgM6D8A2fPj1MfPnj5V6Khz584GhobFRUW4Bb1w/nxWAgMCrl+/Bn3UUq2trQsLC3/99ddpU6deuXpFS0vrY/zp1g2jVRbZHFqgoylgUmipQ01TXePpNnodTI0nNuen0aoTaKk9vMGndP/1tZAxs3elgXL06K/79+sHsq9v+mqZiZwQASJABIhARQQc3Ob6bo9Mvr/VrSLNyo6/Os+VjYT0icD7TCAvP//4sWPLly1/UcA87CCPQlP0p6i4OCEhIezQIVbYG6ooPHdo3x7jFy6cXxUYePbsWb+lfvPnzw8ICLhz57a8K/V7NGVVi0plP0dYWCooZF43pllHm0uj1U2g4R2m8IlGVWXn7l9fRqZ+O8NtxMjZc31fxgnZEoH3nEBV//qS3ftK4HWu237uugDPCS7C6p/z1Xmu/ljJIxF4RwiYNzDv+8knrgMGsOLs7IwkGBXlefPmFcl9ywS35np169aqVVtDQ8PN7cvzFy7EXLx4+MiRoUO/gELHjp1MTM2Ki4ufPcv548+/MLR3796vRnz15ElG2st9L4dsDl2aXywoQdor0GrEvN4Oc5fmFj9dd7vw+jO02TS6dt8GJt+11Kqrix7+M9DY5YvYQ0kZ45M/UMn26JEjSIgAEXizBCr5t5bUiQARIAJEgAhUhYA6XwAs79fAyFBTU0NfX9/M1PTevXt5efn169efPHkS0nEjQ0M9XT1NTU0XF5c6dYzu308uKy1tbdvay2umpWVTWVeV2ZfLoZ8UleaVwINWPT1NEyZLRlsmjTYa2azCBBq28MDY5pWUPilCg4QIEAEiQASIABEgAkSACKggoM4XAMubh/wW8vXXX8+ZPXvY8GFfurl9v3hRVlaWUZ06XT/8kFVGDv3gQfL48eM//99nv/32G8rSlpaWjr16s6NV28rm0CWZBSUp+fClaaqrY818mBFtCD+Nxi5ERQUao7CFBzTgDT7RICECRIAIqCZAo0RAAQF7OwdWFIxVrUvisNpeYywJg40TW0lHNfyGN0aUe2JGsSLlCspGxIZVevWyxFaZ74r6EbBIKtKj8feQQFxc3DI/v4W8zwUmJd1XhwP/dVsRERHZWVkaGhoNGzfKePIkPz8PHm7evAUdNP67daugoFBLU7N2rdrYrbLI5tCCMsGLc0/KSgQaOpq1etTnf76Qn0arTqBhBVt4gB94g88qx0eGRIAIEAEi8D4SsB+6Izoy+X5ictjBMFbQvp8YvX2oA4vDbWt0dCQj6z2t2B6BcyDbI9rukPmAof1Q0ecOE5PvSxyGoQ2Hc91kkunKe64gVHF4Kn/Zz4UTZjnR7Ccj7dx8tkZjyezamVAjd7jZ8VzIK2AtEjg8PZmmgxsmYhZeHmwkA1aGg4wls2vntj2SF9VBHKBoBFyxIWMsEEhi5g4BFng/codPxWGz9m/tlgJ7IwTMG5i7u7v37s3UkmfPnn0hJiYqKqpf//4IBp3GpqZlZWUPkh+cP38+IyMDna1atRRaMf9atGjZUk9Pt6i4ODPrCfqrLHI5tEBQcCW7JJ15sZ2uvbF+ZzO+azaNzjv6OEvyHmj+KNeGFWyxCz/whgYJESACRIAIEAE1CTj4RCaH+bsIFXxM0MqlXyuJFyuhUCSSffwW9zD92JMIk/nBocLPHVq5eAaGRfqWzwJFbuFE4gC/FXsWqBkqHFQovEntfKMPBk5wZv63l5oJXQIOItMVdShUEFi5+IdFzxVfY4j0ym2Yy5LEsABPRR+/FDK2YTJpejlrgf3c6PsHA12E5aMSWAmdAVD2iqW8KbPHzC6/KIwIXSaoDBsqJERAEYFvv/l20uTJvvN9kTFfuHAh++lTozp15s6d9/vvvy9a/L2ZqWl6WvqJEydgGn48vKCgAH+Jt2zZ8seff305bJi2tva9xHthf/+N0SqLghy69FlR/vFUlJBRSDYaaaltUa7QjTQ655ck9otUFM4KfVjBFh7gB94UqlEnEagmAuSGCBCBd4oAstKwCeLsOSExImimt5dIgoIiEhLLrTQhMVEkvE5xD9PP6xW0sFbg0AsOxUrCCWFs9Ve8L3ILJ+Jd5pciz+qHynhQ+6f/9oMMAMnag8ITEyS2yHRRjUaGzSgIEsODRHCCgsM5MkLPdT78crXEUjB0B3NZIt5lwLK2M4PhX9wrECJNV5wNI4EO4+r9goTwYPFxETOE4Y99JF4U/ebPnigOe6Y3phYvTeipKvtX5JH6iMCzZ89KS0tf5Bc8f/4cFeg1q1bdv59cu3YtFJuxRXv9+nXst4L89NOO4KDgzMxMY2Pjpk0tgO78hQsrV67IzWWe8cBu1URBDg1H+SfSCuOz0dA00TWZ0lyrgR7a6gg0oQ8rKMMD/KBBQgSIABEgAkRAPQJDZ0yQ5LtBgxwdPZaGhIaIZKmfh6Oj0HXmmpusoxDsOjk6OjlODRbnYYIIL+xKZGwIqyfeJoR7uzYV8h2GwKErZ+vcn3v2Q13PQ9UNVRyCzC9lu84uLoKEoEEWkrUvdXdydPUOl6gjzZ0AQonBrk2dxvqJ4PgtH+so9JJoWE2Yzi1FYoS6tb+LeCcxyFXEgbUNWQ7/Fq7eQZIs3CWg3OWEyMjOV/rATIQXzN2Xi48LGDYd5MWk8OXL9iIzyYY3Oz/skFBM7cgdAqHnDLm4JR7oNxFgCMTExAzo379Tx45LlizB/ooVKyaM9xzrPgb92D106ND/Phvi5uY2b3p1pw4AABAASURBVN48bNFGD/pZ2bnzp0/69h09+uvZs2Y5O/WeMmlSXFwcO1TlreIcuqyo9NmWhOJkJj3XalzLzMdWr51JhXNAB5rQhyZs4QF+0CYhAkSACBABIqAWAXsbyaMCERv84uVNYkPiY+V7K+iJPzBV6OgeqsAwbvkGSerp0m9oBW5khl9JqKI5EoOnyaw9LnQ1l+QyKolBU5fLLCdkNXc9YNWi/KMpArfpTNrNGro6LZXPHOJClzpyabrzlPKVbAefH3nmHiGy5vEh7k5cBs9MIvMjnR0XObJhC+KWO86MYC1cpih/EIXVoC0RKE/gytUrMu94TkxIOHb0KLblFcV7169dj4iIeMnys9iXQKA4h8Zw6bOi7A13StKYB6M1UY2e3sJ0VisdKwOBBgbLi4YA/Rg1md5CU/Q6PFjBFh7K671NexQLESACRIAIvNUE5BLBl4g2Vjbtk/oKOSJO4ATWNkqfJJaqK2xVZ6iYIHyDXKIpEMQeOC4ptwsE4RsV5cEHj4tryULr1nDDiZ3vFGd2JyHoGwWG7JhAmqZb9R3EQ2E3uC/q3oySCnNeBs9o8n6ks4fP9Ch/b0CiFXJEfCEjtOYedpeM0W8i8PYSUJpDI+SSR/mZfjeK/nsuKBMINDV07YzNFtmZB3cyW9jGxKslK2ibB3VCP0ahA03owwq28FC9oqunp6snfmV19Xomb0SACBCBmkHgnY8y7rYkUxROWC/3xoxqX769nYNNVZ2+qlAT795QFJJ0OkH4kVBFGvH/3RV3W9nwHom2H9RHnAMnHj+goLQvtuGn6cI+g6WV7FaSh8lVmvPC4xwyDensEYcVZ9DQunlXnP1X89UIXJMQgVdHQFUOjVlLswoz/a49Db5Tml2IXYiGnpaOtaFeWxNW0NbQ10I/BDrQhD6ssFuN0sTCYtp3M7/znjPDe57n5GmNGjWuRudquho5ZpydQ7sKlRs1auI+fqJRHWPVmrUNag8fNXr+9349HHvJaH421E2+U0aHdokAESAC7ygBaTVUIPQMDEuM3l6NmbSdm9vcHdu3Mq+QY96qJnpx3gRxjbbyPF9RqAn/Ka+ai4JUkmSLxhRsWltLHo+pyLM0DxZKK9nSR1ZUm3N5cPkQpLNb9d++FfAVyXRJls+bt7wb2qvxBN7FBVSQQzNLLhO8OP0kY8blp+tuF93KERSjKM10S3+Ky9CPUehAE6Vo6VB1tHR0dD/p5xofG/vD0sUBy5c8eZLRXS7vrI55XquPz92GFRYWpjx6+FpnpcmIABEgAm89gVg/J1fes79WzLvnEpOjVb52rcJFMS9WS0y+fzAwwNPFxZl5i1yFJmoovJJQ1Zi3iiqJd8Ufx1RqrygP5pLgis2V+hUNCEFemUiyfJEibYhADSGgRg4tWklZieBFTGam3/VU9wvpky89mXWVFbTRg36MQkekW80bbW1tHV2dokKmEF5SUrLvtz2/7/0Vc6Dcy9ZrUfSdMGUaWyRu0bLltzNnz/ZdMMtnwcc9ekJNV08XCets30Xo9Jg42djEFJ2c2hzfhXJqizi1NvYOXnPmwRD+9XTFj5HI22poaLgOHjLHdxEmderTR6ChgSlUy6GDB7CQ0pJSVo0ttCOYb728jY1lP77ZpetHM+f4zpo3H/67ftwdJhq8Gb/6eoy75yRAQGfffv3hBAFPm+EFn9AkIQICgYAgEIGaRQC5qYWrd3ii+AY/E7xQ6BJwMLkSX+fBGIl/mPey+UtfipyYGB4ewb5ezdV1kKvkA21i5Ur+quZQKzl7zVKXvDEwUUWjZq2Ion3PCaibQ/MxlT4vLn78ghW0+UOvop2fn3fpwoWPu/cYP3lqb+c+hkZGymapbVC7t0vfq/9eWrH0+727d1nb2JgYm3Ts/GGDhg03r1+70m/J45SUdu07QM3R2eXsqX+gti/kt3YdOjZo2AhqpmZm61avRKn7eU5OD8deyEp79nY6e+ok1E6eiDIxM8O8Cm3b2Le1at5857Zg/2VLnj59Wrt2uTdqa2lptW5jB0OYc/JE9JU57C5y3z79BiTevbvC7/uQX38xq1uX7We3MNevVSvs4F8rly2Jijhm365trVq1W7RqLbSyZmdMz0hnZ0RnS9s2v/68EwHHx8V90s8VJXzWCW2JABEgAjWMQFzoWEcnC9dBvLc4CwRC58Aw+TevqV4Z771sicGuTYUWjk5j3T3Yt8LFxlXhLR9y01VbqHKeq7dDWOEn9lpJHn2u3olF3hKDp0leO8i8kVBJe2z5Z6ZFlrQhAm8pgark0K9/KRdjzq8JXBFz7qxlM+HkadPZWqx8GOYNGmppa8fHMS/8uZ907+ft27KfZltaNruXkPA0O6usrOzv/X9GR0VATU9X38TEBBl5M6FQR1u7vnkDqOXn5X3YtVvPXk4CDU1z8wZ169VFnfjmjWsCgSAxISE7KwsNhbYNGzVKT01NSXmEKf67eb2gkHmZCZRZsbJuPnDwkDZ2bdld+a2pmamerl583FWYw0laWipfB6X3c2dO1TYwGPK5W+cuXRG5to6ORVPLjPQ0KMPk3t077IzozMrMxMJhfvPaNZjUq1cfbRIiQASIQE0lEBfPvMW5KTJpribtHLi9Mi+hU/1itWrk8vKhVmMwfFc37ko+plnRJ/akjz4nSj/XyJlXkIIryb/VNedHTG0iUDMI1IwcGiwLCwqv/Hvpp23BJ/850bZ9e9Ri0SkvJcUlpaXFMv2FRUUyPdxucXHxv5cupjx6wPWgkZry6NLFC2ioFoW28ia3/7vlv3xpzPmz8kPq9KAc7j5+Umtbu4cPkv+NiSkqYp5pgSEyaWxlpFj5SmU01dglFSJABIjAW0IAmTTvIWnrSryEzoF733T4EYUlTqlC9ay16qFWz/zyXnifFOwz2E5+XNrTWvLpw8TjB7jPNUrNeV9DI7WRtNz6Sb7DRdLD/lbTnFWmLRGoUQRqQA5dr775NzNnden6EcBqaGg0atS4oLDwxYv8wqJC1I/RU79+PQMDQ4ympT4uLSuxsmLeVNTUstlo9/HmDRokJd1DjVlXTxeaAz/9rG+//lBD4TbzSWZUxHEk0LAtLi6Bmr5+rdOnUKc+nvM8B/pZmVla2lqst0aNGxkbG7NTyNtmpKfXrVe/Xr16ULBsZoVSMRrqS86z58XFxTY2LWECJ2YiP2izYmhgqKOjHRVxDFl4LQMDFNrRn3z/nqmZGSrcWlpa9u06sDMm309CJzxAoVWbNgUFL7Kzmdo5dkmIABEgAjWdQCz31rsKVlKp4aEzJohf/FYpM9XKryZU1XOqGJW+P8RKwVcYcobSFzknHDvI3M8Vj0g/aOgyRdl3oEhtxUbSX6GHxS9/Fig3l2pTiwjUIAI1IIfOSE87d+qko5PzbN8Fs30WfmDR9PiRQ2VlZZcunBdaN5/ts6DfwE+LRNXZvNy8MydP9uztNGve/BGjxqQ8epielnYp5nxuXt533nOh2bRZs7irV/lqE6d+q6ur++xptlTNl/kwYmrK46fZWcionfr0hTfnT/o/zc7GcVVoe/XypQfJ9z0mTYVmPXPzvLw8aKovCB4rsnNwwAKHjxr9LIuZiDN//PjR/aSkUaPd4byNnUNZaSmG/rt58/LFi5+7DfOaPS/3+fO83FxR5w1ojps4GX7ad+gYHRmen1+5SOCEhAgQASLwJgnYD/V1U1wrlRaM797mZXj8YBU8qyBNZ10m+0rfecxa2blt574Bm+1RtlXgWfBSoSqb6JX0x/ptlOSxzoHRCvNgoDgouZqI2FDuWxLjl26QfA2N0DNsu/yDNHxbBfFLv35FZO6gQAVdcDKXvuobIEhqEIHqzqFfzdIvnDsbsNxvW9DmDT+u/nGV/8PkZMxz4/q1wB/8fli6eOPaVetWB8bHXkHntbjY1f4/bA3a4L98yTFRql1YULh39861gf4bflyzce3qlJRHrNqqlcuDNq5f7b/swJ/7kJGzaqv8l29ev+7HgJX3k+5BLfbKlcAVy+Bt66YN24I3cVPI2MJ8/x+hawNWblr/456fd24P3pTz7CnMK5TtWzafjD4BNUy3JmDllk3rN/64hnleRdSJfgjrHMHD+fo1AZvWrWWd3/7vVsAPfiuXLUlMuKulo81qhh34C8vHEhDhf7duoZOECBABIlCjCNj0Yd6/EbnDZ6iDvZ1EhvpujwyTpHiy3zASJ/1+vglhkUjBGSu3ub5sRsZ9B55AKB2FZ7e5O6IPBroIEvhv/5AhpdqzoPKhyvh/fbuhY7nXjyCRvc/gdQMEkbj5bI0WoRCFkxjkKvdtgiFrgrjH0V38k6O3iiHb27lJMAoSI8I5HZEj6SZu+QZJCi9w8Q+LZmZnjpFodgccqe2RzGsHXaylJtQiAjWBQM3IoUESqSQK0mz6iF0VAs3MJ5klJSV8HVRk5W1RaUbqzFfDLjr5PfADb/DJ70QbalBGgxOFU3CjFTYwBSbCdAo1ETyEG2rYsPGXI0aOHT9x0JD/uQ76NDUlhRtFVIiN06QGESACahEgpbeKgFDoMsE/LOygRPwnSF5NlxA0SO7VDfEHjnHpm3BCgMgqwFOSkUmfZBAgjWZH4TnAk3EZ7u24IUH50lV7FtlVLlSRyRvZhHi4zoyQLJXBGwgIIgmc4Gwlfp4FCbTTUu5JaGmc8UsdvaUpstBZDDlM9L5t2CZGeDmukXxJotSMa4W4D/Li7GWIBXhyB5fTpwYRqBEEakwOXSNovrYgU1Iebt6w9uKFs9lZmb/s+gml9Nc2NU1EBIgAEXiFBFD65ZItmWmQqLkKHcs9ZiDWiGW+liVCkiCKO7lfzKg0feS6E8ODvF3dFX5pNqcjYGyDlHiuUqhS16+9FRvi4ejqHaQEb0K4t2tThQk0G2joWEdeHsz2ibaMoaNc6Vo0xNvEh7g7MUk8d7HDG0MzITzYS77+jQGSCgmQwpsjQDn0m2P/cjOj3hx75crJ6BOpj1NezhNZEwEiQATeHgLxS92dLJoKXV0Hec30FovrINF7nT1CFJRIxZHH+nk4MlasCaPPL1cz6aN01NvLdZBFU6exfqGxsA7xwHSMOC5ndtFTXpR7rmKo5d1L9uKWI34mjKbK8tHQsU2FIgWlmW6IO6ug+EqDmSkuVIKXBcVsgRpuHd1FNBglZT9MHgxNV1fGijk0DEahxBC1avHsfPJ8X8xRcBRa4GhyR3bmIFdXxsrRfbmCg1sxE757ahOB102AcujXTZzme5cI0FqIABF4RQRi4+JDQkLFovY3ocTGsSaKvzlFMhoaEhdf2bAltgo8Vy3UygZQvfqS5TC4EH+lnEttK4+RmQhHkzuyIfGxyq+LGGX6IQJvMQHKod/ig0OhEQEiQASIABF4BQTIJREgAi9PgHLol2dIHogAESACRIAIEAEiQATeLwKUQ7/+400zEgEiQASIABEgAkSACNRsApRD1+zjR9ETASJABF4XAZqHCBABIkAEpAQoh5ayoBYoHPrNAAAQAElEQVQRIAJEgAgQASJABIjAu0XgVa2GcuhXRZb8EgEiQASIABEgAkSACLyrBCiHflePLK2LCLwdBCgKIkAEiAARIALvIgHKod/Fo0prIgJEgAgQASJABF6GANkSgYoIUA5dESEaJwJEgAgQASJABIgAESAC5QlQDl2eB+29HQQoCiJABIgAESACRIAIvM0EKId+m48OxUYEiAARIAI1iQDFSgSIwPtDgHLo9+dY00qJABEgAkSACBABIkAEqofAu5RDVw8R8kIEiAARIAJEgAgQASJABFQToBxaNR8aJQJEgAi8agLknwgQASJABGoeAcqha94xo4iJABEgAkSACBABIvCmCbzv82s2s7YmIQJEgAgQASJABIgAESACREB9Apr37t4lIQJEoOYRoL+5RIAIEAEiQASIwJsjQM9yvO93Imj9RIAIEAEiQAReHwGaiQi8KwQoh35XjiStgwgQASJABIgAESACROB1EaAc+nWRfjvmoSiIABEgAkSACBABIkAEXp4A5dAvz5A8EAEiQASIwKslQN6JABEgAm8bAcqh37YjQvEQASJABIgAESACRIAIvO0E1Mmh3/Y1UHxEgAgQASJABIgAESACROB1EqAc+nXSprmIABF4nQRoLiJABIgAESACr4oA5dCviiz5JQJEgAgQASJABIhA5QmQRc0gQDl0zThOFCURIAJEgAgQASJABIjA20OAcui351hQJG8HAYqCCBABIkAEiAARIAIVEaAcuiJCNE4EiAARIAJE4O0nQBESASLweglQDv16edNsRIAIEAEiQASIABEgAjWfAOXQ1XMMyQsRIAJEgAgQASJABIjA+0OAcuj351jTSokAESACsgRonwgQASJABKpGgHLoqnEjKyJABIgAESACRIAIEIE3Q+BtmJVy6LfhKFAMRIAIEAEiQASIABEgAjWJAOXQNeloUaxE4O0gQFEQASJABIgAEXjfCVAO/b6fAeqvv3PnzrZtbNXXh2YVTGBFUikCH3frJrSyqpQJKRMBIkAEqp2Ag9tQN5E42FeXbzvWoZubnUP1uCQvRKA6CdSMHNo/IGD37t1Dh35RnUuvyNc8H589e/fKCDphhy2//6edO0ePHmNgUBtDfFmwYAHUdu362XXAAH5/9bY9J3hiFoRUvW5lvI0cNWr1mjXr1q3r3bu3zJCy3SqY8F11795jx08/bdi0qbKJO9+JOm2cVzi7cI6pVgZhcAZt1Wr80X79++PcwCpkoE2ePBmu1qxd06TJB3z9KrS9Zs5czfxZVS2UWBSIjRMEjwNRhcDIhAi8BQSG7ohOTL6fmBw9lzKwV3847AZP8Q8MYGRG62qazX7QFJHDwIDprarJJbkhAtVIoCo5tKahtnZDfVbQrsZolLlqZmnZvHnzevXqKVN4Ff0fWHxgI/cHnZgLW/6InZ3d1GlTkW107NQRo6ygBNv1o4+ghuTG2aUP2/kqtg3MG2AWhKSm899///3cuXMTJ05SrY/4Dx0+HP3PSfYCICEhIT09IzU1LSMjQ5lhBSbKzJT0G9cxsrKyFjZrZlDbQIlK9XTjvMLZhXNMtTsQBmfQVq3GH71x/ZqRoRFODxcXF64f11q9eveCq/y8/IcPH3D9VWskJd7LyspKfZyWdO9e1TzwrVq0bNGyVSvExgkOAQ4EX4faRKCmEHDwmewiFAUr9FznYydq0YYIEAEiUG0E1M2hNbQE+p3NzHxsG2zvUn9jx7or27KCNnrQj1HoVFtcFTlCIvLhhx9iK6+Iznbt23P9yAZU3Olu17adeQNzTlm+cSHmQqeOyI3FMpmXerJDjj17BPj7p6elI9uYNHkyZmedwMjM1DT1cWp+fr5NSxt7e1U3t5BnQ1hD+S2Ww7nlRhE2gud2lTXgFqJslOuHf2U8WZ0zp09/NuTTEcOHx8XFsT3YIjAVbBWaIGxMBFt5QT9G5fuV9SDmnj17ciaIBPHIK0MBnuX7Ya5Qn9WElTp4WWVl26Sk+zdv3sBoq1atLS2bogHp7eTcoEHDvLz8S/9exi4rmAszsm2ZrcwJjJixUk4ndF9ov08+mTDBMzc3j+uEK4VLZhVUzKWlqamhobH/r7+4c96xZ4+wQ4dYQ9oSgZpLIOF2fNWDJ0siQASIgCICauTQGgL9bnXrrW5vPM1Gp6WRQFtD1o+2BvoxCh1oCuTGZfVfbv/jbt1+3bMnMurEho0bT0T/s/uXX5DewuXGzZsuXrq0ddu2Q4ePrl37I6qngz/99PCRI1AOCQnZ8dNP+/b9wdVfkT8tWbLkzNmzW7dvCws79Meffw2o6uMWyF327t27ZWtwTk5OC5uW/QcMRDAQRKWlrf3vv5eys7PqmtXt1KkTOuVlzpzZJ6Kjd+36GYLGd15erM7vompxYOAqdG7dujUiInLZsmXsEFKo7Tt2HDz4N4I/euxY3Xp12X6ZLVxFRUXBLeTU6dPzfHzYOjHMtXV0PMZ7YAqYKOQJegGBq8zNzQFq8ZIlYIse1KRRmYYTdIIefCIwOMEsnhM8oaDCBBMNGDAAnEEbBw4HYs3aNUj10I94tmzZgh70Y/Sv/RUcCwRzISYG5n//HbZq9eq9e/aOGDECxxdHGfEcOnQYE8EtBA2FMyJ+GAIp9IHXxsYGypzACjEgEuDFGrFS6HOjlW2cPXcuLy+vvrl55y5dWdsO7dvXrl07PS0t5sI5eIZ/+fMQMEH72PHje3/7DSew9yxvVhPxIGaOORzifgLQoQdtCIJXuGRAw98OQMMJg3Xh5Nm0eTOX1sOQlUZNGpeWlhYWFrK7tCUCNZpArN83XuGJCYmJCUHeY0Nq9FIoeCJABN5GAhXk0JqmumY+bYw9m2ua6LLhl70oKbr7vOBqNitolxWUsEPQgSb0YcX2vPxWxgP+1//mm2+Q9CQl3UN57F5iYqtWrZBGcGq4b/70aVZCwl0Ts7oeHuNwm/7Ro4eRkZH16tdr8kETTs3La2bfTz55npt7/NgxFFabNGk8ZcpUZIecAtfQ1dVFQsNKv/79kfBxQ/zGn3/8mZGerq+vZ23F3DtEYmppafnixYu4+Li7d+/CSecPu/D12TYiHzLkMw0NzWhcDfzzj5aWFnaHfuHGjiIF7/pR1/tJ9+/cuaOhqdm9R092aNq0aahqo7wNK9zH79xZgWeEOmjgQD19/fDw8D179iASV1fXbt26R0REpKenFxcXX4i5gPTR0rKpQp73k5OjoiJznj1DOnXixIlzZ8+xIbFbN7cvnV1c4ASr/vvvvzW1tId9Oaxhk8YqTMAWhC0sPgANMAd5BDNjxndwiASxfYcO6RnpRw4fvnbtWuPGTcaN80BgGFImmpqaXbp8eDcxIeleklGdOlOmTm3SpMnp06eePX2KvPyz//0PhvwZo6KikMhixjlz5mJo6tRp3bp1Ky0rQ/Z5N+EuIGtqaaEfwlo1bNgoNi4WcYIbzhOcLRiqmkRFRqSmPq5VS799u7bwgFTYzp65p3zz5o2kpPvwDP+ggbni4uKalD8PTUxMGjVqlJCQ8CD5gULmsIVPTtjgWcgyS2Z1UJzOzMy8cf1GWWlp+/btcZHJ9nNbPT095NBdPuwKMjEXLyLh/mrkSG6UGkSgphGID3F3cnR0cvQLrWmRU7xEgAjUAAKqcmitxrXMfFrrtDAUaAgEpWWF8U8zF8WnTbiYufhaduAtVtBO87yIfoxCB5rQhxVsX8XqkXZ8//3ipUuXzp49Z+H8+SglFhcVNWncBDe42enir8UPHzZs7JgxRoaG5vXNkS8uXbJ09qxZG9ZvyM/LY3Vs29h27tK5oKAwOCh47ty506ZOuX3nNqq5yKtYBf4WN75RiGVl0cKFn/T9hD/Kb6c/yUBuh3QZnV06d65Tpw4q0PHx8Tdv3EKQza2bI7HGEF9+/nnngoULFy1c4OX13XczZjxOSUGy1aKFtCx69MiRr78eNc59LK4Z2CHkSbhsKCkp2bcvFFYYSryXyPfJtpFI6dfSz8rM2hcaGhgQMHfOnOnffrtmzWq0n+fkCMrKYq/Ebli/XhlPpO9hf/+d/+JFUVFxVETErp07WbfstlHjxjo6Oo8ePUL1fdHChT7z5k2dNnXblq0qTJycnOqb10cCPXuWN5jv2rUL1z95ublw6L/SH7XYeXPn+fr6/v7bb/n5L3AsWtu2wZAKuXDhvIf7uJUrVzx58gRp39atW6Z/Oz00NBSoTYyNYcif0XvmzODgYCTEdvb2ffr0wS0CHKmoyMipU6aMG+seE3MBuzCBsFZIZ9GPOH/ZvRuo23foqCynR8qORJa9xMLW2dkZWTL8cJKbm/fv5ctlZWWtRI9zSB7kyDt77lyF52FBQcGmjRu/dHNb5uenkPmxo0e5idBgg2ch85fcs2dPjEIuX7mCvx2jRo3898plbW3tli1lP6VTS78W+o2N65w5c/bevXtmZmYe48bJp9pwRVIhAVIgAkSACBCBd5uA0hxas46OyZTmWub6WH9pdmH2mv+yVt4sSsgVlKGjvJQJ0I9R6EATY7CCLTygXe2CHMvB3n7ZsmXhERGjx3ytraPDJO7a2uxESYn3kLWgXb9+PQylPE6JiYnB7pHDh7kPwwktm9WpY6Knp4sSLO6YHzp81NrKCqmDkZERNGUEmQQK3qwcPnz4xo3rMgrcbh2jOkiVSkpL0WPf1gEOHz54hLkys55kZWcjpe7YsQOG+IJQkXINHz58/4H9J0+daiYUamhoaGmKD0pJcXFaWjr0oZbx5Ak7ZF6/voGB4bOcnKtXrrJDt27cRENGrly+kvkkq0HDBuvWrz906PCoUaMMDQ1ldNhd1TxZHZktLgyQiLdo0WI//hzY7zrQtVatWjI6MruWwmbIUxMTEpC1YwhJ+RdffIHUGe3c3OfNLC1nz5lz7Phx5NEAoqGhqaWB6zYMKpUnGU8wVir+U5ab81y0Kz07ZWZEZpmTk2NkaNCtezdDQyPk01djGYCwYi9y0ICwVra2tjgxIO7jxuFqAaeKuXkDjMpL546dfHx82UssbL1mzrSVy/5xLHKfP2cf52hjawtWyQ+SUZ/GuaH6PETACJudVB3mbPAcZNjCA5bcslVL1knKw0f8hqaW+ExjO7EdMXz4gAH9Px08CJd0Y0Z/DUQo8/dy7IUhEiJABIgAEXgfCNAa1Scg+58oa6mho1lnvJW2BfOytpJH+Zl+1wuuZLNDKrbQgSb0oQNbeIAftKtR7O3t165dg8JY/fr1kd3GxFxE3VEd/0xmplluscXFxY8fo5b68NGjh/fuJd2+fTs1LVXeVVp6GgrerCxevDg6OlpeBz2oMZubmxcWFj548MDFxQXpETo7de6EvGrWrNmIFik1VyzHECuzUU6fM9fBwQHl3jNnzqSlprH9L7/FlcOcObOR+qekPDY2MUZ4fsuWeXpOkPFcNZ4H9u+fN29eVFRUVmY2Vo3KboB/AA6KjHP53QK5B21xXAIDVo36+usPXXgXVQAAEABJREFUPvjgccrj06dPI7uVN6xyj/yM6rjKzsrGWQHB0bxz507ivXu5eUzJXN72fnLy0aNHwJmViIiINLmz6MjhwzjJatXSt23dqoUonUVCjOsi1pua56H6zKu2ZDYYbHESsrFhiytS9OjX1seWhAi8jQTs7RxYqbbgxA4r6U9s9UreoMcuUMWH0itUULYYsWGVXr0ssVXmu6L+V0msorlpnAhUF4FyaSXntFYvc107E+yirpy94U5JagHa6gg0oQ8rKMMD/KBRjdL1o67IRx88eDjOfSzuMt+8cUOZ8/T0DKTXFhZN2Vfz9u/vWq+e+NV4aenpqH2WlJSGhITi1jZk+7ZtG9ZvCA4KVuZNdT9u6KPQa2pqmpqadvLkyY4dOxoaGWVmZSF5YlOr06dOIRgk1kiv+a7ad+iA0vXevXuHfv7594sX5b/I548qbLPB1zE0ZOudyEFbtpa9Iw9D2za2nTp1irty9bMhn/bt4xIXF6erq+vQVvZfePV5wicnnTt3dnBoe/TIEdQsR4z46sGDZKzXzs6OU5BvPH6UgpU2t2mOgDHq9qXb3t9+wwWG68CBFk0twGrB/Plffz3q1KlTpaXSWjI0qywyM1o2E6L8nJ+XH3Mh5vnzHD09PRYg/Ldq3RK3LNCAsFaPUh7hrIDgUmHPnj0B/v7Xrym+/wCwy/z82EssbAMDAthCO1zx5WJMTElJCa5YGjVshJo0KtMYZQ9l+fNQ6XmoDnM2eA4yt+TEBAVP+yAAGcEUf/z5F+4GDP18KDvUqEljNF7kvcCWhAi8RQTsh+6IjmRe/Bx2MIyV+8x7oKO3Dy33b5z9XKhFR0dGR28Vf8qEW4PckIPb1uhoOBE7ZJzDSkXmyriyc/PZGs1MLbYKY9qRO3zKh8FoqvEjG5LEObvAMMQWucON/8+svEKiLAFF0zq4AQu8JSaLPR9kw2ZsK1gv3Nm5bY9kliyxBSgGb8WGsIVIYr5fHcTgj4QIvFECmvKza9bRqdWngYaWoKyoNGd3UnGy+DFieU2FPdCHFWzhAX7gTaFa1TpTUh4XFhaZmpoM/cJtvOf4gQMHaUme4pBxePTYUZQPzUxNFy5a9Pvvv0+ZOkVfT4/VQY326tWrqAt6TvBcsGDBipUrfecvWP7D8mHDhrEK/K1Ncxvu+ybQmOfjw42yQ/uQd/y1H/nH8+fPf929Oy31sX1bB01NzfjYOF9fX+RVkK1bt6Y/yUCi2alLZ84cjcysTC0tLdSnvxz2JerETUQpC/pVCIK/efMmVj1s+LDAwFVBwVuaWzeX13d2dp4wceLUb6b5+Pii/GxhYYEcLjGRyaVKSkth3r1Hd2/vWSp4FhUXC8oE+vp6gz8dMmnSZP4UQ4YM8Rjv4T1r1jfffDN69GhT07oFBQX3799XYYJLi6zsbBBbt37D8uXLPcZ7WltbGxkaJiYk5uXlGRoYoJj91Vdf4VIEM/LnqnJbZkYcaxzx6zduhIWFXbp4CQl9375912/YsGXLli5dPsQuOxFrhWR367ZtuIGwYsUP8+fP9/H1YVN/VqcK2wsxMc+ePbNo2hQXcqhJ4+IKTnAo1T8PlTGHH07Y4DnI3JLDw8M5HRUNxIPSu6mp6eQpk/0DAnb/8kvHDh1znj07EX1ChRUNEYHXTMDBJzI5zN9FyHx6W2ZqK5d+MhUFK6FQJDKKzK6oH6NMG3lhWICz6APhzK74R+gcGBbpqyw7tEcefzBwgrPcd4QKXSb4h1XpW114Idn5Rss7F7oEHESmKwpPoYLAykXl1EzMiWEBnuLXZoscSTZCxjZMJk2XDLK/7edG3z8Y6CKUWbIVAwqGrJLyLTO7/KKgX3ViMCYhAm+QgIIcWq+diVZ95u5tYdzTFzGZagbHPLwxuhn78AasYAtD+IE3NKpLoiIjzp49q6+vP3z48HHjPB6nPs7PV1wkS0xICAjwR5kQUyNxOXXyZPID6fdZrF696vTpU6jmDv70U6SbZWWlf/75J+rBUJYRpBQ2vD8fWHzAKbBDlpZNS0tLMdG8uXND94X26OnYpFFj5JRXrl7hNDF6984dTU1NFGv52djfBw6mp6WjE+ls69a2CerVC9etWxcfH29gYODYy7Fxo0YXzp/nJuIaqKwfP3YMtefP/vfZVyNHIkk9ffr05k0boRAZGZGbm9uqVStU6DMyMpTxvHL58qnTp5B5d+zUERPBkBPkl0gKjY2Nvx49euDAgVpamocOHdr9888qTKKiooI3B2VmZiI97dO3bx0jI3jAQpC3RUVGlZWV9R8wYPqMGUVFRdnZ2dxEL9NQOGNgYCB8rl+/7syZMzgcXbt2tW3T5vz58yW4YMCAQMBaZWVltW3bFsVy/P9769atTRs35vLevixSrNzmzOnTSUlJ2qLrvbjYWM5Y/fNQGXPOFRps8DKQ2SVjVB2BMmjUrl0b50bLli1xLDDvgf371bElHSLwGggggQ6bIM6eExIjgmZ6e4kkKChCvX8+FcTYwicSeaFAkBgeFCzyFhwuKjeIVIUTwuRq2MzA0B1MHs+0RIaSMMITE9g+oWfV0mjWuv/2g8wqJQsM4twKmCwZ1Whk2IwCE7No6iDEzJoKBELPdYq/UIYfs4ChFySynRkM/5yxS8DBHbJFe9EgEugwTy57TggPFsMHeWZc6BLwo8rvEuPPnhguntobU1cLMSYE+iECr52AXA6tIdDvWhcl5LKi0vyT6ahEqhMSEmiT71rWcmpgPM2GSaPLBLCFB/iBN0EFHw+reIYvvvgC6c7mzZuQysyZPRsJNEq87u5j3ceO7dmj+4D+/ZGKTZ44qVPHjuzH1FiPN2/cWLt2bS9HRyQEP/ywXENDGkdaatr0b6cP/nQwbtajAZ1VouyKNWS3rEP45As6MYotv7N7t25jx4xBMoohlBgxXbePP961s9y7LDALTEaO+ApLgBorSD0HDOiPodmzZv3vsyFfjRgBHXYJ3JJZTXZGdgiXB5hu9OivYTVo0MBvvvkGVlBgNdktZkEBtY+LM3Qgzs5O382YgU6MBgcFY71TJk8eM3b0ubNnlfGE5jI/P3iA5jj3sWGHDjn27MGiRgBTJk0aPHgQSw/Lhyb0IWgoNMFQaOjv/fv1mzxpEo4djiA8wA/6/f1XfvbZELjC0Ijhw/v26YOJMB0EDXZGqHGClWK9LAocdyhADcpQwBmC8wTo0IYomxEcZkyfjhiwdkyNtjIrNze3kV99hbo1vPHnxW5lZbyHB8Lu0rnzihUrOFtl5yGWg0VhaVggqwxWICbPXM0lywQPeggGnaxzbotZpk6ZgrMF5wxOsE/69v1l925ulBpE4E0TGDpjgiSBDhrk6OixNCQ0RCRL/TwcHYWuM9co+IR1BUE7T4DPcG/Xpk5j/ZaLvC0f6+jkGsTctROZOveXzSlRA/YXf/VoYrDIUBKGu5OjazCXFM6QNRT5q3jj7OIiSAgaZCFZ4FLGrTd3OwlpLkIW8Kf2Q8xCL4mG1YTpcjPzYhYkBrkKGXp+bNjL4d/C1ZtbsUuA/GWDne96LoGO8IK5+/KlIvIhIN90kFc4cImL+orWx5udH3ZIKKauDmKK5qS+d4jA0KFf7N69e8/evTKC2638VX700Ue4vRweERH9z8l9+/6YMnWqgQHzuT5Wx3XAgF27fj4RHY16E8pDUGb7X2Yrm0NrmelpNWLesVCaVVh0V/HnqGTmYxNorbq6Ag2BtrmeppkOFGALD2jAm5aZ+CEK7FaL4H96pKrKHlFlp+jevce+P/7cvGlTQGDgl8O+XLP2RwsLi2c5ObFx0iogMphjR4+eOnWSNXlTWwQQEYHacOWemcHyK7RCsggdCBoyq0O5EctnO1XwhCE0sWU1+VuYK6QHZWUmMEdSiGOHGdHmhHWFIa6nGhtwKz8j/COG8PBwTI22vCizktd8+R7EoJCkvGc1NV8yeBxBnDM4weQDoB4iUJ0EKuvL3kZSB43Y4Bcvbx0bEi/9911+WFkPsjr3UBnDWL9vpDllP/EnBMQO3KYzKSyzE+HluFzGUBC33HFmBDMoELhMmVvu+Wy2V51tYvA0mQXGha7mAmI8JAZNlZ06ZLUkfRdYtZB5BEUaMxJop6VxjItyP3GhSx25NN15ig//wWuBg8+PkiXD3CNE1jw+xN2Jy+DLuWV3pLO/MmLsRLR9RwnUq1evefPmvGcCxM0GvFdmfdyt2/z5C1ARMzIy0tHRtmxmOWbMmMXfL2GRDP186MxZs2zb2BoaGhrVqdOuXTs/Pz90sqNV3srm0Jp1dTRra8FdSUZBaXYhGqpFmkALBCWP8rNW3ypJZT6ACFt4gC28wScar1mQmP7xx778/PxevXp5e88CL9zjRgkWd9VfcyQ0HREgAkSACFQrAbkc8SW8h2+QTUZFzuIPHENtVdS0tuGlwna+U5xFvYLwmR4hbEtmG3JEXBEWWreSGVJvV2FIsQeOiyvccBK+UVEefPC4OGShdWsocSKNOSHoGwWGYkVpmm7VdxB/yYP7cuV/pea8DF7sTvJLOvurIyaZi36/mwSOHju6aPHihfPnsxLyW0hRUWFxUVF6ega34L59+5o3MM/MyvL18cH9cFSCMNShfXsXFxdLy6bDR4wwMjJCHcrru++WLFmChLCOsfGnQ4bwC9XQV1ckenI5dC1tgRbzzENJivQ5Y806OibTW+jaGUusxL+VJdDssNiDloYmfLJdr3eLjLl3796TJk4EdNyY7t+vH27uv94QaDYiQASIABGoJgJxtyVJpHDC+rluMqXWKk6SePeGYstYbjb+uP2gPuJ8MuKw4gwa2jfvinPZquX6SkKSLl8QfiQU08hJ/H93xX1WNrxCsjTmxOMHFNTvxTYCgTRNF/YZLMXbylq8ZJXmvPA4h0xDOvurI8bMQz/vMAHcN8b95LBDhyAnok/Y29vr6OjeSbj788/Sh2b1dHVB4NnTp2fOnMbd1KR797hPOn3YtSvS6+Li4sNHDkVHR+//669TJ0+WlZU1adL4o48+hlWVRTaHFuhoCpgUWupQ01TXeLqNXgdT44nN+Wm06gRaag9v8Cndr3Rr5+5fX0amfjvDbcTI2XN9X8YJ2RKB95lApf/SvkoD8v0eE5AWSgVCz8CwxOjtL59JJ/wn+2SCSsCtra3E41b9t2/doVimS/JsmXqw2LKiXxWGpCTJVuZXGnNFnqV5MC9y6SM0qs25K4fycUhnf3XEys9IezWcAPLdvp984jpgACvOzs78arGb25fWza0LCwsjjocjV+bWevbcuZxnzywsLJYv/2HixEn9B/TX0ta+cevW2bNnhM2Eenp6BQWFCXfEV5n3kpLgoXbt2s2b23AeqtDQlLUpKpX9HGFhqaCQeWWvZh1tLo1WN4GGd5jCJxpVldEjR5AQASLwBglU9e8u2RGBaiYQ6+fkGiSu8cK1lQuTSSdHR+4o9+JkjLwGEbq4OCsTSZ79GsKozBSJdyv6zKWiPIwThykAABAASURBVJhLgis2Vx1MDSSmekGVHCV1NQl07tjJx8d38ZIlrHjNnMl9n4OBQW3nPi66uroJiQkhIb/xHR7Yv3/Xrl15eXkfd+vmMd6jQYOG8fHxWzYHIc+GvqamZm7u87T0dNbkSXp6cXEJ236ZrWwOXZpfLChB2ivQasS83g6uS3OLn667XXj9GdpsGl27bwOT71pq1WXK5vxnoKHAF7GHkjLGJ3+A2kSACBABIkAEqkQAabSFq3e49PVzAoFQ6BJwMDl6azU93aFuWAmJiRWKur7eD70KcUHh/SBBq1RFQMUXAKMIbSW0QglZpggNd9OnT58wYUKtWrWuxl49cvhwaupje3v7JX5LOncWfy+Hvp6esbH4mWQtHR0NDQ1YvaTI5dBPikrzmNxcq56epgmTJWMCmTTaaGSzChNo2MIDY5tXUvqkCA0SIvDmCNDMRIAIvEME4kLHOjpZuA7yCoqQPCGNTNo5UPG7nF/NwhODpzk6OVYkY5U+M/1qolLtVVjhZxxbWYsffVbtqEqjNZFYlRZKRi9JIC4ubpmf30LJJwgDAwKSku7Dp4oitG0bW2cXZ20dnSuXr4wb6+7r6xsUFJybm4tqtJOT0+PUx0VFhXr6+o2bNIEfiJVQqKurU1hU/CjlEXarLLI5dElmQUlKPtxpmurqWBugwQo/jWZ7VFSgoQBbeEAD3uATDRIiQASIABEgAtVGIC5e/HJi6dMdzoHby7+Hrtomkzi6cVectVecj0pMXup3dRhzMcu/807GvfTRZ94j15x5BUtWkn+ray4TCu0SAVkC8kXojp06enpOEFpZ6erqaWlqw0BPXw+pNhr6+rU0NDTxR1dX99q16zk5z/X09D766COMQv/jbh9ra2tnZmbcuXMbylUW2RxaUCZ4ce5JWYlAQ0ezVo/6Ag2pZ34arTqBhhVs4QF+4A0+pV6oRQSIABEgAkSg2gggk3aSPiRtzX8PXbXNIXUk/dSd/HevSLXerpY0ZmGfwbz3dchH2VryicnE4we4j1pKzVUu2a2f+HtnZNyqaS5jRbtEoDwBA7knoS0tm86b6+M5wXPxosVXLl++eetmWVmZra3ttu07fvzxxymTJ9WuXSvn2bOrsbFnTp8+c/ZsSUlJp06d9h84uGXLFqHQqri4OPpE9Et+E4JcDi0QFFzJLklnXmyna2+s39mMvwo2jc47+ph7DzR/lGvDCrbYhR94Q6NGCgVNBIgAESACNYGA4vfQvZLIQw+LX/78El+h8koCU+FU+j4TRV9hyBlKX+SccOxgLNctuCl5VZ+KJUttpXbiVk0kJg6dfr09BNgidEFBQfgx8es4MjIy8vLykDc/y2E+sDff1+fI4cNFRUXNmzf/uFs3QyOj+/eT16xde2D/fqzCf8UPx44ehbmJ6E9ubm5ISMjq1asx9DKiIIcufVaUfzwVJWQUko1GWmpbSL8pETMhjc75JYn9IhXsygv0YQVbeIAfeJPXoR4iQASIABF4VQTeVb/2Q32VvH/DgfsGw7u3ecnfKwEh/TIRoWfY9qEOiiexc9s+103x0BvojfXbKMn8nQOjFX57IgI+KPkywogN5b4lMX7pBvE3LwpES5ZbAN9WblAgqInEFCyDut4ogR07dnz80UfdPv74p592sIHk5uZNnOg5fvz4qVOmoAe78+fP796t26SJE+fNm+fqOuB/nw3Z/9dfGIJwo9O/nQ7p5ei4KjAQ/S8pCnJoeMw/kVYYn42GpomuyZTmWg300FZHoAl9WEEZHuAHDRIiQASIABEgAi9NwKYP8/6NyB0+Qx3s7SQy1Hd7ZJgk+1Py5SMvPTPfQdzyDZKEVODiHxZdPh63uYgn+f7BQBdrvtGbboeOlXwDOZMH32didpMwdPPZGh2NgNkYE4NcPWQ/CRmyRvrMuYt/cvRWXMyw/N3c5u5gbRMjwqVvHWRdSbY1kpgkePr9OglUci5kxlcuX5YxiomJQck5LTVNpp/dPXXqJIRtv/xWcQ5dVlT6bEtCcXIeJtBqXMvMx1avnQnaqgU60IQ+1GALD/CDNgkRIAJEgAgQgeohIBS6TPAPCzsoEf8JLuJ3SSQEDXo978EIcR/kxSWMMvEEeHLxVM96q8tLiIfrTO41JgzDQAnDwAnOVmKESKCdlnJPQkunjl/q6M2tWCB0nhAg5h8Y4MngT4zwclwj/voKqZW0VSOJScOnFhFQTEBxDg3d0mdF2RvulKQxD0ajrmwyvYXprFY6VgYCDQyWFw0B+jFqMr2Fpuh1eLCCLTyU13vZPV09PV098ev2XtYX2RMBIlCeAO0RgbedQNzB49I8rnywyOFchY7lnkAor1DNe/Eh7k5MSqqk8poQHuwlX82t5hgq7S42xMPR1TtICcOEcG/XpgoTaHai0LGOvCsHtk+0ZQwd5UrXoiHepkYS48VPTSKggIDSHBq6JY/yM/1uFP33nHmxhqaGrp2x2SI78+BOZgvbmHi1ZAVt86BO6MeoQFMDmtCHFWzhobqkiYXFtO9mfuc9Z4b3PM/J0xo1alxdntX3M3LMODuHdhXqN2rUxH38RKM6xqo1axvUHj5q9Pzv/Xo49pLR/Gyom3ynjA7tEgEiQATePwLxS92dLJoKXV0Hec30FovrINemQgvkcPLV07jljhhiRC7DUzHEYQ3xwFyMOC5X+Iw1k5I6Ci0QABfMzEGurkKYOLovD5GPh/OssFFxSKFjmbXAv9JMN8Qdo4wovZyIC5UwlACc6Q2eophDFS6TFyyTB0PT1VVi6zoIu47urCFq1czU6FF2N6DSxCpmwkRHP0TgTRFQlUMjptKswky/a0+D75RmF2IXoqGnpWNtqNfWhBW0NfS10A+BDjShDyvsVpfo6Oh+0s81Pjb2h6WLA5YvefIko7tc3lldc702P5+7DSssLEx59PC1zUgTEQEiQATeDQKxcfEhIaFiiYuvKPN7xYtGAFwwIfGxlU2dX3F0ytzHxkkAhoSCpzI1hf1S27h4hQoVdNZMYhUsiobfSwIV5NAMkzLBi9NPMmZcfrrudtGtHEEx803gTD/3U1yGfoxCB5ooRXMj1dLQ1tbW0dUpKmSS+JKSkn2/7fl976/wjHIvW69F0XfClGlskbhFy5bfzpw923fBLJ8FH/foCTVdPV0krLN9F6HTY+JkYxNTdHJqc3wXyqkt4tTa2Dt4zZkHQ/jX0xU/RiJvq6Gh4Tp4yBzfRZjUqU8fgYYGplAthw4ewEJKS0pZNbbQjmC+9fI2NpZ99LxL149mzvGdNW8+/Hf9uDtMNHgzfvX1GHfPSYCAzr79+sMJAp42wws+ofmKhdwTASJABIgAESACROB9JKBGDi3CUlYieBGTmel3PdX9QvrkS09mXWUFbfSgH6PQEelW8yY/P+/ShQsfd+8xfvLU3s59DI2MlE1Q26B2b5e+V/+9tGLp93t377K2sTExNunY+cMGDRtuXr92pd+Sxykp7dp3gJqjs8vZU/9AbV/Ib+06dGzQsBHUTM3M1q1eiVL385ycHo69kJX27O109tRJqJ08EWViZoZ5Fdq2sW9r1bz5zm3B/suWPH36tHbtcm8D1NLSat3GDoYw5+RJRgbXRu7bp9+AxLt3V/h9H/LrL2Z163JDaMBcv1atsIN/rVy2JCrimH27trVq1W7RqrXQypqdMT0jnZ0RnS1t2/z6804EHB8X90k/V5Tw4YGECBABIkAEyhOgPSJABIjAyxJQN4fmz1P6vLj48QtW0OYPvaL2xZjzawJXxJw7a9lMOHnadLYWKz+XeYOGWtra8XHMnb37Sfd+3r4t+2m2pWWzewkJT7OzysrK/t7/Z3RUBNT0dPVNTEyQkTcTCnW0teubN4Bafl7eh1279ezlJNDQNDdvULdeXdSJb964JhAIEhMSsrOy0FBo27BRo/TU1JSUR5jiv5vXCwqZD2JCmRUr6+YDBw9pY9eW3ZXfmpqZ6unqxcddhTmcpKWl8nVQej935lRtA4Mhn7t17tIVkWvr6Fg0tcxIT4MyTO7dvcPOiM6szEwsHOY3r12DSb169dEmIQJEgAgQASJABIgAEaheAlXJoV82girZFxYUXvn30k/bgk/+c6Jt+/aoxSp0U1JcUlpaLDNUWFQk08PtFhcX/3vpYsqjB1wPGqkpjy5dvICGalFoK29y+79b/suXxpw/Kz+kTg/K4e7jJ7W2tXv4IPnfmJiiIuaZFhgik8ZWRoqVr1RGk3aJABEgAkSACBABIkAEqkygBuTQ9eqbfzNzVpeuH2GRGhoajRo1LigsfPEiv7CoEPVj9NSvX8/AwBCjaamPS8tKrKxs0G5q2Wy0+3jzBg2Sku6hxqyrpwvNgZ9+1rdff6ihcJv5JDMq4jgSaNgWF5dATV+/1ulTqFMfz3meA/2szCwtbS3WW6PGjYyNjeFWoW1GenrdevXr1asHBctmVigVo6G+5DxDYb/YxqYlTODETOQHbVYMDQx1dLSjIo4hC69lYIBCO/qT798zNTNDhVtLS8u+XQd2xuT7SeiEByi0atOmoOBFdjZTO8cuCREgAtVAgFwQASJABIgAEZAQqAE5dEZ62rlTJx2dnGf7Lpjts/ADi6bHjxwqKyu7dOG80Lr5bJ8F/QZ+ylZn83Lzzpw82bO306x580eMGpPy6GF6WtqlmPO5eXnfec+FZtNmzeKuXuWrTZz6ra6u7rOn2VI1X+bDiKkpj59mZyGjdurTF96cP+n/NDsb0BTaXr186UHyfY9JU6FZz9w8Ly8PmuoLgseK7BwcsMDho0Y/y2Im4swfP350Pylp1Gh3OG9j51BWWoqh/27evHzx4uduw7xmz8t9/jwvN1fUeQOa4yZOhp/2HTpGR4bn51cuEjghIQJEgAgQASJABN4pArSYV0OgBuTQWPiFc2cDlvttC9q84cfVP67yf5icjM4b168F/uD3w9LFG9euWrc6MD72CjqvxcWu9v9ha9AG/+VLjolS7cKCwr27d64N9N/w45qNa1enpDxi1VatXB60cf1q/2UH/tyHjJxVW+W/fPP6dT8GrLyfdA9qsVeuBK5YBm9bN23YFryJm0LGFub7/whdG7By0/of9/y8c3vwppxnT2FeoWzfsvlk9AmoYbo1ASu3bFq/8cc1zPMqok70Q1jnCB7O168J2LRuLev89n+3An7wW7lsSWLCXS0dbVYz7MBfWD6WgAj/u3ULnSREgAgQASJABIgAESAC1U6gZuTQWDZSSRSk2fQRuyoEmplPMktKSvg6qMjK26LSjNSZr4ZddPJ74Afe4JPfiTbUoIwGJwqn4EYrbGAKTITpFGoieAg31LBh4y9HjBw7fuKgIf9zHfRpakoKN4qoEBunSY03ToACIAJEgAgQASJABN49AjUmh3730L/MilJSHm7esPbihbPZWZm/7PoJpfSX8Ua2RIAIEAEiQARkCNAuESACqglQDq2az9s7inpz7JUrJ6NPpD5OeXujpMiIABHpeEY/AAAQAElEQVQgAkSACBABIvAuEqAc+u08qhQVESACRIAIEAEiQASIwNtLgHLot/fYUGREgAgQgZpGgOIlAkSACLwvBCiHfl+ONK2TCBABIkAEiAARIAJEQBGBqvRRDl0VamRDBIgAESACRIAIEAEi8D4ToBz6fT76tHYi8HYQoCiIABEgAkSACNQ0ApRD17QjRvESASJABIgAESACbwMBiuH9JkA59Pt9/Gn1RIAIEAEiQASIABEgApUnQDl05ZmRxdtBgKIgAkSACBABIkAEiMCbIkA59JsiT/MSASJABIjA+0iA1kwEiMC7QYBy6HfjONIqiAARIAJEgAgQASJABF4fgfcth359ZGkmIkAEiAARIAJEgAgQgXeVAOXQ7+qRpXURASLwLhGgtRABIkAEiMDbRYBy6LfreFA0RIAIEAEiQASIABF4Vwi8y+ugHPpdPrq0NiJABIgAESACRIAIEIFXQYBy6FdBlXwSgbeDAEVBBIgAESACRIAIvBoClEO/Gq7klQgQASJABIgAEagaAbIiAjWBgGYza2sSIkAEiAARIAJEgAgQASJABNQnoHnv7l0SIsAnQG0iQASIABEgAkSACBAB1QToWY6acLeAYiQCRIAIEIGKCNA4ESACROB1EqAc+nXSprmIABEgAkSACBABIkAE3gUC1ZVDvwssaA1EgAgQASJABIgAESACREAdApRDq0OJdIgAEXhXCdC6iAARIAJEgAhUhQDl0FWhRjZEgAgQASJABIgAEXhzBGjmN0+Acug3fwwoAiJABIgAESACRIAIEIGaRYBy6Jp1vCjat4MARUEEiAARIALlCTi4DXUTiYN9+YGq79mxDt3c7Byq7oQsicCrIkA59KsiS36JABEgAm+awNAd0YnJ9xOTo+dSCiIQEA2B4BWekXaDp/gHBjAyo3U1TWM/aIrIYWDA9FbV5JLcEIFqJEA5dDXCJFdEgAgQgbeIgIPPZBehKB6h5zofO1Hr/d0Qjff32NPKicCrIUA59Kvhqsgr9REBIkAE3hSBhNvxb2rqt3BeovEWHhQKiQjUOAKUQ9e4Q0YBEwEiQATUIhDr941XeGJCYmJCkPfYELVMFCq9G51E4904jrQKIvD2EKAc+u05FhQJESACRKB6CcSHuDs5Ojo5+oVWr9+a6Y1o1MzjRlETgaoSeNV2lEO/asLknwgQASJABIgAESACROBdI0A59Lt2RGk9RODtIEBREAEiQASIABF4lwlQDv0uH11aGxEgAu8qAQd7O0aqc3kih3BbBZ+wEouaxi8xl5ozvIyaeC3qv5P4VS5HHIzy9VSooMxUbKj+MnmOJLa8rko1XyWxSgWiQJm6iIC6BCiHVpcU6REBIkAEXjcB+7k7oiOjGdnqJprbwQ09zCufw8IOMnIf7cgdPkMVv/5ZzlzkQ2Zj5+azNZrxI3IIt0xbuU+etYP9UITHvH8aVmJBPInR25W9jlrduRwQErPqyOjtSpbmtlWERamCgw/LLTKae6lfhTTkl8OgwHKUxCBQdzk8ZsqbsuFJnEvBRu5w47+gUF5BRajSeblTKFns+WAYs0wRyYq/G8XObXskc7ZIbHH0o6O3ulVsyAYgifl+pU821p62ROCtIkA59Ft1ON77YAgAESAC5QlYCYUiYXqRF4YFeIpf+cx0sD9Clwn+Ychj2L3yW5EtPJTv5faYrPFg4ARnK65H3GB9KkuFoWTnuz0yLMzfRci+gBo9UrFy6TNYPqmqzFyxB+4K2IW79FP45Rpu/ZzFS1OsYDe4L1bNCP81dmITBSELwDZZ6XIUxVCZ5UjRqGzxwrPzjZY/LkKXgIO4qBD5UKggsHLBmaD8qDExJyo6heBSyNiGyaTp6OeJ/dzo+wcDXYQyZ4uV0DmQMeRpKmwys8svCqoVnmzQISECbyMByqHfxqPydsbUuXNn2za2lYqtCiaV8k/KIPBxt25CK5n/1NBN8q4RaOETGTaByf4SwoODZnp7QYIiEhIly0QeU+kvIxy6g8kaWQ+J4UEinzO9g8ITE9g+oWeYYp9MAjeBy+UTYRvMxANbfkisE/G2knPFHTwuXppzf7YCL/bD/rJrYc02sFWo0MqaQYXRiMNqvNQPCTTLFgYJiRFivJVaDg4H9CtGhxkqlv7bDzKHWhKJ9IgIBMh0UY1Ghs0oCEBedNSCgsPFuAQCpd+nwz8EAmaZ4iMeDP+SmIRI03coAC4QIIEO8+T+oUngn4SMMQx/7MM0lP3wZ0+szMmmzCH1E4E3T6Bm5ND+AQG7d+8eOvSL1wlsno/Pnr17ZQSdiAFbfv9PO3eOHj3GwKA2hviyYMECqO3a9bPrgAH8/upte07wxCwIqXrdyngbOWrU6jVr1q1b17t3b5khZbtVMOG76t69x46fftqwaVNlE3e+E3XaOK9wduEcU60MwuAM2qrV+KP9+vfHuYFVyECbPHkyXK1Zu6ZJkw/4+lVoe82cuZr5s6raKfn5+SFIBI8DUYXAyKS6CThPQNKUGOHlKnR0X740JDQE4ufh6DjIK0iSPQk9ZyjMfhSHgjzY34UdSgx2beo01k/kMyR0qbuTo2swl0bL+3RjMzzGNjF85iALR9guZ+KBrR9CErrO3PgfM8r9VGGu+P/uis1d+g0Vt7hf9oP6iFNkpkuBgls/8dLCj6iRQg+dMUHsLiFokKOjx1KwFclS8XLW3GTm4X6qsBzOVp2Gs4uLAJFYSCIRHRHvcIkp0lwmXv5R81s+1lHoJdGwmjBd7kTgxSxIDMJZBOfiI74c/i1cvbnzyCVA/OCQZEL8tvNdzyXQcidh00FeTAqPqj80FQpvdn7YOGEqOtkUuqNOIvCWEKhKDq1pqK3dUJ8VtF/DSppZWjZv3rxevXoK53pFnR9YfGAj9wedmA5b/oidnd3UaVORbXTs1BGjrKAE2/Wjj6CG5MbZReX1OWtQ1W0D8waYBSGp6eD3338/d+7cxImTVOsj/kOHD0f/c5K9AEhISEhPz0hNTcvIyFBmWAUTZa7Qb1zHyMrKWtismUFtA+y+OsF5hbML55jqKUAYnEFbtRp/9Mb1a0aGRjg9XPBfomQA11q9eveCq/y8/IcPH0i6q/g7KfFeVlZW6uO0pHv3quhCkdmYMWN79e6NIHEIcCAUqVDf6yeQGDTVIyROZt74EL9vpNnPFOX38WXs3KYzeRjTGeHluDyWafB+4pY7zoxg911kfLptDRTnp4LwmU5jQxR8/WEsElB+nFWaK+SIOACBtY3M094Og/uICqIR4WzWKK9gIxoXCMKPqPFibHuJtiBig5/C5cSX41Ol5bAw1d0mBk+TiSQudDV3mBkvOBlkj1rIasmVj8CqhcyzNNKYkUA7LeUfHcabQBAXutSRS9Odp3APkYtGHXx+lJwtMFd0Ero7cRm8yKL8Rjp7JU+28m5ojwi8bQTUzaE1tAT6nc3MfGwbbO9Sf2PHuivbsoI2etCPUei8tuUhEfnwww+xlZ8Rne3at+f6kQeouNPdrm078wbmnLJ840LMhU4dkRuLZTIv9WSHHHv2CPD3T09LR7YxafJkzM46gZGZqWnq49T8/Hybljb29jL/pLFa4i3ybIh4R+4XlsO55QYRNoLndpU14BaibJTrh39lPFmdM6dPfzbk0xHDh8fFSf/1RWAq2Co0QdiYiPUps0U/RmU6Vewi5p49e3ImiATxyOtDAZ7l+2GuUJ/VhJU6eFllZdukpPs3b97AaKtWrS0tm6IB6e3k3KBBw7y8/Ev/XsYuK5gLM7Jtma3MCYyYsVJOJ3RfaL9PPpkwwTM3N4/rhCuFS2YVVMzFKsD/4E8/1dLSKi4uZnto+4oIVMptQtA3ClIfxkX80g2SdFOo6EFkRkfmx853ijPbFT7TQ3GlNuQIm6AKhNb8J5Ld+okNBeHeYxVbso65bVXnkgYgsyi7wX1FZePEu4fvimrwsquWKAjUepCDC1Qgn3ryxiTNqi5HYq/O7/ANsvkxrGIPHBffHMBO+EYFJ4P0ARihdWsocSKNWflZBGVpmm7VdxDvuoXjieq4spNQwMvg4Yov0tkre7LxvVCbCLyFBNTIoTUE+t3q1lvd3niajU5LI4G2huwytDXQj1HoQFMgNy6r/3L7H3fr9uuePZFRJzZs3Hgi+p/dv/yC9BYuN27edPHSpa3bth06fHTt2h9RPUUecPjIESiHhITs+Omnffv+4OqvyJ+WLFly5uzZrdu3hYUd+uPPvwZU9XEL5C579+7dsjU4JyenhU3L/gMGIhgIotLS1v7330vZ2Vl1zep26tQJnfIyZ87sE9HRu3b9DEHjOy8vVud3UbU4MHAVOrdu3RoREbls2TJ2CCnO9h07Dh78G8EfPXasbr26bL/MFq6ioqLgFnLq9Ol5Pj5snRjm2jo6HuM9MAVMFPIEvYDAVebm5gC1eMkSsEUPatKoTMMJOkEPPhEYnGAWzwmeUFBhgokGDBgAzqCNA4cDsWbtGqR66Ec8W7ZsQQ/6MfrX/gqOBYK5EBMD87//Dlu1evXePXtHjBiB44ujjHgOHTqMieAWgobCGRE/DIEU+sBrY2MDZU5ghRgQCfBijVgp9LnRyjbOnjuXl5dX39y8c5eurG2H9u1r166dnpYWc+EcPMO//HkImKB97Pjxvb/9hhPYe5Y3q4l4EDPHHA5xPwHo0IM2BMErXDKg4W8HoOGEwbpw8mzavJlL62HIF/exYz/4oMmNGzdQ4eb3U/uNEkg8fkBBiVQcEpduIuEtlzyJx2V/SZ+FUJFl3mQT1PKZ5dD+4iJ0YtBqNUq8mLjqc4Ue5rL4cosSP+uccOxgiDitFPYZzH9hhVhBkHi3/DMYiEaRxN2W5KbCCevnVvCKiaovR9HUivsS7zKX3nJj0jiV1delD8BY2fCASGNWeRYJBNI0HUSlZR8JT4FKc1545UKXzl7Zk62cG9p5qwhQMCyBCnJoTVNdM582xp7NNU10WYOyFyVFd58XXM1mBe2yghJ2CDrQhD6s2J5q3+J//W+++QZJT1LSvbBDh+4lJrZq1QppBDeRnZ3d06dZCQl3TczqeniMw236R48eRkZG1qtfr8kHTTg1L6+ZfT/55Hlu7vFjx1BYbdKk8ZQpU5EdcgpcQ1dXFwkNK/3690fCxw3xG3/+8WdGerq+vp61FVMgQWJqaWn54sWLuPi4u3fvwknnD7vw9dk2Ih8y5DMNDc1oXA388w8qf9gd+oX4STak4F0/6no/6f6dO3c0NDW79+jJDk2bNg1VbZS3YYUsp3NnBZ4R6qCBA/X09cPDw/fs2YNIXF1du3XrHhERkZ6ejvrihZgLSB8tLZsq5Hk/OTkqKjLn2bPCwsITJ06cO3uODZjdurl96eziAidY9d9//62ppT3sy2ENmzRWYQK2IGxh8QFogDnII5gZM76DQySI7Tt0SM9IP3L48LVr1xo3bjJunAcCw5Ay0dTU7NLlw7uJCUn3kozq1JkydWqTJk1Oh2kUywAAEABJREFUnz717OlT5OWf/e9/MOTPGBUVhUQWM86ZMxdDU6dO69atW2lZGbLPuwl3AVlTSwv9ENaqYcNGsXGxiBPccJ7gbMFQ1SQqMiI19XGtWvrt27WFB6TCdvbM/203b95ISroPz/APGpgrLi6uSfnz0MTEpFGjRgkJCQ+SHyhkDlv45IQNnoUss2RWB8XpzMzMG9dvlJWWtm/fHheZbD9/O3ToFz169nzy5Mmff/5ZVlbGH6L2GyWQ8J/0DpB8IFy+KyiXPMkrsj2trcXPOgis+m/fukOxTJc8c8wrakofe1AdDzuNaFvluQSCEMnjHOWeeBY/6yzK5yR5mxW/Vi5WECQcO1juGQxROIo2odLHJISegWGJ0duVZ9IvsRxFUyvsq5CtkiRboTN0SmOuyLOEp4B/MabuQZeehJhTKtLZK3mySV1Qiwi8pQRU5dBajWuZ+bTWaWEo0BAISssK459mLopPm3Axc/G17MBbrKCd5nkR/RiFDjShDyvYvooVJyXd//77xUuXLp09e87C+fNRSiwuKmrSuAlucLPTxV+LHz5s2NgxY4wMDc3rmyNfXLpk6exZszas35Cfl8fq2Lax7dylc0FBYXBQ8Ny5c6dNnXL7zm1Uc5FXsQr8LW58oxDLyqKFCz/p+wl/lN9Of5KB3A7pMjq7dO5cp04dVKDj4+Nv3riFIJtbN0dijSG+/PzzzgULFy5auMDL67vvZsx4nJKCZKtFC2lZ9OiRI19/PWqc+1hcM7BDyJNw2VBSUrJvXyisMJR4T3Q3k+9XIEAipV9LPysza19oaGBAwNw5c6Z/++2aNavRfp6TIygri70Su2H9emU8kb6H/f13/osXRUXFURERu3bu5Ltv1Lixjo7Oo0ePUH1ftHChz7x5U6dN3bZlqwoTJyen+ub1kUDPnuUN5rt27cL1T15uLtz6r/RHLXbe3Hm+vr6///Zbfv4LHIvWtm0wpEIuXDjv4T5u5coVyPZKS0u3bt0y/dvpoaGhQG1ibAxD/ozeM2cGBwcjIbazt+/Tpw9uEeBIRUVGTp0yZdxY95iYC9iFCYS1QjqLfsT5y+7dQN2+Q0dlOT1SdiSy7CUWts7OzsiS4YeT3Ny8fy9fRjLaSvQ4h+RBjryz585VeB4WFBRs2rjxSze3ZX5+PObBHPNjR49yE6HBBs9C5i+5Z8+eGIVcvnIFfztGjRr575XL2traLVvyb9FjXIBl/u/z/+nr68NzyqNHTBf9vOMEhC4uuCJWLJI8m4eAy4cS76pV4uWZIi2r3FywvXFXXCF26ScuLQgE4odJkEIzFxWSWrW8guqiKZzzJNbPyZX3tLGVC5NJJ0erfNGboJLoeNO9sWbFR01RHvxSB52/1hpIjB8+tYmAHAGlObRmHR2TKc21zPVhUppdmL3mv6yVN4sScgXylakyAfoxCh1oQh9WsIUHtKtdkGM52NsvW7YsPCJi9JivtXV0mMRdW5udKCnxHrIWtOvXr4ehlMcpMTEx2D1y+DD3YTihZbM6dUz09HRRgo3+5+Shw0etrayQUhgZGUFTRu7dYwreqHlDDh8+fOPGdRkFbreOUR2kSiWlpeixb+sAhw8fPMJcmVlPsrKzkVJ37NgBQ3xBqAYGtYcPH77/wP6Tp041Ewo1NDS0NMUHpaS4OC0tHfpQy3jyhB0yr1/fwMDwWU7O1StX2aFbNxT8X3bl8pXMJ1kNGjZYt379oUOHR40aZWhoCH15Uc1TXh89uDBAIt6iRYv9+HNgv+tA11q1aqFfhVgKmyFPTUxIQNYONSTlX3zxBVJntHNznzeztJw9Z86x48eRRwOIhoamlgau2zCoVJ5kPMFYqfhPWW7Oc9Gu9OyUmfHu3bs5OTlGhgbdunczNDRCPn01lgEIK/YiBw0Ia2Vra4sTA+I+bhyuFnCqmJs3wKi8dO7YycfHl73EwtZr5kxbuewfxyL3+XP2cY42trZglfwgGfVpnBuqz0MEjLDZSdVhzgbPQYYtPGDJLVu1ZJ2kPBSnxWxDU0t8prGj2A4bNtza2vr2f7eDg4OwS/I+EEhITKxQqotDhRNBodxc0gd8rSQfkhvaX/QwCVdjltSqnSWvwJO89i7x+AEmyS7nT8UO0mgLV+/wRF5JQih0CTiYHK30C0QQbYWiYsb3b0hQIS4ovIdYaMk1l4Dsf6LsSjR0NOuMt9K2YF7WVvIoP9PvesGVbHZIxRY60IQ+dGALD/CDdjWKvb392rVrcA+6fv36yG5jYi6i7qiOfyYz0yy32OLi4sePUUt9+OjRw3v3km7fvp2alirvKi09DQVvVhYvXhwdHS2vg56Pu3UzNzcvLCx88OCBi4sL0iN0durcCXnVrFmzES1Saq5YjiFWZqOcPmeug4MDyr1nzpxJS01j+19+iyuHOXNmI+9PSXlsbGKM8PyWLfP0nCDjuWo8D+zfP2/evKioqKzMbKwald0A/wAcFBnn8rsFhYUynTgugQGrRn399QcffPA45fHp06eR3crovMyu/IzqeMvOysZZAcHRvHPnTuK9e7l5TMlc3vZ+cvLRo0fAmZWIiIg0ubPoyOHDOMlq1dK3bd2qhSidRUKM6yLWm5rnofrMq7ZkBIPz09Gxl4aGhlEdw63btnvP8jY2Nsb1w9hx4zwneEKB5C0m0EryOmQB/ytFKg44MXiao5NjRTI2pGJPFWtUZa74A8fYpFbYh33i2Z59hwYS5HjxjCHizz6Kn/ewF7/2LuGYmg9yiN0wv+JCxzo6WbgO8gqKENe/0St0DgyTf9ebQFCV5cDdGxUh/5EXhZFITySFwy/VWROJvdSCyfjdJ1AureSWW6uXua6dCXZRV87ecKcktQBtdQSa0IcVlOEBftCoRun6UVfkow8ePBznPtZj3LibN24oc56enoH02sKiKftq3v79XevVE78aLy09HbXPkpLSkJBQ3NqGbN+2bcP6DcFBwcq8qe7HDX0Uek1NTVNT006ePNmxY0dDI6PMrCwkT2xqdfrUKQSDxBrpNd9V+w4dULreu3fv0M8//37xovwX+fxRhW02+DqGhmy9Ezloy9ayd+RhaNvGtlOnTnFXrn425NO+fVzi4uJ0dXUd2vI+aQ0lgUB9niJ18aZz584ODm2PHjkyYED/ESO+evAgGeu1s7MTDyv69fhRClba3KY5Asa425due3/7DRcYrgMHWjS1AKsF8+d//fWoU6dOlZZKa8nQrLLIzGjZTIjyc35efsyFmOfPc/T09FiA8N+qdUtt3M1ASyBgrR6lPMJZAcGlwp49ewL8/a9fU3z/IS4ubpmfH3uJhW1gQABbaBc5k24uxsSUlJTgiqVRw0aoSaMyjTH2UKp5HqrDnA2eg8wtGWVGTFeh6OBmjoZAU1Pzgw8sbGxsrKysQQnXfs2aNavUG/0E9Of1ExBnlphYvSdluWckhBUmVfDJE6lhHwXfRMhTlDalJpWcS+QiVvypQYGVyNpB/FY7/nO9kscPrEWvwGstftS7ctcSornEm7j4ED8Px6bIpNn0Hd3Ogdsl76h+ueXA1xsQLmYBV85XEoXCE4kzr+BsacVdyJXzrq55OSPaIQI1goCmfJSadXRq9WmgoSUoKyrN2Z1UnCx+jFheU2EP9GEFW3iAH3hTqFa1zpSUx4WFRaamJkO/cBvvOX7gwEFa+I9fka+jx46ifGhmarpw0aLff/99ytQp+np6rCJqtFevXkVdENW1BQsWrFi50nf+guU/LB82bBirwN/aNLfZw/umlXk+PtwoO7Tvjz//+ms/Upznz5//unt3Wupj+7YOSETiY+N8fX2RV0G2bt2a/iQDiWanLp05czQyszK1tLRQ//ty2JeoEzdp0hidqgXB37x5E6seNnxYYOCqoOAtza2by5s4OztPmDhx6jfTfHx8UX62sLBADpcouk1ZUloK8+49unt7z1LBs6i4WFAm0NfXG/zpkEmTJvOnGDJkiMd4D+9Zs7755pvRo0ebmtYtKCi4f/++ChNcWmRlZ4PYuvUbli9f7jHe09ra2sjQEOldXl6eoYEBitlfffUVLkUwI3+uKrdlZsSxxhG/fuNGWFjYpYuXkND37dt3/YYNW7Zs6dLlQ+yyE7FWSHa3btuGGwgrVvwwf/58H18fNvVndaqwvRAT8+zZM4umTXEhh5o0Lq7gBIdS/fNQGXP44YQNnoPMLTlc/BJdTlFxA/EM6N+fe5PjpIkT09LSUC/HCcw+daPYjHpfEwHuWQUF80kySwzxk0vsKhHpR8dUuVVgLDUU9mELwwqUyndJTSo5F+uGe5zDpZ+bQPKStXLfnCKpVSMie8nT0gIVr4Bg/Va4RSbt5BokSaPZBB1GL7kceHj9Io0ZjFQVOwStxVcgAhT6uSdhpOYqj6Cb5HttZBaoprmMFe0SgZpAQEEOrdfORKu+PoIvjHv6IiYTDXWEeXhjdDP24Q1YwRZW8ANvaFSXREVGnD17Vl9ff/jw4ePGeTxOfZyf/0Kh88SEhIAAf5QJMYrE5dTJk8kPpN9nsXr1qtOnT6GaO/jTT5FulpWV/vnnn6gHQ1lGTE1NUZPj5AOLDzgFdsjSsmlpaSkmmjd3bui+0B49HZs0aoyc8srVK5wmRu/euYPEGsVaAwPmCRl26O8DB9PT0tGJdLZ1a9uEBMm/1+ywku26devi4+MNDAwcezk2btTowvnz8oqorB8/dgy158/+99lXI0ciST19+vTmTRuhGRkZkZub26pVK1ToMzIylPG8cvnyqdOnkHl37NQRE8GQE+SXSApxo//r0aMHDhyopaV56NCh3T//rMIkKioqeHNQZmYm0tM+ffvWMTKCBywEeVtUZFRZWVn/AQOmz5hRVFSUnZ3NTfQyDYUzBgYGwuf69evOnDmDw9G1a1fbNm3Onz9fggsGDAgErFVWVlbbtm1RLBcKrW7durVp40akkqLxKm7OnD6dlJSEmi7s42KlrwpQ/zxUxhwOOWGDl4HMLpnToUbNJSB+VkHBAobyvmZvjXqPXUg+iicQyH6FigL//C6podUE+S/D42tybalJJediPXDva3Pu7yZ+TkPmm1MktWokiENbWIusyiXZop4qbWK5t95JzV9yOVJHr7EVyr14ROVRs/OVvDU8odyTMDclbzlUcbZIbeXWVROJyS2COoiAIgJyObSGQL9rXZSQy4pK80+moxKpyEq2Dwm0yXctazk1MJ5mw6TRZQLYlhWVwg+8CSr4eJisN/n9L774AunO5s2bkMrMmT0bCTRKvO7uY93Hju3ZozuKZ0jFJk+chBIav2B288aNtWvX9nJ0RLL4ww/LNTSkcaSlpk3/dvrgTwfjZj0a0Fklyq74U7MO4ZMv6IQOtvzO7t26jR0zBskohlBixHTdPv54185y77LALDAZOeIrLAFqrCD1HDCgP4Zmz5r1v8+GfDViBHTYJXBLZjXZGdkhXB5gutGjv4bVoEEDUQyGFRRYTXaLWVBA7ePiDIJhg84AABAASURBVB2Is7PTdzNmoBOjwUHBWO+UyZPHjB197uxZZTyhuczPDx6gOc59bNihQ449e7CoEcCUSZMGDx7E0sPyoQl9CBoKTTAUGvp7/379Jk+ahGOHIwgP8IN+f/+Vn302BK4wNGL48L59+mAiTAdBg50RapxgpVgviwLHHQpQgzIUcIbgPAE6tCHKZgSHGdOnIwasHVOjrczKzc1t5FdfoW4Nb/x5sVtZGe/hgbC7dO68YsUKzlbZeYjlYFFYGhbIKoMViMkzV3PJMsGDHoJBJ+tc4RZTIwCEgWAUKlDn6ybg4h/tI19EtHPb7i/6lB3CSVT1DmmM80T6jRhCz7DtQ2Uf8xJrwvlc7oUYbJ/UUOAcGD1X4duUHdyG8vulJpWcSzzjEfE3yFhNEX894WGZC4U4capr1Xcy+z4+mSSb9aN0az/U100eLKPuwD59jebd29y170suB85ev8T6bRS/a1t01BQdbhzrgxOYV7MiuogN5b4lMZ73JT7M2QKN8sK3LT8i2quJxESB0+atI/DRRx9t3hy0Z+/eDZs22baxlY/PvIH5ou+/P3b8ePQ/J/ft+2P06DGsjn9AAKxkRJkT1kSdrWwOrWWmp9WIecdCaVZh0V3Fn6OS8csm0Fp1dZEra5vraZrpQAG28IAGvGmZiR+iwG61CPIJpKrKHlFlp+jevce+P/7cvGlTQGDgl8O+XLP2RwsLi2c5ObFx3L+EAmQwx44ePXXqJGvyprYIICICteHKPTOD5VdohWQROhA0ZFaH4iuWz3aq4AlDaGLLavK3MFdID8rKTGCOzAzHDjOizQnrCkNcTzU24FZ+RvhHDOHh4ZgabXlRZiWv+fI9iEEhSXnPamq+zuDlg3wbet7VGKwmHEyO3oqEz8HeDuLmszU6+mCgJINOCPpGwXfXKWMRt3yDJKsSuPiHRUfu8BkKn2Jxm+u7PTL5PpyzdV2el7jljjPFSa1A6BkYFil6mzITD2zd3ObuiE4MC5jcgmchqPJcrBPJpwathKIUT0GNWVLpxJ0jxqSyD3LY9GHev1GeABLr7ZFhkqSyXFL+ksthInz9P6FjeUct7D6zWDfRWcQctXInUmKQq4fMRYogZA33VAvOlnInIXPEcZ4IBIkR4crupNZIYq//GNGMFRD4zstr+Q8/dOrcycbGRtismUFtAxkDoZXV2rU/DnR1NTMzMzCobdnMctLkSYsWL4ZaM0tLWMmIQidQVl9kc2jNujqatbVgX5JRUJpdiIZqkSbQAkHJo/ys1bdKUpkPIMIWHmALb/CJxmsWJKZ//LEvPz+/V69e3t6z2rVrh3vcKMHirvprjoSmIwJEgAi8NIEIL9dg5mURQucJAQfDwhgJnOBsJcoq4TwhaJBjudoh+iqQEPdBXlzWIxS6TPBn3TLbAM8JLhLX8m5CPFylKZVQ9DZlJh4YBgZ4KrSr+lzM7JIUmWkr/oq+EEmtmlFJrMK7qwUCGQJh/hwBsJV5M8nLLYeJ8Q384KjN5N43whzuQNFZxBw16YmEBNpJ0ZVY/FJHb+5kEfBPQvaIJ0Z4Oa65q3xVNZKY8uW8uhHyrIxA586dXVxcdHV1795VeqKhYGptbV1YWPjrr79Omzr1ytUrWlpaH+NPt27IrRfOn8/KooULE0UfD3v27FlamoIXsimLQb5fLoeupS3QYp55KEmRPmesWUfHZHoLXTtjGXtlCTSrJvagpaEJn2zX690iY+7du/ekiRNBbfTor/v364eb+683BJqNCBABIlBNBFAAZl5gLOctMSJoZqUTaJGX+BB3J1fkVUrKhwnhwV7yJUmRpYK3KYv6RZvE8KCNB7hPpIm6BIKqzwUHvBRZyYtHbki+jUWg/tcTwrFI4g4el6aHoh5ug9TQVajo4uSllsO5f82N2BAPR1fvICWLTQj3dm2qMIFmwwwd68i76GL7RFvG0FGudC0a4m1qJDFe/NR88wTuJd1bsWJFfJzsPy5sZCg8d2jfXlNT88KF86sCA8+ePeu31G/+/PkBAQF37txGXTXs0CFW6tWr37hx48LCwqNHjiQl3WfNq7aVzaEFOpoCJoWWetM01TWebqPXwdR4YnN+Gq06gZbawxt8Svcr3dq5+9eXkanfznAbMXL2XN+XcUK27wGBlzrN3nk+lf57SwavggDzAmMh8wLjmd5ejAxybSq0cPRYGiJ5WbLMpEi7ocCI0hSHyascGZ+ujEOJW1ehRVOho/vyEMX/W4mmYYJxgpqrK2vFbF0ZQ6exfqHSZ+ZEuuym6nOFeGAikShJ8qQrVZjyiuaX6sjQiF/qzi5kkIgqsxAvVzFbFQSqvhxROLIbpeFxiqFjmUOJQ6MEgkAQ4o5RRhTl/SI/caGSxYqWKTrorq6DwNbRXfFRE5mxGyYPhqb0iJczRK2amRoKMmV71hjbShOrmAm8kryDBFB1dh0wgBPbNrYxMTGTJ07a/9dfylbbsWMnE1Oz4uLiZ89y/vjzr5iLF/fu3fvViK+ePMlIS5V++QZSbZe+Lnp6egmJCSEhvynzpma/pqxeUans5wgLSwWFzCt7Netoc2m0ugk0vMMUPtGoqoweOYKECBCBN0ugqn99ye4VEIiLDwkJFUm8wlS1KlPGxceKfcJzfKyK1FnOe2wcTMSiluFLzCU3eXV2xErBhoYgSDV9Q7Oq6NScoSK1qoyXP2pKrsGUOJbaxlXOUOyvZhITB0+/XguBsePcFy9Zwsnn//u8wmmNDA31dPVQh3ZxcalTx+j+/eSy0tLWtq29vGZaWjblzN3cvrQSWqEIHXE8PDe3cp9D45xwDdkcujS/WFCCtFeg1Yh5vR30SnOLn667XXj9GdpsGl27bwOT71pq1dVFD/8ZaOzyReyhpIzxyR+gNhEgAkSACBABIkAEiAARUETg3Nlz7HMX7PbylSuKtBT0aWpqPniQPH78+M//99lvv/2GsrSlpaVjr96sKorQzn2Yh6qrpQgNn3I59JOi0rwSDGjV09M0YbJktGXSaKORzSpMoGELD4xtXknpkyI03k2hVREBIkAEiAARIAJEgAhUH4FdO3cunC/+CCAafx88WKHvjCdP8vOZuvLNm7cSE5gPYP9361ZBQaGWpmbtWuLv5ajeIjRCks2hSzILSlLymQFTXR1r6XtD+Gk0RiEqKtAYha2mKZOCwxt8ooeECBABIkAE3hYCFAcRIAJEoOYT6N27t7u7u3kD8/Pnz2dkZGBBrVq1FFpZodGiZUs9Pd2i4uLMrCfYrfYiNHzK5tCCMsGLc0/KSgQaOpq1etQXaEBHLPw0WnUCDSvYwgP8wBt8il3QLyJABIgAESACRIAIEAEiUDUCPKuPu3WbO2/epMmTvb1noTv8eHhBQYFQaLVly5Y//vzry2HDtLW17yXeC/v7b4xWexEaPuVyaIGg4Ep2STrzYjtde2P9zmZQ4oRNo/OOPubeA80N8Ruwgi164Afe0CAhAkSACBABIkAEiAARIALVRSDn2bMXL/JLS0vRgM+fftoRHBScmZlpbGzctKkFes5fuLBy5Yrc3DyuCH33zt2Xfx0HPLOiIIcufVaUfzwVJWQUko1GWmpbiJ8jYQ2QRuf8ksR+kQrbI7OFPqxgCw/wA28yCrRLBIhADSVAYRMBIkAEiAAReIMElixZ0qljxwH9+8fExMTFxXl4eEwY7/n999+zIe3c+dMnffuOHv317FmznJ16T5k0CToYQho9csRXMBw1aiTa6KkWUZBDw2/+ibTC+Gw0NE10TaY012qgh7Y6Ak3owwrK8AA/aJAQASJABIhAVQjQ+3GrQo1siIAsAdp/VwmkpaZduSr71o7r165HRERUY66sjJ7iHLqsqPTZloTiZOYTjlqNa5n52Oq1M1HmguuHDjShjx7YwgP8oE1CBIgAESACRIAIEAEiQATeJQKKc2issPRZUfaGOyVpzIPRqCubTG9hOquVjpWBQAOD5UVDgH6MmkxvoSl6HR6sYAsP5fVedk9XT09Xj3nXx8s6IvvXSoAmIwJEgAgQASJABIjAu0ZAaQ6NhZY8ys/0u1H033NBmUCgqaFrZ2y2yM48uJPZwjYmXi1ZQds8qBP6MQodaEIfVrCFh+qSJhYW076b+Z33nBne8zwnT2vUqHF1eVbfz8gx4+wc2lWo36hRE/fxE43qGKvWrG1Qe/io0fO/9+vh2EtG87OhbvKdMjq0SwSIABEgAq+YALknAkSACKgioCqHhl1pVmGm37WnwXdKswuxC9HQ09KxNtRra8IK2hr6WuiHQAea0IcVdqtLdHR0P+nnGh8b+8PSxQHLlzx5ktFdLu+srrlem5/P3YYVFhamPHr42makiYgAESACRIAIEAEiQASqi0AFOTQzTZngxeknGTMuP113u+hWjqAYRWmmW/pTXIZ+jEIHmihFS4eq3OIZamtr6+jqFBUySXxJScm+3/b8vvdXjKPcy9ZrUfSdMGUaWyRu0bLltzNnz/ZdMMtnwcc9ekJNV08XCets30Xo9Jg42djEFJ2c2hzfhXJqizi1NvYOXnPmwRD+9XTFj5HI22poaLgOHjLHdxEmderTR6ChgSlUy6GDB7CQ0pJSVo0ttCOYb728jY1lHz3v0vWjmXN8Z82bD/9dP+4OEw3ejF99PcbdcxIgoLNvv/5wgoCnzfCCT2iSEAEiQASIABEgAkSACFQ7ATVyaNGcZSWCFzGZmX7XU90vpE++9GTWVVbQRg/6MQodkW41b/Lz8y5duPBx9x7jJ0/t7dzH0MhI2QS1DWr3dul79d9LK5Z+v3f3LmsbGxNjk46dP2zQsOHm9WtX+i15nJLSrn0HqDk6u5w99Q/U9oX81q5DxwYNG0HN1Mxs3eqVKHU/z8np4dgLWWnP3k5nT52E2skTUSZmZphXoW0b+7ZWzZvv3Bbsv2zJ06dPa9cu9zZALS2t1m3sYAhzTp6Ivk2H3UXu26ffgMS7d1f4fR/y6y9mdeuy/ewW5vq1aoUd/GvlsiVREcfs27WtVat2i1athVbW7IzpGensjOhsadvm1593IuD4uLhP+rmihM86oS0RIAKviQBNQwSIABEgAu8HAXVzaD6N0ufFxY9fsII2f+gVtS/GnF8TuCLm3FnLZsLJ06aztVj5ucwbNNTS1o6Pi8XQ/aR7P2/flv0029Ky2b2EhKfZWWVlZX/v/zM6KgJqerr6JiYmyMibCYU62tr1zRtALT8v78Ou3Xr2chJoaJqbN6hbry7qxDdvXBMIBIkJCdlZWWgotG3YqFF6ampKyiNM8d/N6wWFzAcxocyKlXXzgYOHtLFry+7Kb03NTPV09eLjrsIcTtLSUvk6KL2fO3OqtoHBkM/dOnfpisi1dXQsmlpmpKdBGSb37t5hZ0RnVmYmFg7zm9euwaRevfpokxABIkAEiAARIAJEQBUBGqs8gark0JWfpRosCgsKr/x76adtwSf/OdG2fXvUYhU6LSkuKS0tlhkqLCqS6eF2i4v+vzaWAAAEAUlEQVSL/710MeXRA64HjdSUR5cuXkBDtSi0lTe5/d8t/+VLY86flR9SpwflcPfxk1rb2j18kPxvTExREfNMCwyRSWMrI8XKVyqjSbtEgAgQASJABIgAESACVSZQA3LoevXNv5k5q0vXj7BIDQ2NRo0aFxQWvniRX1hUiPoxeurXr2dgYIjRtNTHpWUlVlY2aDe1bDbafbx5gwZJSfdQY9bV04XmwE8/69uvP9RQuM18khkVcRwJNGyLi0ugpq9f6/Qp1KmP5zzPgX5WZpaWthbrrVHjRsbGxnCr0DYjPb1uvfr16tWDgmUzK5SK0VBfcp6hsF9sY9MSJnBiJvKDNiuGBoY6OtpREceQhdcyMEChHf3J9++Zmpmhwq2lpWXfrgM7Y/L9JHTCAxRatWlTUPAiO5upnWOX5C0nQOERASJABIgAESACNYtADcihM9LTzp066ejkPNt3wWyfhR9YND1+5FBZWdmlC+eF1s1n+yzoN/BTtjqbl5t35uTJnr2dZs2bP2LUmJRHD9PT0i7FnM/Ny/vOey40mzZrFnf1Kl9t4tRvdXV1nz3Nlqr5Mh9GTE15/DQ7Cxm1U5++8Ob8Sf+n2dk4tAptr16+9CD5vsekqdCsZ26el5cHTfUFwWNFdg4OWODwUaOfZTETceaPHz+6n5Q0arQ7nLexcygrLcXQfzdvXr548XO3YV6z5+U+f56XmyvqvAHNcRMnw0/7Dh2jI8Pz8ysXCZyQEAEiQASIABFQkwCpEYH3mUANyKFxeC6cOxuw3G9b0OYNP67+cZX/w+RkdN64fi3wB78fli7euHbVutWB8bHMlz1ei4td7f/D1qAN/suXHBOl2oUFhXt371wb6L/hxzUb165OSXkEW6itWrk8aOP61f7LDvy5Dxk5q7bKf/nm9et+DFh5P+ke1GKvXAlcsQzetm7asC14EzeFjC3M9/8RujZg5ab1P+75eef24E05z57CvELZvmXzyegTUMN0awJWbtm0fuOPa5jnVUSd6IewzhE8nK9fE7Bp3VrW+e3/bgX84Ldy2ZLEhLtaOtqsZtiBv7B8LAER/nfrFjpJiAARIAJEgAgQASJABKqdQM3IobFspJIoSLPpI3ZVCDQzn2SWlJTwdVCRlbdFpRmpM18Nu+jk98APvMEnvxNtqEEZDU4UTsGNVtjAFJgI0ynURPAQbqhhw8Zfjhg5dvzEQUP+5zro09SUFG4UUSE2TpMaRIAIEAEiQASIABEgAtVOoMbk0NW+8hrtMCXl4eYNay9eOJudlfnLrp9QSq/Ry6HgiQAReJ8I0FqJABEgAu8CAcqha+pRRL059sqVk9EnUh+n1NQ1UNxEgAgQASJABIgAEagZBGSj/D8AAAD//zUnxLEAAAAGSURBVAMAjUSwBUpteqMAAAAASUVORK5CYII="}}},{"cell_type":"code","source":"from scipy.signal import butter, filtfilt\n\nleads = ['I', 'II', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']\ntemplate_len = 500\nlead_templates = {}\n\nfor lead in leads:\n    signals = []\n    for _, row in train.iterrows():\n        csv_path = os.path.join(TRAIN_DIR, str(row['id']), f\"{row['id']}.csv\")\n        \n        if not os.path.exists(csv_path):\n            continue\n        \n        try:\n            df = pd.read_csv(csv_path)\n            if lead not in df.columns:\n                continue\n            \n            s = df[lead].dropna().values.astype(np.float32)\n            if len(s) < 50:\n                continue\n            \n            s_norm = (s - s.mean()) / (s.std() + 1e-8)\n            s_resamp = np.interp(\n                np.linspace(0, 1, template_len),\n                np.linspace(0, 1, len(s_norm)),\n                s_norm\n            )\n            signals.append(s_resamp)\n        except:\n            continue\n    \n    if signals:\n        lead_templates[lead] = np.mean(signals, axis=0)\n    else:\n        t = np.linspace(0, 1, template_len)\n        lead_templates[lead] = np.sin(2 * np.pi * t)\n\npredictions = {}\nmin_val, max_val = 0.0, 0.1\n\nfor _, row in test.iterrows():\n    base_id = row['id']\n    lead = row['lead']\n    n_rows = row['number_of_rows']\n    fs = row.get('fs', 500)\n    \n    template = lead_templates.get(lead, lead_templates['II']).copy()\n    \n    if len(template) != n_rows:\n        signal = np.interp(\n            np.linspace(0, 1, n_rows),\n            np.linspace(0, 1, len(template)),\n            template\n        )\n    else:\n        signal = template\n    \n    if len(signal) > 10:\n        nyq = 0.5 * fs\n        normal_cutoff = min(15.0 / nyq, 0.99)\n        b, a = butter(2, normal_cutoff, btype='low')\n        signal = filtfilt(b, a, signal)\n    \n    s_min, s_max = signal.min(), signal.max()\n    \n    if s_max - s_min < 1e-8:\n        signal = np.full(n_rows, (min_val + max_val) / 2)\n    else:\n        signal = (signal - s_min) / (s_max - s_min)\n        signal = min_val + signal * (max_val - min_val)\n    \n    predictions[(base_id, lead)] = signal.astype(np.float32)\n\nsubmission_data = []\nfor _, row in test.iterrows():\n    base_id = row['id']\n    lead = row['lead']\n    n_rows = row['number_of_rows']\n    signal = predictions[(base_id, lead)]\n    \n    for i in range(n_rows):\n        submission_data.append({\n            'id': f\"{base_id}_{i}_{lead}\",\n            'value': float(signal[i])\n        })\n\nsubmission = pd.DataFrame(submission_data)\nsubmission.to_csv('submission_baseline.csv', index=False)\nsubmission.head(30)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-24T10:45:35.320086Z","iopub.execute_input":"2025-10-24T10:45:35.320414Z","iopub.status.idle":"2025-10-24T10:46:49.272362Z","shell.execute_reply.started":"2025-10-24T10:45:35.320390Z","shell.execute_reply":"2025-10-24T10:46:49.271452Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"![image.png](attachment:5192e778-91d3-4edf-b55a-4a833b2314ae.png)","metadata":{},"attachments":{"5192e778-91d3-4edf-b55a-4a833b2314ae.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAA7IAAAA9CAYAAACQjlVkAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAACDmSURBVHhe7d0PWFRloj/wr4qI6BAKycqIOak7aklWTGVj3sD+YD6h/YHrKtUV3E25rWIb8LSC3gTtQdtVuq3abtjWotui25r9TMoCfyZrOXRTTJObNiYyhoEaY4gj2n3fc87ADAwwg3+S9ft5nnk8c/7NOe85x+d8ed/znm6DIob8CCIiIiIiIqIuorv2LxEREREREVGXwCBLREREREREXQqDLBEREREREXUpDLJERERERETUpTDIEhERERERUZfCIEtERERERERdCoMsERERERERdSkMskRERERERNSlMMgSERERERFRl8IgS0RERERERF3K5Q2y1/VA98E9mz7drtPGE10FwqNiYB6lfemkS7EOujoZx8fAFKF9ISIiIqKrSrdBEUN+1IYvge7obg5B4KQwBA7sg249tNGuzp/HuWPf4szmajSUnsMl/HGvJb20ARMH21H+xgxkb9ZGXiViM19HcqRO++bOXv4aknKKPM9ztgblW/OR/6YFNm2Uq8kLX0fiTXKZGpT+bhZWfKqOv+JmLsP6BwxN+/KTmbQM7y+KQYijAgVjp2GFNtonl2IdrSRjaeFEGOx7kJ+cjZ+whMT+ZWHNU7dAd2QL4p/L10a2z3luWj94HOmvaSN9NWYe8n5rRrg4k0t+MwcrK7XxmvBnXkbe+HDg21Jkz1mOcm38JTV3Hf75pBH+1cXIfCjtsh6H8Nh5yEiaiMiBOvQ6b0dV+RbkL1uOImW/Y5GVPxNt/JdwceVMRERE1IVdshrZbvfqEfxHE65/+kb0GdRGiJV69EDPQXoEPX0bBvxxJALv7alNuHL0BgMMQw0ID9ZGXEVCI9Rt8/iJCG17nlEmTJ67Ghs2LkNCq1qkRMRGR2rzmjBxaqw2/ieg17ZX2xdfZf39M5R9tg2rpmsjvJKIVdvkchuQpY1B+R58Xe2A3WaFVRvVvkuxDm+Ea8dJj86V0CUUrFe3xSBCo5eazk29NqIzdu+A3V+ux4zYWSZtpJMeydFm5Td6ndlxeUKs9PmXqLKLY3vMevl+Q5q+GmsXJ8I8NAS95PdeITCYE5Hz59VIVK7jUIQr50Mbn4spZyIiIqIu7OKDbGAfBC4RoTRpEHoFaOOEHxvO4dzR73Dmf10+R3/AhXPaDFJAEHRJt+H6JQPhF6iN+0npYRpvElGiPWKeKKM2rBEh0psmiOFRJrRYsk32suWIuv12t0/00wXaVFXzPHGYnVOI8lrAPyIGqdkp7vswfRxGyBqdw1altjbkpocwWZnghVFmr5rOGsW+tV1uRlFmPt5xe/m77rw5fprKAsx+aCyiH5uPd7RRTuFRZu+alLazDrnP5vHtH23jeLPX54NnYn/va7GOCBPM7ZZ1R9vl4fz2qBPH1GsWvLOvShkKv2kK3KJsRBIih8iBKpS/bVFGOXl1fXm4Vj0e7+3ZiL9XHNvklR5aOHR8bCVvtid1kgny0rRtnoO77xmLu8fOUWtig02Y/JScQ5xj97r/PxA1JR/l9XJaLWyfy3+JiIiIrj0X17Q4TATRBSMR2NTs7QIaDx6FfV01HAcvaONa6zYsCL2nDUHfYb3RTRsHezVOLjoMR7X2/TKStXqTh9hh+f29mL1WjIhIxNIVMxEzpLn9nv1QEVb9Zj4KxU1l4qvbkBqlg62sFL0izQjxt+Kd2x9H9vgMrF2YAKNWs2vfXwhLYIJYj8u6RTBIWJSH2bEG6LRaaofNgsJls7Biu/rdlfO3ZEhtGVyd2pznziys/+8pMPSoQlFKHDK15sNJ+TuRMsYf1k3ZsJmzYA6pRekLD2DuJnW6J+a017EgLhIhzj8wOKpgeSsXs/NKla9NZbi5HPr7zQj3FyPPO2Ddmov4+RuVecRakP6XxUgYpZWrvQKFZf5IiDa0uX/muauR8bgIo87fPW9Hxfr5mL7MgFXb5sHUfIhEMN+IqMeyOzh+okw+E2WijZfU3zao4+0WrLh3logLrY8T6qtQ+qfZmPtmkpfrEMbPw6q0BJiUAhFEmdg+eQ1pc/JRoYwQZZK/EFNGh8Df4/mgba/rOl04j33Fh0XQjY9Vy90hzseM5Tj7S5eyri1HQc6M5nOsw+3SI3H5KqSY9ep2ieNdutWGSBm0nOUstVqPHdailZi7oFAJfM7ts266HfEvqLN0SoQoh42yzN3P5fCFG7ApThyJyiLMnjJfRN6Ori9neVag1BYOs9F53ZS0c7xFiJ6+GiXPdrTv7mXYtO8lhXBEif8XnIeiLB+ZT68U29pSDNJfTYGpn3sTaud6arfPx4PzWjdqNuduQt59ejgqCvD4tOUeHyUgIiIi+lfX+RrZnr0R+LxLiG2og/2Vz1C76Fi7IVb68WAd6heV47tXjuBsgzZSF4Z+z+vhd8VbGpuQtWKeEoLsh0pRtLkIpYcd0A2NRep/JWvzqGTNjX+NFdZDVnF7PQV5mWqIddjKUby5FLawiYhpUbNjysxD6iRxs1xvRalYd/HuWrEiExIXrkaiNo9Hva5H7KRYl48XHc98mo1y5WZYD4NZGSMkY9zP5Y13lQgcG1F8QPw+QkRAaefXxU18zlQRYv1qUS4CU9GH5ajtoYfpyQzk3KnNo9DBFBsJ7LfAcqgWjh7+MMQmNc0zeflCNViJUCTXU3o8HFPGu8bBFiIykDJdhlgRXjfnY+XrxbDW62CcuhB5cSKIvF+k1DqLFcJaKrarZK8Y7uj47UWJODZWrQarXEzfUqpGNzdxGUiWx+m8DOwrsfIfYp/99DDPXobUCC/XIZsfL0xUgo66LXIZf4SbU7AsVz0gia+KsDkmBKiR54w4H/bWwl+eD88udq957IAxWqyvvBSWSjvgb8Dk3JeQEFGD0g8tsIlRCIlEwlPO87fFdsl5tO16efkUZY7wZ5chZbwIsRDBVJRt8X5RsvertYXNnOsBbLuLtf3TwTBpnrgWfNl6L1SuQflhOaBHZJxz3XokR6rnj23fRiUYen196YwwG8R5I65da2VNB8dbW8ZNx2XoZBg/BaHHxTWxu0qUpjgUUU9g9kx1mrtiLH36ccQnuD4HbEZkuCx1B6r2eXgyNyIFyXfLmvBasd0MsURERHTt6mSQ7YFez98EXX/t64ljOJH1Jep3tR9gW/px1zGcytqHH05oI/oPEmG2f3Mt7RVhQXbWHGQuy0amuKHMXDAfc9/eo9yA+g8xIUGdSeHYm4/pD8sbzzSsmR6LSJFHYLdg5cMzkL5gDqYnb9Fqt5xikSiCm7+sVcp4HHPFutOTZ2CjnCn4FsTMVufyRDc6ETmLFrt8FiBpvDaxHdYaueVi250JZPY4jJC1m7WHUfwp8M6nXyv7phs1AUnKDB6szUZ6Ri6yX5iBpIz5yMyYgR1aQDber8zRxLp5OuKSZ2F2wgxs0YKHOk8iYm+VBeRAeX6csp65CdOx8aCc1gazAeGydqz2AN7500qseSUNc59LQ/rTMzB3kwUFL86HVdm9s6j5VGzXK7Lmt6PjtxErF+xCzXm5XB2sYvrSN1vXjYn0IOK9UFmKFcvysSZnBubIbU6ehhWVXq7j2YdgkrXzhzciXdmWOYj/vQjjNitq6q9XZinIEfsjtnWR2Kd0eT78x8fq87URRnG2eK/2kxWIe3oOZk95DRZlZ79DsTzHMmYhbo1F3f9grZF1y+2S8ywoVpuZR01GigyIZqM4T0VA3JqG+DnyPI1DrkX5q0EzbT2O3a+JY56m7V+piFT+MNyZ4EUQF2HS7Y8z4nNfW03Bq5D/iXo1NTUvbmpWbIVltSx/X66vKlE+D4hr93G1k7F2j7eygLsOy7CZ/fM/4MEEcU2IMswvU44EDLe5/k/SFiOSXl2MGFEgjopCZHroxMn8jPh/R1zPjootyG2nRQURERHRv7pOBdluDxtw3TCtPV5DLU69eATnOtskuPo0Tr/4FRq0mtnuw25E0MOXrA8q79gdCB0Wi5TCDXh/2078c65WEyV2UWtEqKiyujwvNyxUmcdxuLy5+WdlLg643QSLm3TlblmPmJfEej+Wnw2YMkyO84dugPzXM8dhtXax+bMFpe2FQE2oTkuwSugCUu4coeyD3XoAkMHhVDmsp8SIwBEY57GWSKpCXWAozPF5WP/uNmW7JysBQnA2w1TYUXNQfZZRqfHVQrQ6j0Fsi/in/gDKmm7Iq7C0op1ukd7agQNy20JMSP/7TpS8twHLfil7r+2Al8evXVv3wOoQ/w5NwNpd4rff3YCMuBt96nAp0ahuqe1gUXMz0k0iGD4swtMLWnPryjroQs1I+N0GbPpIbOtO9ybL3qqrdTbfdnLA7uyJWjv2Th63a/seVMnDFWjA6KnR2nlahQObmgO6848eTs71+I+aqZ3L4pNmUgOhOO86fmo0Aaluf5wRn8xkxGhTW7K9tQcVcl8iRmOyrOWPH6mWVYUF+cp15sP1Zbeh3LUpv4/Hu+MyVEdJNbbmBuEFthp1oEdHZ6IMsa8iJUonrv0iLMrwUNvaVBtbhR2rWRtLRERE17ZOJMY+6BMbotWansUPbxzE2Yt9rrX6BL5/45h2/90DAbHh8FOGr4QpyMtfjdRHTDCG+KPukAXF261uN/De06OXW9DTnBeB71gVqpyfw2rzZKszA3pwtmaXUrvY/MlFQYevzElGpNIs0tkJzDyYR6k30LqoZC08iHlkzZK40R8xzrUeqVn4s+vw+kIRMEbrEeSwYc8/30P5FXh2WenY5j/mY81mCyqqz6JXPwOMUVOQ+uradnopvkTH79NsxCfnorBUBP0a8dthBkRGJyMrfwOy3JpTe+FsWxFDj9R1f0bWjBhEDgwSs+3Bjs3l4mhdIW1ul28cJ13OZfFRzuVDVdAiWztkE23XP86Iz/s7WrRicFGZC4vyxxs9jLF6pGjNiis+L3APcZ24vjp9vC9RGbozI/XPWoitLMaK1Pnaq3fcNdXG7i3y+Hw9ERER0bXE9yB7f1jTc7E/fvk1Tqt9//gurD90z7k0Iy49grovtaok3c/Qu0UT1stmqtZE2FaM2RPixM3tHGSWd3xLjoM1avPNoSakOJ+pG58iWyy6sEKppOxRh/LfySbJ6mfpmgLkvzL/Er//UdboPKHc6EIEpMLN4t9nRbiTwdpeAfk8ZlN4+LBC3fZR45Aq/m0p+S7ZzNSB8tfH4sHHpmF2xhrYZO2VT7R9l7VVTSFUj3RjO/WPo8yIve8W+JdnY/rD9+LusXEo2C9/WIcR5jaaZnb2+LUQHhWDhBgDajbMQPxD4rfvmINiGd79RcDx8lws1mrfwoclNtcij1+Mte9uwJrFiYhEEsxGf6C+HGvueQDxT8xC+hs2nNVmvVw8btedI6GX1/F5GyreKoFNSdN6GGKbGwhPvvNGt2dknevBsS1N53J8wnIUvL0SS//Lm3feyibarn+cEZ8XC5prOD1Y8XG5OBMBgzEDo4eIsjtfgT1vORNq568vX493x2WojvKdDLGLkThahNjqUqz8dZrSwVwrrrWxazz1pExERER0bfE5yPYae31TbWz9R3XKkM9EiA3KHo7AyOEIzWwOs46Pvm2qlQ0c63wA9zIr1WqSQo1IfCYBsTOysDaxZSc3Hqz9C4oPi1vsQBFP1n2A9bJZa24MQt2adRag4FN5061HbO4G5KVlID33deQszELO8mVY2s4zr7phCco6XT9rMt2fomya5+8foGTXOqVGR+nBdlkaSmVo1F6jUvv5m8rzmE3hIUN9VQ96GGFKa/0KFetJJebCcHcWkiYlIOXlVcpze74pQNHn8kd0MD2zCWtyFyPvL69jSnvtT++eiQXiGCT+ZhVWpSUj4ZlUmCNkjbIDtopCdR6lfEWwjRXll5bo5fETx0lZTg/TysXImt36aVT9I6lInyH2deE65Ij1JGXOxC395JRaWJXa7Y7XYVtTrL4WxTgFr+cvQ86il7F+USyM4QboA0+j3CXcmzOTETs1BXn/HdNx0+mL5HG7cmOV37WXvYcVIhzll1aogXHSMqx/eTGWvroJqXcp7XabONfjP2YmNuUvRnpaFvIKRVmkLUNey1c+XSqrduCA3PZhZrXH6oMWLG0Kep2/vjo+3u46LsPOSRLlKEOsPL/sZ8Ix+Xcu1/xLzZ3NsTaWiIiIyJ2PQTYI/jdog/Un0LBLG/aFFmJ7a++c7R7aF35h6jB2nUSDvFmUbgjCFenAuHIN1n9YBYe/HuYZGeKm9iHo9lm8qPEoRXbqchTtr4VdxCb9wFDUbV2O4hYLlmakYaUIWw4ZXqYmIOG+SISI4GB5Mxfp7d2QButhGGpw/0S0eILPOc+QEOjO21G7vwhLkx9HtlxvRCJuUZ4VrMWBD1vWlTl7Lxb35be61DBpCta+owRd3agpSFmUgaTI07BUyATmm3fmvYDC/WI5UbaR98XCPKwOW0raeUb2tUys2GyFXfaQLEJe+owYGAJFiN3+B6TlqbPk/79SpXZYNyoGMbHjkODV8StEQVEF7Of9EX5nLCZHt+6WyDI/FwVlYj3BRhGGM5DyiDhOPeyoeOsFpMvabS/WgcqVyFxciIpT/ggZE4PYSWYYdGL7ywqwaJ58plUEL9k77nkdjI+kIEeE9Ui7BZ0oWt+0s12ZKerznLbfi/N0u9j/HjoYzLGIuTUAB4rUTqOauKwnfEwsEqZOgXmo7C1ahKusy1VLmI8d/ysOuNZkv6JsuTqg6ez11fHxbsGLMuwMfbDzTy5ivUNaXO8G7cqMmKf0KM3aWCIiIqJmvr1HdnAY+uUMEbdcwsGvUL3I2d2wl1qEWKW34xYdRfVacCeClQB2GnWZ+3DmiDL68oswISYyBGetRSjdr43rkAhcUX1hKXM+5adHzsZNiI1wfY+sk+yx1QCd3Yqi7W0+FXhVMY6PhcFfvoLHm2DfDtlkONzhw3pEud43WgSStn/bON4MbC9tfr7Sq+Mn1js+HFXb29uOjo6TN+tQm65GhgG15cWwtGoqqv5Gr+q9KBZh6kpqf7sEpRz9UbXZpWw96HA9V1xnry/fl7v69p2IiIjop2REQuY8xEffAr0Iimft8tWbryH390Vt3y+PSkDOb7T37tv3ID85G+EvbcDEwerkVrR5XKvnfAuykXqEPDdI6YipcWc5aledUcdL8r2ymYPw48avcMZD0zxvQqzUY3YkQsf2FkNnYH+pXD5OeJVKwZptyYjU2WEtKcT6krMwPvrvmDgmBP5io1feMwNrtDmJiIiIiIj+9eiR+OpapMpHHFuo/XQ5ZqS06KRTMM9djYypJoQrtaOC3YIV986C4e+fNb8ppSVtHtd2cL41LQ7o5vkdr4F9EJhzE3QGEVZnj0TvW7XxTl6GWHfd0d05/1VpJTLXlMJWr4MhOhnpi1IwWYbYUxUoXJzJEEtERERERP/aZuao/fTAAevm+Yi7fRqWlqgtDkPufBIZccqgi0QkPiJCrHykq0J91NGpMM+lTx/lk42iw9rE2qpWrQV9C7INP8Jj9e25Rlw4pw0HBLmH2U6FWOkCLmjvlr1a2d6cg7h7bkfc02lKYc994nZETRAHr+jKNhclIiIiIiK60lLGjVAfO60sxtIFsilxBQqfK4RF6WglBCOiW3eQaq+0oGDBdEyvcO84uGK7y1tW5KdHDExKDW0tSt/IbvWmi0v3jGzPXghYOBrXDdZ6ZWmoQ926avSc5luIvdzPyL5RsE4bIiLq+p5KnKYNEREREV1Jscj7YDHMISKcli1H9NPNDX9T132GRPm2ksMbEfVYtjqypYUbUBZn8NhsWDZZTl+3CQlyHRWFiJuW26qJsm9BFkHQ5Y9EoOxOuP4YamcdQaM6QdUyzLryqia2D/quvhl95LtQxYwnkg/DWdFLREREREREV4tErNo2T3lFonXT7Yh/QRstJL66TX1utrNBNu5lvL/QrHT+WvrCA5i7SRvvwsfX79TB8Y02GNgfAXdow07nzqLhhb34/ojby1S9b058Rz8EyBArfVPHEEtERERERHSVCwpu3YS48/RInypDrFDxEXI9hFjJxyALnN35nfacbC8ETghShty0DLNePxML+E/4mfa6yPOo3+njq32IiIiIiIjoCimA7bg6FBI2Uh1QmBAZrvZibK+xKv/6JC4DE2STYlkb+1brJsVOPgdZbK1GvfLwLtBt5I3oa1aH3TjD7GdHvQ6xMA9G0EitSbL9W5zZqg4SERERERHR1afoK62T22HRyBmvDoY/mQxTuBxywLZfazA8agpSnkmEKUL92jbvamMl34MsfsAPRbVNtbJ9nhqGXmHKF3cyzOZVeRdiw/rjuqcGNtXGNhTZ3J+9JSIiIiIioquKZXUhLKfEQA89Yl/aifff/QAbnjFBqY+17cDKPDlgQs6SLCTNmIe8JfPkiLZ5WRsrdSLIAj++a8X3B7WmwwEhCH5+MHp6CrPeCOuLvs8PR4DWs/GFg1+j7t0L6hciIiIiIiK6OlUWYPYLBbDYHCLM+iMkPAT+PQD7oSIsTUlDqTKTDTX1Yrpw1v6d8q9nzbWxjoot7dbGSj72WuyiZ28ELouErr/2vaEO9tcqUL/L+xDa7Y6BuG7mYPRqej3PUdSmVaGRvTwRERERERF1HREmxETqYC8vhqVSG9dED1NUX1jKKrTvF6/zQVYKC4JuwUgEqs/yChfQePAo7Ouq4TjYdqDtNiwIvacNQd9hvdFNGwd7NU4uOgyHN02RiYiIiIiI6Jp1cUFWCuyDwEwjdIPky2VdnDuHc9Wn0FivfZcCA+Ef2hs9AtxbNF84egQnc465z0tERERERETkwcUHWU23e/W4btqg5mbC3pDNkdcdRP22n6gtcYAOAY12NLBnKSIiIiIioi7jkgVZVXd0N4cgcFIYAgf2QTftbTpuzp/HuWPf4szmajSUntN6P77C+t+F+F8+hOH+jWiEH3ByH95d9zd8eZW+utb0q+cRWvIitnyljfBKNJ7IGICS3L/hqDamQ34DcOsvkvHwz3Wwbv0t/rJdG+9iUPzziK5+0eM0IiIiIiKiK6FTvRa37QIulH6H07/9AsdnfIrqX5fhu8z/afoc/7UYN6MMJ357FGd+qhALHcZOfQi6suVYkr0IS7OXYH3tjZj44E3a9GvXzU8+g3sat6Kk1cPZREREREREV49LXCPbFYQjOnUWfla2BH/d0aCNa2b61RKMOOCsjTRiYsajwNtabegND+GJxLsQId9y2/0cvtq6Gus/OSkm9MPIX8zEI8N06vtvT5Zh3epNOCq/NC0j+J3BV+97scyIOPwqPgqhF4DGH3bjywYjGj9quQ2C2/r8MOjhZzDttv7y7wk4/cU+NPwc+MCHGtm+/fvh9ImT7mXgWnvdeBx7j/ZD6Deta2T73jYdT04you+FRviJsjnw/it4e5ddTGlvu/wQ+uBMJN0RrkzzazyETX96A19cpTXjRERERER0dbjENbJdgQ07t1Xg+gm/xZxnnsLECbcg1Kvnegfj/sfvAj5ejiWLF2HJHz6B3x0PwCAn3fUoHhl0BG++uABLsxdhXdWNuGfcAG2ZKDRsWSKWWYAlq/eg3/2P4mY/bZmQfVipLLME734/BnGTBosJIjxPjoJ9q7rM0neACL38Eamd9Q1/DI+OseO9XHUb3r4wGIPUhVoIwM/G3IJguUwLMsS684Pp8Thcf+ANtfZ6+W4ED2nqotpFAPqGNGLvejGPLJvNNgy7OxrBclJEHOJctmtdw8Dm7RLTEu4ASv4gpon9eblMh7ip94u1ERERERERte0aDLJAw+61eHnxMqzbYQMMDyAp43k8eoengOYi6GYY+n6Dz7drYe9ECf768t9gFYPDh+px6sDHam0qGnH0Hyvw123H1WUCT+B0v3swdsL9GHtzAM5cGIAIkX7lMqe/P4eR/ybGT7gHwRfOIDRsqFhmKCL8DqHsE6222PoJvqxVB9tbX8CNAxHwdRm+UBZrxLe7KlCjLNTCwAcQ99hjuH+M9r1dNyHi+uPYW3JI/drwCfZWKjvZQgO+LXkP1qBo3BP/FH4x7gYE9OqNvmJKwKjB6OuyXad2f920XXJacOVuWLQa2NMf7UPlwKEYrX4lIiIiIiLy6JoMsopGO2p2b8WW15Zh6fvHMcJZg9iexkac1QZbanTIZrQdsePgJx9jr6d35X67G1tF+PRNO+try7FN+GPWAqwv075fEkZMnJeGR406NH5Thi3bD+GUNkXy69lbG2qt8dwZbYiIiIiIiMg712CQvQvxCxcg/jZnA1Y/DIoYAL96uxK+zjga0S/MqE7qPxQRslpRqvsC1oYbMHqM1ia3fzR+MWc6hovpXx2qQujwe5QaSGV9j6Qi6cHB6jL1/dG7tgQ7P9qKnf//CPz695Z5WFmm73U9sVeOF5+9P+gQ7C8m1H2DYxdcf2cMhoeog+2tr6HqOzQOFvNqi/WNvBGh6uBF+BrHTg7AyLH91K9+t+DnAz20SUY4QoOOY+df38bOXfvQGNJPKwuxXfu/xqmIKIwNk8sFwHDvTU3b1bD/CE67bvOEmxDxnQ0+dc5MRERERETXnGuwsyetY6KHjQgWCbCxewACHBV429nJ0PA4/GrqXSJsNaDx++Oo6d0Px5ydPclOmB6LQjAa4efnh28/eQ1r3j8iJvTDzU/MQtyQnnKKCMUV+Mfra9XX+bgsg+5+OP3FWqz8h6x5dVnmglgG32Drqny1me2IRzHn38cgQCZU+xf40uHS2VOb65MBeh6ejNSJYNuI0/srfO7sycm9sycR2P8zGoYLsrOnKuy1DfDQ2ZP22zeLUC32pcFuR0DvI1indegUPE6Ud/RQEW7t+GrrPgT8m057LZDWEdQYtSMov+7H8fGbq/HxN2LfiIiIiIiI2nBNBlmVHwLCBiDgh+M4ddqX4CSXE2Gx+qSMku4C+iHYz+5xfX79B8Cv7jgaWk7y0yE46BxOnWjZg3IA+ob1xOlqz02W21xfO9vQee3ss6u+A0TAPtHit0WQ1ffHqarj6teIR5EyrSc2uQbsNsuAiIiIiIiotWs4yNIVEXQX4v/zIUTU7oOl2g8jRxkRsH8tXn7H1+eBiYiIiIiIVAyydPn56RB68y0YHtKIbz/fDStrXomIiIiI6CIwyBIREREREVGXcu2+foeIiIiIiIi6JAZZIiIiIiIi6lIYZImIiIiIiKhLYZAlIiIiIiKiLoVBloiIiIiIiLoUBlkiIiIiIiLqUhhkiYiIiIiIqEthkCUiIiIiIqIuhUGWiIiIiIiIuhQGWSIiIiIiIupSGGSJiIiIiIioCwH+D9C6yRtlqex5AAAAAElFTkSuQmCC"}}},{"cell_type":"markdown","source":"## Thank you for the wonderful work: [AmbrosM](https://www.kaggle.com/code/ambrosm/ecg-original-explained-baseline)\n\nImprovements:\n\n1. A more accurate definition of outliers is using the median of neighboring points instead of a simple zeroing.\n\n2. Adaptive thresholds - automatic adjustment to signal characteristics\n\n3. Improved noise reduction - better quality of the original image\n\n4. Increased reliability - better handling of borderline cases\n\n5. Improved performance - code runs faster without visualizations","metadata":{}},{"cell_type":"code","source":"import scipy.signal\nimport scipy.optimize\n\nfrom tqdm import tqdm\nfrom glob import glob\nfrom typing import Tuple\nfrom scipy.signal import medfilt\n\nLEADS = ['I', 'II', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']\nMAX_TIME_SHIFT = 0.2\nPERFECT_SCORE = 384\n\nOUTLIER_LOW_THRESHOLD = -1.6\nOUTLIER_HIGH_THRESHOLD = 0.85\nMARKER_ARTIFACT_THRESHOLD = 0.2\nMEDIAN_FILTER_SIZE = 5\nIMAGE_BINARIZATION_THRESHOLD = 160\nSCALING_FACTOR = 80\nTAIL_CORRECTION_FACTOR = 2.0\n\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef compute_power(label: np.ndarray, prediction: np.ndarray) -> Tuple[float, float]:\n    if label.ndim != 1 or prediction.ndim != 1:\n        raise ParticipantVisibleError('Inputs must be 1-dimensional arrays.')\n    finite_mask = np.isfinite(prediction)\n    if not np.any(finite_mask):\n        raise ParticipantVisibleError(\"The 'prediction' array contains no finite values (all NaN or inf).\")\n\n    prediction[~np.isfinite(prediction)] = 0\n    noise = label - prediction\n    p_signal = np.sum(label**2)\n    p_noise = np.sum(noise**2)\n    return p_signal, p_noise\n\ndef compute_snr(signal: float, noise: float) -> float:\n    if noise == 0:\n        snr = PERFECT_SCORE\n    elif signal == 0:\n        snr = 0\n    else:\n        snr = min((signal / noise), PERFECT_SCORE)\n    return snr\n\ndef align_signals(label: np.ndarray, pred: np.ndarray, max_shift: float = float('inf')) -> np.ndarray:\n    if np.any(~np.isfinite(label)):\n        raise ParticipantVisibleError('values in label should all be finite')\n    if np.sum(np.isfinite(pred)) == 0:\n        raise ParticipantVisibleError('prediction can not all be infinite')\n\n    label_arr = np.asarray(label, dtype=np.float64)\n    pred_arr = np.asarray(pred, dtype=np.float64)\n\n    label_mean = np.mean(label_arr)\n    pred_mean = np.mean(pred_arr)\n\n    label_arr_centered = label_arr - label_mean\n    pred_arr_centered = pred_arr - pred_mean\n\n    correlation = scipy.signal.correlate(label_arr_centered, pred_arr_centered, mode='full')\n    n_label = np.size(label_arr)\n    n_pred = np.size(pred_arr)\n    lags = scipy.signal.correlation_lags(n_label, n_pred, mode='full')\n    valid_lags_mask = (lags >= -max_shift) & (lags <= max_shift)\n\n    max_correlation = np.nanmax(correlation[valid_lags_mask])\n    all_max_indices = np.flatnonzero(correlation == max_correlation)\n    best_idx = min(all_max_indices, key=lambda i: abs(lags[i]))\n    time_shift = lags[best_idx]\n    start_padding_len = max(time_shift, 0)\n    pred_slice_start = max(-time_shift, 0)\n    pred_slice_end = min(n_label - time_shift, n_pred)\n    end_padding_len = max(n_label - n_pred - time_shift, 0)\n    aligned_pred = np.concatenate((np.full(start_padding_len, np.nan), pred_arr[pred_slice_start:pred_slice_end], np.full(end_padding_len, np.nan)))\n\n    def objective_func(v_shift):\n        return np.nansum((label_arr - (aligned_pred - v_shift)) ** 2)\n\n    if np.any(np.isfinite(label_arr) & np.isfinite(aligned_pred)):\n        results = scipy.optimize.minimize_scalar(objective_func, method='Brent')\n        vertical_shift = results.x\n        aligned_pred -= vertical_shift\n    return aligned_pred\n\ndef _calculate_image_score(group: pd.DataFrame) -> float:\n    unique_fs_values = group['fs'].unique()\n    if len(unique_fs_values) != 1:\n        raise ParticipantVisibleError('Sampling frequency should be consistent across each ecg')\n    sampling_frequency = unique_fs_values[0]\n    if sampling_frequency != int(len(group[group['lead'] == 'II']) / 10):\n        raise ParticipantVisibleError('The sequence_length should be sampling frequency * 10s')\n    sum_signal = 0\n    sum_noise = 0\n    for lead in LEADS:\n        sub = group[group['lead'] == lead]\n        label = sub['value_true'].values\n        pred = sub['value_pred'].values\n        aligned_pred = align_signals(label, pred, int(sampling_frequency * MAX_TIME_SHIFT))\n        p_signal, p_noise = compute_power(label, aligned_pred)\n        sum_signal += p_signal\n        sum_noise += p_noise\n    return compute_snr(sum_signal, sum_noise)\n\ndef score(solution: pd.DataFrame, submission: pd.DataFrame, row_id_column_name: str) -> float:\n    for df in [solution, submission]:\n        if row_id_column_name not in df.columns:\n            raise ParticipantVisibleError(f\"'{row_id_column_name}' column not found in DataFrame.\")\n        if df['value'].isna().any():\n            raise ParticipantVisibleError('NaN exists in solution/submission')\n        if not np.isfinite(df['value']).all():\n            raise ParticipantVisibleError('Infinity exists in solution/submission')\n\n    submission = submission[['id', 'value']]\n    merged_df = pd.merge(solution, submission, on=row_id_column_name, suffixes=('_true', '_pred'))\n    merged_df['image_id'] = merged_df[row_id_column_name].str.split('_').str[0]\n    merged_df['row_id'] = merged_df[row_id_column_name].str.split('_').str[1].astype('int64')\n    merged_df['lead'] = merged_df[row_id_column_name].str.split('_').str[2]\n    merged_df.sort_values(by=['image_id', 'row_id', 'lead'], inplace=True)\n    image_scores = merged_df.groupby('image_id').apply(_calculate_image_score, include_groups=False)\n    return max(float(10 * np.log10(image_scores.mean())), -PERFECT_SCORE)\n\ntrain = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/train.csv')\ntest = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/test.csv')\n\ndef fit_mean_model(train):\n    mean_dict = defaultdict(list)\n    for idx, row in tqdm(train.iterrows(), total=len(train)):\n        labels = pd.read_csv(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}.csv')\n        for lead in labels.columns:\n            values = labels[lead]\n            values = values[~values.isna()]\n            mean_dict[lead].append(values)\n    \n    for lead in mean_dict.keys():\n        mean_dict[lead] = [\n            np.interp(np.linspace(0, len(values)-1, 20000), np.arange(len(values)), values)\n            for values in mean_dict[lead]\n        ]\n        mean_dict[lead] = np.stack(mean_dict[lead])\n\n    return mean_dict\n\ndef validate_mean_model(val, mean_dict):\n    snr_list = []\n    for idx, row in tqdm(val.iterrows(), total=len(val)):\n        labels = pd.read_csv(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}.csv')\n        sum_signal = 0\n        sum_noise = 0\n        for lead in labels.columns:\n            label = labels[lead]\n            label = label[~ label.isna()]\n            pred = mean_dict[lead].mean(axis=0)\n            pred = np.interp(np.linspace(0, 1, len(label)), np.linspace(0, 1, len(pred)), pred)\n            assert len(label) == len(pred)\n    \n            aligned_pred = align_signals(label, pred, int(row.fs * MAX_TIME_SHIFT))\n            p_signal, p_noise = compute_power(label, aligned_pred)\n            sum_signal += p_signal\n            sum_noise += p_noise\n    \n        snr = compute_snr(sum_signal, sum_noise)\n        snr_list.append(snr)\n    \n    snr = np.array(snr_list).mean()\n    val_score = max(float(10 * np.log10(snr)), -PERFECT_SCORE)\n    print(f\"# Validation SNR for mean prediction: {snr:.2f} {val_score=:.2f}\")\n\ntrain_test_split_loc = 780\nmean_dict = fit_mean_model(train.iloc[:train_test_split_loc])\nvalidate_mean_model(train.iloc[train_test_split_loc:], mean_dict)\nmean_dict = fit_mean_model(train)\n\nclass MarkerFinder:\n    \n    def __init__(self):\n        ima = np.max([\n            cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/4292118763/4292118763-0001.png'),\n            cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/4289880010/4289880010-0001.png'),\n            cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/4284351157/4284351157-0001.png'),\n        ], axis=0)\n\n        absolute_points = np.zeros((17, 2), dtype=int)\n        for i in range(3):\n            absolute_points[5 * i] = np.array([707 + 284 * i, 118])\n            for j in range(1, 5):\n                absolute_points[5 * i + j] = np.array([707 + 284 * i, 118 + 492 * j])\n        absolute_points[15] = np.array([1535, 118])\n        absolute_points[16] = np.array([1535, 118 + 492 * 4])\n\n        template_positions = [None] * 17\n        template_points = [None] * 17\n        \n        for i in range(len(absolute_points)):\n            if absolute_points[i][1] < 118 + 492 * 4:\n                if i % 5 == 0:\n                    template_positions[i] = (absolute_points[i][0] - 87, absolute_points[i][1] - 50)\n                else:\n                    template_positions[i] = (absolute_points[i][0] - 37, absolute_points[i][1] - 13)\n                template_points[i] = np.array([\n                    absolute_points[i][0] - template_positions[i][0],\n                    absolute_points[i][1] - template_positions[i][1]\n                ])\n\n        template_sizes = np.array([(105, 60)] * 17)\n\n        self.templates = [None] * 17\n        for i in range(len(template_positions)):\n            if template_positions[i] is not None:\n                template = ima[\n                    template_positions[i][0]:template_positions[i][0] + template_sizes[i][0],\n                    template_positions[i][1]:template_positions[i][1] + template_sizes[i][1]\n                ]\n                self.templates[i] = template\n\n        self.template_positions = template_positions\n        self.template_sizes = template_sizes\n        self.template_points = template_points\n        \n    def find_markers(self, ima):\n        markers = [None] * 17\n        \n        for j in range(len(self.templates)):\n            if self.templates[j] is not None:\n                t = self.template_positions[j][0] - 100\n                l = max(self.template_positions[j][1] - 100, 0)\n                search_range = ima[\n                    t:self.template_positions[j][0] + 100 + self.template_sizes[j][0],\n                    l:self.template_positions[j][1] + 250 + self.template_sizes[j][1]\n                ]\n                res = cv2.matchTemplate(search_range, self.templates[j], cv2.TM_CCOEFF)\n                _, max_val, _, max_loc = cv2.minMaxLoc(res)\n    \n                top_left = max_loc\n                markers[j] = np.array((\n                    t + top_left[1] + self.template_points[j][0], \n                    l + top_left[0] + self.template_points[j][1]\n                ))\n\n        for i in range(3):\n            if markers[5 * i + 3] is not None and markers[5 * i + 2] is not None:\n                m = markers[5 * i + 3] * 2 - markers[5 * i + 2]\n                markers[5 * i + 4] = m\n\n        if markers[14] is not None and markers[9] is not None:\n            markers[16] = ((markers[14] * (284 + 260) - markers[9] * 260) / 284).astype(int)\n\n        return markers\n        \n    @staticmethod\n    def lead_info(lead):\n        begin, end = {\n            'I': (0, 1),\n            'II-subset': (5, 6),\n            'III': (10, 11),\n            'aVR': (1, 2),\n            'aVL': (6, 7),\n            'aVF': (11, 12),\n            'V1': (2, 3),\n            'V2': (7, 8),\n            'V3': (12, 13),\n            'V4': (3, 4),\n            'V5': (8, 9),\n            'V6': (13, 14),\n            'II': (15, 16),\n        }[lead]\n        return begin // 5, begin, end\n\nmf = MarkerFinder()\n\ndef find_line_by_topdown_sweep(ima):\n    top = np.argmin(ima, axis=0)\n    median_top = int(np.median(top))\n    top[top == 0] = median_top\n    top[top > median_top + 300] = median_top\n    \n    strip_width = 64\n    for strip_left in range(0, ima.shape[1], strip_width):\n        median_top_strip = int(np.median(top[strip_left:strip_left + strip_width]))\n        if median_top_strip > median_top + 300:\n            median_top_strip = median_top\n        strip = ima[median_top_strip + 80:, strip_left:strip_left + strip_width]\n        all_white = strip.all(axis=1)\n        if all_white.size > 0:\n            first_white_row = np.argmax(all_white)\n            if first_white_row > 0 or all_white[0]:\n                first_white_row += median_top_strip + 80\n                mask = top > first_white_row\n                mask[:strip_left] = False\n                mask[strip_left + strip_width:] = False\n                top[mask] = median_top_strip\n\n    mask = np.tile(np.arange(len(ima)).reshape(-1, 1), reps=(1, ima.shape[1]))\n    mask = mask >= top\n    ima &= mask\n\n    bottom = np.argmax(ima, axis=0)\n\n    bottomx = np.maximum(bottom, np.median(top) + 100)\n    mask = np.tile(np.arange(len(ima)).reshape(-1, 1), reps=(1, ima.shape[1]))\n    mask = mask < bottomx\n    ima |= mask\n    ima[:, :-1] |= mask[:, 1:]\n    ima[:, 1:] |= mask[:, :-1]\n\n    return top, bottom\n\ndef get_lead_from_top_bottom(tops, bottoms, lead, number_of_rows, markers):\n    line, begin, end = mf.lead_info(lead)\n    top = tops[line]\n    bottom = bottoms[line]\n    begin, end = markers[begin], markers[end]\n    baseline = np.linspace(begin[0], end[0], end[1] - begin[1])\n\n    pred0 = (top[begin[1]:end[1]] + bottom[begin[1]:end[1]]) / 2\n    if len(pred0) < len(baseline):\n        baseline = baseline[:len(pred0)]\n    pred = baseline - pred0\n\n    pred /= SCALING_FACTOR\n\n    if lead in ['aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']:\n        pred[:4] = np.where(pred[:4] > MARKER_ARTIFACT_THRESHOLD, pred[4], pred[:4])\n    if lead in ['I', 'II-subset', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3']:\n        pred[-5:] = np.where(pred[-5:] > MARKER_ARTIFACT_THRESHOLD, pred[-6], pred[-5:])\n    if lead in ['I', 'II-subset', 'III', 'II']:\n        pred[:2] = pred[2]\n\n    pred = np.interp(np.linspace(0, 1, number_of_rows),\n                     np.linspace(0, 1, len(pred)),\n                     pred)\n    \n    outlier_mask = (pred < OUTLIER_LOW_THRESHOLD) | (pred > OUTLIER_HIGH_THRESHOLD)\n    if np.any(outlier_mask):\n        for i in np.where(outlier_mask)[0]:\n            start_idx = max(0, i - 3)\n            end_idx = min(len(pred), i + 4)\n            neighbors = pred[start_idx:end_idx]\n            valid_neighbors = neighbors[(neighbors >= OUTLIER_LOW_THRESHOLD) & (neighbors <= OUTLIER_HIGH_THRESHOLD)]\n            if len(valid_neighbors) > 0:\n                pred[i] = np.median(valid_neighbors)\n            else:\n                pred[i] = 0\n\n    pred = medfilt(pred, kernel_size=MEDIAN_FILTER_SIZE)\n\n    if lead in ['II']:\n        n_tail = number_of_rows // 48\n        tail_values = pred[-n_tail:]\n        tail_median = np.median(tail_values)\n        tail_std = np.std(tail_values)\n        threshold = min(OUTLIER_HIGH_THRESHOLD, tail_median + TAIL_CORRECTION_FACTOR * tail_std)\n        pred[-n_tail:] = np.where(np.abs(pred[-n_tail:]) <= threshold, pred[-n_tail:], tail_median)\n        \n    if lead in ['V4', 'V5', 'V6']:\n        n_tail = number_of_rows // 12\n        tail_values = pred[-n_tail:]\n        tail_median = np.median(tail_values)\n        tail_std = np.std(tail_values)\n        threshold = min(OUTLIER_HIGH_THRESHOLD, tail_median + TAIL_CORRECTION_FACTOR * tail_std)\n        pred[-n_tail:] = np.where(np.abs(pred[-n_tail:]) <= threshold, pred[-n_tail:], tail_median)\n        \n    return pred\n\ndef convert_scanned_color(ima, markers, n_timesteps):\n    crop_top = 400\n    ima = ima[crop_top:, :, 2] > IMAGE_BINARIZATION_THRESHOLD\n\n    iima = ima.astype(np.uint8)\n    ima = (iima[:-2, :-2] + iima[:-2, 1:-1] + iima[:-2, 2:]\n           + iima[1:-1, :-2] + iima[1:-1, 1:-1] + iima[1:-1, 2:]\n           + iima[2:, :-2] + iima[2:, 1:-1] + iima[2:, 2:]) >= 7\n    \n    tops, bottoms = [], []\n    for i in range(4):\n        top, bottom = find_line_by_topdown_sweep(ima)\n        tops.append(top)\n        bottoms.append(bottom)\n\n    tops = [t + crop_top for t in tops]\n    bottoms = [b + crop_top for b in bottoms]\n\n    n_timesteps['II-subset'] = n_timesteps['I']\n    preds = {}\n    for lead in LEADS + ['II-subset']:\n        pred = get_lead_from_top_bottom(tops, bottoms, lead, n_timesteps[lead], markers)\n        preds[lead] = pred\n\n    preds['II'][:len(preds['II-subset'])] = (preds['II'][:len(preds['II-subset'])] + preds['II-subset']) / 2\n    del preds['II-subset']\n\n    apply_einthoven(preds)\n\n    return preds\n\ndef apply_einthoven(preds):\n    residual = preds['I'] + preds['III'] - preds['II'][:len(preds['III'])]\n    correction = residual / 3\n    preds['I'] -= correction\n    preds['III'] -= correction\n    preds['II'][:len(preds['III'])] += correction\n    \n    residual = preds['aVR'] + preds['aVL'] + preds['aVF']\n    correction = residual / 3\n    preds['aVR'] -= correction\n    preds['aVL'] -= correction\n    preds['aVF'] -= correction\n\n    residual = 2 * preds['aVR'] - 2 * preds['aVF'] + 3 * preds['II'][len(preds['I']):len(preds['I']) + len(preds['aVR'])]\n    correction = residual / 17\n    preds['aVR'] -= 2 * correction\n    preds['aVF'] += 2 * correction\n    preds['II'][len(preds['I']):len(preds['I']) + len(preds['aVR'])] -= 3 * correction\n\ndef is_color_image(ima):\n    return ima.std(axis=2).mean() != 0\n\nsubmission_data = []\nold_id = None\n\nfor idx, row in tqdm(test.iterrows(), total=len(test)):\n    if row.id != old_id:\n        path = f\"/kaggle/input/physionet-ecg-image-digitization/test/{row.id}.png\"\n        ima = cv2.imread(path)\n        shape = ima.shape\n        good_shape = shape[0] == 1652\n        \n        if good_shape and is_color_image(ima):\n            markers = mf.find_markers(ima)\n            n_timesteps = {lead: row.fs * 10 if lead == 'II' else row.fs * 10 // 4 for lead in LEADS}\n            preds = convert_scanned_color(ima, markers, n_timesteps)\n        else:\n            preds = None\n            \n        old_id = row.id\n\n    if preds is not None:\n        pred = preds[row.lead]\n    else:\n        pred = mean_dict[row.lead].mean(axis=0)\n        pred = np.interp(np.linspace(0, 1, row.number_of_rows),\n                         np.linspace(0, 1, len(pred)),\n                         pred)\n\n    for timestep in range(row.number_of_rows):\n        signal_id = f\"{row.id}_{timestep}_{row.lead}\"\n        submission_data.append({\n            'id': signal_id,\n            'value': pred[timestep]\n        })\n\nsubmission = pd.DataFrame(submission_data)\nsubmission.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}