{"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 - phat ML\n\nThe [CMI 2023 - Sleepy Times notebook](https://www.kaggle.com/code/dan3dewey/cmi-2023-sleepy-times) outputs a \"phat_frame\" of predicted Train events which also includes additional columns: the actual event steps, predicted night duration, and various features from the fitting of the event location. The phat_frame.csv that was created is uploaded to the dataset, CMI 2023 - Sleepy Times - phat data . This notebook looks at that phat data to better see what kind of errors are produced.\n\nUsing the phat data, an ML model is trained to predict the level of error for an onset/wakeup event from its fit features. The dill'ed models are output and then saved to the dataset so they can be used by the main notebook to assign scores to the predicted events for submission.\n","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# v1   Initial version: doesn't look like errors depend on duration.\n# v2   Added \"above_rat\" feature and eye-balled score heirarchy for onsets and wakeups.\n# v3   Created the scaled error, -1 to 1, from error in steps.\n# v4   Using v63 of phat frame with more features :)\n# v5   Added side-by-side onset/wakeup plots of sclerr vs features; added scatter matrix plots.\n# v6   Fit abs(sclerr) w/min leaf= 30 and then assign score as 1 - yh\n# v7   Fit sclerr w/min leaf= 20 and then assign score as 1 - abs(yh)\n# v8   Better to fit abs(sclerr) -- how bad is the error --\n#      rather than sclerr -- how bad is the error and in which direction.\n#      Added histograms of the errors.\n# v9   Generate v9 of the dill'ed models using all 6 ratio features from phat v67.\n# v10  Read in phat_frame_v68 which has scores set by the v9 models.\n#      Add scatter plots of score vs error from the phat file.\n# v12  Add two additional features, avem and avep, for improved modeling.\n# v13  Read in the v69 phat output with newest score values.\n# v14  Read in v70 phat frame with p/m 14 ratios added. Create model.\n# v15  Arg: I was writing out onset_tree to the wakeup dill  :(  Try again.\n# v16  Correctly wrote out the models. Use them in making phat frame v72.\n# v17  Read in phat frame v72 to check that its scores agree with expectations - OK.\n# v18  Use Random Forest w/ 10 x 30 min leaf. (LB 503)\n# v19  \"  w/10 x 60 min leaf. (LB 499)\n# v20  \"  w/10 x 15 min leaf.  (LB 495  )  Map out LB vs min leaf a little...\n# v21  Use phat_frame_v76, made with v18 model, and re-run to be up to date.\n# v22  Looked at how score is used to calculate the precision-recall AP.\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-11-01T01:00:35.096745Z","iopub.execute_input":"2023-11-01T01:00:35.097101Z","iopub.status.idle":"2023-11-01T01:00:35.106088Z","shell.execute_reply.started":"2023-11-01T01:00:35.097073Z","shell.execute_reply":"2023-11-01T01:00:35.104502Z"},"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\n##from scipy.optimize import leastsq, curve_fit\n# ML models\nfrom sklearn.tree import DecisionTreeRegressor\nfrom sklearn.ensemble import RandomForestRegressor\n\n# To save/restore trained ML models\nimport dill\n","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-11-01T01:00:38.288983Z","iopub.execute_input":"2023-11-01T01:00:38.289371Z","iopub.status.idle":"2023-11-01T01:00:39.142916Z","shell.execute_reply.started":"2023-11-01T01:00:38.289340Z","shell.execute_reply":"2023-11-01T01:00:39.141425Z"},"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}\".format(sum(AP_values)/len(AP_values)),\n             \"     (score not used)\\n\")\n    # Return the array of values for each tolerance level\n    return AP_values","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:00:40.495007Z","iopub.execute_input":"2023-11-01T01:00:40.495572Z","iopub.status.idle":"2023-11-01T01:00:40.507658Z","shell.execute_reply.started":"2023-11-01T01:00:40.495539Z","shell.execute_reply":"2023-11-01T01:00:40.506090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ \n\n### Look at the Phat Frame\n\nA previously output phat_frame_vNN.csv was uploded to the kaggle dataset `cmi-2023-sleepy-times-phat-data`, attached to this notebook.","metadata":{}},{"cell_type":"code","source":"# Look at previously saved phat_frame\n# _v70: has amplpp, dcval, rampsteps, duration, rat[m/p][14,25,60,90]\n# _v71: updated scores from v70.\n# _v72: _v71 with the updated scores.\n# _v76: updated with the v18 model scores.\nphat_frame = pd.read_csv(\"../input/cmi-2023-sleepy-times-phat-data/\" +\n                        \"phat_frame_v76.csv\")\n\n# Set the score values to NaN when the step values are NaN\nphat_frame.score = phat_frame.score + 0*phat_frame.step\n\n# Add the event errors\nphat_frame[\"error\"] = phat_frame.step - phat_frame.truth_step\n# Add scaled errors going from -1 to 1\nerr_scale = 120.0  # steps\nphat_frame[\"sclerr\"] = 2.0/(1.0+np.exp(-1.0*phat_frame.error/err_scale)) - 1.0\n\nphat_len = len(phat_frame)\nphat_frame","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:00:44.295555Z","iopub.execute_input":"2023-11-01T01:00:44.295973Z","iopub.status.idle":"2023-11-01T01:00:44.361297Z","shell.execute_reply.started":"2023-11-01T01:00:44.295942Z","shell.execute_reply":"2023-11-01T01:00:44.359555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Add combined features: average plus and minus ratio values\nphat_frame[\"avem\"] = (20*phat_frame[\"ratm25\"] + 35*phat_frame[\"ratm60\"] +\n                      30*phat_frame[\"ratm90\"]) / 85.0\nphat_frame[\"avep\"] = (20*phat_frame[\"ratp25\"] + 35*phat_frame[\"ratp60\"] +\n                      30*phat_frame[\"ratp90\"]) / 85.0","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:00:49.686673Z","iopub.execute_input":"2023-11-01T01:00:49.687196Z","iopub.status.idle":"2023-11-01T01:00:49.697069Z","shell.execute_reply.started":"2023-11-01T01:00:49.687159Z","shell.execute_reply":"2023-11-01T01:00:49.695955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get sub frames with just onsets and just wakeups\nphat_onsets = phat_frame[phat_frame.event == \"onset\"].copy()\nonset_len = len(phat_onsets)\nphat_wakeups = phat_frame[phat_frame.event == \"wakeup\"].copy()\nwakeup_len = len(phat_wakeups)\n\n# Calculate the onset and wakeup APs (ignores score)\nprint(\"Onsets:\")\ntol_APs = calc_AP_values(phat_onsets.step,\n                          phat_onsets.truth_step, OUTPUT=True)\n\nprint(\"Wakeups:\")\ntol_APs = calc_AP_values(phat_wakeups.step,\n                          phat_wakeups.truth_step, OUTPUT=True)","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:00:55.766331Z","iopub.execute_input":"2023-11-01T01:00:55.766680Z","iopub.status.idle":"2023-11-01T01:00:55.809456Z","shell.execute_reply.started":"2023-11-01T01:00:55.766654Z","shell.execute_reply":"2023-11-01T01:00:55.808149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show the relationship between the score and error in the phat file:\nif True:\n    # Show the assigned scores vs errors in the phat file\n    max_range = 90  # minutes\n    fig, (axon, axwk) = plt.subplots(1, 2, figsize=(8,4))\n    \n    axon.plot([0.,0.],[0.,1.],c=\"orange\")\n    axon.plot(phat_onsets.error/12, phat_onsets.score, '.b',\n            markersize=5, alpha=0.10)\n    axon.set_xlim(-1.0*max_range,max_range)\n    axon.set_ylim(0.0,1.05)\n    axon.set_xlabel(\"Error (min)\")\n    axon.set_ylabel(\"Score value\")\n    axon.set_title(\"Onsets\")\n    #\n    axwk.plot([0.,0.],[0.,1.],c=\"orange\")\n    axwk.plot(phat_wakeups.error/12, phat_wakeups.score, '.b',\n            markersize=5, alpha=0.10)\n    axwk.set_xlim(-1.0*max_range,max_range)\n    axwk.set_ylim(0.0,1.05)\n    axwk.set_xlabel(\"Error (min)\")\n    axwk.set_ylabel(\"Score value\")\n    axwk.set_title(\"Wakeups\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:00:57.269338Z","iopub.execute_input":"2023-11-01T01:00:57.269957Z","iopub.status.idle":"2023-11-01T01:00:57.624147Z","shell.execute_reply.started":"2023-11-01T01:00:57.269926Z","shell.execute_reply":"2023-11-01T01:00:57.622623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Histograms of the errors\nfig, (axon, axwk) = plt.subplots(1, 2, figsize=(10,4))\n    \nmax_range = 90  # minutes\nbin_width =  2 # minutes\naxon.hist((phat_onsets.error)/12, bins=int(2.0*max_range/bin_width),\n              range=(-1*max_range,max_range), log=True)\naxon.set_title(\"Onset Errors\")\naxon.set_xlabel(\"Error (minutes)\")\n\naxwk.hist((phat_wakeups.error)/12, bins=int(2.0*max_range/bin_width),\n              range=(-1*max_range,max_range), log=True)\naxwk.set_title(\"Wakeup Errors\")\naxwk.set_xlabel(\"Error (minutes)\")\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:01:28.233112Z","iopub.execute_input":"2023-11-01T01:01:28.233506Z","iopub.status.idle":"2023-11-01T01:01:29.277665Z","shell.execute_reply.started":"2023-11-01T01:01:28.233477Z","shell.execute_reply":"2023-11-01T01:01:29.276104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Histograms of the scores\nfig, (axon, axwk) = plt.subplots(1, 2, figsize=(10,4))\n    \naxon.hist(phat_onsets.score, bins=50,\n              range=(0.0,1.0), log=True)\naxon.set_title(\"Onset Scores\")\naxon.set_xlabel(\"Score values\")\n\naxwk.hist(phat_wakeups.score, bins=50,\n              range=(0.0,1.0), log=True)\naxwk.set_title(\"Wakeup Scores\")\naxwk.set_xlabel(\"Score values\")\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:01:49.076051Z","iopub.execute_input":"2023-11-01T01:01:49.076440Z","iopub.status.idle":"2023-11-01T01:01:50.016922Z","shell.execute_reply.started":"2023-11-01T01:01:49.076409Z","shell.execute_reply":"2023-11-01T01:01:50.014554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Demo: how the precision-recall AP is calculated when score thresholds are used.\n# Do it for a given tolerance level and for either onsets or wakeups:\ndfphat = phat_wakeups\ntol = 12\n\n# The precision-recall ignoring the score:\nnum_GTruth = sum(dfphat.truth_step > 0)\nnum_Predicted = sum(dfphat.step > 0)\nnum_TP = sum(np.abs(dfphat.error) < tol)\nprint(\"Precision-Recall AP ignoring score = {:.4f}\".\n              format((num_TP/num_Predicted)*(num_TP/num_GTruth)))\n\n# Set the score thresholds based on quantiles of the scores\nthreshs = np.nanquantile(dfphat.score, [0.8, 0.65, 0.40, 0.0])\nprint(\"Score thresholds:\",threshs)\n# Keep track\ncumu_recall = 0.0\ncumu_pr = 0.0\nfor thresh in threshs:\n    # Use only predictions with score above the threshold value\n    num_Predicted = sum(dfphat[dfphat.score > thresh].step > 0)\n    num_TP = sum(np.abs(dfphat[dfphat.score > thresh].error) < tol)\n    prec = (num_TP/num_Predicted)\n    d_recall = (num_TP/num_GTruth) - cumu_recall\n    cumu_pr += prec*d_recall\n    cumu_recall += d_recall\n    print(\"{:.3f}:  {}  {}  {}    {:.4f} * {:.4f} = {:.4f}   {:.4f}\".format\n          (thresh, num_GTruth, num_Predicted, num_TP,\n           prec, d_recall, prec*d_recall, cumu_pr))\n","metadata":{"execution":{"iopub.status.busy":"2023-11-01T01:32:08.390558Z","iopub.execute_input":"2023-11-01T01:32:08.390988Z","iopub.status.idle":"2023-11-01T01:32:08.416329Z","shell.execute_reply.started":"2023-11-01T01:32:08.390960Z","shell.execute_reply":"2023-11-01T01:32:08.414579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Show Scaled Error vs the Features","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"# Show the tolerence error values and their scaled values\nerr_tols = np.array([12, 36, 60, 90, 120, 150, 180, 240, 300, 360])\ntols_scaled = 2.0/(1.0+np.exp(-1.0*err_tols/err_scale)) - 1.0\nprint(\"Relation between actual error in steps [minutes] and the scaled error\")\nfor itol in range(len(err_tols)):\n    print(\"  {} [{}] steps [min]  -->  {:.2f} scaled error\".\n          format(err_tols[itol], err_tols[itol]/12, tols_scaled[itol]))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show the Wakeup errors, duration, and Onset errors for all nights\n\n# Show the Wakeup errors\nplt.figure(figsize=(9,3))\nplt.plot(range(1,wakeup_len+1), phat_wakeups.sclerr,\n         \"b.\", markersize=2)\nplt.plot([0,wakeup_len+2],2*[0.0],c='orange')\n# 120 and 360 levels\nplt.plot([0,wakeup_len+2],2*[0.46],c='orange')\nplt.plot([0,wakeup_len+2],2*[-0.46],c='orange')\nplt.plot([0,wakeup_len+2],2*[0.91],c='orange')\nplt.plot([0,wakeup_len+2],2*[-0.91],c='orange')\nplt.ylim(-1.05,1.05)\nplt.xlabel(\"Wakeup number\")\nplt.ylabel(\"Scaled Error (orange: +/- 10, 30 min)\")\nplt.title(\"Wakeup Scaled Errors\")\nplt.show()\n\n# Show the Sleep Duration\none_hr = 60*12  # steps\nplt.figure(figsize=(9,3))\nplt.plot(range(1,wakeup_len+1), phat_wakeups.duration/one_hr,\n         \"b.\", markersize=2)\nfor this_hr in [4,8,12]:\n    plt.plot([0,wakeup_len+2],2*[this_hr],c='lime')\nplt.ylim(0,18)\nplt.xlabel(\"Onset-Wakeup number\")\nplt.title(\"Sleep Duration (hours)\")\nplt.show()\n\n# Show the Onset errors\nplt.figure(figsize=(9,3))\nplt.plot(range(1,onset_len+1), phat_onsets.sclerr,\n         \"b.\", markersize=2)\nplt.plot([0,onset_len+2],2*[0.0],c='orange')\n# 120 and 360 levels\nplt.plot([0,onset_len+2],2*[0.46],c='orange')\nplt.plot([0,onset_len+2],2*[-0.46],c='orange')\nplt.plot([0,onset_len+2],2*[0.91],c='orange')\nplt.plot([0,onset_len+2],2*[-0.91],c='orange')\nplt.ylim(-1.05,1.05)\nplt.xlabel(\"Onset number\")\nplt.ylabel(\"Scaled Error (orange: +/- 10, 30 min)\")\nplt.title(\"Onset Scaled Errors\")\nplt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plots of the errors vs fit features\n# duration is not included for now since it requires onset and wakeup values for a night.\nfit_feats = ['amplpp', 'dcval', 'rampsteps',\n                'ratm90', 'ratm60', 'ratm25', 'ratm14',\n                'ratp14', 'ratp25', 'ratp60', 'ratp90',\n                'avem', 'avep']\n\nfor this_feat in fit_feats:\n    \n    # Onset and Wakeup side-by-side\n    fig, (axon, axwk) = plt.subplots(1, 2, figsize=(7,3))\n    \n    if ('ratm' in this_feat) or ('ratp' in this_feat) or ('ave' in this_feat):\n        # use squareroot scaling\n        axon.plot(np.sqrt(phat_onsets[this_feat]), phat_onsets.sclerr,\n                 \"b.\", markersize=2, alpha=0.2)\n        axon.set_title(\"Onset Errors vs sqrt(\"+this_feat+\")\")\n        feat_minmax = [np.sqrt(phat_onsets[this_feat].quantile(0.005)),\n                   np.sqrt(phat_onsets[this_feat].quantile(0.995))]\n    else:\n        axon.plot(phat_onsets[this_feat], phat_onsets.sclerr,\n                 \"b.\", markersize=2, alpha=0.2)\n        axon.set_title(\"Onset Errors vs \"+this_feat)\n        feat_minmax = [phat_onsets[this_feat].quantile(0.005),\n                   phat_onsets[this_feat].quantile(0.995)]\n    axon.plot(feat_minmax,2*[0.91],c='orange')\n    axon.plot(feat_minmax,2*[0.46],c='orange')\n    axon.plot(feat_minmax,2*[0.0],c='orange')\n    axon.plot(feat_minmax,2*[-0.91],c='orange')\n    axon.plot(feat_minmax,2*[-0.46],c='orange')\n    axon.set_xlim(feat_minmax[0], feat_minmax[1])\n    axon.set_ylim(-1.05,1.05)\n    axon.set_ylabel(\"Scaled Error (lines: +/- 10, 30 min)\")\n    \n\n    if ('ratm' in this_feat) or ('ratp' in this_feat) or ('ave' in this_feat):\n        # use squareroot scaling\n        axwk.plot(np.sqrt(phat_wakeups[this_feat]), phat_wakeups.sclerr,\n                 \"b.\", markersize=2, alpha=0.2)\n        axwk.set_title(\"Wakeup Errors vs sqrt(\"+this_feat+\")\")\n        feat_minmax = [np.sqrt(phat_wakeups[this_feat].quantile(0.005)),\n                   np.sqrt(phat_wakeups[this_feat].quantile(0.995))]\n    else:\n        axwk.plot(phat_wakeups[this_feat], phat_wakeups.sclerr,\n                 \"b.\", markersize=2, alpha=0.2)\n        axwk.set_title(\"Wakeup Errors vs \"+this_feat)\n        feat_minmax = [phat_wakeups[this_feat].quantile(0.005),\n                   phat_wakeups[this_feat].quantile(0.995)]\n    axwk.plot(feat_minmax,2*[0.91],c='orange')\n    axwk.plot(feat_minmax,2*[0.46],c='orange')\n    axwk.plot(feat_minmax,2*[0.0],c='orange')\n    axwk.plot(feat_minmax,2*[-0.91],c='orange')\n    axwk.plot(feat_minmax,2*[-0.46],c='orange')\n    axwk.set_xlim(feat_minmax[0], feat_minmax[1])\n    axwk.set_ylim(-1.05,1.05)\n    ##axwk.set_ylabel(\"Scaled Error (lines: +/- 10, 30 min)\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Onsets Scatter-Matrix Plots - color-coded by scaled error\npd.plotting.scatter_matrix(phat_onsets.dropna(axis=0)[fit_feats],\n                           figsize=(9,9), s=20, alpha=0.3,\n                           cmap='plasma_r', # plasma_r\n                           c=np.clip(np.abs(phat_onsets.dropna(axis=0).sclerr),0.5,1.0))\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Wakeup Scatter-Matrix Plots - color-coded by scaled error\npd.plotting.scatter_matrix(phat_wakeups.dropna(axis=0)[fit_feats],\n                           figsize=(9,9), s=20, alpha=0.3,\n                           cmap='plasma_r', # plasma_r\n                           c=np.clip(np.abs(phat_wakeups.dropna(axis=0).sclerr),0.5,1.0))\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fit the abs(sclerr) using the fit features","metadata":{}},{"cell_type":"code","source":"# Predict yh using y = abs(sclerr) and X = fit features.\n# Can then assign score as 1 - yh .\nX = phat_onsets.dropna(axis=0)[fit_feats].copy()\ny = np.abs(phat_onsets.dropna(axis=0)['sclerr'].copy())\n\n# Simple tree\n##tree_regress = DecisionTreeRegressor(max_depth=None, min_samples_leaf=30)\n\n# Random Forest\n# class sklearn.ensemble.RandomForestRegressor(n_estimators=100, *, \n# criterion='squared_error', max_depth=None, min_samples_split=2, min_samples_leaf=1, \n# min_weight_fraction_leaf=0.0, max_features=1.0, max_leaf_nodes=None,\n# min_impurity_decrease=0.0, bootstrap=True, oob_score=False, n_jobs=None,\n# random_state=None, verbose=0, warm_start=False, ccp_alpha=0.0,\n# max_samples=None)[source]¶\ntree_regress = RandomForestRegressor(min_samples_leaf=30, n_estimators=10)\n\nonset_tree = tree_regress.fit(X,y)\n\nprint(\"X and y have lengths of\",len(y))\nprint(\"Feature importances:\\n\", fit_feats,\"\\n\",\n      onset_tree.feature_importances_)\n\nyh = onset_tree.predict(X)\n\n# Assign score as 1 - yh\nplt.plot(y, 1 - yh,'b.',alpha=0.2)\nplt.plot(2*[0.91], [0,0.93],c='orange')\nplt.plot(2*[0.46], [0,0.93], c='orange')\nplt.xlim(-0.02,1.02)\nplt.ylim(-0.02,1.02)\nplt.text(0.46-0.05,0.97,\"10 min\",c='orange')\nplt.text(0.91-0.05,0.97,\"30 min\",c='orange')\nplt.title(\"Onsets: Assigned Score vs abs(Scaled Error)\")\nplt.xlabel(\"abs( Scaled Error )\")\nplt.ylabel(\"Assigned Score = 1 - y_hat\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Predict yh using y = abs(sclerr) and X = fit features.\n# Can then assign score as 1 - yh .\nX = phat_wakeups.dropna(axis=0)[fit_feats].copy()\ny = np.abs(phat_wakeups.dropna(axis=0)['sclerr'].copy())\n\n# Simple tree\n##tree_regress = DecisionTreeRegressor(max_depth=None, min_samples_leaf=30)\n\n# Random Forest\n# class sklearn.ensemble.RandomForestRegressor(n_estimators=100, *, \n# criterion='squared_error', max_depth=None, min_samples_split=2, min_samples_leaf=1, \n# min_weight_fraction_leaf=0.0, max_features=1.0, max_leaf_nodes=None,\n# min_impurity_decrease=0.0, bootstrap=True, oob_score=False, n_jobs=None,\n# random_state=None, verbose=0, warm_start=False, ccp_alpha=0.0,\n# max_samples=None)[source]¶\ntree_regress = RandomForestRegressor(min_samples_leaf=30, n_estimators=10)\n\nwakeup_tree = tree_regress.fit(X,y)\n\nprint(\"X and y have lengths of\",len(y))\nprint(\"Feature importances:\\n\", fit_feats,\"\\n\",\n      wakeup_tree.feature_importances_)\n\nyh = wakeup_tree.predict(X)\n\n# Assign score as 1 - yh\nplt.plot(y, 1 - yh,'b.',alpha=0.2)\nplt.plot(2*[0.91], [0,0.93],c='orange')\nplt.plot(2*[0.46], [0,0.93], c='orange')\nplt.xlim(-0.02,1.02)\nplt.ylim(-0.02,1.02)\nplt.text(0.46-0.05,0.97,\"10 min\",c='orange')\nplt.text(0.91-0.05,0.97,\"30 min\",c='orange')\nplt.title(\"Wakeups: Assigned Score vs abs(Scaled Error)\")\nplt.xlabel(\"abs( Scaled Error )\")\nplt.ylabel(\"Assigned Score = 1 - y_hat\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Saving and restoring the models to dill files for use in other notebooks","metadata":{}},{"cell_type":"code","source":"# Save the models to byte dill files\ndill.dump(onset_tree, open(\"onset_tree.dill\", 'wb'), recurse=True)\ndill.dump(wakeup_tree, open(\"wakeup_tree.dill\", 'wb'), recurse=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Restoring","metadata":{}},{"cell_type":"code","source":"# These are needed in the notebook using the dill'ed models.\n##import dill\n##import pandas as pd\n\n# Load a model back in:\nsave_tree = dill.load(open(\"/kaggle/working/onset_tree.dill\", 'rb'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use the model to predict and assign a score to the fit\n\nthis_X = pd.DataFrame(data={'amplpp':[5.61],\n                    'dcval':[-6.08],\n                    'rampsteps':[1752.5],\n                    'ratm90':[0.000],\n                    'ratm60':[0.511],\n                    'ratm25':[0.121],\n                    'ratm14':[0.100],\n                    'ratp14':[0.010],\n                    'ratp25':[0.016],\n                    'ratp60':[0.003],\n                    'ratp90':[0.096],\n                    'avem':[0.288],\n                    'avep':[0.041]\n                    }) \nassigned_score = (1 - save_tree.predict(this_X))[0]\nassigned_score","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# --- The End ---","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}