')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load all dependencies you need\n from coffee import ***** ","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport random\nimport pandas as pd\nimport pydicom\nimport os\nimport matplotlib.pyplot as plt\nfrom timeit import timeit\nfrom tqdm import tqdm\nfrom PIL import Image\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold, GroupKFold, StratifiedKFold\n\n#color\nfrom colorama import Fore, Back, Style\n\nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as Layers\nimport tensorflow.keras.models as Models\nimport warnings\nwarnings.filterwarnings('ignore') #Ignore \"future\" warnings and Data-Frame-Slicing warnings.\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's start seeding everything to make results somewhat reproducible. Anyway, in keras it is quite hard to get 100% reproducible results.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed): \n random.seed(seed)\n os.environ['PYTHONHASHSEED'] = str(seed)\n np.random.seed(seed)\n tf.random.set_seed(seed)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# First glimpse at the data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = '../input/osic-pulmonary-fibrosis-progression'\n\ntrain_df = pd.read_csv(f'{ROOT}/train.csv')\nprint(f'Train data has {train_df.shape[0]} rows and {train_df.shape[1]} columnns and looks like this:')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.sample(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To get an idea of the meaning of the weeks column, let's check some of their values.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_unique_df = train_df.drop_duplicates(subset = ['Patient'], keep = 'first')\ntrain_unique_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can clearly see, some patients took the first measure of FVC before and some after their baseline CT images.\n> the relative number of weeks pre/post the baseline CT (may be negative)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# CHECK FOR DUPLICATES & DEAL WITH THEM\n# keep = False: All duplicates will be shown\ndupRows_df = train_df[train_df.duplicated(subset = ['Patient', 'Weeks'], keep = False )]\ndupRows_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can eighter drop all except the first or last duplicate or average them.As there are only a few duplicates, we can drop them without a bad conciousness for loosing to much data for our first apporach.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.drop_duplicates(subset=['Patient','Weeks'], keep = False, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'So there are {dupRows_df.shape[0]} (= {dupRows_df.shape[0] / train_df.shape[0] * 100:.2f}%) duplicates.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = pd.read_csv(f'{ROOT}/test.csv')\nprint(f'Test data has {test_df.shape[0]} rows and {test_df.shape[1]} columnns, has no duplicates and looks like this:')\ntest_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data Wrangling","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Getting the format right","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In this section we are going to do all the Data-Wrangling and pre-processing. For this we are going to define some functions and transformations, which then are applied to the data.\nIt's good practice to concatinate all tabular data (train, test, submission), to ensure all data get's the same & correct treatment.\nIf you don't do that, you need to be careful with some steps, e.g.: \n* Standardization or Normalization (e.g. MinMax Scaling) in ```test_df``` will not have the same range of values (e.g. min/max values) and therefore scaling than in ```train_df```.\n* The categorical features might have different categories in ```test_df``` than in ```train_df``` (e.g. ```test_df``` only contains male, Ex-smokers).\n\nSo let's concatinate all our data first and then start with the transformations.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## CHECK SUBMISSION FORMAT\nsub_df = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\n\nprint(f\"The sample submission contains: {sub_df.shape[0]} rows and {sub_df.shape[1]} columns.\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"sub_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need to get this into a format which we can easily use for making our predictions, so let's split up the ```Patient_Week``` column into a ```Patient``` and a ```Weeks``` column to align it with the train & test data-format.\nThen we merge our info from the test data to the ```submission_df```: that's the fastest way of getting the correct format for predictions and submissions later on.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# split Patient_Week Column and re-arrage columns\nsub_df[['Patient','Weeks']] = sub_df.Patient_Week.str.split(\"_\",expand = True)\nsub_df = sub_df[['Patient','Weeks','Confidence', 'Patient_Week']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df = sub_df.merge(test_df.drop('Weeks', axis = 1), on = \"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# introduce a column to indicate the source (train/test) for the data\ntrain_df['Source'] = 'train'\nsub_df['Source'] = 'test'\n\ndata_df = train_df.append([sub_df])\ndata_df.reset_index(inplace = True)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\nThe first big challenge is data wrangling: \nWe could see that some patients take FVE measurements only after their baseline CT-Images, and some took measurements before that.\nSo let's first find out what the actual baseline-week and baseline-FVC for each Patient is. \nWe start with the baseline week:\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_week(df):\n # make a copy to not change original df \n _df = df.copy()\n # ensure all Weeks values are INT and not accidentaly saved as string\n _df['Weeks'] = _df['Weeks'].astype(int)\n # as test data is containing all weeks, \n _df.loc[_df.Source == 'test','min_week'] = np.nan\n _df[\"min_week\"] = _df.groupby('Patient')['Weeks'].transform('min')\n _df['baselined_week'] = _df['Weeks'] - _df['min_week']\n \n return _df ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df = get_baseline_week(data_df)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What we can see here, is that the Patient with ID ending on \"430\" had his first FVC measure 4 weeks before the first (baseline) CT images ( = \"Weeks\" column -4) were taken. Then the patient took the next FVC measurement 9 weeks later. \nIn the next step we need to baseline the FVC values. Note, that the **BASELINE-FVC it not the minimum FVC**, but the first measurement, meaning the measurement taken in the \"min_week\" or ```baselined_week = 0```.\n\nFor getting the baselined FVC I first wrote the following straightforward function:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_FVC_old(df):\n # copy the DF to not in-place change the original one\n _df = df.copy()\n # get only the rows containing the baseline (= min_weeks) and therefore the baseline FVC\n baseline = _df.loc[_df.Weeks == _df.min_week]\n baseline = baseline[['Patient','FVC']].copy()\n baseline.columns = ['Patient','base_FVC'] \n \n # fill the df with the baseline FVC values\n for idx in _df.index:\n patient_id = _df.at[idx,'Patient']\n _df.at[idx,'base_FVC'] = baseline.loc[baseline.Patient == patient_id, 'base_FVC'].iloc[0]\n _df.drop(['min_week'], axis = 1)\n \n return _df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This apporach works fine, but as it contains a lot of look-ups, its slow and didn't feel right. \nBtw: there is an even worse approach: Using ```for row in df.iterrows()``` is roughly 8 times slower than using ```for idx in df.index```. \nSo I looked up how other people solved it and I found a rough equivalent to the following function:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_FVC(df):\n # same as above\n _df = df.copy()\n base = _df.loc[_df.Weeks == _df.min_week]\n base = base[['Patient','FVC']].copy()\n base.columns = ['Patient','base_FVC']\n \n # add a row which contains the cumulated sum of rows for each patient\n base['nb'] = 1\n base['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\n \n # drop all except the first row for each patient (= unique rows!), containing the min_week\n base = base[base.nb == 1]\n base.drop('nb', axis = 1, inplace = True)\n \n # merge the rows containing the base_FVC on the original _df\n _df = _df.merge(base, on = 'Patient', how = 'left') \n _df.drop(['min_week'], axis = 1)\n \n return _df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The second apporach is using ```transform```, which is not as known as ```apply```, but faster for basic-operations not involving multiple columns of a dataframe. Here is an interesting post about it for those, who want to learn more:\n[Apply vs transform.](https://stackoverflow.com/questions/27517425/apply-vs-transform-on-a-group-object)\n\n\nI wanted to know how much this speeds up the processing, you can find the results in the following:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def old_baseline_FVC():\n return get_baseline_FVC_old(data_df)\n pass\n\ndef new_baseline_FVC():\n return get_baseline_FVC(data_df)\n \n\nduration_old = timeit(old_baseline_FVC, number = 3)\nduration_new = timeit(new_baseline_FVC, number = 3)\n\nprint(f\"Taking the old, non-vectorized version took {duration_old / 3:.2f} sec, while the vectorized version only took {duration_new / 3:.3f} sec. That's {duration_old/duration_new:.0f} times faster!\" )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df = get_baseline_FVC(data_df)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Not the data has the format we need to work with.\nIn the next section, we go trough two possibilities on how to normalize, standardize and prepare the data for the neural Network.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Preparing the data for the Neural Network","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### The first apporach is using sklearn, as it is super famous and used frequently.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# import the necessary Encoders & Transformers\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler\nfrom sklearn.compose import ColumnTransformer\n\n# define which attributes shall not be transformed, are numeric or categorical\nno_transform_attribs = ['Patient', 'Weeks', 'min_week']\nnum_attribs = ['FVC', 'Percent', 'Age', 'baselined_week', 'base_FVC']\ncat_attribs = ['Sex', 'SmokingStatus']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To get a ColumnTransformer who outputs the whole DataFrame compatible format, we need a class that takes attributes which we dont want to change, and simply passes them through.\nThis class needs a ```fit``` and a ```transform``` method, so that the ColumnTransformer itself can use ```fit_transform``` like for the numerical and categorical attributes.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.base import BaseEstimator, TransformerMixin\n\nclass NoTransformer(BaseEstimator, TransformerMixin):\n \"\"\"Passes through data without any change and is compatible with ColumnTransformer class\"\"\"\n def fit(self, X, y=None):\n return self\n\n def transform(self, X):\n assert isinstance(X, pd.DataFrame)\n return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## GET TRANSFORMED DATAFRAME\n\n# create an instance of the ColumnTransformer\ndatawrangler = ColumnTransformer(([\n # the No-Transformer does not change the data and is applied to all no_transform_attribs \n ('original', NoTransformer(), no_transform_attribs),\n # Apply MinMax to the numerical attributes, here you can change to e.g. StdScaler() \n ('MinMax', MinMaxScaler(), num_attribs),\n # OneHotEncoder all categorical attributes. \n ('cat_encoder', OneHotEncoder(), cat_attribs),\n ]))\n\ntransformed_data_series = []\ntransformed_data_series = datawrangler.fit_transform(data_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Okay, now we encoded all the data and want to have a look at it. But wait, we only get back a series element? How to put that in a Dataframe again?\nSadly, not as easy as I had hoped. We need to use ```pd.DataFrame(data, columns=column_names)``` function, for which we need the column names. But we used OneHot-Encoding. So the number of columns now depends on how many different values/categories a categorical value has, because for each unique value we get a separate column: e.g. If for the column \"SmokingStatus\" we only have the values \"Smoker\" and \"Never-Smoked\", we have two resulting columns if there are additional possible values like \"Ex_Smoker\" we get more columns. And we also should not get things wrong and mix the columns up. The following code is getting our data back to a Dataframe and preserving the correct order.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# get column names for non-categorical data\nnew_col_names = no_transform_attribs + num_attribs\n\n# extract possible values from the fitted transformer\ncategorical_values = [s for s in datawrangler.named_transformers_[\"cat_encoder\"].get_feature_names()]\nnew_col_names += categorical_values\n\n# create Dataframe based on the extracted Column-Names\ntrain_sklearn_df = pd.DataFrame(transformed_data_series, columns=new_col_names)\ntrain_sklearn_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Okay wow, that was not as easy as I expected. The upside is, it's easy to add some additional features and pump them through the Pipeline or to change the Pipeline itself (e.g. exchanging MiNmaxScaler with StdScaler, RobustScaler etc).\nThe downside is: it's not intuitive (at least not for me) and takes quite some lines of code.\nDoes anybody have an idea on HOW to improve this? ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### The 2nd approach is doing all the legwork ourselfs.\nThe good thing is: we dont need a NoTransformer here, as we simply can work in the DataFrame itself and not change any data which we want to preserve.\nDownside is, we need to implement the MinMaxScaler by hand. Make sure to not call it MinMaxScaler and shadow the already important MinMaxScaler from Sklearn!\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def own_MinMaxColumnScaler(df, columns):\n \"\"\"Adds columns with scaled numeric values to range [0, 1]\n using the formula X_scld = (X - X.min) / (X.max - X.min)\"\"\"\n for col in columns:\n new_col_name = col + '_scld'\n col_min = df[col].min()\n col_max = df[col].max() \n df[new_col_name] = (df[col] - col_min) / ( col_max - col_min )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def own_OneHotColumnCreator(df, columns):\n \"\"\"OneHot Encodes categorical features. Adds a column for each unique value per column\"\"\"\n for col in cat_attribs:\n for value in df[col].unique():\n df[value] = (df[col] == value).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## APPLY DEFINED TRANSFORMATIONS\nown_MinMaxColumnScaler(data_df, num_attribs)\nown_OneHotColumnCreator(data_df, cat_attribs)\n\ndata_df[data_df.Source != \"train\"].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get back original data split\ntrain_df = data_df.loc[data_df.Source == 'train']\nsub = data_df.loc[data_df.Source == 'test']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Okay, so the second apporach (using our own implementation) was more straightforward and less code. Downside: if you want to replace the MinMaxScaler with another scaling method (RobustScaler, StdScaler), you need to implement it first.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Model & Loss\nIn this section we are going to define the loss & a first model.\nFirst we are taking care of the loss. We are trying to minimize the following:\n\n\n\nThe global minimum of this function is achieved for delta = 0 and sigma = 70, which [results in a loss of roughly -4.59.](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/168469).\n\nGetting our model to predict the Confidence and FVC values (which is what we need!) is not working fine so far, as you can read [here](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/167764).\nCurrently the way to go seems to be pinball loss. \n","attachments":{"image.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAVsAAACDCAYAAAAj1bgbAAAgAElEQVR4Ae1972vbTNfm+wfqX9AnE7i9BX+xeT6IQEWWioDeQFAhIhQHIroEFSJK1l26zr73um8fh+2jwo322eDsFkHXheBQMIQVBEHhWs7MSBrZcn7btZMplMj6MXPmGumaM+ecOfMvUP8UAgoBhYBCYO4I/Mvca1AVKAQUAgoBhQAU2aqXQCGgEFAILAABRbYLAFlVoRBQCCgEFNmqd0AhoBBQCCwAAUW2CwBZVaEQUAgoBBTZqndAIXArBBIMoxjjX7e6eaE3JecDDC7ShdapKrs7Aops746ZeuLZITBCf7sO+8/Rcrb8agC/2YJ/pgh3OTuIS6XIdpl7R8m2FAiMPprQd/oYz0Ga0UcDWtOBdxggmPx/PEBCdf4aITxwYNou2jsmDNtHeDEhTBygpTnoX06cVz+XBgFFtkvTFUqQpUTgZw+WZqLzYx7SpYj2NGha1X9daNJj9Ld0aJtdjIQJY/SnBU2z0Pspy0Rl6ai/HUDptzIuy3OsyHZ5+kJJsoQIxIcNaC+7mI8BYYTuuouQqa9F40efbLTeRlyrjQM0NA3WJ0mvHvdhaxrq7+LiIQBp5EHXbPSlW0s3qB+/FQFFtr8VflX5ciMQI2hoMD4M5yPmrwGCnX6ZyMkc0PQxEAQcv28wzdc7lUUYwCdt+I8AJbpNI7Q1Dc7JBHvLj6rj34aAItuFQp9g+KkNq1mrnDqWtJeFyrUslY0QvrFhGjXU7Q7iiwGC1xbct22YL0y0P4+QXITwdyy4Ow7MJj83OW0eRwEcw4SzR3bOAN3jAMEXmTATDA4d2LttOC9NtD9FCA8dOFTmRhthNj0/78LQNLhfJ2sQeP0aIzp0YLzQK/vTP7sjrmnMHF1BzqCFmaGSbDUX4ZVcB2nKGrT9SD6pjpcEAUW2C+uIEXpke1uz4H8KEZ34MHUN2t889KIIEYUVzfimFybivSoaomsbMNZv+99Gd4b9k6bBrXcxkr/anLyaPiLh8BkcCLvmywCxUNz4ubJTiKbgum6jJxxIyYnDytIPBqJ1KeJ3Rh5ZkH512fXGuxjxUYsd54PemQ9Na6AgPwkgFgGgQTccdE4ihMcuWqRtbncQUX+eDZHcMUwsft+Cthty8wGraozeJm93NdlaeTu5ZAnCHQ3aq15ZW5bEVoe/DwFFtgvBnj7wFjSJBKha7ugok0WlOJcD+Js1aPq0fa/y/kc5OUK4b6Oh6fCixYwCg0MDwVmKwQFpig3434qGDN4S6Vjo5l74TOuTHEXMI6+BiDP7l3zhZJprp1cRvO2CjFg0gCDUceTD3enmZD7+JBxReZ1ZqcJpJU336QqTcXJqnz1y098khKtPatF3JVuAD0AesqHlpmrV9cUhoMh2EVgzj3aZBKha/qFLZHGNLOxeyVGTfG2jprXQ+X7NQw+9dOZDn/J6P7TQm54XU2HdlwhDnFvvoDAGDODTzCA/lyLaJ5LWIWuBnLhnYUykqUHTPQwqtNCZZMtw0eB8lj1RgvxLct/U1uI618DlwYSuZQOKVmoTIGy2U2YERbYFost3pMh2AX2SaVdeyZQmpnyNCSdHpTycFGqHhcaGXymSZL4aJ48BlQmuUrjHPSk87TSdzlsnzhWmAADCS184r7gzSysRUBVJS+IKh5K23Zem7sX1WWTLnVaTxCjqn1FWUWrVkXgXKuJkr3WQVbw7SrOtwnc5zimyXUA/jD/bPC5Sno4Kbdf8OBFU9GuE/p4N56ALb6sF72tCMT1o59P5FIN3NoymnnvJ4yMTrRc1GK9duDs+2lsG2l/HQNyB2ayjZjhwd1z4ezaMvZAH54tAeXs/gLdpIvjGqS2NO3A22vD3LNTWNOj7UUF6lVjd1WZroXuNNp4Ke60lrdbKzuWmAJqyM1ODie45mWMctL/8k2u6r6QwLYExI+m4A+tthPQqRpewiMaA0FCt46wPRuhte4gypxOz2danbLbxO9Kgy1P19MxnJpf2X/kQwdBKz0N0DgP04vL5MpQxl32iTHZPVejXjJkSUGWzTTAgB+FRpOy4ZdAX/kuR7SIgP+/C1OrwTsUH94s7y/StXh6ozsXgtsAaEdxlH44mCPWbz5w+ffKSn/mo73bR3dGgM2fKAL4RoEPkI8wMTCNrdvDfDwwER2QK4KQE9pGS6YGmpzpqwhlD9zOt+Wcftl7jcl6FcDUddmmqPH+wuCZnlBYR8HOybVvYMsmEkAzgNV30x8KM8LdMExcOSYpR/XMEipelaT8f+DS03g8QvRVhVWLGQc4140BaFHBBCxom7ahZPKvknEpouayG1oFY8SXBxImZHKGZXNLF7DDT5qvIFsI+LA0izNY/Yf/nRQlNvmJWQAsn5MEqq1r9XRwCimwXhHVyGsBpUjiSC3vDgX8ST3urmRZTKzmGSLzSdD5NkFyGTNN1SesFmROG6L7SOGHSDJtiM9cCnCUJhscWO2YGCFF+cELkzx1f6WUE3zCYZsuf83nsJrt3zjbhKezH6G/rqG3LsacJwl0+MMgW0vHXNgv9sjYddL6J8IQkRseuw9h24Ww4CKIBeq/r0JstmOQUI7vszxBtw4Tz2oR9GCI8slFv2mjvWHAOo4lEM9w00JqKs01ZCJ/ZtNBm4WUuOtG4egaQjDCIOnAqiVQAQM4ximSYZe8VsxBjw0F714a57aP/vSKWdkac7TiO0H9rsEFnCnJ1YmEIKLJdGNQ3V8Q00qmIA2GvfRdieDritkWa/tJ9P1P+gQszAydfrvVxrZc7WPgxkJX/7z0Xmm7BP+EhSmM2bZafA9cAdR9RklSTyM3NeRJ3sBVk65Jp4j6tohnFZhEBUVVE+nOI4eV1poaqp8rn+Aoyp3IF2eCAZizl+9WvxSKgyHaxeF9bG7NN6kKzJJ31LID75xe2isna92C94TGYpIHqOz309g1uT2RaaINHJtCHrXNNFeCaWeOIfPhEpjqMwxhjCoeS6sHlEIPzCxajqbOwqRThrgZ9u4POrg1mvrhW8id8kZleymaNu7aW4mfLkQt3LeE293PTUKMqN8JlCLfpI66IuLhNyeqex0FAke3j4PhIpYzQf92C9a6H/nsX7kEfwysgvRggigYYCccN2R11Mkkccvsi01jX6jDfBPC3LbS/CIcPI4oa6httBAcOrL2QO0nSEXqvW3COQoTHHtz9HobkhzsLYBptBO8dONsWak0DpiD4R2rgShbDFhvcK8oAwI8OzCnb/BxgoBhjvT2VZ4HMTOEbA/5phdlhDmKoImcjoMh2NjYrcqVsKpCF5pryIhdCyLU/pWPKZ1tbXpsnrWYzFKEu+xunyHbZe+hG+bipwMzDl4oHmMNLWghRXFFHd0dA7NTwMLPq3au9xRNsp4ZzpbneAqrfeosi298K/0MrFzGUIul0eF6Ul5x1i2TUmVmhuKyOFAIKgQUjoMh2wYCr6hQCCoHniYAi2+fZ76rVCgGFwIIRUGS7YMBVdQoBhcDzRECR7fPsd9VqhYBCYMEIKLJdMOCqOoWAQuB5IqDI9nn2u2q1QkAhsGAEFNkuGHBVnUJAIfA8EVBk+zz7XbVaIaAQWDACimwXDLiqTiGgEHieCCiyfZ79rlqtEFAILBgBRbYLBlxVpxBQCDxPBBTZPs9+V61WCCgEFoyAItsFA66qUwgoBJ4nAopsn2e/q1YrBBQCC0ZAke2CAVfVKQQUAs8TAUW2z7PfVasVAgqBBSOgyHZOgNO+YYM8e36K8bcY4yXM8j+n5qtinwkC6c8Y8c8VfbF/xoiqtoSfU989AtmmiPZ1aH94GKwo5o+NbXrmwzB8DPKdSmhnWwu9i+mahsc2jPUWapoGTdPYJovGugH+vzivverhywcTRrPG7svvfR9LhcboGMV1/YWB9tdciOK+yxj9QxdWswa6xzBaaNkBBpdA8peH1m642tuXX0XwCMMXOsdKrws8Ba7ZeU2Ddxqjs2Ggtcbx17QaWuu0O3EBF+IOjPy6jrrRRngpXWeHKcZRF+1tA3Wdl9Fqmmh/GiK9GqK7ZaLzY/KZ1f9Nm41af4oNRleuOSP0tmqwPy1G/oeT7XkXpiCK1QU9e0vGCPcsGDt9jLNTd/17FaGtt/gW4/mzs8mW3zKAxzC00PuZP8QPaDO/pgbt7YD/Hvdhs3tdhGK33fITMfw/Gmh/qWpBguGxg5qmw3jTQ3xZjI7J9w4sXYeuabAq9jMr17Eav0Z/Woxs9QOBnST26JMNXSsGQNqxmAYwbdZA881HvdFGWAXrOIRn6NDWLATRCGm2Zfgvep8a0HUifQ9Rdl6SY9UP70+2KcZnHThrFe88gORbF65twtlrwzZMuMcxJtWG0VcfzoYNd8+Fadjwv06Q5q8RwgMHpu2ivWPCsH2EkwoPbfOuNya+1/n0ygPJNsXgoIGWYbCPVPvDh6xnzUfkOZaaEPAaWnvRPcmW46HvhRMvxg1ke96FQR/6egfDiubRxo2NwwzZAfyMmCdfHABEMNP1U6EJov0WNK0F92TipRR1Dg6IFHT4ZxVCrNwpvuswEaj7tRhU8mZQX2sO+pmGeuZzsn3V49u95zfSwQi9Tb1ylpD+6MLWNegbAeJJNqBH2XbyGrTN3j3fqZIgS/fjzmQ7DtFeN9BqtlDXaTZRDHhZ49JvPlpaA14kAL0awGtoaBwM8hkXDY40WHazb+CiB4sUhU/ZaDhGf0uHttnFSAxyfPCdJvfRRxM0c6z+KjKpHv73YWTLXiQH/Z98h1d6se3PWWMfLtzKlXDZh6PpaP81+XFfT7a5VrUfVTaZXpJi1jBC7xWf8k6RIo3Say76FV0Qvyei1dDKSXu6KvpwNG2Wxjx9/3KfieGLjzn/IEsC02zCQ67zio9V0/zinLg/+eKiVjXbYVqRBk2vxpw/zgdH40PVMFoSaCV/3Jls81bSN1FFtgn624RpuR/id6QI2OLdFnxTGsCIXDXkCl8coFEiXwBiVlh/lykuQhg2O6/Dl01HuZyPd/Agso0PG2gIwZMTh2sGq67dPgBb+iiryep6sh285eRZjMqkoTpof+Ej++BAg5fzcIJwh99f1thIk5thf8q0K72NqNL0wBvNyHa7P6GVPwCQ3/loRp7yR5sO4G90+OyLMHnVLbQZpukSrhODDZmF1uwKezufxdAAZn68Ticisl3MNPV3wP3oZJv1w4SmyRUBDc5JAggizU1rouH0nWhanZkEaDZIfeOdyqiIWeEfwcQMfITuuob6NYqIXMp9j+9PtkyLcwotil7kBicBBsh9JbrXcyniIxum0ULtZRvhjyF6ezacfQ+OUYd9OECSxOju2rB3yTHUgkPnsrqu6H4L1ksDdbuDOFNML0K0bROGXoP1IcbotAN3y4W3Z6HVdBCc5SUAENPWSlPAdWSbzQqMwoHCpk3S70xO8Ze/VJqk7QLpqYeGNGWSH6FBkV68bGCUr5WO0wRJkjW+dGXlfvCBr2yDZdPIWTZZVJlnUgzeNqpt2Oz9p/dd+gZmoJQmSWHHnXHPqp5+dLLNTGozyLb1YYj0rzZX7DI/hgAv+y7cr/8P0R7nokqynRxQs293Qxp859Ah9yZbptW+LWwoJFs2+miNyZGjWnLuic887zf/tY9nTMVoGrDdx/iHsH3qNnrnnDQKmVz0hX2Hn+MjICPJ/Rrck3HeibyDaDrTgP+NPjjeca2DCAmz/4ipfGmEFNOiSs3wGrLNtE6NvNwcA+YZ1z0MZjhU8jZlL1saw2/O8nZnZK5PjPLVffJUznL7swZtrcUjEYw6d/7N9JyL/pO0IWY7fNnBsKIfcjLfnL+tb5n7hN7FwsR1F0kzvCdstjNs5/I7Lx/LNWZka/35v4WJYoZmW2En5pqwZFaSC36k4/uRLbNVTYBEAqURPGYnm+GUeCShJ4sZfXLgfh5hxGyOGhzJbsyN4pKxndwdH43CEURhQpsdxL+EZppNtfPzgqwaPuL8o+PTDnJ6FBNIoRllBFgScjbZ5qO0pHENj4zZXnGCORvZBbGTPbYlOQ/KVWfRC5IzqHTDU/wh+keTZgdseir9nmp24VDjM7MYQbMF/6xa0+c2RE1yXE4VuIATCQbHLosGMucdQXIewttqQNM9RBIkT4VsOYFXcNoj9uK9yJa8d/pelHsGZXlGxyZX8dfnq5LLdfLj7GORvY3ZOdkOl42o8jkAYlqoT5KlmNaUwoeENlp2etyPbLOPVtYOaJSVf0+19ZvPoz9o2kNafcND5ridujfXnG8YteOARWFI3xGS8xija2y8U3U94MRdZjkzZzhZ/Vl4nDw7ILKVf2f3Sn+zviDione88TYqTE3SfXRYaFHFcDtxC4sA6W+b6J5PX3m0M6x/pQiS7x1m9nKFvX9mPb/GiO8Y0M9mCyWnFJ/NXvuuzhQg+w4nCO4GMwJ9c7myMfGtZn1ydzNCNiufkGWm7Pe7cHeyZXGk17xA7Dqfdk975e8n5O2eEt7nks1UnJNfkMzWJmmSVD4fJMhsUK4tmy7Kzig+Ck5qSfch20wDm3CgXCW41nSaE6gF61V1SFLRisyMMDG4FDcwUgh3WuW2M01QhyerMaVnlvdH/jGWTDrpjfbofHr6yoKltxHKJvmJ5mbvBdkQZ/6jAYzMWzNvePgF1lbdlWS9uZ1UK5s2v7yLQsTf1dZRub2Prtk+2EHGv6VrHWQVZs6l1GyJlKrjOIsXJwsz0m7Qbu+izdCKKutjuaOLGgFkGqg82olzsgbKQeXhWempD/MDhYGM0H2ZxbmO0Nv2hNc+04zlEU/cy15UujdbTSRG6tIHnkk4w4yQvVh3DnjPFkFo0KtCkrJqxd+sP+QBQ76FAvyNA8lhSBdPPRbSxAlnxFZAGc06TNuBs9lC7XWI0VkAa8NAY00Q03kPznqGBwARVG7vB/A2TQTfZL1ZluBxj7MP7c6LM6jNLIZZZzb8a6XKwr5Kg7v0BC1GMeQohgThGwP1tRqsnTbctx5cw0L3e4roLV8ZaG7acLYtNIwgN1mxwH3bQ7BvwTyM+WySOW5tePsujDUdmiDN0WcHZrOB2hspzjuJ0Xltwz0K0DZMdL+TjNznUCJO4UB2DgK4Lymck+4jx7MD842P9maNLYaZHHynyZZMGwGCo0gysUm45IczNFuI0K+JEDyutWbOSKFAyEoURHmZua8q9EsoKVWOYv7OyLO/27Yjb9CNB3fTbIVNtrXtIzgMZv/fs1BnL20d3uliPrBM05AJhZ8ra42800gToEB1sdJLdAKRcnLqwcwdf5m9VnL4TdzbksiOOdL+VrUwYQbZZh93JUFf13dCLnIEZkHd192e0pJUilM0yxEU6RjROwvm2xDj3B49RvSBPjgNWtOBd9hB9HcPxvsugoaG1vsI4b4BY+/f8HHLRv8LxTPymQ7DOw+1ooFKR22Xf/j0UdbmHFrDIcjikMv9fh08+TXxgepbsi0+vzp1MP7qokVRHjt9jKTXPIl7jLC68jT9Zw/2TgfBlpab4Ni7+K+7+I87XXR3aOAMEX20YJAPQdjm9TWXLw2m924tQMycobowMw3R+ZuGGoVfUvlGB12KPGmKdzC7l+y53wLUtBp2//O/Idi30dB0mG8CBJ+IwEmB0GEc8fhTMhc4J2OwhQNrfBl++tWFptuChAsopsg2M+HMWkySPzqLbAG+qKFezKrEogbZL1G9qEGXlt6KRQ1SeB/z31R+M0KpkgfOW7cjb9CNB3ciW64VchMB1wBucbw132lU1sL4sAad8hFINsbhUQu6pCXQvWncgf3ChGlbaOcrqWhZpYla04S13SlWAo37cPQanNLaaX5v66VVvpcm47eMs03/8piHfHI9vjG18ixr3eRfPup2JpcnTt5W+p0gPm7DIq885Qp4ZcHeCdCvXPZUMWVkcau1Us6ANElL01H6SPVs2swCxbkZIr2M4BuUb0BipJJsj/FjiK5NeR541AG9nyzvw7qJzkQM+8zakgG6h53pJZ0zHwDSixDBNs+twN6fVw7ax9F00qFfNL2P4GmcyGhFH8VLc7xoppKdzyrjMyh9P0KaJogODBiHMcb0jmWkJ0w9LuW/oPIpB8O6hkxj5YPfNEFOmh7SyIMuBsz0Rw/2i2KhEiPyzOyQkXgmoohAmrTZjuMI/bfGDL9DjA7lrZD6KYsY6UgmPLZcd8OAvdfmS3KPIkkh4AIwrd+gJb0u7A0H/smwbGMXMytjw0F714a57aMvD4B5O/j7XqzS5Beub0f+8K0P7kS2ty71ud4o7MHTtuoZmu2y4sQ+4lrJhss/3LInmplf8o87I+gY6ZUYeHQL/kmE6GyIsTQILmuz5y4XM2u10GFT+QFb4caIkWnUE34QQaTWuz6iaIChSBnHtOGM9Jij1ELvXMTxsllXC504ZXG9pXulxtGUWZeUIDaFXnPRjSJE8Yj7C1hZOhiRg6/O0in0ccKZMKXZinoGB7XVCTVkisG0v4aa8pjtUGQrvYQPPxSriiacb2D2JNnu+/Ca5llCrvlcjtDb4bMFMpHkWmteOXcKsiXabKDR4R0FMHdDXDANTMqVcTmUUk7mBTyrg3zmkwCUGa6xxm2jbHpLJgIZDUG2hcM2wfB0hH/SyihBlMMPBrS1AJ33BhsYWfnNAP2PFigagYU4SqQ6+tyGH/3Itd808tmsjZGtdB+lB43/L+WO0PmAe8WPnQ+0qKc8U60kW7JnN+VQSblhy3fMciOUwjiFjI/cDkW2j933zHEyaS9cLc02/RbAXLPQ3ncQnHKXfPI9wuBi0gxAZgQT9S0P/mtynNXQeulwc0E6Qu91C85RiPDYg7vfw/Aa7/5jd8MylsfMLC9aIIdhe9NhGijJSeQ2nVc1ZUu2W9sdhF+68HY99GgKfNGH88KG/8GDs23DWGvBpFWPVA4lcGkYMMmGTAVfDRC8pIxZPXT3XbhH5AQdo/+6DnOnDedNjydpITPPSxPepxD9ozbcQ5qypxgcmjD2AgRZ3xrmVDKeabJNmSPQF+/NMvZDSSb2vVYtCHr8diiyLSH/OD9YPtumJ9mPV4tsHwcFVUoZgczMck1ETfmBlfg1TbYrIbYQkvLZ1id8MvOTX5HtnLBl2sqPTJWjnRoGC1sgMKcmqWIfgoAwC0yGTj2kyGV4lnZqmJ7xLINkt5CBdmqI5xkFXZZBkW0ZD/VLITAHBEYI81DJrrSDxxyqUkUuLQKKbJe2a5RgCgGFwFNCQJHtU+pN1RaFgEJgaRFQZLu0XaMEUwgoBJ4SAopsn1JvqrYoBBQCS4uAItul7RolmEJAIfCUEFBk+5R6U7VFIaAQWFoEFNkubdcowRQCCoGnhIAi26fUm6otCgGFwNIioMh2abtGCaYQUAg8JQQU2T6l3lRtUQgoBJYWAUW2S9s1SjCFgELgKSGgyPYp9aZqi0JAIbC0CCiyXdquUYIpBBQCTwkBRbZPqTdVWxQCCoGlRUCR7dJ2zWoJlnyf2GxvtcRX0ioE5o6AItu5Q/wMKqD9x5oB4nw79GfQZtVEhcAdEVBke0fA1O3TCMSHDfzrfwoQ5AmyZxx/pD2w1D+FwPNEQJHt8+z3x2v1o+xAmiLa16H94WEwuafk40mqSlII/FYEFNn+VvhXv/L4fQvWn2wvV96YJEZ3x4a5Xoeu0y6uHUQ3bfN03oWpadA0rVzW6sOjWqAQyBFQZJtDoQ7ujMBVhPZaG9GVeDKN4Rs2eudCPU0G8JtEohZ6F7NKTzE4aKBlGNCJcP/w2bbcs+5W5xUCq4qAIttV7bklkHt0bML8WGi1yReXaaeNgwEyawBtdU0aq/Z2UC3xzx4szUH/Z4ygwbVb+/NNqnB1UeqsQmCZEVBku8y9s8yypRG8tTZCyeOVfHWZdqrvhIUj7MznZPuqh4KWi4aRc63xLmYnkhOH36u02wIgdfRkEFBk+2S6crENIY21dchJUq45TRIkmVoLINN29SrNlkLGSKvNFNl0AF9ot86JxOJyBepYIbCiCCiyXdGOm7vYV0P0jyNkPFiqj5GiRJKli/KPMXqbZBpoIZjmZTCt9m1hcqAnc7NDI1C2WxlKdbzyCCiyXfkunEMDLvpwmWOrDv/bdPk03c+m/tNXizPjzzZ0rQbnU4UBgULG9ArHGZkndG67db9KKnJRrDpSCKwkAopsV7Lb5ij0VQT/dYDwzzaPDtjqT2i35MiqIMlJkS56sPUWvK+VujFGH03oe1HuSJMfJ8cbc6qtdyvtvPK96lghsCoIKLJdlZ5auJwjdF+Shmmge15UTjbY1sTUv7gqji76cF4Y8E+F3TUdI/42LoiVQsZ0s1RuqQx2nWu37b+UdlvCRv1YWQQU2a5s181fcCJWFvuaa6BEwNeQJIl0NYBvSERL57750PejXGDSXPU9KWIhv1Ic0GIJpd0WeORHl0MML/Nf6mCFEFBku0KdtXhRs9hXm0UMpH+1UcuJt0qaEXpbOuqv2vClPAneVgPGhyF/QNhkW9v+9bkU9izU2aqyOrzTp6LdphidDTC6d6BFgv72f4C9MyP3hIR5+KOqf1bxXILR6QCjbOHMKjZByKzIdoU7bxGiZ/bTxmGPabWdaz5icogxbVQsvZWPs1CuPNqg4h75/tLxlN14ES0XdYxDtJs10S4PM5Zm3EKgFIN3BoyDByTjiQM0tidt6LeoesVvSc98tJo+BitOuM+PbK8SJPfotDTuwFm3EHx7KlrWLb9AKTpA33l+HzpHaYTeK7Ih359saVagPygNZYJwp1WKDhl99eFsmDDWNNSaDvyTp5tTmJmVtvvFYplbvr7LdNszI1sxLb5Hpw2PWtDIqfN9mbpvMbJw+2l1rOxiJPjdtWTxwvckWxaXrKP99d72AyAO0NosVuHRLMI4iJCIHMKjPynMTkPjJufl74byvvUzp+lqm5SeF9nSNEzTCvvhfTv+uT13OUD342AiBOw5gfAwsmXLkHUpYc+doUsR7dVQRGZktnQLvZ9ZYQP4LD75FmF52SMr9Vek4S1ccVMAAAy+SURBVPydJqUH4vWsyHb00YCm6fBOH4iaevyZIfAQsqXpvwZtNyxC3+6KHqWgfCnHHA/RYYtO5NlGJqP2ZN9vMsVotLx7RaMxloRsU8RHNkyjhdrLNsIfQ/T2bDj7HhyjDvtwgITypO7asHddWM0WHDo3+dKOI3R2bFhbLtqvTZh7IUZIEO4bMNYNtNZ47GatacDY6CCGVO9mB9FZD23bhvvaZs9SOH5yGsDZNGG8MNGuCNBPL0JWp73rwcvrnBRM/V5tBDIiE2YE4TSrvahDf9VBeOLD3fXQ3qyx97f3Q7LrpxHaD5pNkUZXmzZB/EqRyEkoIGm7pXSWCYaf2rByJx//BjIHpPWpetHJUvbXeReGpsErogiXUsxZQi0H2dLITV7WHxxMTS9yoube64aLvniJ+Ll6ab39iILttRb8KKPgEbrrGvwz0fQkhEsecMnuhazeU5GZijyelzGCLAfr/+rCZE4NXpY2kbkqr5MF7ycId7nXulJzvorgrXPSJ+K/+b9X5Imd1Xvq/IIQmCBbVmvmNNNgHY+E1ioIT175xlJIasiiMe4sMD1/C8caeezJRFZ2YvJQPG3Ngv8pRHTiwyRTw9889KIIURRjLI0Ld5Zt0Q+Ib7iUrH7RMjygvqUg29EnB+7nEUYi96kj5TMd/Un5UBvwchIFMnNATqQi038pUJ6Wi665RUYpCqwnB4KUqYrCmqjjsnr5WvwhejsO3OMYZ8c22l/GIFKmEVXO3cocFlSeSA8IJIgOLDgHIUZLvPFhptE8lb8PePfv8GgV2WbnXIR5dEt2TnKkiRSTlQMwkyDF8KSLKLe9ymLxxOry9yBfzY8paXtTg274GGS6Bs3a3pFT1y4lbuff081T8fRc5MdodiAipPPq5ncwQrhvo0GmvqhqFBjAJ4WpKoPc/IR6tJKXgmx5a8gJQFMc2eifnat6oYtz8fsGi4Osb7YRHHpwty1Y2z763/M3Lyfo6Zc+q2P2C0jZqcq7DQhjvdYoadeP1is3FjRE176Ndkz32OhKsbH/+Mc/8JT+Xw/V/XEql1tBorjluevI9tcI/V2+Uq6eD9pSzZSCsuHfsGsxEXIL+lYXw+J1B4RGXSgDvFyuqMjfmFRf6ZB/F7Tyj9NegvBNDdq8yfeMlKJZ8imyLXXR/X/E3Ju6Lo+k4txmr/CEsxyossMhe+lndRBJJJwUWkHQhZxFHRW5qYAszrQUTC6e0TxES6zFFm1URw9DIHvHJI31wWSbIjpwEHztoc2iCPgqPVlOGuSvnzLzhRL17V4+m0p+RGw5b5ZHuGzfFN/BrdJX0juuwzmRbLoUo16lcMpCP/CYDQYzCV2R7QPhFY+LqXopybQ4ly/1zPOd6iwMJj31YX74J/f2VgWcpwlSRoYDeDT9yOJrz3twsgD9ijrkBqWRx8wPZGKg+uxjmlQJsp2w4dJz6QW55Cr+3cNmK++CUFGiOrUwBB5Atrew2War9AxpiyGwXYu9a1dNjT7ZkIkWGKO/ZaF7AfDVfBNhYEKWkjmMMCQNmxzSB114W5SpLRGmM5EHg/Jd2AZauoEOizMnLddAfa0Ga6cN960H17B4DHrcgdmso2Y4cHdc+Hs2jL2QK0u/RggPHNj7AbxNM18gxBYMbbTh71morWksj0Ylp8+w2abzHgEe6T1bGjNCNhLLOUz5ufJUfXBApgYXYTJCb5OHvvCEKWXNIPneg0ubD5JTTbxkXEsgp4EB/0xMjti+WbPDwcr1FS8vC/Sf2L5lHPkwS3azR+olVcxvRuABZAs+0Ley3BBVLWG5fem9biMzVVIKyilSlJ5lS1i1FpwDOU+CCzObvTE/hrQI4Bd3lulbhRbMiyOC1lEjcwGbNeosDp2R9RrffHNwUId73IWr6XCJiH/2YO90EGxpeZpM9p3shvifBwaCIzIFCKJm314Lne9kltBR2+UJiMjJXSP/yc8+bL3G819ckRNbx8w96IRiVHBEZs7TUBqoJJyW6XBpyDY+rHEDf+5sAGjVlm4EJZsVjYL2CxOmbaF9kk38U8THLqx1B22y2W7acI9CKTtSisGhibphwdl0EEh7a1O92jUB57y+FowNuT6uDdAobWw48N66sF85aNPOBpVD8jJ1uZLlTghQmNd6CzWRy4HCBjtfK859qTjHEq/fLs6W+wU0TjSTuxZPCZyZscphXMzpKZkIWNhi04Sz58LeoOW8cb7iLC+SLfSplZYBA8JeK4iRbXX0tQ1NJyWH3n0KO4vgaZmZgbdR3+5jlCQYHlvQ1sROG6L84IS2q+eOr/QyYpnhaOk787cIUienc0MjYs6lKx1UxdnGR7Qrcw0uObKX/N/SkO3vwylFqgjy98H/DGrmG1kWWmtlk0VEjdYI0COt9jpNuLKA+51kYZQZieZFcHutfTxAHHMSGxzo0HdDjK/Ex8K0zIwY+eq11hGZ2MpEnZX/7z0Xmm7BP4kQnQ0xZkoVnzFQucThTJvWfURJUrEARGixJd9JLvBKHCiyXYluUkKuNAJiI8ti+lvVmmJKTBpkvglm1a2PeI5pizo3F1Cx6VkA97/8Vzadb+97gvQpfliH+6kHj2aa+UaeXNNlMb5rtB09lcBjjRuMeIlMdRiHMcZkrpPqweUQg/ML5m/RWSRGinBXg77dQWfXFmVJDWW5ERornW5Tka3Un+pQITAvBJhfQZriV9ZDyWY0rXLX4sr7H+XkCP3XLVjveui/d+Ee9DG8ApLvtOghyyLG7bqtlw4C4etgmu6LFnN2tTcddGKh8TIbbQ31jTaCAwcWW8VJLD5C73ULzlGI8NiDu99joWpE7qbRRvDegbNtgcw05pvpxPLkIynF0T9K2xdbiCLbxeKtanu2CPAwreu3FEow+NjFYOnNj3xFJTcblDuUa8rCtlu+dO9fpDkbT8DxrMj23q/AE3lwHCHYNlAXO9qWV5b5D0iW/UTwedRm0E4NPA72UYtddGEsBKt6lRdzeJWS5jxUOLFTQ2U85UPLXuzzimwXi/dS1cbDh3QYrzvoRyG6YjWTc0RTyAiDH0/gDV8qxJ+CMCOE+fY7XWl5MJCcdYutjr5kkUJPoc2P0wZFto+D4+qVwuIbNbRK27TwmNC6lD+iumEpRp8p8Y+Gqqlk9TMPP0sfs0s7/j6q5vRwuVQJCoHbIKDI9jYoPcF7yMHBcoPK9kGRDlA/uMVOW+xeaSr5vQNDp3jHeWrDIlToNvI9wT5TTVptBBTZrnb/3VN6kQrwlZyQmqJ2+E4Wt0oHSFnUdEcKUZrMr3pP0a57TBA8W8V03X3qmkJgCRFQZLuEnTJ/kcQKpFKqOp7Or2o1XfKtA8d20Tlsw9jospR7LGGImM6PPjswmw3UspCdyxBto47amgV3z4W378LY7GL465o19RS7+a0L13bgH7o8vzEDYoRwz4a9T+vva8Uy0PmDpGpQCDwqAopsHxXOVSlMBNBLiXSSU58nX2eJ0It2pN98tHQL3XOALW1eayNMKO5SQ42C0WmdvNFBl9JQimxN40+0XDqArfGEQQDP1uT+t+OZa+pTWkFFiU4oYp7uZ1ozz8mqCzmHH1rQsqWdhYjqSCGwEggosl2JbpqDkFe09ZCJ1mabbSFEiUWin5PrlsWa98lN9uTpPK2TvxqyXTFyZxml4vvLg5aZGUSokPP30Yw19T183te54+tXiuGflM2qj7F4jicmEbKIpZ1zQEQVqRCYKwKKbOcK76oXzjXSnESz5jB7rY3uaYyYlmhmmZ3iVKS0FLtpZHlJWUJosY6+ck393xE0NNR2uyzkLL4Qa+NZ4u1s/T03fVDmtkRKVpSJpP4qBJYdAUW2y95Dv1U+vlKoSHk3Qv+Nj//xmda5t+HtmyxDE0uF2QzQ/2iJaASR6WqHll1yW3BNJBDhaTMn19Rzh11RDwX/xxj/H3LYidSZPzowtBqCDwGMql0NfitOqnKFwM0IKLK9GaNnfQdbu/7SRfdTF96ui47Y3HJIGwaKbYeYXbdhwNzpg4eyk81VR71pwztsw9ruIBba6Kw19QnlAn7pofe1j84bV6TBpHX7ddgHHXivHdhGDa2XtrDrPutuUY1fQQQU2a5gpy29yCVTgSzt7DX18l3qWCHwFBFQZPsUe/U3t4mZCnQv33UgF+eaNfX5PepAIfBEEVBk+0Q79rc16zws1scfD6T92Gavqf9tsqqKFQILRECR7QLBVlUpBBQCzxcBRbbPt+9VyxUCCoEFIqDIdoFgq6oUAgqB54uAItvn2/eq5QoBhcACEVBku0CwVVUKAYXA80VAke3z7XvVcoWAQmCBCCiyXSDYqiqFgELg+SKgyPb59r1quUJAIbBABBTZLhBsVZVCQCHwfBH4/yPe5+FhkkCiAAAAAElFTkSuQmCC"}},"execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# CONFIG Section ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In this section you can configure the following:\n* Features used for training\n* Basic training setup: BATCH_SIZE and EPOCHS,\n* Configuration for the loss function\n* Optimizers, Learning-Rate-Schedulers incl. Learning Rate start- & endpoint\n* Custom Logging Callback\n* Checkpoint-Saving Callback\n\nThe Learning-Rate scheduler below is inspired by Chris great [Melanoma-detection notebook](https://www.kaggle.com/cdeotte/triple-stratified-kfold-with-tfrecords). \nFeel free to experiment with the scheduler and it's max/min and decay values.\n\n**Ever wondered why lr_max is scaled by BATCH_SIZE and therefore bigger for larger batches?** The reason for this is the following: the larger the BATCH_SIZE, the more averaged & smoothened a step of gradient decent is and the bigger our confidence in the *direction* of the step is. As there is less \"randomness\" in a huge averaged batch (compared with for example Stochastic Gradient Decent (=SGD) with batch size = 1) and our confidence in the direction is higher, the learning rate can be bigger to advance fast to the optimum.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"######## CONFIG ########\n# be careful, the resulsts are VERY SEED-DEPENDEND!\nseed_everything(1989)\n\n\n### Features: choose which features you want to use\nfeatures_list = ['baselined_week_scld', 'Percent_scld', 'Age_scld', 'base_FVC_scld', 'Male', 'Female', 'Ex-smoker', 'Never smoked', 'Currently smokes']\n\n### Basics for training:\n\nEPOCHS = 1500\nBATCH_SIZE = 256\n\n\n### LOSS; set tradeoff btw. Pinball-loss and adding score\n_lambda = 0.8 # 0.8 default\n\n\n### Optimizers\nADAM = tf.keras.optimizers.Adam(lr = 0.1,\n beta_1 = 0.9, \n beta_2 = 0.999\n )\nSGD = tf.keras.optimizers.SGD()\n\n# choose ADAM or SGD\noptimizer = ADAM\n\n### Learning Rate Scheduler\ndef get_lr_callback(batch_size = 64, plot = False):\n \"\"\"Returns a lr_scheduler callback which is used for training.\n Feel free to change the values below!\n \"\"\"\n lr_start = 0.00001\n lr_max = 0.00001 * BATCH_SIZE # higher batch size --> higher lr\n lr_min = 0.000001\n # 30% of all epochs are used for ramping up the LR and then declining starts\n lr_ramp_ep = EPOCHS * 0.3\n lr_sus_ep = 0\n lr_decay = 0.991\n\n def lr_scheduler(epoch):\n if epoch < lr_ramp_ep:\n lr = (lr_max - lr_start) / lr_ramp_ep * epoch + lr_start\n\n elif epoch < lr_ramp_ep + lr_sus_ep:\n lr = lr_max\n\n else:\n lr = (lr_max - lr_min) * lr_decay**(epoch - lr_ramp_ep - lr_sus_ep) + lr_min\n\n return lr\n \n if plot == False:\n # get the Keras-required callback with our LR for training\n lr_callback = tf.keras.callbacks.LearningRateScheduler(lr_scheduler,verbose = False)\n return lr_callback \n \n else: \n return lr_scheduler\n \n# plot & check the LR-Scheulder for sanity-check\nlr_scheduler_plot = get_lr_callback(batch_size = 64, plot = True)\nrng = [i for i in range(EPOCHS)]\ny = [lr_scheduler_plot(x) for x in rng]\nplt.plot(rng, y)\nprint(f\"Learning rate schedule: {y[0]:.3f} to {max(y):.3f} to {y[-1]:.3f}\")\n\n\n# logging & saving\nLOGGING = True\n\n# defining custom callbacks\nclass LogPrintingCallback(tf.keras.callbacks.Callback):\n \n def on_train_begin(self, logs = None):\n print(\"Training started\")\n # self.val_loss = [] not used for now\n self.val_score = [] \n \n def on_epoch_end(self, epoch, logs = None):\n # self.val_loss.append(logs['val_loss']) not used for now\n self.val_score.append(logs['val_score'])\n if epoch % 50 == 0 or epoch == (EPOCHS -1 ):\n print(f\"The average val-loss for epoch {epoch} is {logs['val_loss']:.2f}\"\n f\" and the score is {logs['val_score']}\")\n \n def on_train_end(self, lowest_val_loss, logs = None):\n # get index of best epoch\n best_epoch = np.argmin(self.val_score)\n # get score in best epoch\n best_score = self.val_score[best_epoch]\n print(f\"Stop training, best model was found and saved in epoch {best_epoch + 1} with score: {best_score}.\"\n f\" Final results in this fold (last epoch):\") \n \n \ndef get_checkpont_saver_callback(fold):\n checkpt_saver = tf.keras.callbacks.ModelCheckpoint(\n 'fold-%i.h5'%fold,\n monitor = 'score',\n verbose = 0,\n save_best_only = True,\n save_weights_only = True,\n mode = 'min',\n save_freq = 'epoch')\n \n return checkpt_saver","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Loss Function","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# create constants for the loss function\nC1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n\n# define competition metric\ndef score(y_true, y_pred):\n \"\"\"Calculate the competition metric\"\"\"\n tf.dtypes.cast(y_true, tf.float32)\n tf.dtypes.cast(y_pred, tf.float32)\n sigma = y_pred[:, 2] - y_pred[:, 0]\n fvc_pred = y_pred[:, 1]\n \n sigma_clip = tf.maximum(sigma, C1)\n delta = tf.abs(y_true[:, 0] - fvc_pred)\n delta = tf.minimum(delta, C2)\n sq2 = tf.sqrt( tf.dtypes.cast(2, dtype = tf.float32) )\n metric = (delta / sigma_clip) * sq2 + tf.math.log(sigma_clip * sq2)\n return K.mean(metric)\n\n# define pinball loss\ndef qloss(y_true, y_pred):\n \"\"\"Calculate Pinball loss\"\"\"\n # IMPORTANT: define quartiles, feel free to change here!\n qs = [0.2, 0.50, 0.8]\n q = tf.constant(np.array([qs]), dtype = tf.float32)\n e = y_true - y_pred\n v = tf.maximum(q * e, (q-1) * e)\n return K.mean(v)\n\n# combine competition metric and pinball loss to a joint loss function\ndef mloss(_lambda):\n \"\"\"Combine Score and qloss\"\"\"\n def loss(y_true, y_pred):\n return _lambda * qloss(y_true, y_pred) + (1 - _lambda) * score(y_true, y_pred)\n return loss","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Neural Network Model\nIn this section we build an initial neural Network. The code of this section is derived from [Ulrich's](https://www.kaggle.com/ulrich07) great [notebook](https://www.kaggle.com/ulrich07/osic-multiple-quantile-regression-starter), which also inspired me to change my loss to the above coded version. Please support the original Notebook creators! The chosen quartiles are simply derived by testing; using 0.25 and 0.75 leads to worse results.\n\nFor the architecture: It's good practice to use numbers of units following the schema 2^x, with x element of N (= resulting in 1, 2, 4, 8, 16, 32, 64, 128,..). \nWe are going to use dropout for regularization and not a too broad and deep network, as the training data is very limited.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model():\n \"Creates and returns a model\"\n inp = Layers.Input((len(features_list),), name = \"Patient\")\n x = Layers.Dense(128, activation = \"relu\", name = \"d1\")(inp)\n x = Layers.Dropout(0.25)(x)\n x = Layers.Dense(128, activation = \"relu\", name = \"d2\")(x)\n x = Layers.Dropout(0.2)(x)\n # predicting the \n p1 = Layers.Dense(3, activation = \"relu\", name = \"p1\")(x)\n # quantile adjusting p1 predictions\n p2 = Layers.Dense(3, activation = \"relu\", name = \"p2\")(x)\n preds = Layers.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis = 1), \n name = \"preds\")([p1, p2])\n \n model = Models.Model(inp, preds, name = \"NeuralNet\")\n model.compile(loss = mloss(_lambda), optimizer = optimizer, metrics = [score])\n \n return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create neural Network\nneuralNet = get_model()\nneuralNet.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## GET TRAINING DATA AND TARGET VALUE\n\n# get target value\ny = train_df['FVC'].values.astype(float)\n\n\n# get training & test data\nX_train = train_df[features_list].values\nX_test = sub[features_list].values\n\n# instantiate target arrays\ntrain_preds = np.zeros((X_train.shape[0], 3))\ntest_preds = np.zeros((X_test.shape[0], 3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the following we want to create leak-free folds to get a robust **cross-validation strategy** in order to evaluate all our models & our training. The idea is to avoid having the same patient (= PatientID) in training- and in validation-Data, as this might lead to evaluate a higher CV-score for a model which is luckily learning/memorizing the data for a particular patientID which is also frequently occuring in the validation-data. \n\n\nThe idea on how to do that is coming from @PAB97 [Pierre's great notebook (CHECK IT OUT!)](https://www.kaggle.com/rftexas/osic-eda-leak-free-kfold-cv-lgb-baseline#kln-440)\nPlease note, that we still don't use propoer stratification based on 'Age', 'Sex', 'SmokingStatus'.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## Non-Stratified GroupKFold-split (can be further enhanced with stratification!)\n\"\"\"K-fold variant with non-overlapping groups.\nThe same group will not appear in two different folds: in this case we dont want to have overlapping patientIDs in TRAIN and VAL-Data!\nThe folds are approximately balanced in the sense that the number of distinct groups is approximately the same in each fold.\"\"\"\n\nNFOLDS = 6\ngkf = GroupKFold(n_splits = NFOLDS)\n# extract Patient IDs for ensuring \ngroups = train_df['Patient'].values\n\nfold = 0\nfor train_idx, val_idx in gkf.split(X_train, y, groups = groups):\n fold += 1\n print(f\"FOLD {fold}:\")\n \n # callbacks: logging & model saving with checkpoints each fold\n callbacks = [get_lr_callback(BATCH_SIZE)]\n if LOGGING == True:\n callbacks += [get_checkpont_saver_callback(fold), \n LogPrintingCallback()]\n\n # build and train model\n model = get_model()\n model.fit(X_train[train_idx], y[train_idx], \n batch_size = BATCH_SIZE, \n epochs = EPOCHS, \n validation_data = (X_train[val_idx], y[val_idx]), \n callbacks = callbacks,\n verbose = 0) \n \n # evaluate\n print(\"Train:\", model.evaluate(X_train[train_idx], y[train_idx], verbose = 0, batch_size = BATCH_SIZE))\n print(\"Val:\", model.evaluate(X_train[val_idx], y[val_idx], verbose = 0, batch_size = BATCH_SIZE))\n \n ## Load best model to make pred\n model.load_weights('fold-%i.h5'%fold)\n train_preds[val_idx] = model.predict(X_train[val_idx],\n batch_size = BATCH_SIZE,\n verbose = 0)\n \n # predict on test set and average the predictions over all folds\n print(\"Predicting Test...\")\n test_preds += model.predict(X_test, batch_size = BATCH_SIZE, verbose = 0) / NFOLDS","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the next section we are going to use the ```train_preds``` to calculate the optimized sigma, which is a measure for certainty or rather uncertainty. We can do that, as we have both: the model's estimate and the real data. We subtract the lower quartile from the upper quartile (defined in the loss function) and average it.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## FIND OPTIMIZED STANDARD-DEVIATION\nsigma_opt = mean_absolute_error(y, train_preds[:,1])\nsigma_uncertain = train_preds[:,2] - train_preds[:,0]\nsigma_mean = np.mean(sigma_uncertain)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## PREPARE SUBMISSION FILE WITH OUR PREDICTIONS\nsub['FVC1'] = test_preds[:, 1]\nsub['Confidence1'] = test_preds[:,2] - test_preds[:,0]\n\n# get rid of unused data and show some non-empty data\nsubmission = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()\nsubmission.loc[~submission.FVC1.isnull()].head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.loc[~submission.FVC1.isnull(),'FVC'] = submission.loc[~submission.FVC1.isnull(),'FVC1']\n\nif sigma_mean < 70:\n submission['Confidence'] = sigma_opt\nelse:\n submission.loc[~submission.FVC1.isnull(),'Confidence'] = submission.loc[~submission.FVC1.isnull(),'Confidence1']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Okay, we made it! Let's finally check our stats and submit it!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.describe().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the following last step we overwrite our predictions with the known data from the orginal submission file to not waste known data.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"org_test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\nfor i in range(len(org_test)):\n submission.loc[submission['Patient_Week']==org_test.Patient[i]+'_'+str(org_test.Weeks[i]), 'FVC'] = org_test.FVC[i]\n submission.loc[submission['Patient_Week']==org_test.Patient[i]+'_'+str(org_test.Weeks[i]), 'Confidence'] = 70","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index = False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Thanks a lot for reading, I hope you could gain as much insights from reading this as I got from writing it. If you liked it, an upvote is highly appreciated. If you are interested in more content like this, feel free to follow me! ;)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}