{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## CMI 2023 - Sleepy Times\nSummary of the current algorithm:
\n - Use sliding 1 minute windows to create std(anglez) and median(enmo)\n - Create EAprod as the log() of the product of these two.\n - Create LIDS using usual 10-minute, 30-minute windows, as in the [Feature Eng ideas and LightGBM - CMI](https://www.kaggle.com/code/renatoreggiani/feature-eng-ideas-and-lightgbm-cmi) notebook. Here I've flipped its direction and changed the scale to have: 0.0=asleep, 1.0=awake. Also create LIDSprod as the log() of the product of LIDS, angleStd, and enmoMed with 0.0001 min value.\n - For each series:\n - - Fit its LIDSprod waveform to a one-day-period cosine wave to get an average middle-of-night offset amount.\n - - Validate nights by i) checking for repeat (\"filler\") data and ii) checking for not-worn, i.e., when max(anglezStd) is less than ~ 25 degrees.\n - - Estimate the average time asleep over the valid and not-too-short nights.\n - - Use the middle-of-night offset and average time asleep to calculate approx onset and wakeup times.\n - - Use these to select many-hour regions to fit for onset using a step down function and for wakeup with a step up function. Fit twice, re-centering for the second fit. The first fit is to a very smoothed times series, the second much less smoothed. Update the approx onsets and wakeups with these more accurate fit step locations and save ancillary information (\"features\") about the fit.\n - - Assign `score` values based on ML models applied to the onset and wakeup fit features. Models are created in the [CMI 2023 - Sleepy Times - phat ML](https://www.kaggle.com/code/dan3dewey/cmi-2023-sleepy-times-phat-ml) notebook using previously output event data (next step). \n - Assemble the results from all the series, and output a \"phat_frame\" which includes the ground truth and the fit features for each event.\n - Drop the extra columns and any events (rows) with NaN step values to form and output the submission.csv \n \nUsually each version of the notebook is run on all 277 test=train series\nand prints a single line summary for each series, including AP metrics (onset, wakeup, and average; score is ignored) and other values. This is useful to find hard-to-predict series. A set of plots are also made from the series' values.\n\nThe \"phat_frame\" csv is uploaded to a dataset so that the notebook CMI 2023 - Sleepy Times - phat ML can examine the features and also fit and save models which can predict the magnitude of each event's error from the fit features. One minus this predicted scaled error is used to assign the score in the steps above.\n \nRegarding some aspects of Average Precision in this competition see the notebook [CMI 2023 - AP metric and score](https://www.kaggle.com/code/dan3dewey/cmi-2023-precision-vs-recall)","metadata":{"papermill":{"duration":0.008947,"end_time":"2023-09-07T16:55:33.307435","exception":false,"start_time":"2023-09-07T16:55:33.298488","status":"completed"},"tags":[]}},{"cell_type":"code","source":"### Versions\n# v 1 Start, use pyarrow to read parquet format\n# . . . See v21 for these intermediate version notes.\n# v21 Based on discussion .../discussion/438438 :\n# Add simple detection of nights when the device is \"not worn\",\n# 0.043 reject the night if max(angleStd) during it is < 30 degrees.\n# v22 Included \"not worn\" detection in excluding nights.\n# 0.037 Params: 12, -5.5, 25, 0.5\n# v23 Get approx onsets, wakeups as in v22 (Params: 12, -5.5, 25, 0.50)\n# then focus in on the regions and fit a step-down or step-up function\n# 0.296 to determine each onset and wakeup.\n# . . .\n# v51 Use quantile(0.70) for EAprodMed. Use mean, quantile for LIDSprod\n# LB / AP and fit onset three times changing smoothing and scale between them.\n# xxx/420 and fit wakeups twice with smooth/scale changes also.\n# v52 432/458 Same as v51, median offsets adjusted -47.7, +45.6 steps\n# v53 xxx/457 Same as v52, but only do the first two onset fits (is 3rd doing much? - No)\n# v54 434/460 Same as v53, set wakeup fits(2) to the same params as the onset ones.\n# . . .\n# v73 503/460 Use v18 model: random forest 10 x 30 min leaf.\n# v74 499/460 Use v19 model: random forest 10 x 60 min leaf.\n# v75 495/460 Use v20 model: random forest 10 x 15 min leaf. Tune min leaf :)\n# v76 Repeat v18 model running on all test to get full phat_frame.\n","metadata":{"papermill":{"duration":0.018555,"end_time":"2023-09-07T16:55:33.334438","exception":false,"start_time":"2023-09-07T16:55:33.315883","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:27:50.256495Z","iopub.execute_input":"2023-10-30T19:27:50.256993Z","iopub.status.idle":"2023-10-30T19:27:50.288709Z","shell.execute_reply.started":"2023-10-30T19:27:50.256950Z","shell.execute_reply":"2023-10-30T19:27:50.287474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Observations:\n\n# - Train has 488 row groups and there are 277 series_ids in events:\n# --> each row group is 262185 rows = 3x5x7x11x227 (?!?) ~ 15.17 days\n# So, a given series can extend over 1(?), 2, 3, 4(?) row groups.\n# 488 x 15.17 days / 277 series = 26.7 days per series on average.\n# - The Train Events list has a 3 double onsets and a missing onset;\n# added code to correct these after reading the events file.\n# Corrected in a second update ~ 9/20.\n# - Train Events has 7254 pairs of onset/wakeup.\n# 4791 of onsets and 4794 of wakeups are real values (rest are NaN.)\n# - Looks like row_id in submission does not have to be 0,1,2,...\n# (In dropping the NaNs I removed row_ids and didn't renumber.)\n# - The sklearn.metrics.average_precision_score docs gives:\n# AP = Sum (R_n - R_(n-1)) * P_n and n is over the SCORE.\n# - Can detect \"not worn\" by checking if max(anglezStd) is less than ~ 25 degrees.)\n# - Understood \"score\"'s role in the metric:\n# It is the value that is thresholded to create a Precision-Recall curve and AP\n# --> for each event-type X tolerance combination; so 20 precision-recall curves\n# and 20 AP values which are then averaged.\n# Using a constant score=1.0 makes these P-R AP values easy to calculate:\n# AP for event X tolerance = Precision * Recall =\n# (TPs / Number of actual events) * (TPs / Number predicted events)\n# Score can be useful to indicate events that may be innaccurate,\n# e.g., if there is activity after (before) a predicted onset (wakeup).\n# - Noticed that 5 series' events do not have a night = 1 entry.\n# - Detect identical repeated \"filler\" data in the series:\n# check for the same (within a tolerance) anglez 24 hours before and after.","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:50.291382Z","iopub.execute_input":"2023-10-30T19:27:50.292247Z","iopub.status.idle":"2023-10-30T19:27:50.299862Z","shell.execute_reply.started":"2023-10-30T19:27:50.292205Z","shell.execute_reply":"2023-10-30T19:27:50.298332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Things to use","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n# Use this for parquet: https://arrow.apache.org/docs/python/parquet.html \nimport pyarrow\nimport pyarrow.parquet as pq\nimport pyarrow.dataset as pads\n\n# Fitting routines\nfrom scipy.optimize import leastsq, curve_fit\n\n# To load pre-trained models\nimport dill","metadata":{"papermill":{"duration":1.141685,"end_time":"2023-09-07T16:55:34.484112","exception":false,"start_time":"2023-09-07T16:55:33.342427","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:27:50.301888Z","iopub.execute_input":"2023-10-30T19:27:50.302361Z","iopub.status.idle":"2023-10-30T19:27:51.377068Z","shell.execute_reply.started":"2023-10-30T19:27:50.302321Z","shell.execute_reply":"2023-10-30T19:27:51.376103Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_ids_parquet(parquet_of_file):\n '''\n Return a list of all series ids in the train/test file.\n '''\n parquet_groups = parquet_of_file.metadata.num_row_groups\n # Go through all the row groups and assemble a list of series_id s\n all_ids = []\n for igroup in range(0,parquet_groups):\n this_group = parquet_of_file.read_row_group(igroup,\n ['series_id']).to_pandas()\n this_ids = this_group.series_id.dropna().unique()\n for serid in this_ids:\n if not(serid in all_ids):\n all_ids.append(serid)\n return all_ids","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.379821Z","iopub.execute_input":"2023-10-30T19:27:51.380683Z","iopub.status.idle":"2023-10-30T19:27:51.387708Z","shell.execute_reply.started":"2023-10-30T19:27:51.380644Z","shell.execute_reply":"2023-10-30T19:27:51.386748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_series_pads(a_series_id, pads_of_file):\n '''\n Read in the data for this series - omit the timestamp (need it?)\n '''\n # *** Parameter *** ***PARAMETER***\n enmo_clip_level = 1.0\n # *** ***\n this_series = pads_of_file.to_table(\n columns=['series_id','step','anglez','enmo'],\n filter=(pads.field('series_id') == a_series_id)).to_pandas()\n # Don't expect NaNs in anglez and enmo (train doesn't have them), but incase:\n this_series.anglez = this_series.anglez.fillna(0.0)\n # clip enmo to 1.0\n this_series.enmo = np.clip(this_series.enmo.fillna(0.0), 0.0, 1.)\n # If there remains an NaN (in series_id or step) drop that row\n return this_series.dropna(axis=0)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.389563Z","iopub.execute_input":"2023-10-30T19:27:51.389966Z","iopub.status.idle":"2023-10-30T19:27:51.403712Z","shell.execute_reply.started":"2023-10-30T19:27:51.389935Z","shell.execute_reply":"2023-10-30T19:27:51.402233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_series_events(a_series_id, all_series_events):\n the_id_events = all_series_events[(\n all_series_events.series_id == a_series_id)].copy()\n # Reset the index\n the_id_events = the_id_events.reset_index(drop=True)\n\n # Is 2 the number of the first night in the events?\n if (the_id_events.loc[0,'night'] == 2):\n ##print('Oops, starts at night 2, adding a night 1 with NaNs.\\n')\n # Duplicate night 2 after itself\n the_id_events.loc[1.25] = the_id_events.loc[0] \n the_id_events.loc[1.50] = the_id_events.loc[1]\n # Sort and Reset the index\n the_id_events = the_id_events.sort_index().reset_index(drop=True)\n # Make the first pair be night 1 with NaNs\n the_id_events.loc[0] = [this_id, 1, 'onset', np.NaN, np.NaN]\n the_id_events.loc[1] = [this_id, 1, 'wakeup', np.NaN, np.NaN]\n \n # Include the onset and wakeup values\n the_id_onsets = 1.0 * np.array(the_id_events[\n the_id_events['event'] == 'onset'].step)\n the_id_wakeups = 1.0 * np.array(the_id_events[\n the_id_events['event'] == 'wakeup'].step)\n return the_id_events, the_id_onsets, the_id_wakeups","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.405536Z","iopub.execute_input":"2023-10-30T19:27:51.406271Z","iopub.status.idle":"2023-10-30T19:27:51.422107Z","shell.execute_reply.started":"2023-10-30T19:27:51.406238Z","shell.execute_reply":"2023-10-30T19:27:51.421234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_features(this_series, wind_width):\n # *** Parameters *** ***PARAMETERs***\n lids_offset = 0.02\n lids_sum_width = 10*12 # 10 minutes as steps\n lids_smooth_width = 30*12 # 30 minutes as steps\n # *** ***\n \n # Calculate the std of anglez in the window\n this_series['anglezStd'] = this_series.anglez.rolling(\n wind_width+1, center=True, min_periods=3\n ).std().fillna(method='bfill')\n # and median of enmo in the window\n this_series['enmoMed'] = this_series.enmo.rolling(\n wind_width+1, center=True, min_periods=1).median()\n # Create the log of the product of enmoMed and anglezStd\n this_series['EAprod'] = np.log( this_series['enmoMed'] * \n this_series['anglezStd'] + 0.0001)\n # Create LIDS too\n this_series['LIDS'] = np.clip(this_series.enmo - lids_offset,0.0,np.inf).rolling(\n lids_sum_width+1, center=True, min_periods=1).sum()\n this_series['LIDS'] = (this_series.LIDS/(1.0 + this_series.LIDS)).rolling(\n lids_smooth_width+1, center=True, min_periods=1).mean()\n # LIDSprod is like EAprod\n this_series['LIDSprod'] = np.log( this_series['enmoMed'] *\n this_series['anglezStd'] *\n this_series['LIDS'] + 0.0001)\n # Get the series down-sampled to once every window width\n pruned_series = this_series.loc[\n range(int(wind_width/2), len(this_series), wind_width)].copy()\n #\n # Additional PRUNED features available:\n #\n # EAprod rolling \"median\" (quantile 0.70), width is in pruned (1 minute) units\n pruned_series[\"EAprodMed\"] = pruned_series.EAprod.rolling(\n 60+1, center=True, min_periods=1).quantile(0.70)\n # Pruned step is an unsigned 32 bit integer, convert to float\n pruned_series.step = 1.0 * pruned_series.step\n return pruned_series","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.425414Z","iopub.execute_input":"2023-10-30T19:27:51.425877Z","iopub.status.idle":"2023-10-30T19:27:51.440005Z","shell.execute_reply.started":"2023-10-30T19:27:51.425843Z","shell.execute_reply":"2023-10-30T19:27:51.438864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_func_sleep(x, dcval, ampl, phase):\n '''\n # Fit Function used to find the average phase of the middle-sleep times\n # by fitting a one-day period cosine wave to a feature, e.g. EAprod.\n '''\n # Fit function is \n # EAprod_hat = dcval + ampl *\n # ( -1*cosine(omega*(x - phase)) )\n # where\n # x is the input value of steps\n # dcval is offset (<0)\n # ampl (>0) is amplitude of the cosine curve\n # phase (>0) is the location in steps of the first night (minimum)\n # phase=0 gives x=0 at the middle of a daily sleep time\n omega = 2.0*np.pi/(12*60*24) # radians per step (2pi/day)\n out = dcval - ampl*np.cos(omega * (x - phase))\n return out\n\ndef fit_func_onset(x, dcval, amplpp, phase, slope):\n '''\n # Fit Function used to fit an Onset.\n '''\n # Fit function is \n ##slope = (8/(15*12))\n out = dcval - np.clip(slope*(x - phase), -amplpp/2, amplpp/2)\n return out\n\ndef fit_func_wakeup(x, dcval, amplpp, phase, slope):\n '''\n # Fit Function used to fit a Wakeup\n '''\n # Fit function is \n ##slope = (8/(15*12))\n out = dcval + np.clip(slope*(x - phase), -amplpp/2, amplpp/2)\n return out","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.441541Z","iopub.execute_input":"2023-10-30T19:27:51.441957Z","iopub.status.idle":"2023-10-30T19:27:51.460073Z","shell.execute_reply.started":"2023-10-30T19:27:51.441912Z","shell.execute_reply":"2023-10-30T19:27:51.458656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_first_night_phase(df_series, SHOW_PLOT=True):\n '''\n Find the location in steps of the average middle-of-nights\n by fitting the EAprod with a cos curve.\n Returns the phase in steps to the first middle of night.\n '''\n # Tried different features, EAprodMed looks good (small dips at filler nights.)\n fit_feature = \"EAprodMed\"\n try:\n dcval, ampl, phase = curve_fit( fit_func_sleep, \n df_series.step, \n df_series[fit_feature],\n bounds=([-9.0, 0.0, 0.0 ],\n [ 0.0, 6.0, 24*60*12 ] ) )[0]\n except:\n # if it fails use a default value for phase\n dcval=-5.0; ampl=4.0; phase = 10 * 60*12\n if SHOW_PLOT:\n print(\"\\n Fit dcval, ampl & phase: {:.3f}, {:.3f}, {:.0f}\\n\".\n format(dcval, ampl, phase, phase/(60*12)))\n print(\"The night phase is fit as {:.2f} hours from the start.\\n\".\n format(phase/(60*12)))\n plt.figure(figsize=(10,4))\n plt.plot(df_series.step, df_series[fit_feature], '.',\n markersize=2,c='lightblue')\n select_asleep = df_series.asleep > 0.0\n plt.plot(df_series[select_asleep].step, df_series[select_asleep][fit_feature], '.',\n markersize=2,c='magenta')\n plt.plot(df_series.step,\n fit_func_sleep(df_series.step, dcval, ampl, phase),c='red')\n plt.title(fit_feature + \" - \" + this_id)\n plt.savefig(this_id+\"_phase_fit_plot.png\")\n plt.show()\n return phase","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.461929Z","iopub.execute_input":"2023-10-30T19:27:51.462705Z","iopub.status.idle":"2023-10-30T19:27:51.478103Z","shell.execute_reply.started":"2023-10-30T19:27:51.462661Z","shell.execute_reply":"2023-10-30T19:27:51.476747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def setup_prelim_events(phase, this_duration, OUTPUT=True):\n # Assume these are available: this_series, pruned_series\n '''\n Setup onset and wakeup lists with 1.0s if the night is OK\n and NaNs if the night is not valid because:\n - it's from an identical filler waveforms\n - it's not worn\n '''\n # *** Parameter *** ***PARAMETER***\n notworn_angle = 25.0 # max of anglezStd when not worn\n # ***\n #\n predict_onsets = []\n predict_wakeups = []\n max_step_in_series = pruned_series.step.max()\n for inight in range(1,int(this_duration+0.5)+1):\n day_offset_steps = (inight-1)*(24*60*12)\n steps_6hrs = 6*60*12\n early_onset = round((day_offset_steps + phase - steps_6hrs), 1)\n middle_of_night = round((day_offset_steps + phase), 1)\n late_wakeup = round((day_offset_steps + phase + steps_6hrs), 1)\n # Check if too large for the series\n if late_wakeup > max_step_in_series:\n if OUTPUT:\n print(\"Night {} late_wakeup beyond series steps, not used.\".format(inight))\n break;\n # Check if the anglez value is identical 24 hours away\n night_is_filler = False\n # Check this with following night\n if int(middle_of_night)+24*60*12 < max_step_in_series:\n if np.isclose(this_series.loc[int(middle_of_night),'anglez'],\n this_series.loc[int(middle_of_night)+24*60*12,'anglez'],\n rtol=0.0, atol=1e-6):\n night_is_filler = True\n # Check this with previous night\n if int(middle_of_night)-24*60*12 > 0:\n if np.isclose(this_series.loc[int(middle_of_night),'anglez'],\n this_series.loc[int(middle_of_night)-24*60*12,'anglez'],\n rtol=0.0, atol=1e-6):\n night_is_filler = True\n # max anglezStd to check if not worn\n max_anglezStd = np.max(pruned_series[\n (pruned_series.step > early_onset) &\n (pruned_series.step < late_wakeup)].anglezStd)\n # Is this not a filler or not-worn?\n predict_onset = 1.0 ; predict_wakeup = 1.0 # OK values\n if night_is_filler:\n if OUTPUT:\n print(\"Night {} is not included - it's repeating filler.\".format(inight))\n predict_onset = np.nan ; predict_wakeup = np.nan\n elif (max_anglezStd < notworn_angle):\n if OUTPUT:\n print(\"Night {} is not included - not worn.\".format(inight))\n predict_onset = np.nan ; predict_wakeup = np.nan\n predict_onsets.append(predict_onset)\n predict_wakeups.append(predict_wakeup)\n # all done, return the lists\n return predict_onsets, predict_wakeups\n","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.483921Z","iopub.execute_input":"2023-10-30T19:27:51.484711Z","iopub.status.idle":"2023-10-30T19:27:51.505172Z","shell.execute_reply.started":"2023-10-30T19:27:51.484666Z","shell.execute_reply":"2023-10-30T19:27:51.503567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def assign_ave_events(phase, predict_onsets, predict_wakeups, OUTPUT=True):\n '''\n Find the average sleep-time per valid night and use that and the phase\n to assign onset and wakeup values.\n '''\n # *** Parameters *** ***PARAMETERs*** \n ave_feature = 'EAprodMed'\n feature_low = -6.0 # threshold level for enmoMed \n short_night_hours = 3.0 # reject nights shorter than this, hours\n # *** ***\n\n # Estimate average sleep per night, ignore nights that\n # - have little sleep time\n steps_asleep = 0\n this_nights_asleep = 0\n for inight in range(1,len(predict_onsets)+1):\n if not np.isnan(predict_onsets[inight-1]):\n day_offset_steps = (inight-1)*(24*60*12)\n steps_6hrs = 6*60*12\n early_onset = round((day_offset_steps + phase - steps_6hrs), 1)\n late_wakeup = round((day_offset_steps + phase + steps_6hrs), 1)\n # Rough amount of sleep\n steps_featlow = sum( (pruned_series.step > early_onset) &\n (pruned_series.step < late_wakeup) &\n (pruned_series[ave_feature] <= feature_low) ) * wind_width\n if (steps_featlow < short_night_hours * 60*12):\n if OUTPUT:\n print(\"Night {} is not included - too short.\".format(inight))\n predict_onsets[inight-1] = np.nan\n predict_wakeups[inight-1] = np.nan\n else:\n # include this night in sleep time estimate\n this_nights_asleep += 1\n steps_asleep += steps_featlow\n # Done with all nights\n if this_nights_asleep == 0:\n # default if none were included\n ave_steps_asleep = 5*60*12\n else:\n ave_steps_asleep = steps_asleep / this_nights_asleep\n if OUTPUT:\n print(\"\\nApprox average sleep time per night over\",this_nights_asleep,\n \"nights is {:.3f}\".format(ave_steps_asleep/(60*12)),\n \"hours ({:.0f} steps).\".format(ave_steps_asleep))\n \n # Generate onset and wakeup values\n # based on the phase and ave_steps_asleep.\n # Skip nights that are already assigned NaNs.\n half_steps_asleep = ave_steps_asleep/2 \n\n for inight in range(1,len(predict_onsets)+1):\n day_offset_steps = (inight-1)*(24*60*12)\n # Approx/average events for this night\n predict_onset = round((day_offset_steps + phase - half_steps_asleep),1)\n predict_wakeup = round((day_offset_steps + phase + half_steps_asleep),1)\n if not np.isnan(predict_onsets[inight-1]):\n predict_onsets[inight-1] = predict_onset\n predict_wakeups[inight-1] = predict_wakeup\n # Done with the nights\n predict_onsets = 1.0 * np.array(predict_onsets) # float\n predict_wakeups = 1.0 * np.array(predict_wakeups) # float\n \n return predict_onsets, predict_wakeups, this_nights_asleep","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.508366Z","iopub.execute_input":"2023-10-30T19:27:51.508876Z","iopub.status.idle":"2023-10-30T19:27:51.532792Z","shell.execute_reply.started":"2023-10-30T19:27:51.508829Z","shell.execute_reply":"2023-10-30T19:27:51.531062Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def pad_same_length(arry1, arry2, OUTPUT=True):\n '''\n Add NaNs as appropriate so that the two arrays have the same length.\n '''\n len1 = len(arry1)\n len2 = len(arry2)\n diff = abs(len1 - len2)\n if len2 > len1 :\n arry1 = np.append(arry1, diff*[np.NaN])\n if OUTPUT:\n print(\"\\nSecond has more ({}) than First ({}): padded First.\".\n format(len2,len1))\n if len1 > len2 :\n arry2 = np.append(arry2, diff*[np.NaN])\n if OUTPUT:\n print(\"\\nFirst has more ({}) than Second ({}): padded Second.\".\n format(len2,len1))\n return arry1, arry2","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.535089Z","iopub.execute_input":"2023-10-30T19:27:51.535638Z","iopub.status.idle":"2023-10-30T19:27:51.552952Z","shell.execute_reply.started":"2023-10-30T19:27:51.535594Z","shell.execute_reply":"2023-10-30T19:27:51.551769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# These two functions are defined further below just before they are used.\n\n# def refine_onset_values()\n\n# def refine_wakeup_values()","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.554972Z","iopub.execute_input":"2023-10-30T19:27:51.555657Z","iopub.status.idle":"2023-10-30T19:27:51.570670Z","shell.execute_reply.started":"2023-10-30T19:27:51.555610Z","shell.execute_reply":"2023-10-30T19:27:51.569224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_AP_values(vals_pred, vals_true, OUTPUT=False):\n '''\n # Calculate the AP for already matched values, which may contain NaNs.\n # This assumes that all of the score values are the same, e.g., 1.0.\n '''\n # These countings will handle NaN values OK\n num_input = len(vals_pred) # equals len(vals_true) also.\n num_pred_pos = max([sum(abs(vals_pred) > -1), 1])\n num_actual_pos = max([sum(abs(vals_true) > -1), 1])\n tol_levels = [12,36,60,90,120,150,180,240,300,360]\n prec_values = np.zeros(len(tol_levels))\n recall_values = np.zeros(len(tol_levels))\n AP_values = np.zeros(len(tol_levels))\n for itol, tol_level in enumerate(tol_levels):\n num_TP = sum(abs(vals_pred - vals_true) < tol_level)\n prec = num_TP/num_pred_pos\n recall = num_TP/num_actual_pos\n # Calc the average-precision assuming all scores are the same, e.g., 1.0,\n # so there is only one point in the P-R curve (plus the recall=0, prec=1 point.)\n AP_values[itol] = prec * recall\n prec_values[itol] = prec\n recall_values[itol] = recall\n # Show output info:\n if OUTPUT:\n print(\"calc_AP_values: {} of {} predicted pos, {} of {} are actually pos\".\n format(num_pred_pos, num_input, num_actual_pos, num_input))\n ##print(\"Precision values for each tolerance width:\\n {}\".format(prec_values))\n print(\" Recall values for each tolerance width:\\n {}\".format(recall_values))\n print(\"Recall average: {:.5f}\".format(sum(recall_values)/len(recall_values)))\n print(\" AP values for each tolerance width:\\n {}\".format(AP_values))\n print(\"AP average: {:.5f}\\n\".format(sum(AP_values)/len(AP_values)))\n # Return the array of values for each tolerance level\n return AP_values","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.573237Z","iopub.execute_input":"2023-10-30T19:27:51.574712Z","iopub.status.idle":"2023-10-30T19:27:51.609422Z","shell.execute_reply.started":"2023-10-30T19:27:51.574550Z","shell.execute_reply":"2023-10-30T19:27:51.607825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The Input Files","metadata":{}},{"cell_type":"code","source":"# Show the input file names\nimport os\nfor dirname, _, filenames in os.walk('../input'):\n for filename in filenames:\n print(os.path.join(dirname, filename))\n \n# file name variables (use\"/\"\")\nsample_file = \"../input/child-mind-institute-detect-sleep-states/sample_submission.csv\"\ntest_file = \"../input/child-mind-institute-detect-sleep-states/test_series.parquet\"\ntrain_events_file = \"../input/child-mind-institute-detect-sleep-states/train_events.csv\"\ntrain_file = \"../input/child-mind-institute-detect-sleep-states/train_series.parquet\"","metadata":{"papermill":{"duration":0.03276,"end_time":"2023-09-07T16:55:34.525576","exception":false,"start_time":"2023-09-07T16:55:34.492816","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:27:51.617151Z","iopub.execute_input":"2023-10-30T19:27:51.621505Z","iopub.status.idle":"2023-10-30T19:27:51.648046Z","shell.execute_reply.started":"2023-10-30T19:27:51.621423Z","shell.execute_reply":"2023-10-30T19:27:51.646479Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Read the Test file - is it Real?","metadata":{}},{"cell_type":"code","source":"# Ways to access the test series file\nparquet_test = pq.ParquetFile(test_file)\npads_test = pads.dataset(test_file)\n\n# Info about the test file of series\nprint(parquet_test.schema_arrow)\ntest_groups = parquet_test.metadata.num_row_groups\ntest_rows = parquet_test.metadata.num_rows\nprint(\"\\nThe provided Test has\", str(test_groups), \"row groups,\\nwith\",\n test_rows,\"rows.\\n\")\n\n# Is this the real test data?\nif test_rows > 1000:\n REAL_TEST = True\n print(\"This is the real test data! Will never see this output :)\")\nelse:\n # Not the real test data - use train instead.\n REAL_TEST = False\n print(\"This is useless test data...\")\n print(\"Use the train data as test.\")\n # redefine these using train file\n parquet_test = pq.ParquetFile(train_file)\n pads_test = pads.dataset(train_file)\n test_groups = parquet_test.metadata.num_row_groups\n test_rows = parquet_test.metadata.num_rows\n print(\"Now Test has\", str(test_groups), \"row groups,\\nwith\",\n test_rows,\"rows.\")\n","metadata":{"papermill":{"duration":0.241078,"end_time":"2023-09-07T16:55:34.775530","exception":false,"start_time":"2023-09-07T16:55:34.534452","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:27:51.650453Z","iopub.execute_input":"2023-10-30T19:27:51.657930Z","iopub.status.idle":"2023-10-30T19:27:51.732890Z","shell.execute_reply.started":"2023-10-30T19:27:51.657855Z","shell.execute_reply":"2023-10-30T19:27:51.731764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Read the Ground-truth Events file\nOf course this information won't be available for the real test data.","metadata":{}},{"cell_type":"code","source":"# Read the events file\ntrain_events = pd.read_csv(train_events_file)\nprint(\"\\nContains\", len(train_events.series_id.unique()),\n \"unique series_id values.\\n\")\n\n# Train Events has 14508/2 = 7254 pairs of onset/wakeup.\n# 4791 of onsets and 4794 of wakeups are real values (rest are NaN.)\n# --> 66% of events have non-NaN values.\ntrain_events.info()","metadata":{"papermill":{"duration":0.07404,"end_time":"2023-09-07T16:55:34.952633","exception":false,"start_time":"2023-09-07T16:55:34.878593","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:27:51.734416Z","iopub.execute_input":"2023-10-30T19:27:51.735376Z","iopub.status.idle":"2023-10-30T19:27:51.825031Z","shell.execute_reply.started":"2023-10-30T19:27:51.735336Z","shell.execute_reply":"2023-10-30T19:27:51.823965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Versions v42 and earlier had a check here to see that\n# the onsets and wakeups strictly alternate - they do, so leave it out.\n\n# Some series do not have night=1 but do have night=2:\nprint(sum(train_events.night == 1),\n sum(train_events.night == 2),\n sum(train_events.night == 3),\n sum(train_events.night == 4),\n sum(train_events.night == 5))\n# A night 1 will be added when individual series are extracted.","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:27:51.826498Z","iopub.execute_input":"2023-10-30T19:27:51.827642Z","iopub.status.idle":"2023-10-30T19:27:51.854192Z","shell.execute_reply.started":"2023-10-30T19:27:51.827599Z","shell.execute_reply":"2023-10-30T19:27:51.852900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ \n\n## Train: EDA and New Features","metadata":{}},{"cell_type":"markdown","source":"### Read the Train Series file to Get All series_id values","metadata":{}},{"cell_type":"code","source":"# Ways to access the train series file\nparquet_train = pq.ParquetFile(train_file)\npads_train = pads.dataset(train_file)\n\n# Info about the train file of series\nprint(parquet_train.schema_arrow)\ntrain_rows = parquet_train.metadata.num_rows\nprint(\"\\nTrain has\", train_rows, \"rows.\\n . . . working . . .\")\n\n# Get a list of all the series_id in Train\ntrain_ids = get_ids_parquet(parquet_train)\n \n# Gives 277 ids from the whole file - agrees with the events data. \nprint(\"Number of series ids:\",len(train_ids))","metadata":{"papermill":{"duration":0.039418,"end_time":"2023-09-07T16:55:34.869403","exception":false,"start_time":"2023-09-07T16:55:34.829985","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:43:33.893608Z","iopub.execute_input":"2023-10-30T19:43:33.894032Z","iopub.status.idle":"2023-10-30T19:44:09.832936Z","shell.execute_reply.started":"2023-10-30T19:43:33.893995Z","shell.execute_reply":"2023-10-30T19:44:09.832076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Read in one of the series","metadata":{}},{"cell_type":"code","source":"# The list of train series ids is train_ids, pick one by index:\n##this_id = train_ids[276]\n\n# OR\n\n# Various ones Useful for algorithm checking:\n\n# weird:\n##this_id = '31011ade7c0a' # [55] AP(v63) = 0.227 32/11 nights, 12-33 OK,but not given???\n\n# These gave double_scalar Warning for ave_ratio calc - added check for empty select.\n##this_id = 'a4e48102f402' # [169] Night 24 multiple empty selects\n##this_id = 'a9a2f7fac455' # [174] Night 37 ends with non-selects\n##this_id = 'fa149c3c4bde' # [272] # Night 24 ends early,no steps in p90 wakeup\n \n# medium APs:\n##this_id = 'e69aff66e0cb' # [251] AP(v63) = 0.481 21/19 nights\n##this_id = '0d0ad1e77851' # [12] AP(v63) = 0.597 20/24 nights\n##this_id = 'c6788e579967' # [202] AP(v63) = 0.432 23/23 nights\n##this_id = '1f96b9668bdf' # [39] AP(v43) = 0.392 15/34 nights, added Night 1\nthis_id = '1087d7b0ff2e' # [19] AP(v63) = 0.305 24/24 nights\n# low APs:\n##this_id = '0402a003dae9' # [2] AP(v63) = 0.269 18/24 nights\n##this_id = '154fe824ed87' # [27] AP(v63) = 0.455 30/32 nights\n##this_id = 'b750c8c1556c' # [185] AP(v63) = 0.289 15/22 nights\n##this_id = 'fe90110788d2' # [276] AP(v63) = 0.349 33/34 nights\n# really low:\n##this_id = 'b7fc34995d0f' # [186] AP(v63) = 0.147 3/34 nights, identical after 5\n\nprint(\"Looking at series id:\", this_id)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:09.834516Z","iopub.execute_input":"2023-10-30T19:44:09.835036Z","iopub.status.idle":"2023-10-30T19:44:09.842632Z","shell.execute_reply.started":"2023-10-30T19:44:09.835006Z","shell.execute_reply":"2023-10-30T19:44:09.841065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"this_series = read_series_pads(this_id,pads_train)\n# Length in days:\nthis_duration = (this_series.step.max() - this_series.step.min())/(12*60*24)\n\nthis_series","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:09.843881Z","iopub.execute_input":"2023-10-30T19:44:09.844184Z","iopub.status.idle":"2023-10-30T19:44:13.759700Z","shell.execute_reply.started":"2023-10-30T19:44:09.844158Z","shell.execute_reply":"2023-10-30T19:44:13.758411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Get the known Events for this series and add Asleep column\nUse these to add an \"asleep\" column to series to help visualize, develop, and test algorithms.","metadata":{}},{"cell_type":"code","source":"# Find the events and event values for this series\n(this_events, actual_onsets,\n actual_wakeups) = get_series_events(this_id, train_events)\n# Doing this now so that \"asleep\" column can be made.\n# The known event values will be used further below when fitting, etc.\n\n# Summary of duration and nights\nthis_num_nights = this_events[\"night\"].max()\nthis_num_rated = int(sum(this_events.step > 0.0)/2)\nprint(\"\\nThis series covers\", this_duration, \"days.\\n\\n\" +\n \"The Events File lists\", this_num_nights, \"nights;\", \n this_num_rated, \"of them have onsets/wakeups identified.\\n\")\n\n# Show the events\nthis_events","metadata":{"papermill":{"duration":0.048026,"end_time":"2023-09-07T16:55:35.309860","exception":false,"start_time":"2023-09-07T16:55:35.261834","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:44:13.763389Z","iopub.execute_input":"2023-10-30T19:44:13.763868Z","iopub.status.idle":"2023-10-30T19:44:13.799776Z","shell.execute_reply.started":"2023-10-30T19:44:13.763834Z","shell.execute_reply":"2023-10-30T19:44:13.798238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Add a column \"asleep\": 0 = awake, 1 = asleep - can only do this for train\n# Creating \"asleep\" to use it for color-coding only.\nthis_series[\"asleep\"] = 0\nminstep = this_series.step.min()\nmaxstep = this_series.step.max()\n# Go through all the events and apply each (in range) to future steps\nfor ievent in range(len(this_events)):\n evtype = this_events.iloc[ievent,2] # *** Uses actual column numbers ***\n evstep = this_events.iloc[ievent,3] #\n ##print(evstep, evtype)\n # skip NaNs and ones out of range\n if evstep == evstep and evstep >= minstep and evstep <= maxstep:\n evstep = int(evstep)\n # if onset then this and future are asleep\n ##print(evstep, evtype)\n atandabove = this_series.step >= evstep\n if evtype == \"onset\":\n this_series.loc[atandabove,['asleep']] = 1\n else:\n this_series.loc[atandabove,['asleep']] = 0","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:13.801180Z","iopub.execute_input":"2023-10-30T19:44:13.801561Z","iopub.status.idle":"2023-10-30T19:44:13.934509Z","shell.execute_reply.started":"2023-10-30T19:44:13.801529Z","shell.execute_reply":"2023-10-30T19:44:13.933292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show the anglez and enmo values vs step, color coded by asleep.\nthis_series.plot.scatter('step', 'anglez', s=2, figsize=(10,3), title=this_id,\n c='asleep', colormap='cool', colorbar=False)\nplt.savefig(this_id+\"_anglez.png\")\nplt.show()\n\n# Note enmo has been clipped at 1.0\nthis_series.plot.scatter('step', 'enmo', s=2, figsize=(10,3), title=this_id, \n c='asleep', colormap='cool', colorbar=False)\nplt.savefig(this_id+\"_enmo.png\")\nplt.show()","metadata":{"papermill":{"duration":1.242122,"end_time":"2023-09-07T16:55:37.281914","exception":false,"start_time":"2023-09-07T16:55:36.039792","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:44:13.936151Z","iopub.execute_input":"2023-10-30T19:44:13.937166Z","iopub.status.idle":"2023-10-30T19:44:47.526780Z","shell.execute_reply.started":"2023-10-30T19:44:13.937121Z","shell.execute_reply":"2023-10-30T19:44:47.525858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Add New Features and Down-sample too\n#### Added: anglezStd, enmoMed, EAprod, LIDS, LIDSprod\nThe first three are created using a sliding window of given width. LIDS is from the Random Forest paper in Nature, but flipped to have asleep is 0.0 and awake is near 1.0 .","metadata":{}},{"cell_type":"code","source":"# Adds features to this_series and also create a pruned_series\n# *** Parameter ***\nwind_width = 12 # one minute ***PARAMETER***\n# *** ***\npruned_series = add_features(this_series, wind_width)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:47.528108Z","iopub.execute_input":"2023-10-30T19:44:47.528784Z","iopub.status.idle":"2023-10-30T19:44:47.778131Z","shell.execute_reply.started":"2023-10-30T19:44:47.528749Z","shell.execute_reply":"2023-10-30T19:44:47.776658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Look at the new features in a selected range of days","metadata":{}},{"cell_type":"code","source":"# Show values down-sampled to once per - rolling window time\n\n# Can select a range of times to narrow in on\nday_start = 7\nday_end = 17\n\nstart_step = int(day_start * 24*60*12)\nend_step = int(min([ day_end * 24*60*12, pruned_series.step.max()]))\nsubselect = ((pruned_series.step > start_step) & (pruned_series.step < end_step))\n\n# Plot of enmoMed\npruned_series.loc[subselect,].plot.scatter('step', 'anglezStd', s=2,\n figsize=(10,3), title=this_id,\n c='asleep', colormap='cool', colorbar=False)\nplt.show()\n\n# Plot of anglezStd\npruned_series.loc[subselect,].plot.scatter('step', 'enmoMed', s=2,\n figsize=(10,3), title=this_id,\n c='asleep', colormap='cool', colorbar=False,\n ylim=(0.0,pruned_series.enmoMed.max()/2))\nplt.show()\n\n# Plot of EAprod \npruned_series.loc[subselect,].plot.scatter('step', 'EAprod', s=2,\n figsize=(10,3), title=this_id,\n c='asleep', colormap='cool', colorbar=False)\nplt.show()\n\n# Plot of EAprodMed\npruned_series.loc[subselect,].plot.scatter('step', 'EAprodMed', s=2,\n figsize=(10,3), title=this_id,\n c='asleep', colormap='cool', colorbar=False)\nplt.show()\n\n# Plot of LIDS \npruned_series.loc[subselect,].plot.scatter('step', 'LIDS', s=2,\n figsize=(10,3), title=this_id,\n c='asleep', colormap='cool', colorbar=False)\nplt.show()\n# Plot of LIDSprod\npruned_series.loc[subselect,].plot.scatter('step', 'LIDSprod', s=2,\n figsize=(10,3),title=this_id, \n c='asleep', colormap='cool', colorbar=False)\nplt.show()","metadata":{"papermill":{"duration":0.703342,"end_time":"2023-09-07T16:55:38.180450","exception":false,"start_time":"2023-09-07T16:55:37.477108","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-10-30T19:44:47.779931Z","iopub.execute_input":"2023-10-30T19:44:47.780285Z","iopub.status.idle":"2023-10-30T19:44:50.968264Z","shell.execute_reply.started":"2023-10-30T19:44:47.780252Z","shell.execute_reply":"2023-10-30T19:44:50.966960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Experiment with other features based on these\nif True:\n \n if False:\n # EAprodMed using quantile instead of median\n pruned_series[\"test_feature\"] = pruned_series.EAprod.rolling(\n 60+1, center=True, min_periods=1).quantile(0.70)\n \n # For fitting coarse onset/wakeup location: mean of a quantile\n if True:\n pruned_series[\"test_feature\"] = pruned_series['LIDSprod'].rolling(\n 30+1, center=True, min_periods=1).quantile(0.80)\n pruned_series[\"test_feature\"] = pruned_series[\"test_feature\"].rolling(\n 120+1, center=True, min_periods=1).mean()\n \n \n pruned_series.loc[subselect,].plot.scatter('step', 'test_feature', s=2,\n figsize=(10,3),title=this_id, \n c='asleep', colormap='cool', colorbar=False)\n plt.savefig(this_id+\"_coarse_feature.png\")\n plt.show()\n\n # Histogram of the test_feature pruned values in this series\n pruned_series['test_feature'].plot.hist(figsize=(6,3),bins=100,\n title=\"test_feature - \"+this_id)\n plt.ylim(0,2000)\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:50.970304Z","iopub.execute_input":"2023-10-30T19:44:50.970747Z","iopub.status.idle":"2023-10-30T19:44:52.479488Z","shell.execute_reply.started":"2023-10-30T19:44:50.970709Z","shell.execute_reply":"2023-10-30T19:44:52.478348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Histograms of the EAprod and LIDSprod pruned values in this series\npruned_series['EAprodMed'].plot.hist(figsize=(6,3),bins=100,\n title=\"EAprod - \"+this_id)\nplt.ylim(0,2000)\nplt.show()\n\npruned_series['LIDSprod'].plot.hist(figsize=(6,3),bins=100,\n title=\"LIDSprod - \"+this_id)\nplt.ylim(0,2000)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:52.484516Z","iopub.execute_input":"2023-10-30T19:44:52.485598Z","iopub.status.idle":"2023-10-30T19:44:53.571349Z","shell.execute_reply.started":"2023-10-30T19:44:52.485550Z","shell.execute_reply":"2023-10-30T19:44:53.570280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Algorithm to determine onset and wakeup times","metadata":{}},{"cell_type":"markdown","source":"### Predict approximate (average) onsets and wakeups","metadata":{}},{"cell_type":"code","source":"# Find the location of the first/average middle-night in steps. Uses EAprod Med.\nphase = fit_first_night_phase(pruned_series, SHOW_PLOT=True)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:53.572866Z","iopub.execute_input":"2023-10-30T19:44:53.573190Z","iopub.status.idle":"2023-10-30T19:44:54.354202Z","shell.execute_reply.started":"2023-10-30T19:44:53.573162Z","shell.execute_reply":"2023-10-30T19:44:54.352918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create initial lists of event values for this series:\npredict_onsets, predict_wakeups = setup_prelim_events(\n phase, this_duration, OUTPUT=True)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:54.355749Z","iopub.execute_input":"2023-10-30T19:44:54.356109Z","iopub.status.idle":"2023-10-30T19:44:54.398800Z","shell.execute_reply.started":"2023-10-30T19:44:54.356078Z","shell.execute_reply":"2023-10-30T19:44:54.397454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Assign values to the valid nights based on average time asleep and phase.\npredict_onsets, predict_wakeups, this_nights_asleep = assign_ave_events(phase,\n predict_onsets, predict_wakeups, OUTPUT=True)\n\n# Get the known ave sleep per valid night\nprint(\"Actual average sleep per night over\",this_num_rated,\n \"rated nights is {:.3f}\".format(\n this_series['asleep'].sum()/(this_num_rated * 60*12)),\n \"hours.\\n\")","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:54.400278Z","iopub.execute_input":"2023-10-30T19:44:54.400994Z","iopub.status.idle":"2023-10-30T19:44:54.542145Z","shell.execute_reply.started":"2023-10-30T19:44:54.400949Z","shell.execute_reply":"2023-10-30T19:44:54.541049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Already read in this series' known/actual events values --\n# both predict and known come in pairs of onset, wakeup.\n# Pad as needed to have the same lengths\npredict_onsets, actual_onsets = pad_same_length(predict_onsets, actual_onsets)\npredict_wakeups, actual_wakeups = pad_same_length(predict_wakeups, actual_wakeups)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:54.543378Z","iopub.execute_input":"2023-10-30T19:44:54.543711Z","iopub.status.idle":"2023-10-30T19:44:54.551847Z","shell.execute_reply.started":"2023-10-30T19:44:54.543683Z","shell.execute_reply":"2023-10-30T19:44:54.550299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fit each onset and wakeup region","metadata":{}},{"cell_type":"code","source":"def get_ave_ratio(the_section, low_min, high_min, the_feature, dcval,amplpp,phase):\n '''\n Measure the average value of the feature between low_min and high_min (in minutes)\n and return it's ratio to amplpp.\n '''\n select_vals = the_section[(the_section.step > phase+low_min*12) &\n (the_section.step < phase+high_min*12)][the_feature]\n if len(select_vals) == 0: # catch empty selects\n ave_ratio = 0.0\n else:\n ave_ratio = sum(np.clip(select_vals - (dcval - amplpp/2),\n 0, amplpp)) / (amplpp*len(select_vals))\n return ave_ratio","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:54.553376Z","iopub.execute_input":"2023-10-30T19:44:54.554052Z","iopub.status.idle":"2023-10-30T19:44:54.563117Z","shell.execute_reply.started":"2023-10-30T19:44:54.554017Z","shell.execute_reply":"2023-10-30T19:44:54.561868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def refine_onset_values(predict_onsets, known_onsets=[]):\n '''\n Go through the onset regions and refine the onset locations by fitting.\n Also returns fit features.\n '''\n nfits = 2\n # Set up any features to save - dictionary of arrays\n default_nans = len(predict_onsets) * [np.NaN]\n fit_feats = {\n 'amplpp': np.array(default_nans),\n 'dcval': np.array(default_nans),\n 'rampsteps': np.array(default_nans),\n 'ratm90': np.array(default_nans),\n 'ratm60': np.array(default_nans),\n 'ratm25': np.array(default_nans),\n 'ratm14': np.array(default_nans),\n 'ratp14': np.array(default_nans),\n 'ratp25': np.array(default_nans),\n 'ratp60': np.array(default_nans),\n 'ratp90': np.array(default_nans)\n }\n for ionset in range(len(predict_onsets)):\n # continue if not NaN\n this_onset = predict_onsets[ionset]\n if this_onset == this_onset:\n for ifit in range(nfits):\n # Set the parameters for this fit\n fit_factor = [1.0, 0.7, 0.7][ifit]\n reg_left = 6 * fit_factor\n reg_right = 6 * fit_factor\n fit_left = 5 * fit_factor\n fit_right = 5 * fit_factor\n slope_low = [0.05*(8/(15*12)), 0.7*(8/(15*12)), 0.7*(8/(15*12))][ifit]\n slope_high = (8/(15*12))\n quant_width = [30, 12, 12][ifit] # even value, 1 will be added\n mean_width = [120, 4, 0][ifit] # even value, 1 will be added\n\n # Select a range around the current location\n select_series = ((pruned_series.step > \n this_onset - reg_left*60*12) &\n (pruned_series.step <\n this_onset + reg_right*60*12))\n # extract that section\n this_section = pruned_series[select_series].copy()\n # Create the feature to fit\n this_section[\"fit_feature\"] = this_section['LIDSprod'].rolling(\n quant_width+1, center=True, min_periods=1).quantile(0.80)\n this_section[\"fit_feature\"] = this_section[\"fit_feature\"].rolling(\n mean_width+1, center=True, min_periods=1).mean()\n # Fit the Onset: Note: amplpp lower is 4.0\n lower_bounds = [-7.0, 4.0, this_onset - fit_left*60*12, slope_low]\n upper_bounds = [ 0.0, 12.0, this_onset + fit_right*60*12, slope_high]\n try:\n dcval, amplpp, phase, slope = curve_fit( fit_func_onset, \n this_section.step, this_section.fit_feature,\n bounds=(lower_bounds, upper_bounds) )[0]\n except:\n # signal a problem fitting with the dcval and amplpp\n dcval=0.0; amplpp=2.0; phase = this_onset; slope=8/(15*12)\n break # stop any further fitting of this event\n # get the first-fit \"ramp steps\" = amplpp / slope\n if ifit == 0:\n rampsteps = amplpp / slope\n # do a repeat fit centering the transition\n this_onset = phase\n # Done fitting this night.\n #\n # Get/make the feature values for this fit...\n ratm90 = 1.0 - get_ave_ratio(this_section, -90, -60, 'fit_feature', dcval,amplpp,phase)\n ratm60 = 1.0 - get_ave_ratio(this_section, -60, -25, 'fit_feature', dcval,amplpp,phase)\n ratm25 = 1.0 - get_ave_ratio(this_section, -25, -14, 'fit_feature', dcval,amplpp,phase)\n ratm14 = 1.0 - get_ave_ratio(this_section, -14, -4, 'fit_feature', dcval,amplpp,phase)\n ratp14 = get_ave_ratio(this_section, 4, 14, 'fit_feature', dcval,amplpp,phase)\n ratp25 = get_ave_ratio(this_section, 14, 25, 'fit_feature', dcval,amplpp,phase)\n ratp60 = get_ave_ratio(this_section, 25, 60, 'fit_feature', dcval,amplpp,phase)\n ratp90 = get_ave_ratio(this_section, 60, 90, 'fit_feature', dcval,amplpp,phase)\n #\n # Update this onset value in the array.\n # Small, constant adjustment based on previous version(s) results.\n # To get to 0, add on the negative of the previous Overall.\n onset_correction = 37.4 - 47.7 - 10.2 + 7.0 # in steps\n new_onset = round(phase + onset_correction, 1)\n # If the amplpp is at lower value, 4.00, then NaN the event\n if amplpp < 4.01:\n predict_onsets[ionset] = np.NaN\n else:\n predict_onsets[ionset] = new_onset\n # Save the features values for this night's fit\n fit_feats['amplpp'][ionset] = amplpp\n fit_feats['dcval'][ionset] = dcval\n fit_feats['rampsteps'][ionset] = rampsteps\n fit_feats['ratm90'][ionset] = ratm90\n fit_feats['ratm60'][ionset] = ratm60\n fit_feats['ratm25'][ionset] = ratm25\n fit_feats['ratm14'][ionset] = ratm14\n fit_feats['ratp14'][ionset] = ratp14\n fit_feats['ratp25'][ionset] = ratp25\n fit_feats['ratp60'][ionset] = ratp60\n fit_feats['ratp90'][ionset] = ratp90\n #\n # Make a plot with comparison to actual onsets, if given\n if len(known_onsets) > 0:\n print(\"\\nsteps error = {:.0f} amplpp = {:.2f}\".\n format(predict_onsets[ionset] - known_onsets[ionset], amplpp) +\n \" dcval = {:.2f} rampsteps = {:.1f}\".\n format(dcval, rampsteps))\n print(\"\\n ratios: m90= {:.3f}, m60= {:.3f}, m25= {:.3f}, m14= {:.3f},\".\n format(ratm90, ratm60, ratm25, ratm14),\n \"p14= {:.3f}, p25= {:.3f}, p60= {:.3f}, p90= {:.3f}\\n\".\n format(ratp14, ratp25, ratp60, ratp90))\n # Show time axis as minutes from the mid-transition (phase) value\n # Adjust the event locations too.\n min_from_mid = (this_section.step - phase)/12\n known_from_mid = (known_onsets[ionset]-phase)/12\n new_from_mid = (new_onset - phase)/12\n plt.figure(figsize=(6,3))\n plt.plot(min_from_mid, this_section.fit_feature, '.',\n markersize=2,c='gray')\n # Plot the known event location\n plt.plot(2*[known_from_mid],[dcval,2],c=\"green\")\n # and the fit (using the unadjusted phase)\n plt.plot(min_from_mid,\n fit_func_onset(this_section.step, dcval, amplpp, phase,\n slope), c='#8888FF')\n # and the new onset location\n plt.plot(2*[new_from_mid],[-12,dcval],c=\"blue\")\n plt.title(\"Night \"+str(ionset+1)+\" Onset \"+this_id)\n plt.xlabel(\"Time from mid-fit (minutes)\")\n plt.savefig(this_id+\"_onset-\"+str(ionset+1)+\".png\")\n plt.show()\n return predict_onsets, fit_feats","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:54.565710Z","iopub.execute_input":"2023-10-30T19:44:54.566064Z","iopub.status.idle":"2023-10-30T19:44:54.603026Z","shell.execute_reply.started":"2023-10-30T19:44:54.566034Z","shell.execute_reply":"2023-10-30T19:44:54.600944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Go through the onsets and fit the onset transition. Uses mean quantile of LIDSprod.\npredict_onsets, onset_feats = refine_onset_values(predict_onsets, \n known_onsets=actual_onsets)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:44:54.604750Z","iopub.execute_input":"2023-10-30T19:44:54.605144Z","iopub.status.idle":"2023-10-30T19:45:08.900081Z","shell.execute_reply.started":"2023-10-30T19:44:54.605111Z","shell.execute_reply":"2023-10-30T19:45:08.898875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def refine_wakeup_values(predict_wakeups, known_wakeups=[]):\n '''\n Go through the wakeup regions and refine the wakeup locations by fitting.\n Also returns fit features.\n '''\n nfits = 2\n # Set up any features to save - dictionary of arrays\n default_nans = len(predict_wakeups) * [np.NaN]\n fit_feats = {\n 'amplpp': np.array(default_nans),\n 'dcval': np.array(default_nans),\n 'rampsteps': np.array(default_nans),\n 'ratm90': np.array(default_nans),\n 'ratm60': np.array(default_nans),\n 'ratm25': np.array(default_nans),\n 'ratm14': np.array(default_nans),\n 'ratp14': np.array(default_nans),\n 'ratp25': np.array(default_nans),\n 'ratp60': np.array(default_nans),\n 'ratp90': np.array(default_nans)\n }\n for iwake in range(len(predict_wakeups)):\n # continue if not NaN\n this_wakeup = predict_wakeups[iwake]\n if this_wakeup == this_wakeup:\n for ifit in range(nfits):\n # Set the parameters for this fit\n fit_factor = [1.0, 0.7, 0.7][ifit]\n reg_left = 6 * fit_factor\n reg_right = 6 * fit_factor\n fit_left = 5 * fit_factor\n fit_right = 5 * fit_factor\n slope_low = [0.05*(8/(15*12)), 0.7*(8/(15*12)), 0.7*(8/(15*12))][ifit]\n slope_high = (8/(15*12))\n quant_width = [30, 12, 12][ifit] # even value, 1 will be added\n mean_width = [120, 4, 0][ifit] # even value, 1 will be added\n\n # Select a range around the current location\n select_series = ((pruned_series.step > \n this_wakeup - reg_left*60*12) &\n (pruned_series.step <\n this_wakeup + reg_right*60*12))\n # extract that section\n this_section = pruned_series[select_series].copy()\n # Create the feature to fit\n this_section[\"fit_feature\"] = this_section['LIDSprod'].rolling(\n quant_width+1, center=True, min_periods=1).quantile(0.80)\n this_section[\"fit_feature\"] = this_section[\"fit_feature\"].rolling(\n mean_width+1, center=True, min_periods=1).mean()\n # Fit the Wakeup amplpp lower is 4.0\n lower_bounds = [-7.0, 4.0, this_wakeup - fit_left*60*12, slope_low]\n upper_bounds = [ 0.0, 12.0, this_wakeup + fit_right*60*12, slope_high]\n try:\n dcval, amplpp, phase, slope = curve_fit( fit_func_wakeup, \n this_section.step, this_section.fit_feature,\n bounds=(lower_bounds, upper_bounds) )[0]\n except:\n # signal a problem fitting with the dcval and amplpp\n dcval=0.0; amplpp=2.0; phase = this_wakeup, slope=8/(15*12)\n break # stop any further fitting of this event\n # get the first-fit \"ramp steps\" = amplpp / slope\n if ifit == 0:\n rampsteps = amplpp / slope\n # do a repeat fit centering the transition\n this_wakeup = phase\n # Done fitting this night.\n #\n # Get/make the feature values for this fit...\n ratm90 = get_ave_ratio(this_section, -90, -60, 'fit_feature', dcval,amplpp,phase)\n ratm60 = get_ave_ratio(this_section, -60, -25, 'fit_feature', dcval,amplpp,phase)\n ratm25 = get_ave_ratio(this_section, -25, -14, 'fit_feature', dcval,amplpp,phase)\n ratm14 = get_ave_ratio(this_section, -14, -4, 'fit_feature', dcval,amplpp,phase)\n ratp14 = 1.0 - get_ave_ratio(this_section, 4, 14, 'fit_feature', dcval,amplpp,phase)\n ratp25 = 1.0 - get_ave_ratio(this_section, 14, 25, 'fit_feature', dcval,amplpp,phase)\n ratp60 = 1.0 - get_ave_ratio(this_section, 25, 60, 'fit_feature', dcval,amplpp,phase)\n ratp90 = 1.0 - get_ave_ratio(this_section, 60, 90, 'fit_feature', dcval,amplpp,phase)\n #\n # Update this wakeup value in the array\n # Small, constant adjustment based on previous version(s) results.\n # To get to 0, add on the negative of the previous Overall Median.\n wakeup_correction = -40.8 + 45.6 + 0.5 - 3.0 # in steps\n new_wakeup = round(phase + wakeup_correction, 1)\n # If the amplpp is at lower value, 4.00, then NaN the event\n if amplpp < 4.01:\n predict_wakeups[iwake] = np.NaN\n else:\n predict_wakeups[iwake] = new_wakeup\n # Save the features values for this night's fit\n fit_feats['amplpp'][iwake] = amplpp\n fit_feats['dcval'][iwake] = dcval\n fit_feats['rampsteps'][iwake] = rampsteps\n fit_feats['ratm90'][iwake] = ratm90\n fit_feats['ratm60'][iwake] = ratm60\n fit_feats['ratm25'][iwake] = ratm25\n fit_feats['ratm14'][iwake] = ratm14\n fit_feats['ratp14'][iwake] = ratp14\n fit_feats['ratp25'][iwake] = ratp25\n fit_feats['ratp60'][iwake] = ratp60\n fit_feats['ratp90'][iwake] = ratp90\n #\n # Make a plot with comparison to actual wakeups, if given\n if len(known_wakeups) > 0:\n print(\"\\nsteps error = {:.0f} amplpp = {:.2f}\".\n format(predict_wakeups[iwake] - known_wakeups[iwake], amplpp) +\n \" dcval = {:.2f} rampsteps = {:.1f}\".\n format(dcval, rampsteps))\n print(\"\\n ratios: m90= {:.3f}, m60= {:.3f}, m25= {:.3f}, m14= {:.3f},\".\n format(ratm90, ratm60, ratm25, ratm14),\n \"p14= {:.3f}, p25= {:.3f}, p60= {:.3f}, p90= {:.3f}\\n\".\n format(ratp14, ratp25, ratp60, ratp90))\n # Show time axis as minutes from the mid-transition (phase) value\n # Adjust the event locations too.\n min_from_mid = (this_section.step - phase)/12\n known_from_mid = (known_wakeups[iwake]-phase)/12\n new_from_mid = (new_wakeup - phase)/12\n plt.figure(figsize=(6,3))\n plt.plot(min_from_mid, this_section.fit_feature, '.',\n markersize=2,c='gray')\n # Plot the known event location\n plt.plot(2*[known_from_mid], [dcval,2],c=\"green\")\n # and the fit (using the unadjusted phase)\n plt.plot(min_from_mid,\n fit_func_wakeup(this_section.step, dcval, amplpp, phase,\n slope), c='#8888FF')\n # and the fit location\n plt.plot(2*[new_from_mid],[-12,dcval],c=\"blue\")\n plt.title(\"Night \"+str(iwake+1)+\" Wakeup \"+this_id)\n plt.xlabel(\"Time from mid-fit (minutes)\")\n plt.savefig(this_id+\"_wakeup-\"+str(iwake+1)+\".png\")\n plt.show()\n # end for ifit\n # end if\n # end for iwake\n return predict_wakeups, fit_feats","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:45:08.902020Z","iopub.execute_input":"2023-10-30T19:45:08.902497Z","iopub.status.idle":"2023-10-30T19:45:08.938453Z","shell.execute_reply.started":"2023-10-30T19:45:08.902444Z","shell.execute_reply":"2023-10-30T19:45:08.937135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Go through the wakeup regions and fit the wakeup transition\npredict_wakeups, wakeup_feats = refine_wakeup_values(predict_wakeups, \n known_wakeups=actual_wakeups)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:45:08.939888Z","iopub.execute_input":"2023-10-30T19:45:08.940228Z","iopub.status.idle":"2023-10-30T19:45:23.840725Z","shell.execute_reply.started":"2023-10-30T19:45:08.940200Z","shell.execute_reply":"2023-10-30T19:45:23.839459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The errors\nonset_errs = predict_onsets - actual_onsets\nwakeup_errs = predict_wakeups - actual_wakeups\n\n# Calculate the average error in the onsets and wakeups\nmed_err_onsets = np.nanmedian(onset_errs)\nmed_err_wakeups = np.nanmedian(wakeup_errs)\n \n# Show the Wakeup errors - x axis is night\nplt.figure(figsize=(5,2))\nplt.plot(range(1,len(wakeup_errs)+1), wakeup_errs, \"b.\")\nplt.plot([0,len(wakeup_errs)+2],[360.0,360.0],c='orange')\nplt.plot([0,len(wakeup_errs)+2],[-360.0,-360.0],c='orange')\nplt.xlabel(\"Night number\")\nplt.title(\"Wakeup errors\")\nplt.show()\n\nprint(\" \")\nprint(\"Median error in Wakeup : {:.0f}\\n\\n\".\n format(med_err_wakeups))\n\n# Show the Onset errors - x axis is night\nplt.figure(figsize=(5,2))\nplt.plot(range(1,len(onset_errs)+1), onset_errs,\"b.\")\nplt.plot([0,len(onset_errs)+2],[360.0,360.0],c='orange')\nplt.plot([0,len(onset_errs)+2],[-360.0,-360.0],c='orange')\nplt.xlabel(\"Night number\")\nplt.title(\"Onset errors\")\nplt.show()\n\nprint(\" \")\nprint(\"Median error in Onset : {:.0f}\\n\".\n format(med_err_onsets))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:45:23.842350Z","iopub.execute_input":"2023-10-30T19:45:23.842767Z","iopub.status.idle":"2023-10-30T19:45:24.314508Z","shell.execute_reply.started":"2023-10-30T19:45:23.842735Z","shell.execute_reply":"2023-10-30T19:45:24.313340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate the APs and average them over tolerances and event type\n# This ignores the score values.\n\nprint(\"Series id\", this_id, \"Wakeup:\")\nAP_valsW = calc_AP_values(predict_wakeups, actual_wakeups, OUTPUT=True)\nwakeup_AP = sum(AP_valsW)/len(AP_valsW)\n\nprint(\"Series id\", this_id, \"Onset:\")\nAP_valsO = calc_AP_values(predict_onsets, actual_onsets, OUTPUT=True)\nonset_AP = sum(AP_valsO)/len(AP_valsO)\n\noverall_AP = (wakeup_AP + onset_AP) / 2.0\nprint(\"Overall AP value = {:.6f}\".format(overall_AP))","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:45:24.316084Z","iopub.execute_input":"2023-10-30T19:45:24.316848Z","iopub.status.idle":"2023-10-30T19:45:24.330229Z","shell.execute_reply.started":"2023-10-30T19:45:24.316806Z","shell.execute_reply":"2023-10-30T19:45:24.328733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use pre-trained ML modela to assign scores based on event features.\n# The models were trained in the 'phat ML notebook, using data from 'phat data.\n# Using _v18: \nonset_tree = dill.load(open(\"/kaggle/input/cmi-2023-sleepy-times-phat-data/\"+\n \"onset_tree_v18.dill\", 'rb'))\nwakeup_tree = dill.load(open(\"/kaggle/input/cmi-2023-sleepy-times-phat-data/\"+\n \"wakeup_tree_v18.dill\", 'rb'))\n\ndf_onset_feats = pd.DataFrame(onset_feats)\n# Add combined features: average plus and minus ratio values\ndf_onset_feats[\"avem\"] = (20*df_onset_feats[\"ratm25\"] + 35*df_onset_feats[\"ratm60\"] +\n 30*df_onset_feats[\"ratm90\"]) / 85.0\ndf_onset_feats[\"avep\"] = (20*df_onset_feats[\"ratp25\"] + 35*df_onset_feats[\"ratp60\"] +\n 30*df_onset_feats[\"ratp90\"]) / 85.0\n# scores via ML\nonset_scores = 1.0 - onset_tree.predict(df_onset_feats.fillna(0.0))\n##print(onset_scores)\n\ndf_wakeup_feats = pd.DataFrame(wakeup_feats)\n# Add combined features: average plus and minus ratio values\ndf_wakeup_feats[\"avem\"] = (20*df_wakeup_feats[\"ratm25\"] + 35*df_wakeup_feats[\"ratm60\"] +\n 30*df_wakeup_feats[\"ratm90\"]) / 85.0\ndf_wakeup_feats[\"avep\"] = (20*df_wakeup_feats[\"ratp25\"] + 35*df_wakeup_feats[\"ratp60\"] +\n 30*df_wakeup_feats[\"ratp90\"]) / 85.0\n# scores via ML\nwakeup_scores = 1.0 - wakeup_tree.predict(df_wakeup_feats.fillna(0.0))\n##print(wakeup_scores)\n\n# Show the scores vs errors\nfig, (axon, axwk) = plt.subplots(1, 2, figsize=(9,3))\naxon.plot(onset_errs/12, onset_scores, '.b')\naxon.set_xlim(-90,90)\naxon.set_ylim(0.0,1.0)\naxon.set_xlabel(\"Error (min)\")\naxon.set_ylabel(\"Assigned Score value\")\naxon.set_title(\"Onsets\")\naxwk.plot(wakeup_errs/12, wakeup_scores, '.b')\naxwk.set_xlim(-90,90)\naxwk.set_ylim(0.0,1.0)\naxwk.set_xlabel(\"Error (min)\")\naxwk.set_ylabel(\"Assigned Score value\")\naxwk.set_title(\"Wakeups\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-30T19:45:24.332021Z","iopub.execute_input":"2023-10-30T19:45:24.332605Z","iopub.status.idle":"2023-10-30T19:45:24.780151Z","shell.execute_reply.started":"2023-10-30T19:45:24.332558Z","shell.execute_reply":"2023-10-30T19:45:24.778891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ \n\n## Make Test Predictions","metadata":{"papermill":{"duration":0.020431,"end_time":"2023-09-07T16:55:38.614072","exception":false,"start_time":"2023-09-07T16:55:38.593641","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Test has these already assigned\n# parquet_test = pq.ParquetFile(...)\n# pads_test = pads.dataset(...)\n\ntest_rows = parquet_test.metadata.num_rows\nprint(\"Test has\", test_rows, \"rows.\\n . . . working . . .\")\n\n# Get all the test series ids\ntest_ids = get_ids_parquet(parquet_test)\nprint(\"Number of series ids:\",len(test_ids))","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:46:17.688009Z","iopub.execute_input":"2023-10-29T21:46:17.688479Z","iopub.status.idle":"2023-10-29T21:46:43.249991Z","shell.execute_reply.started":"2023-10-29T21:46:17.688426Z","shell.execute_reply":"2023-10-29T21:46:43.248055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Big loop over the test series-es,\n# Do all the analysis steps above without plots, etc.\n# and if REAL_TEST is true then there is no Events for comparison.\n\nif REAL_TEST:\n max_id_to_do = len(test_ids) \nelse:\n max_id_to_do = len(test_ids) # Do all test=train series\n \n# Keep track of the event results for submission.\n# Can also include other per-event information and, for not-real, the ground truth.\nsubmission_frame = []\n# Keep track of summary information by Series - median errors, AP values, etc.\nall_onset_meds = []\nall_wakeup_meds = []\nall_onset_APs = []\nall_wakeup_APs = []\nall_rated_ratios = []\nfor series_number, this_id in enumerate(test_ids[0:max_id_to_do]):\n #\n # Read in the data for this series\n this_series = read_series_pads(this_id,pads_test)\n # Length in days:\n this_duration = (this_series.step.max() - this_series.step.min())/(24*60*12)\n #\n # Add other features to the series and create a down-sampled version\n wind_width = 12 # ***PARAMETER***\n pruned_series = add_features(this_series,wind_width)\n #\n # Find the location of the first/average middle-night in steps.\n phase = fit_first_night_phase(pruned_series, SHOW_PLOT=False)\n #\n # Create initial event values lists for this series:\n predict_onsets, predict_wakeups = setup_prelim_events(\n phase, this_duration, OUTPUT=False)\n #\n # Assign values to the valid nights based on average time asleep and phase.\n predict_onsets, predict_wakeups, this_nights_asleep = \\\n assign_ave_events(phase, predict_onsets, predict_wakeups, OUTPUT=False)\n #\n # Get the known events for this series if Test is Train.\n # and do padding to same lengths before fitting.\n if not(REAL_TEST):\n (this_events, actual_onsets,\n actual_wakeups) = get_series_events(this_id, train_events)\n # number of listed and rated nights:\n this_num_nights = this_events[\"night\"].max()\n this_nights_rated = int(sum(this_events.step > 0.0)/2)\n #\n # Pad as needed to have predicted and actual the same lengths\n predict_onsets, actual_onsets = pad_same_length(\n predict_onsets, actual_onsets, OUTPUT=False)\n predict_wakeups, actual_wakeups = pad_same_length(\n predict_wakeups, actual_wakeups, OUTPUT=False)\n #\n # Go through the onset regions and fit the onset transition.\n # Also returns fit features dict.\n predict_onsets, onset_feats = refine_onset_values(predict_onsets)\n #\n # Go through the wakeup regions and fit the wakeup transition.\n # Also returns fit features dict.\n predict_wakeups, wakeup_feats = refine_wakeup_values(predict_wakeups)\n #\n # Features are onset_feats and wakeup_feats and also:\n predict_durs = predict_wakeups - predict_onsets\n # \n # Assign score values based on the fit features (tree models were read-in above)\n df_onset_feats = pd.DataFrame(onset_feats)\n # Add combined features: average plus and minus ratio values\n df_onset_feats[\"avem\"] = (20*df_onset_feats[\"ratm25\"] + 35*df_onset_feats[\"ratm60\"] +\n 30*df_onset_feats[\"ratm90\"]) / 85.0\n df_onset_feats[\"avep\"] = (20*df_onset_feats[\"ratp25\"] + 35*df_onset_feats[\"ratp60\"] +\n 30*df_onset_feats[\"ratp90\"]) / 85.0\n onset_scores = 1.0 - onset_tree.predict(df_onset_feats.fillna(0.0))\n df_wakeup_feats = pd.DataFrame(wakeup_feats)\n # Add combined features: average plus and minus ratio values\n df_wakeup_feats[\"avem\"] = (20*df_wakeup_feats[\"ratm25\"] + 35*df_wakeup_feats[\"ratm60\"] +\n 30*df_wakeup_feats[\"ratm90\"]) / 85.0\n df_wakeup_feats[\"avep\"] = (20*df_wakeup_feats[\"ratp25\"] + 35*df_wakeup_feats[\"ratp60\"] +\n 30*df_wakeup_feats[\"ratp90\"]) / 85.0\n wakeup_scores = 1.0 - wakeup_tree.predict(df_wakeup_feats.fillna(0.0))\n #\n # Final event values and features have been determined.\n # \n # Monitoring summary info for each Series if Test is Train\n if not(REAL_TEST):\n # Show a one-line summary for each series (if not real test) \n # Calculate the AP values (simple routine ignores \"score\")\n ave_prec_onset = sum(calc_AP_values(\n predict_onsets, actual_onsets))/10.0\n ave_prec_wakeup = sum(calc_AP_values(\n predict_wakeups, actual_wakeups))/10.0\n ave_prec = (ave_prec_onset + ave_prec_wakeup)/2.0\n # The onset/wakeup errors\n onset_errs = predict_onsets - actual_onsets\n wakeup_errs = predict_wakeups - actual_wakeups\n #\n # Calculate the MEDIAN error in the series' onsets and wakeups\n # -- ignore outliers: ones with abs(err) exceeding some value.\n outlier_value = 200 # steps\n use_these = (np.abs(onset_errs) < outlier_value)\n if sum(use_these) > 0:\n med_err_onsets = np.nanmedian(onset_errs[use_these])\n else:\n med_err_onsets = np.nan\n use_these = (np.abs(wakeup_errs) < outlier_value)\n if sum(use_these) > 0:\n med_err_wakeups = np.nanmedian(wakeup_errs[use_these])\n else:\n med_err_wakeups = np.nan\n # \n # Summarize this series results in one line\n print(\"[{}] {} : Nights = {:.1f}, {}, {}, {}\".\n format(series_number, this_id, this_duration,\n this_num_nights, this_nights_rated, this_nights_asleep),\n \" medErrs= {:.0f}, {:.0f}, APs= {:.3f}, {:.3f}, {:.3f}\".\n format(med_err_onsets, med_err_wakeups,\n ave_prec_onset, ave_prec_wakeup, ave_prec))\n # Add to the list of all series info\n all_onset_meds.append(med_err_onsets)\n all_wakeup_meds.append(med_err_wakeups)\n all_onset_APs.append(ave_prec_onset)\n all_wakeup_APs.append(ave_prec_wakeup)\n if this_num_nights > 0:\n all_rated_ratios.append(this_nights_rated / this_num_nights)\n else:\n all_rated_ratios.append(0.0)\n #\n # Create/assemble a dataframe of these results in submission format\n # Required for submission:\n row_id=[]; series_id=[]; step=[]; event=[]; score=[]\n # Additional values - \"features\", mostly based on the event fitting\n duration=[] # Predicted sleep duration of the night = wakeup - onset (steps)\n dcval = []; amplpp = []; rampsteps = []\n ratm90 = []; ratm60 = []; ratm25 = []; ratm14 = []\n ratp14 = []; ratp25 = []; ratp60 = []; ratp90 = []\n # Additional value - The Target step value, only non-NaN for test is train\n truth_step=[]\n phat_cols = ['duration',\n 'dcval', 'amplpp', 'rampsteps',\n 'ratm90', 'ratm60', 'ratm25', 'ratm14',\n 'ratp14', 'ratp25', 'ratp60', 'ratp90',\n 'truth_step']\n #\n # Alternate the onset and wakeup events\n for inight in range(len(predict_onsets)):\n # Onset entry in frame\n row_id.append(0) # dummy value\n series_id.append(this_id)\n step.append(predict_onsets[inight])\n event.append('onset')\n score.append(onset_scores[inight])\n # Features\n duration.append(predict_durs[inight])\n amplpp.append(onset_feats['amplpp'][inight])\n dcval.append(onset_feats['dcval'][inight])\n rampsteps.append(onset_feats['rampsteps'][inight])\n ratm90.append(onset_feats['ratm90'][inight])\n ratm60.append(onset_feats['ratm60'][inight])\n ratm25.append(onset_feats['ratm25'][inight])\n ratm14.append(onset_feats['ratm14'][inight])\n ratp14.append(onset_feats['ratp14'][inight])\n ratp25.append(onset_feats['ratp25'][inight])\n ratp60.append(onset_feats['ratp60'][inight])\n ratp90.append(onset_feats['ratp90'][inight])\n # Target, if known\n if not(REAL_TEST):\n # include the known onset value too\n truth_step.append(actual_onsets[inight])\n else:\n truth_step.append(np.NaN)\n # Wakeup entry in frame\n row_id.append(0) # dummy value\n series_id.append(this_id)\n step.append(predict_wakeups[inight])\n event.append('wakeup')\n score.append(wakeup_scores[inight])\n # Features\n duration.append(predict_durs[inight])\n amplpp.append(wakeup_feats['amplpp'][inight])\n dcval.append(wakeup_feats['dcval'][inight])\n rampsteps.append(wakeup_feats['rampsteps'][inight])\n ratm90.append(wakeup_feats['ratm90'][inight])\n ratm60.append(wakeup_feats['ratm60'][inight])\n ratm25.append(wakeup_feats['ratm25'][inight])\n ratm14.append(wakeup_feats['ratm14'][inight])\n ratp14.append(wakeup_feats['ratp14'][inight])\n ratp25.append(wakeup_feats['ratp25'][inight])\n ratp60.append(wakeup_feats['ratp60'][inight])\n ratp90.append(wakeup_feats['ratp90'][inight])\n # Target, if known\n if not(REAL_TEST):\n # include the known wakeup value too\n truth_step.append(actual_wakeups[inight])\n else:\n truth_step.append(np.NaN)\n # Put these into a data frame\n this_submit = pd.DataFrame(data={'row_id':row_id, 'series_id':series_id, \n 'step':step, 'event':event, 'score':score,\n 'duration':duration,\n 'amplpp':amplpp,\n 'dcval':dcval,\n 'rampsteps':rampsteps,\n 'ratm90':ratm90,\n 'ratm60':ratm60,\n 'ratm25':ratm25,\n 'ratm14':ratm14,\n 'ratp14':ratp14,\n 'ratp25':ratp25,\n 'ratp60':ratp60,\n 'ratp90':ratp90,\n 'truth_step':truth_step}) \n # If this is the first series make it the submission frame\n if len(submission_frame) == 0:\n submission_frame = this_submit.copy()\n else:\n submission_frame = pd.concat([submission_frame, this_submit])\n #\n # Go to next series and repeat.","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:48:00.156280Z","iopub.execute_input":"2023-10-29T21:48:00.157237Z","iopub.status.idle":"2023-10-29T21:48:57.894608Z","shell.execute_reply.started":"2023-10-29T21:48:00.157192Z","shell.execute_reply":"2023-10-29T21:48:57.893032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not(REAL_TEST):\n \n # Save the current \"phat\" submission frame and output it as a csv file.\n # Reset the index and fill in the row_id\n phat_frame = submission_frame.copy()\n phat_frame = phat_frame.reset_index(drop=True)\n phat_frame.row_id = phat_frame.index\n \n phat_frame.to_csv(\"phat_frame.csv\", header=True, \n index=False, na_rep='',\n float_format='%.3f')\n \n print(phat_frame[0:30])","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:03.810669Z","iopub.execute_input":"2023-10-29T21:49:03.811061Z","iopub.status.idle":"2023-10-29T21:49:03.844120Z","shell.execute_reply.started":"2023-10-29T21:49:03.811028Z","shell.execute_reply":"2023-10-29T21:49:03.843011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### If not REAL_TEST, Look at some of the per-series values that were saved.\n\nThe full \"phat_frame\" per-event values are explored in a separate notebook, CMI 2023 - Sleepy Times - phat ML , using a phat_frame.csv uploaded to the dataset CMI 2023 - Sleepy Times - phat data .","metadata":{}},{"cell_type":"code","source":"if not(REAL_TEST):\n # Approx average AP\n ave_overall_AP = ( (np.nansum(all_onset_APs) + np.nansum(all_wakeup_APs))/\n (len(all_onset_APs) + len(all_wakeup_APs)) )\n print(\"\\n\\tAverage AP over all series = {} \\n\".format(ave_overall_AP))\n \n series_APs = (np.array(all_onset_APs)+np.array(all_wakeup_APs))/2.0\n \n plt.figure(figsize=(8,2.5))\n plt.plot(series_APs,'.')\n plt.title(\"AP vs Series\")\n plt.xlabel(\"Series number\")\n plt.ylim(-0.02,0.8); \n plt.ylabel(\"AP\")\n plt.savefig(\"AP_vs_series.png\")\n plt.show()\n \n plt.figure(figsize=(8,2.5))\n plt.plot(np.array(all_rated_ratios),'.')\n plt.title(\"Ratio of Valid Nights / Labeled Nights\")\n plt.xlabel(\"Series number\")\n plt.ylim(-0.02,1.02); \n plt.ylabel(\"Ratio (Valid/Labeled)\")\n plt.savefig(\"ratio_valid_labeled.png\")\n plt.show()\n \n plt.figure(figsize=(5,5))\n plt.plot(all_rated_ratios, series_APs,'.')\n plt.title(\"AP vs Ratio Valid/Labeled\")\n plt.ylabel(\"Series AP\")\n plt.ylim(-0.02,0.8); \n plt.xlabel(\"Ratio (Valid/Labeled)\")\n plt.xlim(-0.02,1.02);\n plt.savefig(\"AP_vs_ratio.png\")\n plt.show()\n ","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:29.578620Z","iopub.execute_input":"2023-10-29T21:49:29.579024Z","iopub.status.idle":"2023-10-29T21:49:30.554744Z","shell.execute_reply.started":"2023-10-29T21:49:29.578993Z","shell.execute_reply":"2023-10-29T21:49:30.553657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not(REAL_TEST):\n \n plt.figure(figsize=(8,2.5))\n plt.plot(np.array(all_onset_APs),'.')\n plt.title(\"Onset Average Precisions\")\n plt.xlabel(\"Series number\")\n plt.ylim(-0.02,0.8); plt.ylabel(\"Onset AP\")\n plt.savefig(\"onset_ave_APs.png\")\n plt.show()\n\n plt.figure(figsize=(6,3))\n plt.hist(np.array(all_onset_APs),bins=20)\n plt.xlim(0.0,1.0); plt.xlabel(\"Onset AP\")\n plt.ylabel(\"Frequency\")\n plt.title(\"Histogram: Onset Average Precisions\")\n plt.savefig(\"onset_ave_histogram.png\")\n plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:31.650845Z","iopub.execute_input":"2023-10-29T21:49:31.651908Z","iopub.status.idle":"2023-10-29T21:49:32.201738Z","shell.execute_reply.started":"2023-10-29T21:49:31.651868Z","shell.execute_reply":"2023-10-29T21:49:32.200587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not(REAL_TEST):\n \n plt.figure(figsize=(8,2.5))\n plt.plot(np.array(all_wakeup_APs),'.')\n plt.title(\"Wakeup Average Precisions\")\n plt.xlabel(\"Series number\")\n plt.ylim(-0.02,0.8); plt.ylabel(\"Wakeup AP\")\n plt.savefig(\"wakeup_ave_APs.png\")\n plt.show()\n\n plt.figure(figsize=(6,3))\n plt.hist(np.array(all_wakeup_APs),bins=20)\n plt.xlim(0.0,1.0); plt.xlabel(\"Wakeup AP\")\n plt.ylabel(\"Frequency\")\n plt.title(\"Histogram: Wakeup Average Precisions\")\n plt.savefig(\"wakeup_ave_histogram.png\")\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:33.249797Z","iopub.execute_input":"2023-10-29T21:49:33.250206Z","iopub.status.idle":"2023-10-29T21:49:33.823105Z","shell.execute_reply.started":"2023-10-29T21:49:33.250173Z","shell.execute_reply":"2023-10-29T21:49:33.821943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not(REAL_TEST):\n \n plt.figure(figsize=(5,5))\n # y = x line and the average AP line\n plt.plot([0.0,2.0*ave_overall_AP],[2.0*ave_overall_AP,0.0],c='lightgreen')\n plt.plot([0,0.8],[0.0,0.8],c=\"orange\")\n plt.plot(np.array(all_onset_APs), np.array(all_wakeup_APs),'.',c='blue')\n plt.title(\"Wakeup APs vs Onset APs\\nAve combined AP = {:.3f} (light green)\".\n format(ave_overall_AP))\n plt.xlim(0.0,0.8); plt.xlabel(\"Onset AP\")\n plt.ylim(0.0,0.8); plt.ylabel(\"Wakeup AP\")\n plt.savefig(\"wakeup_vs_onset_aveAP.png\")\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:37.603517Z","iopub.execute_input":"2023-10-29T21:49:37.603908Z","iopub.status.idle":"2023-10-29T21:49:37.954244Z","shell.execute_reply.started":"2023-10-29T21:49:37.603878Z","shell.execute_reply":"2023-10-29T21:49:37.953219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not(REAL_TEST):\n \n plt.figure(figsize=(5,5))\n plt.plot([-180.0,180.8],[0,0],c=\"orange\")\n plt.plot([0,0],[-180.0,180.8],c=\"orange\")\n plt.plot(np.array(all_onset_meds), np.array(all_wakeup_meds),'.',c='blue')\n plt.title(\"Wakeup Median Error vs Onset Median Error\")\n plt.xlim(-180.0,180.8); plt.xlabel(\"Median Onset Error (step)\")\n plt.ylim(-180.0,180.8); plt.ylabel(\"Median Wakeup Error (step)\")\n plt.savefig(\"wakeup_vs_onset_medErr.png\")\n plt.show()\n\n print(\"\\nOverall Median Errors: Onset = {:.1f}, Wakeup = {:.1f}\".\n format(np.nanmedian(all_onset_meds), np.nanmedian(all_wakeup_meds)))","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:40.821535Z","iopub.execute_input":"2023-10-29T21:49:40.822278Z","iopub.status.idle":"2023-10-29T21:49:41.170028Z","shell.execute_reply.started":"2023-10-29T21:49:40.822236Z","shell.execute_reply":"2023-10-29T21:49:41.168941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Make the submission csv","metadata":{}},{"cell_type":"code","source":"# For submission, remove the \"phat\" columns\nsubmission_frame = submission_frame.drop(phat_cols,axis=1)\n\n# Then remove rows with NA in step - show the numbers:\nlen_all = len(submission_frame)\nsubmission_frame = submission_frame.dropna(axis=0)\nlen_no_NaN = len(submission_frame)\nprint(\"\\n\\nTotal events = {}, Number without NaNs = {} --> {:.1f} % non-NaN.\\n\\n\".\n format(len_all, len_no_NaN, 100.0*len_no_NaN/len_all))\n# For comparison:\n# Train Events has 14508/2 = 7254 pairs of onset/wakeup.\n# 4791 of onsets and 4794 of wakeups are real values (rest are NaN.)\n# --> 66.1 % of events have non-NaN values.\n\n# Reset the index and fill in the row_id\nsubmission_frame = submission_frame.reset_index(drop=True)\nsubmission_frame.row_id = submission_frame.index\n\nsubmission_frame","metadata":{"execution":{"iopub.status.busy":"2023-10-29T21:49:45.762905Z","iopub.execute_input":"2023-10-29T21:49:45.763793Z","iopub.status.idle":"2023-10-29T21:49:45.786357Z","shell.execute_reply.started":"2023-10-29T21:49:45.763753Z","shell.execute_reply":"2023-10-29T21:49:45.785342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Output in submission file format:\n# row_id,series_id,step,event,score\nsubmission_frame.to_csv(\"submission.csv\", header=True, \n index=False, na_rep='',\n float_format='%.3f')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"##!head submission.csv\n# Look at some random place too\n##!head -50 submission.csv | tail\n##!tail submission.csv","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# --- The End ---","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}