{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"},{"sourceId":7863631,"sourceType":"datasetVersion","datasetId":4528641}],"dockerImageVersionId":30635,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"papermill":{"default_parameters":{},"duration":83.251521,"end_time":"2024-02-01T20:26:06.174423","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-02-01T20:24:42.922902","version":"2.4.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## HMS 2024 - Brain addicted to tinkering${^{\\scriptsize\\dagger}}$","metadata":{"papermill":{"duration":0.011066,"end_time":"2024-02-01T20:24:47.517269","exception":false,"start_time":"2024-02-01T20:24:47.506203","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"Classifying different types of <b>Harmful Brain Activity</b>, <b>HBA</b>, using EEG recordings and spectra.\n\nA previous notebook, [HMS 2024 - Brain Clusters](https://www.kaggle.com/code/dan3dewey/hms-2024-brain-clusters), looked at general properties of the training targets (the HBA votes/probabilities) including the k-means clustering of the 6-dimension probability vectors of the HBA samples. \n\nThis notebook makes nicer looking vote-cluster plots (~ 1/3 of the way down) and continues by looking at the spectrograms: the spectrum at the event and the flare-rate spectra (fft of a spectral frequency's variation within +/1 128 bins from the event.). Simple features from these are used to classify which cluster each sample belongs to.\nThe features and classifier (random forest) are primative and, unsurprisingly, give mediocre results. Still, it was fun tinkering with the code :)\n\n*($\\scriptsize\\dagger$) to repair, adjust, or work with something in an unskilled or experimental manner*","metadata":{"papermill":{"duration":0.010218,"end_time":"2024-02-01T20:24:47.538042","exception":false,"start_time":"2024-02-01T20:24:47.527824","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"#### Key Hyper-Hyper parameters","metadata":{}},{"cell_type":"code","source":"# Target classes creation:  (will supplant the pre-processed value)\nNUM_CLUSTS = 9      # 6 to 10\n# The clusters are generated from the nominal training samples, defined below.\n\n\n# Use raw features previously saved:\nUSE_PREPROC = True   # Read-in saved meta and features frames\n#  --- OR ---\n# Generate raw features in this code, values below affect the features created\n# Training and Validation rows to use: defined in code below.\n# Number of training/validation samples, downselect factor from nominal sets\nTRAIN_DOWNSEL = 1  # large values for code test; set to 1 to output all.\nVALID_DOWNSEL = 1  #  \"\n# Middle 8 second spectra features - set smooth width which is also the downsample factor\nSMOOTH_WIDTH = 5    # Odd>1: 3,5,7,9,...\n# Flare spectra: do FFT of the amplitude vs time for select spectral frequencies\n# Use:  7(2.0Hz), 15(3.5Hz), 23(5.0Hz)m 60(12Hz), averaged +/-1, +/-2\n#    freqbins = [7, 15, 23, 60]\n# From these ffts make features from these bins\n#    fftfeatbins = [4, 9, 16, 25, 36, 49]\n\n# Train using the combined train and validation events (for best submission?)\nTRAIN_ALL = True\n\n# Use added LR proba features?  (added after reading in pre-proc.)\nLR_REGU = \"l1\"\nUSE_LR1 = True\nLR1_C = 0.05     # smaller --> fewer non-zero coeff.s\nUSE_LR2 = True\nLR2_C = 0.03\n#\nLR_BLUR = 0.10\n\n# Random Forest model parameters set below, e.g.:\n#   n_estimators=200, min_samples_leaf=3, max_features=0.20, max_samples=0.85\n\n# Used the \"tamed cluster\" scheme based on validation best_centers?\n# Otherwise, use the model proba values to weight the cluster centers.\nUSE_TAMED = False\n\n# - - -\n# Directory prefix for the data\nabove_dir = \"../input/hms-harmful-brain-activity-classification/\"\n# and for pre-processed inputs\nabove_dir_preproc = \"../input/hms-2024-brain-data/\"\n# or, offline, use my local directory\n##above_dir = \"D:/Kaggle/input/hms-harmful-brain-activity-classification/\"","metadata":{"papermill":{"duration":0.021426,"end_time":"2024-02-01T20:24:50.700256","exception":false,"start_time":"2024-02-01T20:24:50.678830","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:47.249468Z","iopub.execute_input":"2024-03-19T16:41:47.249916Z","iopub.status.idle":"2024-03-19T16:41:47.266422Z","shell.execute_reply.started":"2024-03-19T16:41:47.249864Z","shell.execute_reply":"2024-03-19T16:41:47.265168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n    \n## Things To Use","metadata":{"papermill":{"duration":0.01061,"end_time":"2024-02-01T20:24:47.591629","exception":false,"start_time":"2024-02-01T20:24:47.581019","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n# For k-means\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import silhouette_score\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# Simple random forest model:\nfrom sklearn.ensemble import RandomForestClassifier\n\n# Logistic Regression (to add features)\nfrom sklearn.linear_model import LogisticRegression","metadata":{"papermill":{"duration":3.065432,"end_time":"2024-02-01T20:24:50.667510","exception":false,"start_time":"2024-02-01T20:24:47.602078","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:47.268529Z","iopub.execute_input":"2024-03-19T16:41:47.268865Z","iopub.status.idle":"2024-03-19T16:41:48.429756Z","shell.execute_reply.started":"2024-03-19T16:41:47.268835Z","shell.execute_reply":"2024-03-19T16:41:48.428674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Functions, Etc.","metadata":{"papermill":{"duration":0.010495,"end_time":"2024-02-01T20:24:50.721613","exception":false,"start_time":"2024-02-01T20:24:50.711118","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Define these since they get used a lot\nHBA_number = 6\nHBA_names = [\"seizure\", \"lpd\", \"gpd\", \"lrda\", \"grda\", \"other\"]\nHBA_expert_names = [\"Seizure\",\"LPD\",\"GPD\",\"LRDA\",\"GRDA\",\"Other\"]\niHBA_of_expert = {\"Seizure\":0,\"LPD\":1,\"GPD\":2,\"LRDA\":3,\"GRDA\":4,\"Other\":5}\nHBA_votes = [\"seizure_vote\", \"lpd_vote\", \"gpd_vote\",\n             \"lrda_vote\", \"grda_vote\", \"other_vote\"]\nHBA_probs = [\"seizure_prob\", \"lpd_prob\", \"gpd_prob\",\n             \"lrda_prob\", \"grda_prob\", \"other_prob\"]\nthe4chains = [\"LL\",\"RL\",\"LP\",\"RP\"]\n\n# Output format for arrays\nnp.set_printoptions(precision=6, suppress=True)","metadata":{"papermill":{"duration":0.023834,"end_time":"2024-02-01T20:24:50.756102","exception":false,"start_time":"2024-02-01T20:24:50.732268","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:48.431168Z","iopub.execute_input":"2024-03-19T16:41:48.431639Z","iopub.status.idle":"2024-03-19T16:41:48.442373Z","shell.execute_reply.started":"2024-03-19T16:41:48.431606Z","shell.execute_reply":"2024-03-19T16:41:48.441556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def kld_score(solution, submission):\n    '''\n    Calculate the average KL divergence score.\n    Ignores the \"row id\" assumed in the first column.\n    '''\n    sumsum = 0.0\n    # Go through the probabilities\n    for prob_col in solution.columns.values:\n        sumsum += np.nansum(-1.0*solution[prob_col] *\n                        np.log(submission[prob_col] / solution[prob_col]))\n    return sumsum/(len(solution))","metadata":{"papermill":{"duration":0.02337,"end_time":"2024-02-01T20:24:50.790522","exception":false,"start_time":"2024-02-01T20:24:50.767152","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:48.443571Z","iopub.execute_input":"2024-03-19T16:41:48.444119Z","iopub.status.idle":"2024-03-19T16:41:48.453730Z","shell.execute_reply.started":"2024-03-19T16:41:48.444084Z","shell.execute_reply":"2024-03-19T16:41:48.452800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_hms_meta():\n    '''\n    Read in the train.csv and test.csv files.\n    Add total_vote, _prob columns, and vote entropy to train_meta.\n    Add extra cols to test to allow the same processing as train:\n        eeg[spectro]_sub_id, eeg[spectro]_label_offset_seconds, label_id\n    Make various plots of the train_meta values.\n    '''\n\n    # Read the test meta data\n    test_meta = pd.read_csv(above_dir+\"test.csv\")\n    test_meta_len = len(test_meta)\n    print(\"Test has length\", test_meta_len)\n    # Add columns to allow similar train/test processing:\n    test_meta[\"eeg_sub_id\"] = 0\n    test_meta[\"eeg_label_offset_seconds\"] = 0.0\n    test_meta[\"spectrogram_sub_id\"] = 0\n    test_meta[\"spectrogram_label_offset_seconds\"] = 0.0\n    test_meta[\"label_id\"] = test_meta.eeg_id\n    # Can decide to replace the not-real test with training data instead\n    if test_meta_len > 1:\n        REAL_TEST = True\n    else:\n        REAL_TEST = False\n        print(\"  --> not the real LB test data.\\n\")\n        # Replace the test_meta?\n        pass\n  \n    # Read the train meta data\n    train_meta = pd.read_csv(above_dir+\"train.csv\")\n    train_meta_len = len(train_meta)\n    print(\"Train has length\", train_meta_len, \" with:\")\n    \n    # Add a total_vote column\n    train_meta[\"total_vote\"] = ( train_meta[\"seizure_vote\"] +\n                    train_meta[\"lpd_vote\"] + train_meta[\"gpd_vote\"] +\n                    train_meta[\"lrda_vote\"] + train_meta[\"grda_vote\"] +\n                                    train_meta[\"other_vote\"] )\n    # Add a max_vote column (i.e. number of votes in the expert consensus)\n    train_meta[\"max_vote\"] = np.max(np.array([train_meta[\"seizure_vote\"] ,\n                    train_meta[\"lpd_vote\"] , train_meta[\"gpd_vote\"] ,\n                    train_meta[\"lrda_vote\"] , train_meta[\"grda_vote\"] ,\n                    train_meta[\"other_vote\"]]), axis=0)\n    \n    # Show various unique numbers\n    for this_col in [\"label_id\",\"eeg_id\",\"spectrogram_id\",\n                     \"patient_id\",\"total_vote\"]:\n        print(\"   \", len(train_meta[this_col].unique()),\n            \"unique \"+this_col+\" values.\")\n    \n    # Look at the total votes values:  1 to 28 with missing 8 and 9\n    print(\"\\nHistogram of the total votes\")\n    plt.figure(figsize=(6,3))\n    plt.hist(train_meta[\"total_vote\"],bins=55,log=True)\n    plt.title(\"Histogram of Total Votes\")\n    plt.show()\n    \n    # Distribution of the \"Expert consensus\" for less/more than 9 votes\n    allvt = train_meta.expert_consensus.value_counts()\n    less9 = train_meta[train_meta[\"total_vote\"] < 9].expert_consensus.value_counts()\n    more9 = train_meta[train_meta[\"total_vote\"] > 9].expert_consensus.value_counts()\n    vnot3 = train_meta[train_meta[\"total_vote\"] != 3].expert_consensus.value_counts()\n    ##print(\"Counts for votes less than 9:\\n\",less9[allvt.index])\n    ##print(\"\\nCounts for votes more than 9:\\n\",more9[allvt.index])\n    ##print(\"\\nCounts for votes not equal to 3:\\n\",vnot3[allvt.index])\n    \n    # Create _prob values from the _vote values\n    ##print(\"\\nHistograms of the probabilites of the different HBAs:\")\n    ##print(\"   (note that the large Prob=0 bin is not included.)\")\n    for col_pre in HBA_names:\n        train_meta[col_pre + \"_prob\"] = (train_meta[col_pre + \"_vote\"] / \n                                     train_meta[\"total_vote\"] )\n        # Show the probability histogram for each type\n        if False:\n            plt.figure(figsize=(6,1.5))\n            plt.hist(train_meta[col_pre + \"_prob\"],bins=20,range=(0.02,1))\n            plt.ylim(0,len(train_meta)/5)\n            plt.title(\"Histogram of   \"+col_pre+\"_prob\")\n            plt.show()\n        \n    # Calculate the entropy for each row (~ amount of vote variation)\n    print(\"Calculating voting entropy values ...\")\n    def calc_entropy(row):\n        the_probs = np.clip(row[16:21+1].values.astype(float), 1.e-8,1.0)\n        return np.nansum(the_probs * -1*np.log(the_probs))\n    # Add an entropy column\n    train_meta[\"entropy\"] = train_meta.apply(calc_entropy, axis=1)\n    \n    return train_meta, test_meta\n","metadata":{"papermill":{"duration":0.037337,"end_time":"2024-02-01T20:24:50.838932","exception":false,"start_time":"2024-02-01T20:24:50.801595","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:48.457409Z","iopub.execute_input":"2024-03-19T16:41:48.458081Z","iopub.status.idle":"2024-03-19T16:41:48.476363Z","shell.execute_reply.started":"2024-03-19T16:41:48.458044Z","shell.execute_reply":"2024-03-19T16:41:48.475434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prob_prob_scatter(name1, name2, probs2plot, clust_ids, iclust_order=[0]):\n    '''\n    Make a prob1 vs prob2 scatter plot.\n    - name1, name2 are 2 of the 6 HBA expert_consensus labels.\n    - probs2plot is 6-column dataframe, e.g., train_meta[HBA_probs]\n    - clust_ids is an array of cluster id integers, e.g., train_meta[\"clust_id\"]\n    Include an x at the cluster centers in the chosen axes.\n    Use sqrt scaling to emphasize lower values.\n    External: HBA_probs, iHBA_of_expert[ ], clust_centers\n    '''\n    # Color-code the vectors by their k-means label, this is in HBA order\n    hba_clrs = [\"orange\",\"blue\",\"red\",\"black\",\"green\",\"purple\",\n              \"green\",\"red\",\"blue\",\"orange\"]  # up to 10 clusters\n    if len(iclust_order) > 2:\n        # adjust the colors to the cluster order\n        kmclrs = hba_clrs.copy()\n        for iord, iclust in enumerate(iclust_order):\n            kmclrs[iclust] = hba_clrs[iord]\n   \n    clstclrs = []\n    for ilab in clust_ids:\n        clstclrs.append(kmclrs[ilab])\n        \n    ixax = iHBA_of_expert[name1]\n    iyax = iHBA_of_expert[name2]\n    lenprob = len(probs2plot)\n    plt.figure(figsize=(5,5))\n    plt.scatter(np.sqrt(probs2plot[HBA_probs[ixax]]) + \n                0.04*(np.random.rand(lenprob)-0.5),\n            np.sqrt(probs2plot[HBA_probs[iyax]]) + \n                0.04*(np.random.rand(lenprob)-0.5),\n           s=3, c=clstclrs, alpha=0.02)\n    # Add the centers\n    for iclust in range(0,len(clust_centers)):\n        plt.plot(np.sqrt([clust_centers[iclust,ixax]]),\n                 np.sqrt([clust_centers[iclust,iyax]]),\n                 c=kmclrs[iclust],marker=\"x\",markersize=15)\n    plt.xlabel(\"sqrt( \"+name1+\" )\")\n    plt.ylabel(\"sqrt( \"+name2+\" )\")\n    plt.show()\n    return kmclrs  # returns cluster colors appropriate for iclust","metadata":{"papermill":{"duration":0.029957,"end_time":"2024-02-01T20:24:50.879699","exception":false,"start_time":"2024-02-01T20:24:50.849742","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:48.477998Z","iopub.execute_input":"2024-03-19T16:41:48.479136Z","iopub.status.idle":"2024-03-19T16:41:48.493843Z","shell.execute_reply.started":"2024-03-19T16:41:48.479091Z","shell.execute_reply":"2024-03-19T16:41:48.492892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def assemble_features(meta_frame, traintest=\"train\", smooth_width=5,\n                     ):\n    '''\n    Create a dataframe of spectrogram features from the meta_frame rows.\n    Will include clust_id (i.e, the y) if it is in the input meta_frame.\n    Assumes these are available: above_dir, the4chains\n    '''\n    # Make a quick figure when doing validation events\n    if traintest == \"validation\":\n        plt.figure(figsize=(9,7))\n    \n    # The frequencies and general trend (for normalization)\n    freqs = np.array(range(100))*0.19525 + 0.59\n    spect_trend = 150.0/(1.0**2.3+freqs**(2.3))\n    # Indices to use when down-selecting after the smoothing, includes first bin\n    baseinds = np.insert(np.arange(int((smooth_width-1)/2), 100, smooth_width),\n                             0, 0)  # bin 0 put at front\n    # Extend the freqs and trend to 4 copies to do all four \"chains\" at once\n    freqs4 = np.array(4*list(freqs))\n    spect_trend4 = np.array(4*list(spect_trend))\n    \n    # Setup for fft: 256 2s bins gives a time range of +/- 4.27 minutes\n    n_fft_bins = 256\n    # Doing the FFT of the time series of EEG amplitude at 4 freq.s\n    # Use:  7(2.0Hz), 15(3.5Hz), 23(5.0Hz)m 60(12Hz), averaged +/-1, +/-2\n    freqbins = [7, 15, 23, 60]\n    freqclrs = [\"red\",\"blue\",\"green\",\"yellow\"]  # for quick plots.\n    # Apodizing window for this\n    apod_wind = np.blackman(n_fft_bins)\n    # Which bins of the fft to use for features\n    fftfeatbins = [4, 9, 16, 25, 36, 49]\n    \n    # Go through the rows and extract the features from spectrogram\n    feats_frame=[]\n    # Save last spectro_id_str, skip a read if it's the same\n    last_spectro_id_str = \"starting\"\n    # show status about every 1/15 of total in units of 100s\n    print_every_nth = max([100, 100*int(0.5+len(meta_frame.index)/(100.0*15))])\n    # plot about 100 rows, each row plots 4x4x6 = 96 spectra\n    plot_every_nth = max([1, int(0.5 + len(meta_frame.index)/100 )])\n    for irow in meta_frame.index:\n        this_row = meta_frame.loc[irow]\n        # Get this row's spectrogram\n        spectro_id_str = str(int(this_row.spectrogram_id)) # make sure int\n        if spectro_id_str != last_spectro_id_str:\n            if traintest != 'test':\n                spectro_file = above_dir+\"train_spectrograms/\"+spectro_id_str+\".parquet\"\n            else:\n                spectro_file = above_dir+\"test_spectrograms/\"+spectro_id_str+\".parquet\"\n            pads_spectro = pads.dataset(spectro_file)\n            this_spectro = pads_spectro.to_table().to_pandas()\n        last_spectro_id_str = spectro_id_str\n        # the time (row) offset of the event center\n        loc_offset = int(this_row.spectrogram_label_offset_seconds/2)\n        fftlocbeg = int(loc_offset + 149 - (n_fft_bins/2 -1))\n        fftlocend = int(loc_offset + 150 + (n_fft_bins/2 -1))\n        \n        # - - -\n        # The middle 8 seconds spectra\n        # Skip first column, the rest include 100 each for \"LL\",\"RL\",\"LP\",\"RP\"\n        middle8s = ((this_spectro.iloc[loc_offset+148, 1:] +\n                     this_spectro.iloc[loc_offset+149, 1:] +\n                     this_spectro.iloc[loc_offset+150, 1:] +\n                     this_spectro.iloc[loc_offset+151, 1:])/\n                        (4.0*spect_trend4)) # 4 copies of trend\n        # Clip and clean the values of 0s, etc.\n        middle8s = np.clip(middle8s,0.001,1000.0)\n        middle8s = middle8s.replace([np.nan, -np.inf, np.inf], 0.001)\n        # Get means/medians for overall and the individual \"LL\",\"RL\",\"LP\",\"RP\"\n        spect_mean = np.mean(middle8s)  # over all 4 spectra\n        spect_median = np.median(middle8s)\n        the4means=[]\n        the4medians = []\n        for ispec in range(4):\n            ibeg = ([0,100,200,300])[ispec]\n            iend = ibeg+100\n            the4means.append(np.mean(middle8s[ibeg:iend]))\n            the4medians.append(np.median(middle8s[ibeg:iend]))\n        # Scale all the spectra relative to the 4-spectra mean\n        # and take the log10\n        middle8spre = np.log10(middle8s / spect_mean)\n        # smooth the logged values with smooth_width\n        middle8s = middle8spre.rolling(smooth_width, min_periods=smooth_width, \n                                        center=True, closed=None).mean()\n        # Set the values of the first few in each spectrum to pre-smoothed values\n        for ioff in range(0, 399, 100):\n            for ibin in range(int((smooth_width-1)/2)):  # assumes width is odd\n                middle8s[ibin+ioff] = middle8spre[ibin+ioff]\n        # downselect based on the smoothing width; includes the first value too\n        select_inds = np.concatenate((baseinds,100+baseinds,\n                                      200+baseinds,300+baseinds))\n        freqs4ds = freqs4[select_inds]\n        middle8sds = middle8s[select_inds]\n        # Put the middle8s ds values in a dataframe\n        middle_feats = middle8sds.to_frame().T\n        \n        # - - -\n        # The Flare-Rate Spectra: FFT of the time series of EEG amplitude at a freq.\n        flarecols = []\n        flarevals = []\n        for ispec in range(4):\n            for ifreq, freqbin in enumerate(freqbins):\n                sum_spect_trend = sum(spect_trend[freqbin-4:freqbin+4+1:2])\n                ifreqoff = freqbin + 1 + ispec*100  # column\n                amplvstime = (this_spectro.iloc[fftlocbeg:fftlocend+1, ifreqoff-4] +\n                          this_spectro.iloc[fftlocbeg:fftlocend+1, ifreqoff-2] +\n                          this_spectro.iloc[fftlocbeg:fftlocend+1, ifreqoff] +\n                          this_spectro.iloc[fftlocbeg:fftlocend+1, ifreqoff+2] +\n                          this_spectro.iloc[fftlocbeg:fftlocend+1, ifreqoff+4]\n                         )/sum_spect_trend\n                # Clean the values of 0s, etc. - NaNs etc are 0.001\n                amplvstime = np.clip(amplvstime,0.001,1000.0)\n                amplvstime = amplvstime.replace([np.nan, -np.inf, np.inf], 0.001)\n                # Normalize to the middle bins, 127, 128 from start\n                amplvstime = 2.0*amplvstime/(amplvstime[fftlocbeg+127]+\n                                             amplvstime[fftlocbeg+128])\n                # Apply a further clip and the apodizing window\n                amplvstime = apod_wind * np.clip(amplvstime, 0.0, 10.0)\n                # Apply some smoothing, reduce high freq.s\n                for iroll in range(2):  # 2 x 2 smoothing gives 1/4, 1/2, 1/4\n                    amplvstime = amplvstime.rolling(2, min_periods=1, \n                                        center=False, closed=None).mean()\n                # Take the fft, it's magnitude, use first half\n                amplfft = np.abs(np.fft.fft(amplvstime))[0:int(n_fft_bins/2)]\n                # Smooth it and take log10(1+...)]\n                amplfft = np.log10(1.0 + np.array(pd.Series(amplfft).rolling(5,\n                            min_periods=1, center=True, closed=None).mean()))\n                # Make feature names and get the values to use as features\n                for fftfeatbin in fftfeatbins:   \n                    flarecols.append(\"fft-\"+the4chains[ispec]+\"{:.1f}-\".format(freqs[freqbin])+\n                              str(fftfeatbin))\n                    flarevals.append(amplfft[fftfeatbin])\n                # Plot the fft spectra for some rows, colored by cluster\n                if (len(feats_frame) % plot_every_nth == 0) and traintest == 'validation':\n                    plt.plot(np.sqrt(range(int(n_fft_bins/2))), amplfft,\n                            c=kmclrs[this_row[\"clust_id\"]],lw=2,alpha=0.01)\n        \n        # Put the flarecols and flarevals into a one-row dataframe\n        flare_feats = pd.DataFrame([flarevals], columns=flarecols)\n        # Combine the two sets of features\n        these_feats = pd.concat([middle_feats, flare_feats], axis=1)\n                         \n        # - - - \n        # Include log10()s of means/medians\n        spect_mean = np.log10(spect_mean)\n        spect_median = np.log10(spect_median)\n        the4means = np.log10(the4means)\n        the4medians = np.log10(the4medians)\n        these_feats[\"Mean\"] = spect_mean\n        these_feats[\"Median\"] = spect_median\n        for ispec in range(4):\n            these_feats[the4chains[ispec]+\"mean\"] = the4means[ispec]\n            these_feats[the4chains[ispec]+\"median\"] = the4medians[ispec]\n        # Include the clust_id (the y) if present in meta_frame\n        if \"clust_id\" in meta_frame.columns:\n            these_feats[\"clust_id\"] = this_row.clust_id\n        # If this is the first row make it the features frame\n        if len(feats_frame) == 0:\n            feats_frame = these_feats.copy()\n        else:\n            feats_frame = pd.concat([feats_frame, these_feats])\n        # Show progress\n        if len(feats_frame) % print_every_nth == 0:\n            print(\"... {} done...\".format(len(feats_frame)))\n    # All done with desired rows\n    if traintest == \"validation\":\n        plt.ylim(-0.05,2.25)\n        plt.xlim(0.0,11.5)\n        plt.xlabel(\"sqrt[  FFT bin number ]\")\n        plt.ylabel(\"log10[1+ FFT amplitude ]\")\n        plt.title(\"Flare-Rate Spectra (colored by cluster)\")\n\n    # Return the features (includes clust_id if given)\n    return feats_frame.reset_index().drop(columns=[\"index\"])\n   ","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:41:48.495411Z","iopub.execute_input":"2024-03-19T16:41:48.495725Z","iopub.status.idle":"2024-03-19T16:41:48.536102Z","shell.execute_reply.started":"2024-03-19T16:41:48.495697Z","shell.execute_reply":"2024-03-19T16:41:48.534942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def find_best_tamed_kl():\n    '''\n    Adjust the taming fraction for each cluster center to optimize KL\n    Assumed inputs in environment:\n        submission, pred_ids, solution\n    Assumed useful values available:\n        clust_centers, NUM_CLUSTS, HBA_number, HBA_votes\n    '''\n    # Adjust, tame, the cluster_center values?\n    mean_all_probs = np.array([0.208319, 0.132120, 0.128532,\n                           0.138913, 0.179294, 0.212822])\n    # Start all with the mean values, i.e., tamed_frac=0.\n    tamed_fracs = 0.0 * np.ones(NUM_CLUSTS)\n    tamed_centers = clust_centers.copy()\n    for iclust in range(NUM_CLUSTS):\n        tamed_centers[iclust,:] = mean_all_probs\n    # Go through the clusters\n    for iclust in range(NUM_CLUSTS):\n        last_kl = 10.0\n        # Find the best overall tame frac for this cluster\n        for this_frac in np.arange(0.03, 1.00, 0.05):  # 0.03--0.98\n            tamed_fracs[iclust] = this_frac\n            this_cent = (tamed_fracs[iclust]*clust_centers[iclust,:] +\n                         (1.0-tamed_fracs[iclust])*mean_all_probs)\n            tamed_centers[iclust,:] = this_cent\n\n            # Put the probs into the submission column by column\n            for iprob in range(HBA_number):\n                this_col_probs = tamed_centers[: , iprob]\n                submission[HBA_votes[iprob]] = this_col_probs[pred_ids]\n            # Evaluate the KL divergence using predicted cluster ids\n            this_kl = kld_score(solution, submission)\n            if (this_kl < last_kl):\n                best_fracs = tamed_fracs.copy()\n                best_centers = tamed_centers.copy()\n                last_kl = this_kl\n            else:\n                # restore the best one\n                tamed_fracs[iclust] = best_fracs[iclust]\n                tamed_centers[iclust, :] = best_centers[iclust, :]\n                break\n                # go to the next cluster\n        #\n    # done with all clusters\n    return best_fracs, best_centers","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:41:48.537824Z","iopub.execute_input":"2024-03-19T16:41:48.538242Z","iopub.status.idle":"2024-03-19T16:41:48.551261Z","shell.execute_reply.started":"2024-03-19T16:41:48.538209Z","shell.execute_reply":"2024-03-19T16:41:48.550164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n\n## Get (and look at) the Meta Data","metadata":{"papermill":{"duration":0.010211,"end_time":"2024-02-01T20:24:50.922057","exception":false,"start_time":"2024-02-01T20:24:50.911846","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Read in the meta data, routine also looks at the train values\ntrain_meta, test_meta = read_hms_meta()","metadata":{"papermill":{"duration":14.696415,"end_time":"2024-02-01T20:25:05.629228","exception":false,"start_time":"2024-02-01T20:24:50.932813","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:41:48.552840Z","iopub.execute_input":"2024-03-19T16:41:48.553276Z","iopub.status.idle":"2024-03-19T16:42:02.533327Z","shell.execute_reply.started":"2024-03-19T16:41:48.553232Z","shell.execute_reply":"2024-03-19T16:42:02.532136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at the values of the spectro sub id\n##print(\"Spectrogram_sub_id unique values:\",\n##            np.sort(train_meta.spectrogram_sub_id.unique()))\n# Show a histogram\nplt.figure(figsize=(5,2))\nplt.hist(train_meta[\"spectrogram_sub_id\"],bins=55,log=True)\nplt.title(\"Histogram of spectrogram_sub_id\")\nplt.show()\n\n# Look at the values of the eeg sub id\n##print(\"eeg_sub_id unique values:\",np.sort(train_meta.eeg_sub_id.unique()))\n# Show a histogram\nplt.figure(figsize=(5,2))\nplt.hist(train_meta[\"eeg_sub_id\"],bins=55,log=True)\nplt.title(\"Histogram of eeg_sub_id\")\nplt.show()","metadata":{"papermill":{"duration":1.531268,"end_time":"2024-02-01T20:25:07.172003","exception":false,"start_time":"2024-02-01T20:25:05.640735","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:42:02.535309Z","iopub.execute_input":"2024-03-19T16:42:02.535774Z","iopub.status.idle":"2024-03-19T16:42:03.702934Z","shell.execute_reply.started":"2024-03-19T16:42:02.535728Z","shell.execute_reply":"2024-03-19T16:42:03.701625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n\n## Clustering of HBA Probability Vectors","metadata":{"papermill":{"duration":0.012,"end_time":"2024-02-01T20:25:07.220923","exception":false,"start_time":"2024-02-01T20:25:07.208923","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Choose the number of clusters\nnum_clusts = NUM_CLUSTS   # 6 to 10\n\n# Choose which rows to use for clustering\n# Use all the HBA samples, except the ends of the very long ones \nclust_rows_bool =(train_meta.eeg_sub_id < 200)","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:42:03.704649Z","iopub.execute_input":"2024-03-19T16:42:03.705761Z","iopub.status.idle":"2024-03-19T16:42:03.712339Z","shell.execute_reply.started":"2024-03-19T16:42:03.705710Z","shell.execute_reply":"2024-03-19T16:42:03.710961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get the vectors and run k-means clustering on them\nprob_vectors = train_meta.loc[clust_rows_bool, HBA_probs]\n\n# Look at the votes distribution of the selected prob_vectors\nprint(\"\\nUsing {} HBA samples for clustering.\".format(len(prob_vectors)))\nprint(\"These include {} unique eeg_ids\".format(\n                train_meta.loc[clust_rows_bool, \"eeg_id\"].nunique()),\n                \"and {} unique patient ids.\".format(\n                train_meta.loc[clust_rows_bool, \"patient_id\"].nunique()))\n##plt.figure(figsize=(6,3))\n##plt.hist(train_meta.loc[clust_rows_bool,\"total_vote\"],bins=55,log=True)\n##plt.title(\"Histogram of total_vote in HBA samples clustered\")\n##plt.show()\n    \n# Find the centers for a specified number of clusters \nprob_array = np.array(prob_vectors)\nkmeans = KMeans(n_clusters = num_clusts, init = 'k-means++', n_init = 10,\n                    max_iter = 300, random_state=None)\nkmeans.fit(prob_array)\n\n# Get the cluster center coord.s/probs\nclust_centers = kmeans.cluster_centers_\n# They are very close to normalized, make them closer  :)\nfor iclust in range(NUM_CLUSTS):\n    clust_centers[iclust ,:] = (clust_centers[iclust ,:]/\n                            np.sum(clust_centers[iclust ,:]))\nprint(\"cluster centers:\")\nprint(clust_centers)\n\n# Map iclust to standard order of the HBAs followed by hybrids\niclust_of_order = []\n# The first dim_centers (6) are the close-to-unit vectors\nfor icol in range(HBA_number):\n    iclust_of_order.append(np.argmax(clust_centers[:,icol]))\n# Add the remaining num_clusts - HBA_number (3) in order by max component\n# index ordered by largest to smallest\nclust_by_max = np.argsort(-1*np.max(clust_centers,axis=1))\n# Add the last ones:\nfor iord in range(HBA_number,num_clusts):\n    iclust_of_order.append(clust_by_max[iord])    \n    \n# Make some names for the clusters in standard order\nkmnames = HBA_expert_names.copy()\nfor ihyb in range(1,(num_clusts - HBA_number)+1):\n    kmnames.append(\"Hybrid-\"+str(ihyb))\n","metadata":{"papermill":{"duration":2.481754,"end_time":"2024-02-01T20:25:09.714848","exception":false,"start_time":"2024-02-01T20:25:07.233094","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:42:03.714310Z","iopub.execute_input":"2024-03-19T16:42:03.714811Z","iopub.status.idle":"2024-03-19T16:42:05.758917Z","shell.execute_reply.started":"2024-03-19T16:42:03.714763Z","shell.execute_reply":"2024-03-19T16:42:05.757865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Not sure what is the best way to view the cluster results.\n# Try scatter plots of pairs of the HBA probs, use sqrt scaling.\n\n# Use all the data: assign the clusters to all of the train_meta\ntrain_meta[\"clust_id\"] = kmeans.predict(np.array(train_meta[HBA_probs]))\n\nall_probs = train_meta[HBA_probs]\nall_ids = train_meta[\"clust_id\"]\n\nif True:\n    kmclrs = prob_prob_scatter(\"Seizure\",\"GPD\", all_probs, all_ids, iclust_of_order)\n    \n    kmclrs = prob_prob_scatter(\"LPD\",\"GRDA\", all_probs, all_ids, iclust_of_order)\n\n    kmclrs = prob_prob_scatter(\"LRDA\",\"Other\", all_probs, all_ids, iclust_of_order)\n\n# List the centers in order with other information\n# Counts in the clusters\nclust_counts = train_meta.clust_id.value_counts()\n##print(\"   Seizure\",\"  LPD\",\"     GPD\",\"     LRDA\",\"    GRDA\",\"    Other\",\n##     \" Clust name\",\"  color\",\" Number in clust\")\n##for iord, iclust in enumerate(iclust_of_order):\n##    print(clust_centers[iclust,:],\" {} {:10} {:6}  {:5}\".format(\n##        iclust, kmnames[iord],kmclrs[iclust],clust_counts[iclust]))\n    \n# And/or show them graphically\n# Also create the reverse mapping from iclust to the usual order\niorder_of_clust = num_clusts*[-1]\nfor iord, iclust in enumerate(iclust_of_order):\n    ##print(clust_centers[iclust,:],\" {} {:10} {:6}  {:5}\".format(\n    ##    iclust, kmnames[iord],kmclrs[iclust],clust_counts[iclust]))\n    iorder_of_clust[iclust] = iord\n    if iord == 0: print(\"The close-to-unit-vector cluster centers:\")\n    if iord == 6: print(\"The Hybrid cluster centers:\")\n    plt.figure(figsize=(5,1))\n    plt.bar(HBA_expert_names, clust_centers[iclust,:],color=kmclrs[iclust])\n    plt.ylim(-0.01,1.01)\n    plt.title(kmnames[iord]+\"  kmclust={} has {} samples\".format(\n        iclust,clust_counts[iclust]), size='medium')\n    if iord < 5: plt.xticks([])\n    plt.show()\n    ","metadata":{"papermill":{"duration":6.6687,"end_time":"2024-02-01T20:25:16.398495","exception":false,"start_time":"2024-02-01T20:25:09.729795","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:42:05.760852Z","iopub.execute_input":"2024-03-19T16:42:05.761264Z","iopub.status.idle":"2024-03-19T16:42:14.242195Z","shell.execute_reply.started":"2024-03-19T16:42:05.761231Z","shell.execute_reply":"2024-03-19T16:42:14.241342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#     Cluster variables summary\n# clust_centers[iclust, HBA probs]  iclust is in k-means order\n# train_meta[\"clust_id\"] is k-means iclust \n# kmclrs are cluster colors in k-means order, iclust\n# iclust_of_order  gives the iclust of the standard order\n# iorder_of_clust  gives the standard order of iclust\n# kmnames are in standard order","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:42:14.247041Z","iopub.execute_input":"2024-03-19T16:42:14.247588Z","iopub.status.idle":"2024-03-19T16:42:14.252077Z","shell.execute_reply.started":"2024-03-19T16:42:14.247553Z","shell.execute_reply":"2024-03-19T16:42:14.251028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### KL score using cluster centers\nCalculate the KL score if all samples have their cluster correctly identified. This KL error is the error due to using a single probabilty set for each cluster type with no tuning to the individual sample.","metadata":{"papermill":{"duration":0.016719,"end_time":"2024-02-01T20:25:16.437745","exception":false,"start_time":"2024-02-01T20:25:16.421026","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Create a 'solution' of actual probabilities from the training data\nsolution_train = train_meta[[\"eeg_id\"] + HBA_votes]\n# In solution_train replace the _votes values with probabilities\nfor col_pre in HBA_names:\n    solution_train.loc[:, col_pre + \"_vote\"] = train_meta[col_pre + \"_prob\"]\n\n# Create a 'submission' of the cluster probabilities\nsubmission_train = solution_train.copy()\n# Use clust_ids to select the correct clust_probs vector\nclust_ids = train_meta[\"clust_id\"]\nif False:\n    # How bad could it be? (with random clusters assigned) Ans: Bad ~ 3.0\n    clust_ids = np.random.choice(9, size=len(train_meta),\n                             replace=True, p=None)\n# Put the probs into the submission\nfor iprob in range(HBA_number):\n    this_col_probs = clust_centers[: , iprob]\n    submission_train[HBA_votes[iprob]] = this_col_probs[clust_ids]\n\n# Evaluate the KL divergence when predictions are the cluster probs\nprint(\"Score if HBA samples are correctly assigned cluster prob.s:\",\n      np.round(kld_score(solution_train, submission_train),4))","metadata":{"papermill":{"duration":0.077321,"end_time":"2024-02-01T20:25:16.531975","exception":false,"start_time":"2024-02-01T20:25:16.454654","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:42:14.253633Z","iopub.execute_input":"2024-03-19T16:42:14.253982Z","iopub.status.idle":"2024-03-19T16:42:14.309602Z","shell.execute_reply.started":"2024-03-19T16:42:14.253951Z","shell.execute_reply":"2024-03-19T16:42:14.308437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n\n## Look at the Spectrograms","metadata":{"papermill":{"duration":0.018243,"end_time":"2024-02-01T20:25:16.635727","exception":false,"start_time":"2024-02-01T20:25:16.617484","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Choose a smoothing window\nsmooth_width = SMOOTH_WIDTH   # Here smooth_width is just for looking,\n                              # can be a value different from SMOOTH_WIDTH\n\n# Choose what rows to show spectra from\nspectro_meta = train_meta[clust_rows_bool].copy()\n# Downsample them by a factor\nevery_nth = 97\n\n# Frequencies of the 100 bins are LL_0.59 LL_0.78 LL_0.98 . . . LL_19.73 LL_19.92\n# Roughly freq.s = 0.59 + 0.19525*ibin\nfreqs = np.array(range(100))*0.19525 + 0.59\n# Down selected freqs\n# Frequencies when downselected, it also includes the first freq bin. \ndownsel_freqs = np.insert(freqs[int((smooth_width-1)/2):100:smooth_width], 0, freqs[0])\n\n# The signal amplitudes vary with frequency - divide out a general trend.\n# Scaled it so the meadian(all_me[di]ans) is near 1.0\n# Power-law with blowup tamed near 0.\nspect_trend = 150.0/(1.0**2.3+freqs**(2.3))\n# Exponential with \"slope\" by eye, lacks increase near 0.\n##spect_trend = 5.0*np.exp(-1.0*freqs/5.0)\n\n# The FFT of the time series of EEG amplitude at some frequencies\n# Use: 10(2.5Hz), 41(8.5Hz), 60(12Hz), 84(17Hz), averaged +/-2, +/-4\n# or    7(2.0Hz)  18(4.1Hz)\n# Use:  7(2.0Hz), 15(3.5Hz), 23(5.0Hz)m 60(12Hz), averaged +/-1, +/-2\nfreqbins = [7, 15, 23, 60]\n# Colors for the 4 frequencies to be fft'ed\nfreqclrs = [\"red\",\"blue\",\"green\",\"yellow\"]\n# 256 2s bins gives a time range of +/- 4.27 minutes\nn_fft_bins = 256\n# Apodizing window for this\napod_wind = np.blackman(n_fft_bins)","metadata":{"papermill":{"duration":0.081997,"end_time":"2024-02-01T20:25:16.735270","exception":false,"start_time":"2024-02-01T20:25:16.653273","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:42:14.311211Z","iopub.execute_input":"2024-03-19T16:42:14.312019Z","iopub.status.idle":"2024-03-19T16:42:14.378435Z","shell.execute_reply.started":"2024-03-19T16:42:14.311971Z","shell.execute_reply":"2024-03-19T16:42:14.376927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plots from the Spectrograms\n# Select plots to be:\n#  - the EEG Frequency Spectra averaged from the middle 8 seconds of the event\n\n#  - the Flare-Rate Spectra of the EEG amplitude vs time at selected frequencies\n\n# Make lists of these values for each spectrogram\nall_medians = []\nall_means = []\nall_clusts = []\n\n# Setup multiple figures: one row for each cluster with spectra and ffts\nfig = plt.figure(figsize=(10,3.5*NUM_CLUSTS))\ngs = fig.add_gridspec(NUM_CLUSTS, 2, hspace=0.05, wspace=0.25)\npltaxs = gs.subplots(sharex='col')\n\n# Colors for the 4 chains [\"LL\",\"RL\",\"LP\",\"RP\"], and freq.a\nchainclrs = [\"red\",\"blue\",\"green\",\"yellow\"]\n\n# Go though the desired rows; monitor the progress     \nprint_every_nth = max([100, 100*int(0.5+(len(spectro_meta)/every_nth)/(100.0*15))])\n\n##print(\"Plotting the 4 spectra from each cluster type\")\n##iloc_rows = [97, 2003, 3003, 7001, 20112, 131, 45, 73, 227]\nprint(\"Plotting {} x 4 processed spectra\".format(int(len(spectro_meta)/every_nth)))\niloc_rows = range(0,len(spectro_meta),every_nth)\n\n# Adjust alpha depending on number of rows\nplt_alpha = min([0.25*(200/len(iloc_rows)), 0.30])\nfor ilocrow in iloc_rows:\n    this_row = spectro_meta.iloc[ilocrow]\n    this_clust = this_row.clust_id   # the k-means iclust\n    # Use the standard cluster order to index the plots\n    this_std_clust = iorder_of_clust[this_row.clust_id]\n    if True:\n        # Get this row's spectrogram\n        spectro_id_str = str(this_row.spectrogram_id)\n        spectro_file = above_dir + \"train_spectrograms/\"+spectro_id_str+\".parquet\"\n        pads_spectro = pads.dataset(spectro_file)\n        this_spectro = pads_spectro.to_table().to_pandas()\n        # offset into the spectro, loc_offset+149,150 are center of event in time\n        loc_offset = int(this_row.spectrogram_label_offset_seconds/2) # 2s bins\n        locbeg = int(loc_offset + 149 - (n_fft_bins/2 -1))\n        locend = int(loc_offset + 150 + (n_fft_bins/2 -1))\n        # Go through the 4 spectra, variable\n        for ispec in range(4):\n            \n            # - - -\n            # The middle 8s spectra and means, medians\n            # The columns for this spectrum\n            ibeg = ([1,101,201,301])[ispec]\n            iend = ibeg+100\n            # The average of middle four spectra, relative to the trend\n            middle8s = ((this_spectro.iloc[loc_offset+148, ibeg:iend] +\n                        this_spectro.iloc[loc_offset+149, ibeg:iend] +\n                        this_spectro.iloc[loc_offset+150, ibeg:iend] +\n                        this_spectro.iloc[loc_offset+151, ibeg:iend])\n                            /(4.0*spect_trend))\n            # Clean the values of 0s, etc. - NaNs etc are 0.001\n            middle8s = np.clip(middle8s,0.001,1000.0)\n            middle8s = middle8s.replace([np.nan, -np.inf, np.inf], 0.001)\n            # The log10(mean) and log10(median) -- mean/median before the log\n            spect_mean = np.log10(np.mean(middle8s))\n            all_means.append(spect_mean)\n            spect_median = np.log10(np.median(middle8s))\n            all_medians.append(spect_median)\n            all_clusts.append(this_clust)\n            # Take the log10 of the trended spectral values\n            middle8s = np.log10(middle8s)\n            # - - - All spectro values are log10() after this point.\n            # Further ratio it to its mean (or median) - subtract logs\n            middle8spre = middle8s - spect_mean\n            # smooth it with smooth_width (5 gives ~ 1 Hz)\n            middle8s = middle8spre.rolling(smooth_width, min_periods=smooth_width, \n                                        center=True, closed=None).mean()\n            # Set the values of the first few to unsmoothed\n            for ibin in range(int((smooth_width-1)/2)):  # assumes width is odd\n                middle8s[ibin] = middle8spre[ibin]\n            # plot this middle 8s spectrum\n            # colored by cluster, c=kmclrs[this_clust], or by chain, c=chainclrs[ispec]\n            pltaxs[this_std_clust,0].plot((freqs), middle8s,\n                               c=chainclrs[ispec],lw=2,alpha=plt_alpha)\n            \n            # - - -\n            # The Flare-Rate Spectra:\n            for ifreq, freqbin in enumerate(freqbins):\n                sum_spect_trend = sum(spect_trend[freqbin-2:freqbin+2+1:1])\n                ifreqoff = freqbin + 1 + ispec*100  # column for frequency/chain\n                amplvstime = (this_spectro.iloc[locbeg:locend+1, ifreqoff-2] +\n                          this_spectro.iloc[locbeg:locend+1, ifreqoff-1] +\n                          this_spectro.iloc[locbeg:locend+1, ifreqoff] +\n                          this_spectro.iloc[locbeg:locend+1, ifreqoff+1] +\n                          this_spectro.iloc[locbeg:locend+1, ifreqoff+2]\n                         )/(sum_spect_trend)\n                # Clean the values of 0s, etc. - NaNs etc are 0.001\n                amplvstime = np.clip(amplvstime,0.001,1000.0)\n                amplvstime = amplvstime.replace([np.nan, -np.inf, np.inf], 0.001)\n                # Normalize to the middle bins, 127, 128 from start\n                amplvstime = 2.0*amplvstime/(amplvstime[locbeg+127]+amplvstime[locbeg+128])\n                # Apply a further clip and the apodizing window\n                amplvstime = apod_wind * np.clip(amplvstime, 0.0, 10.0)\n                # Apply some smoothing, reduce high freq.s\n                for iroll in range(2):  # 0 0 1/4 1/2 1/4 0 0 smoothing\n                    amplvstime = amplvstime.rolling(2, min_periods=1, \n                                        center=False, closed=None).mean()\n                # Diagnostic: plot the amplitude vs time\n                ##pltaxs[this_std_clust, 1].plot(range(n_fft_bins),\n                ##                        (0.0+amplvstime),\n                ##                        c=freqclrs[ifreq],lw=2,alpha=plt_alpha/4.0)\n                # Take the fft, it's magnitude, use first half\n                amplfft = np.abs(np.fft.fft(amplvstime))[0:int(n_fft_bins/2)]\n                # Smooth it and take log10(1+...)]\n                amplfft = np.log10(1.0 + np.array(pd.Series(amplfft).rolling(5,\n                            min_periods=1, center=True, closed=None).mean()))\n                # plot the amplfft vs flare-rate\n                # colored by cluster, kmclrs[this_clust], chain,chainclrs[ispec],\n                #         or by freqbin c=freqclrs[ifreq]\n                pltaxs[this_std_clust, 1].plot(np.sqrt(range(int(n_fft_bins/2))),\n                        amplfft, c=freqclrs[ifreq],lw=2,alpha=plt_alpha/3.0)\n            \n    # Show progress\n    if (ilocrow/every_nth + 1) % print_every_nth == 0:\n        print(\"... {} done...\".format(int(ilocrow/every_nth)+1))\n        \nprint(\"\\n\\n Colors are: LL-red, RL-blue, LP-green, RP-yellow \"+\n             \"         f_Hz-color: 2.0-red, 3.5-blue, 5.0-green, 12.0-yell.\")\n\nfor iax in range(NUM_CLUSTS):\n    # Label the plot\n    if True:\n        # Middle8s plots\n        pltaxs[iax,0].text(12.5, 0.8, kmnames[iax], fontsize=14)\n        # Show the 1 line and locations of every smooth_width frequency for reference\n        pltaxs[iax,0].plot([0.0,20.0],[0.0,0.0],c='black',lw=3,alpha=0.2)\n        pltaxs[iax,0].plot(downsel_freqs, len(downsel_freqs)*[0.0],'.',c=\"gray\")\n        pltaxs[iax,0].set_ylim(-1.5,1.0)  # log already taken\n        pltaxs[iax,0].set_xlim(0.0,20.5)     \n        #\n        # Flare-Rate plots\n        # Show the log=0 line\n        pltaxs[iax,1].plot([0.0,len(amplvstime)],[0.0,0.0],c='black',lw=3,alpha=0.2)\n        pltaxs[iax,1].set_ylim(-0.2,2.0)  # log already taken\n        # x-axis, usual sqrt(fftbin #)\n        pltaxs[iax,1].set_xlim(0.0,11.5)\n        # x-axis Diagnostic pre-fft time series plot\n        ##pltaxs[iax,1].set_xlim(0.0,256.0)  # x-axis is time bin, 0 to 255        \n        \npltaxs[NUM_CLUSTS-1,0].set_ylabel(25*\" \"+\"log10[ Spectra /RefSpectrum /Mean\"+\n                                            \" & smoothed]\")\npltaxs[NUM_CLUSTS-1,0].set_xlabel(\"Frequency (Hz)\")\npltaxs[NUM_CLUSTS-1,1].set_ylabel(\"log10[1+ FFT amplitude ]\")\npltaxs[NUM_CLUSTS-1,1].set_xlabel(\"sqrt[  FFT bin number ]\")\npltaxs[0,0].set_title(\"Middle 4x2s Spectra (n_sm={}, colored by LL,RL,LP,RP)\".format(\n                            smooth_width), fontsize=11)\n#  colored by: freq. fft'ed   or   LL,RL,LP,RP\npltaxs[0,1].set_title(\"Flare-Rate Spectra (colored by freq. FFT'ed)\".format(\n                            freqs[ifreq]), fontsize=11)\nplt.show()","metadata":{"papermill":{"duration":46.794626,"end_time":"2024-02-01T20:26:03.615216","exception":false,"start_time":"2024-02-01T20:25:16.820590","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:42:14.380469Z","iopub.execute_input":"2024-03-19T16:42:14.381371Z","iopub.status.idle":"2024-03-19T16:45:46.976733Z","shell.execute_reply.started":"2024-03-19T16:42:14.381319Z","shell.execute_reply":"2024-03-19T16:45:46.974894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Mean and Median (before the divide by mean or median)\n\nprint(\"\\nMedian of the Means(below): {:.4f}\".format(np.median(all_means)))\nplt.figure(figsize=(6,2))\nplt.hist(np.clip(all_means,-2.0,2.0),bins=100)\nplt.xlim(-2.05,2.05)\nplt.xlabel(\"Mean\")\nplt.title(\"Histogram of the Means of the log(Spectra/Ref-spect)\")\nplt.show()\n\nprint(\"\\nMedian of the Medians(below): {:.4f}\".format(np.median(all_medians)))\nplt.figure(figsize=(6,2))\nplt.hist(np.clip(all_medians,-2.0,2.0),bins=100)\nplt.xlim(-2.05,2.05)\nplt.xlabel(\"log10( Median )\")\nplt.title(\"Histogram of the Medians of the log(Spectra/Ref-spect)\")\nplt.show()\n\n# Scatter plot show the means are greater than or equal to the medians: right skewed.\nmmclrs = []\nfor ilab in all_clusts:\n    mmclrs.append(kmclrs[ilab])\nplt.figure(figsize=(3,3))\nplt.scatter(all_medians, all_means, s=2,c=mmclrs,alpha=0.3)\nplt.xlim(-1.,1); plt.ylim(-1.,1)\nplt.xlabel(\"Median\")\nplt.ylabel(\"Mean\")\nplt.show()","metadata":{"papermill":{"duration":1.189304,"end_time":"2024-02-01T20:26:04.833590","exception":false,"start_time":"2024-02-01T20:26:03.644286","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:45:46.978313Z","iopub.execute_input":"2024-03-19T16:45:46.978711Z","iopub.status.idle":"2024-03-19T16:45:48.018216Z","shell.execute_reply.started":"2024-03-19T16:45:46.978676Z","shell.execute_reply":"2024-03-19T16:45:48.017024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n\n## Generate the Train and Validation Features\nFeatures are made from the spectrograms by assemble_features(). \n    \nThe extracted raw features are saved to csv files so they can be put in a dataset and read-in instead of extracting features each time. Currently features from v47 are read-in.","metadata":{}},{"cell_type":"code","source":"# Select rows for Training features\n# Use one (or more) per eeg by selecting some of the sub_id values.\ntrain_rows_bool = (((train_meta.eeg_sub_id < 50+1) &\n                     (train_meta.eeg_sub_id % 4 == 2)) |  # 2, 6, 10, ..., 50\n                    # Use about 3/4 of the eeg_sub_id=0 ones\n                    (train_meta.eeg_sub_id == 0) &\n                      ((train_meta.eeg_id%23)%8 > 1)) # include odd and even eeg_ids\n# Has 33152:  20583 non-zero eeg_sub_ids and 12569 ones that are 0.\nprint(\"Number of Training rows:\",sum(train_rows_bool))\n\n# Choose ones to use for Validation, select similar-to-training ones without overlap\n# Use one (or more) per eeg by selecting some of the (non-zero) sub_id values.\nvalid_rows_bool = (((train_meta.eeg_sub_id < 51+1) &\n                    (train_meta.eeg_sub_id % 8 == 3)) |  # 3, 11, ..., 51\n                    # Use the other about 1/4 of the eeg_sub_id=0 ones\n                    (train_meta.eeg_sub_id == 0) &\n                      ((train_meta.eeg_id%23)%8 < 2)) # include odd and even eeg_ids\n# Has 15960:  11440 non-zero eeg_sub_ids and 4520 ones that are 0.\nprint(\"Number of Validation rows:\",sum(valid_rows_bool))","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:48.019959Z","iopub.execute_input":"2024-03-19T16:45:48.020854Z","iopub.status.idle":"2024-03-19T16:45:48.075200Z","shell.execute_reply.started":"2024-03-19T16:45:48.020819Z","shell.execute_reply":"2024-03-19T16:45:48.074017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Get the Training Features","metadata":{}},{"cell_type":"code","source":"# Get train meta and features dataframes\nif USE_PREPROC:\n    # Restore previously saved data frames \"/kaggle/input/\n    # v47\n    # v62 Feat.s: Middle 8s spectra(84), Flare FFTs (96), Me[di]ans and clust_id (11)\n    # v66\n    Xy_train_meta = pd.read_csv(above_dir_preproc+\"Xy_train_meta_v66.csv\")\n    Xy_train_feats = pd.read_csv(above_dir_preproc+\"Xy_train_feats_v66.csv\")\n    # Assign the current clusters instead of ones in the file\n    Xy_train_meta[\"clust_id\"] = kmeans.predict(np.array(Xy_train_meta[HBA_probs]))\n    Xy_train_feats[\"clust_id\"] = Xy_train_meta[\"clust_id\"]\n    # Set other saved values:\n    # These are the same: HBA_number, HBA_votes\n    SMOOTH_WIDTH = 5 \nelse:\n    # Get the meta frame and its features\n    # Use the cluster selection, downsampled\n    Xy_train_meta = (train_meta[train_rows_bool])[::TRAIN_DOWNSEL].copy()\n    Xy_train_meta = Xy_train_meta.reset_index().drop(columns=[\"index\"])\n    print(\"Number of samples used for training =\",len(Xy_train_meta))\n\n    # Get its features\n    Xy_train_feats = assemble_features(Xy_train_meta, traintest=\"train\",\n                            smooth_width=SMOOTH_WIDTH)\n    # Save the dataframes\n    Xy_train_meta.to_csv(\"Xy_train_meta.csv\", header=True,\n                    index=False, float_format='%.6f')\n    Xy_train_feats.to_csv(\"Xy_train_feats.csv\", header=True,\n                    index=False, float_format='%.6f')","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:48.076830Z","iopub.execute_input":"2024-03-19T16:45:48.077189Z","iopub.status.idle":"2024-03-19T16:45:49.531804Z","shell.execute_reply.started":"2024-03-19T16:45:48.077157Z","shell.execute_reply":"2024-03-19T16:45:49.530579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The training features\nXy_train_feats","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:49.533465Z","iopub.execute_input":"2024-03-19T16:45:49.534148Z","iopub.status.idle":"2024-03-19T16:45:49.577139Z","shell.execute_reply.started":"2024-03-19T16:45:49.534103Z","shell.execute_reply":"2024-03-19T16:45:49.575833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Get the Validation Features","metadata":{}},{"cell_type":"code","source":"# Get validation meta and features dataframes\nif USE_PREPROC:\n    # Restore previously saved data frames \"/kaggle/input/\n    Xy_valid_meta = pd.read_csv(above_dir_preproc+\"Xy_valid_meta_v66.csv\")\n    Xy_valid_feats = pd.read_csv(above_dir_preproc+\"Xy_valid_feats_v66.csv\")\n    # Assign the current clusters instead of ones in the file\n    Xy_valid_meta[\"clust_id\"] = kmeans.predict(np.array(Xy_valid_meta[HBA_probs]))\n    Xy_valid_feats[\"clust_id\"] = Xy_valid_meta[\"clust_id\"]\nelse:\n    # Get the meta frame and its features\n    # Use the validation selection, downsampled\n    Xy_valid_meta = (train_meta[valid_rows_bool])[::VALID_DOWNSEL].copy()\n    Xy_valid_meta = Xy_valid_meta.reset_index().drop(columns=[\"index\"])\n    print(\"Number of samples used for Validation =\",len(Xy_valid_meta))\n\n    # Get its features\n    Xy_valid_feats = assemble_features(Xy_valid_meta, traintest=\"validation\",\n                            smooth_width=SMOOTH_WIDTH)\n    # Save the dataframes\n    Xy_valid_meta.to_csv(\"Xy_valid_meta.csv\", header=True,\n                    index=False, float_format='%.6f')\n    Xy_valid_feats.to_csv(\"Xy_valid_feats.csv\", header=True,\n                    index=False, float_format='%.6f')","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:49.578664Z","iopub.execute_input":"2024-03-19T16:45:49.579053Z","iopub.status.idle":"2024-03-19T16:45:50.304701Z","shell.execute_reply.started":"2024-03-19T16:45:49.579019Z","shell.execute_reply":"2024-03-19T16:45:50.303731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The 191 features:  \n#  middle: 84=4x21(0-83), flare: 96=4x4x6(84-179), mean,clust: 11=10+1 (180-190) \nprint(Xy_valid_feats.shape)\nXy_valid_feats.iloc[-5: , 80:90]\n##Xy_valid_feats.iloc[-5: , 175:]","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:50.310097Z","iopub.execute_input":"2024-03-19T16:45:50.313438Z","iopub.status.idle":"2024-03-19T16:45:50.336173Z","shell.execute_reply.started":"2024-03-19T16:45:50.313381Z","shell.execute_reply":"2024-03-19T16:45:50.334983Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Combine train and validation to train all?","metadata":{}},{"cell_type":"code","source":"# Use all the samples for training?\nif TRAIN_ALL:\n    print(\"Using all {} features for training.\".format(\n        len(Xy_train_feats) + len(Xy_valid_feats)))\n    Xy_train_feats = pd.concat([Xy_train_feats, Xy_valid_feats])\n    Xy_train_meta = pd.concat([Xy_train_meta, Xy_valid_meta])","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:50.337950Z","iopub.execute_input":"2024-03-19T16:45:50.339131Z","iopub.status.idle":"2024-03-19T16:45:50.384813Z","shell.execute_reply.started":"2024-03-19T16:45:50.339091Z","shell.execute_reply":"2024-03-19T16:45:50.383725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n\n## Add Logistic Regression proba values to the Features\n    \nUse the proba outputs of LR models to provide additional features. This creates linear combinations of the features that are correlated with the targets which may be useful, or more efficient, for the following random forest model. ","metadata":{}},{"cell_type":"code","source":"# USE_LR1, USE_LR2 are set at the very beginning\n#\n# LogisticRegression(penalty='l2', *, dual=False, tol=0.0001, C=1.0, fit_intercept=True,\n#              intercept_scaling=1, class_weight=None, random_state=None, solver='lbfgs',\n#              max_iter=100, multi_class='auto', verbose=0, warm_start=False,\n#              n_jobs=None, l1_ratio=None)\n\n# Train the LR model\n\nX = Xy_train_feats.drop(columns=[\"clust_id\"])\ny = Xy_train_feats.clust_id\n\n# Don't include the last 10 columns - the mean/medians\nXlr = X.drop(columns=X.columns[-10:])\n\nif USE_LR1:\n    # ***** Use the \"middle8s\" features.\n    Xlr1 = Xlr.iloc[:, 0:83+1]\n    lrmodel1 = LogisticRegression(penalty=LR_REGU, C=LR1_C+1.0, solver='saga', max_iter=1500,\n                                multi_class='multinomial', # ovr or multinomial\n                                n_jobs=-1).fit(Xlr1,y)\n    # Cute figure of the coefficients. Means and medians are on the right side.\n    plt.figure(figsize=(8,4))\n    plt.plot(lrmodel1.coef_.T,'-',alpha=0.5)\n    plt.plot(lrmodel1.coef_.T,'.',alpha=1.0)\n    plt.title(\"LR1 coefficients, C=\"+str(LR1_C)+\" (colored by cluster)\")\n    plt.show()\n    \n    # Fraction of non-zero coef.s (over all clusters):\n    thecoefs = lrmodel1.coef_.T.flatten()\n    print(\"Non-zero coef.s fraction: {:.3f}\".format(sum(thecoefs != 0)/len(thecoefs)))\n    print(\"\\nLR1 model score for X,y = {:.1f}%\\n\".format(100*lrmodel1.score(Xlr1, y)))\n\nif USE_LR2:\n    # ***** Use the Flare fft features.\n    Xlr2 = Xlr.iloc[:, 84:179+1]\n    lrmodel2 = LogisticRegression(penalty=LR_REGU, C=LR2_C+1.0, solver='saga', max_iter=1500,\n                                multi_class='multinomial', # ovr or multinomial\n                                n_jobs=-1).fit(Xlr2,y)\n    # Cute figure of the coefficients. Means and medians are on the right side.\n    plt.figure(figsize=(8,4))\n    plt.plot(lrmodel2.coef_.T,'-',alpha=0.5)\n    plt.plot(lrmodel2.coef_.T,'.',alpha=1.0)\n    plt.title(\"LR2 coefficients, C=\"+str(LR2_C)+\" (colored by cluster)\")\n    plt.show()\n    \n    # Fraction of non-zero coef.s:\n    thecoefs = lrmodel2.coef_.T.flatten()\n    print(\"Non-zero coef.s fraction: {:.3f}\".format(sum(thecoefs != 0)/len(thecoefs)))\n    print(\"\\nLR2 model score for X,y = {:.1f}%\\n\".format(100*lrmodel2.score(Xlr2, y)))\n    ","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:45:50.386240Z","iopub.execute_input":"2024-03-19T16:45:50.386852Z","iopub.status.idle":"2024-03-19T16:51:50.788330Z","shell.execute_reply.started":"2024-03-19T16:45:50.386815Z","shell.execute_reply":"2024-03-19T16:51:50.786451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Include a blur so they don't overpower the other features\n# LR_BLUR set at beginning\n\n# Train: Add these values as new columns in Xy_train_wLRfeats\nXy_train_wLRfeats = Xy_train_feats.copy()\nX = Xy_train_feats.drop(columns=[\"clust_id\"])\nXlr = X.drop(columns=X.columns[-10:])\nif USE_LR1:\n    Xlr1 = Xlr.iloc[:, 0:83+1]\n    lrprobas = lrmodel1.predict_proba(Xlr1)\n    for iadd in range(NUM_CLUSTS):\n        Xy_train_wLRfeats[\"lrMid\"+str(iadd)] = (lrprobas[: , iadd] +\n                            LR_BLUR*(np.random.rand(len(lrprobas))-0.5))\n    # Show hist of max(proba) values\n    plt.figure(figsize=(7,2.5))\n    plt.hist(np.max(lrmodel1.predict_proba(Xlr1),axis=1),bins=50)\n    plt.xlim(-0.5*LR_BLUR, 1.01)\n    plt.title(\"Training max LR1 proba values (pre-blur)\")\n    plt.show()\nif USE_LR2:\n    Xlr2 = Xlr.iloc[:, 84:179+1]\n    lrprobas = lrmodel2.predict_proba(Xlr2)\n    for iadd in range(NUM_CLUSTS):\n        Xy_train_wLRfeats[\"lrFFT\"+str(iadd)] = (lrprobas[: , iadd] +\n                            LR_BLUR*(np.random.rand(len(lrprobas))-0.5))\n    # Show hists of max(proba) values\n    plt.figure(figsize=(7,2.5))\n    plt.hist(np.max(lrmodel2.predict_proba(Xlr2),axis=1),bins=50)\n    plt.xlim(-0.5*LR_BLUR, 1.01)\n    plt.title(\"Training max LR2 proba values (pre-blur)\")\n    plt.show()\n\n\n# Validation: Add these values as new columns in Xy_valid_wLRfeats\nXy_valid_wLRfeats = Xy_valid_feats.copy()\nX = Xy_valid_feats.drop(columns=[\"clust_id\"])\nXlr = X.drop(columns=X.columns[-10:])\nif USE_LR1:\n    Xlr1 = Xlr.iloc[:, 0:83+1]\n    lrprobas = lrmodel1.predict_proba(Xlr1)\n    for iadd in range(NUM_CLUSTS):\n        Xy_valid_wLRfeats[\"lrMid\"+str(iadd)] = (lrprobas[: , iadd] +\n                            LR_BLUR*(np.random.rand(len(lrprobas))-0.5))\n    # Show hists of max(proba) values\n    plt.figure(figsize=(7,2.5))\n    plt.hist(np.max(lrmodel1.predict_proba(Xlr1),axis=1),bins=50)\n    plt.xlim(-0.5*LR_BLUR, 1.01)\n    plt.title(\"Validation max LR1 proba values (pre-blur)\")\n    plt.show()\nif USE_LR2:\n    Xlr2 = Xlr.iloc[:, 84:179+1]\n    lrprobas = lrmodel2.predict_proba(Xlr2)\n    for iadd in range(NUM_CLUSTS):\n        Xy_valid_wLRfeats[\"lrFFT\"+str(iadd)] = (lrprobas[: , iadd] +\n                            LR_BLUR*(np.random.rand(len(lrprobas))-0.5))\n    # Show hists of max(proba) values\n    plt.figure(figsize=(7,2.5))\n    plt.hist(np.max(lrmodel2.predict_proba(Xlr2),axis=1),bins=50)\n    plt.xlim(-0.5*LR_BLUR, 1.01)\n    plt.title(\"Validation max LR2 proba values (pre-blur)\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:51:50.791285Z","iopub.execute_input":"2024-03-19T16:51:50.791967Z","iopub.status.idle":"2024-03-19T16:51:52.880651Z","shell.execute_reply.started":"2024-03-19T16:51:50.791903Z","shell.execute_reply":"2024-03-19T16:51:52.879369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n\n## Train a Random Forest model","metadata":{}},{"cell_type":"code","source":"# RandomForestClassifier(n_estimators=100, *, criterion='gini', max_depth=None,\n#    min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0,\n#    max_features='sqrt', # fraction of feat.s to consider when looking for best split\n#                         # Note: the search for a split does not stop until\n#                         # at least one valid partition of the node samples is found,\n#                         # even if it requires inspecting more than max_features features.\n#    max_leaf_nodes=None, min_impurity_decrease=0.0,\n#    bootstrap=True, oob_score=False, n_jobs=None, random_state=None, verbose=0,\n#    warm_start=False, class_weight=None, ccp_alpha=0.0,\n#    max_samples=None, # fraction of samples from X used to train each estimator\n#    monotonic_cst=None)\n\nX = Xy_train_wLRfeats.drop(columns=[\"clust_id\"])\ny = Xy_train_wLRfeats.clust_id\n\nave_oob = []\n##nfits=3    # Do it more than once and average?  No :)\n##for ifit in range(nfits):\nrfmodel = RandomForestClassifier(\n                    n_estimators=200,\n                    min_samples_leaf=3,\n                    max_features=0.2,\n                    max_samples=0.85,\n                    oob_score=True, class_weight=\"balanced_subsample\",\n                    n_jobs=-1, verbose=0, random_state=None).fit(X,y)\nave_oob.append(rfmodel.oob_score_)\n    \n# Show the feature importances\nsort_inds = (rfmodel.feature_importances_.argsort())\nplt.figure(figsize=(4,7))\nplt.barh(rfmodel.feature_names_in_[sort_inds], rfmodel.feature_importances_[sort_inds])\nplt.ylim(len(sort_inds)-40,len(sort_inds)+0.2)\nplt.title(\"Feature Importances (top 40)\")\nplt.show()\n\n##print(\"\\nRF model ave OOB score = {:.1f}% +/- {:.1f}\".format(\n##        100*np.mean(ave_oob),100*np.std(ave_oob)))\nprint(\"\\nRF model ave OOB score = {:.1f}%\".format(100*np.mean(ave_oob)))\n\nprint(\"\\nRF model score for X,y = {:.1f}%\\n\".format(100*rfmodel.score(X, y)))","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:51:52.882604Z","iopub.execute_input":"2024-03-19T16:51:52.883485Z","iopub.status.idle":"2024-03-19T16:55:20.577769Z","shell.execute_reply.started":"2024-03-19T16:51:52.883434Z","shell.execute_reply":"2024-03-19T16:55:20.576601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate the KL score for the train_meta\n\n# Make and add the predicted clusters to Xy_train_meta\nXy_train_meta[\"pred_id\"] = rfmodel.predict(Xy_train_wLRfeats.drop(columns=[\"clust_id\"]))\n\n# How high is the probability of the assigned (maximum prob) class for each X?\nmaxprobs_train = np.max(rfmodel.predict_proba(X),axis=1)\nplt.figure(figsize=(7,2.5))\nplt.hist(maxprobs_train, bins=50)\nplt.xlim(0.,1.)\nplt.title(\"Histogram of max(proba) for model-training samples\")\nplt.show()\n\n# Create a 'solution' of actual probabilities from the Xy_train_meta\nsolution = Xy_train_meta[[\"eeg_id\"] + HBA_votes]\n# In solution replace the _votes values with the probs from meta\nfor col_pre in HBA_names:\n    solution.loc[:, col_pre + \"_vote\"] = Xy_train_meta[col_pre + \"_prob\"]\n\n# Start a submission dataframe using the solution format\nsubmission = solution.copy()\n\nif USE_TAMED:\n    # Create a 'submission' from the predicted cluster ids via clust_centers\n    pred_ids = Xy_train_meta[\"pred_id\"] # putting clust_id gives the ideal value\n    # Calculate the KL score, first doing \"taming\" adjustments on centers\n    best_fracs, best_centers = find_best_tamed_kl()\n    # Using the best_centers\n    # Put the probs into the submission column by column\n    for iprob in range(HBA_number):\n        this_col_probs = best_centers[: , iprob]\n        submission[HBA_votes[iprob]] = this_col_probs[pred_ids]\n    print(\"Tamed fractions:\\n\", best_fracs, \"\\nand centers:\\n\", best_centers)\n\nelse:\n    # Use the proba values to weight the cluster centers as prediced probs\n    probas = rfmodel.predict_proba(Xy_train_wLRfeats.drop(columns=[\"clust_id\"]))\n    pred_probs = probas @ clust_centers  # matrix multiply\n    # Put the probs into the submission column by column\n    for iprob in range(HBA_number):\n        submission[HBA_votes[iprob]] = pred_probs[: , iprob]\n    print(\"Predicted probabilites are proba-weighted cluster centers.\")\n    \n# Evaluate the KL divergence\nthis_kl = kld_score(solution, submission)\nprint(\"\\nKL from predicted centers: {:.4f}\".format(this_kl))","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:20.579536Z","iopub.execute_input":"2024-03-19T16:55:20.580464Z","iopub.status.idle":"2024-03-19T16:55:25.248317Z","shell.execute_reply.started":"2024-03-19T16:55:20.580429Z","shell.execute_reply":"2024-03-19T16:55:25.246842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at where the errors are:\npd.crosstab(Xy_train_meta[\"pred_id\"], Xy_train_meta[\"clust_id\"])","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:25.250665Z","iopub.execute_input":"2024-03-19T16:55:25.251649Z","iopub.status.idle":"2024-03-19T16:55:25.317316Z","shell.execute_reply.started":"2024-03-19T16:55:25.251589Z","shell.execute_reply":"2024-03-19T16:55:25.315821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Apply Model to Validation Data","metadata":{}},{"cell_type":"code","source":"# Calculate the KL score for the valid_meta\n\n# Make and add the predictions to Xy_valid_meta\nXy_valid_meta[\"pred_id\"] = rfmodel.predict(Xy_valid_wLRfeats.drop(columns=[\"clust_id\"]))\n# For reference, randomly assign pred_ids to see zero-accuracy result\n##Xy_valid_meta[\"pred_id\"] = np.random.choice(num_clusts,\n##                        size=len(Xy_valid_meta),replace=True, p=None)\n\n# How high is the probability of the assigned (maximum prob) class for each X?\nmaxprobs_valid = np.max(rfmodel.predict_proba(\n                    Xy_valid_wLRfeats.drop(columns=[\"clust_id\"])),axis=1)\nplt.figure(figsize=(7,2.5))\nplt.hist(maxprobs_valid, bins=50)\nplt.xlim(0.,1.)\nplt.title(\"Histogram of max(proba) for validation samples\")\nplt.show()\n\n# Create a 'solution' of actual probabilities from the Xy_valid_meta\nsolution = Xy_valid_meta[[\"eeg_id\"] + HBA_votes]\n# In solution replace the _votes values with probabilities\nfor col_pre in HBA_names:\n    solution.loc[:, col_pre + \"_vote\"] = Xy_valid_meta[col_pre + \"_prob\"]\n\n# Start a submission dataframe using the solution format\nsubmission = solution.copy()\n\nif USE_TAMED:\n    # Create a 'submission' from the predicted cluster ids via clust_centers\n    pred_ids = Xy_valid_meta[\"pred_id\"]\n    best_fracs, best_centers = find_best_tamed_kl()\n    # Using the best_centers\n    # Put the probs into the submission column by column\n    for iprob in range(HBA_number):\n        this_col_probs = best_centers[: , iprob]\n        submission[HBA_votes[iprob]] = this_col_probs[pred_ids]\n    print(\"Tamed fractions:\\n\", best_fracs, \"\\nand centers:\\n\", best_centers)\n\nelse:\n    # Use the proba values to weight the cluster centers as prediced probs\n    probas = rfmodel.predict_proba(Xy_valid_wLRfeats.drop(columns=[\"clust_id\"]))\n    pred_probs = probas @ clust_centers  # matrix multiply\n    # Put the probs into the submission column by column\n    for iprob in range(HBA_number):\n        submission[HBA_votes[iprob]] = pred_probs[: , iprob]\n    print(\"Predicted probabilites are proba-weighted cluster centers.\")\n    \n# Evaluate the KL divergence\nthis_kl = kld_score(solution, submission)\nprint(\"\\nKL from predicted centers: {:.4f}\".format(this_kl))","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:25.319779Z","iopub.execute_input":"2024-03-19T16:55:25.320738Z","iopub.status.idle":"2024-03-19T16:55:27.230894Z","shell.execute_reply.started":"2024-03-19T16:55:25.320666Z","shell.execute_reply":"2024-03-19T16:55:27.229418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at where the errors are:\npd.crosstab(Xy_valid_meta[\"pred_id\"], Xy_valid_meta[\"clust_id\"])","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:27.233048Z","iopub.execute_input":"2024-03-19T16:55:27.233972Z","iopub.status.idle":"2024-03-19T16:55:27.285174Z","shell.execute_reply.started":"2024-03-19T16:55:27.233923Z","shell.execute_reply":"2024-03-19T16:55:27.283909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Using v47 saved features: Smooth=5, Ntrain=27010, Nvalid=9582\n#\n#                                                         -----  Random Forest  -----\n#  Valid   Train           Clusts  LR1_C  LR2_C  LR_BLUR  n_est  leaf>  feats%  samp%\n# 0.9291  0.6056            10     0.9    N/A     0.10    100     9      0.20   0.80\n# 0.9331  0.6032             \"\n# 0.9307  0.6039             9         Trying different number of clusters\n# 0.9312  0.6125             \"\n# 0.9271  0.6063             8\n# 0.9292  0.6089             \"\n# 0.9292  0.6161             7\n# 0.9224  0.6116             \"\n# 0.9281  0.6162  v49(0.90)  \"\n# 0.9263  0.6204             6\n# 0.9267  0.6201             \"\n# 0.9304  0.6104             8     3.00         Trying different LR C values\n# 0.9307  0.6090             \"     0.25\n# 0.9227  0.6027             \"     0.01\n# 0.9023  0.5917             \"     0.002  only lr0 in top 40 (at #1)\n# 0.9017  0.5881  v50(0.94)  \"     0.0005 all coef.s = 0\n# 0.9240  0.6013  v52(0.92)  9*    0.1   *All rows eeg_id<200 used for clusters\n#  Valid   Train           Clusts  LR1_C  LR2_C  LR_BLUR  n_est  leaf>  feats%  samp%\n# 0.9269  0.6027             9*    0.1    0.1     0.10    100     9      0.20   0.80\n# 0.8901  0.3664  v53(0.92)  9*    0.1    0.1     0.10    300     5      0.20   0.90\n# 0.8418  0.3621  v54(0.95)  9*       No LR feat.s        300     5      0.20   0.90\n# 0.9011  0.3682  v55(0.90)  9*   1.0,L2 1.0,L2   0.10    300     5      0.20   0.90\n# 0.8889  0.3670  v56(0.93)  9*  LRs split,C=1,L2 0.10    300     5      0.20   0.90\n\n# Using v62 saved features: Smooth=5, Ntrain=27010, Nvalid=15960\n#\n#  Valid   Train           Clusts         LR1_C LR2_C  BLUR   n_est  leaf>  feats%  samp%\n# 0.8676  0.3479  v62(-.--)  9   LRs(84,96)L2:1.0,1.0  0.10    300     5     0.20   0.90\n# 0.8663  0.3416  v63(0.90)  9\n# 0.8769  0.3938  v64(0.90)  6\n# 0.8645  0.3423             9   LRs(84,96)L1:.15,.10  0.10    300     5     0.20   0.90\n# 0.8451  0.3388  v65(-.--)  9                .05 .03\n\n# Using v66 saved features: Smooth=5, Ntrain=33152, Nvalid=9582\n#\n#  Valid   Train           Clusts         LR1_C LR2_C  BLUR   n_est  leaf>  feats%  samp%\n# 0.7792  0.3437  v67(0.96)  9   LRs(84,96)L1:.05,.03  0.10    300     5     0.20   0.90\n# 0.8200  0.4576  v68(0.96)  9   LRs(84,96)L1:.05,.03  0.10    300     7     0.30   0.80\n# 0.7866  0.6201  v69(0.72)  Same as v68 but using weighted centers (instead of 'tamed')\n# 0.7559  0.6072  v70(0.72)  v69 without LRs  \n# 0.7252  0.5346  v71(0.72)  v69 without LRs and RF params back to: 300, 5, 0.20, 0.90\n#                 v72(0.71)  v69 (with LRs) RF params back to: 300, 5, 0.20, 0.90\n# 0.7444  0.5104                 (with LRs) RF params: 200, 4, 0.20, 0.85 even min leaf is not good\n# 0.7280  0.4601  v73(0.71)      (with LRs) RF params: 200, 3, 0.20, 0.85\n# 0.5165  0.4390  v74(0.??)  v73 using combined data to train (validation KL not meaningful)","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:27.287682Z","iopub.execute_input":"2024-03-19T16:55:27.288692Z","iopub.status.idle":"2024-03-19T16:55:27.299011Z","shell.execute_reply.started":"2024-03-19T16:55:27.288643Z","shell.execute_reply":"2024-03-19T16:55:27.297329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Scanning RF hyper parameters 'by hand'","metadata":{}},{"cell_type":"code","source":"# Do training fit and validation KL in a loop to chek hyper params.\nif False:\n    # Training X,y\n    X = Xy_train_wLRfeats.drop(columns=[\"clust_id\"])\n    y = Xy_train_wLRfeats.clust_id\n    # Create a 'solution' of actual probabilities from the Xy_valid_meta\n    solution = Xy_valid_meta[[\"eeg_id\"] + HBA_votes]\n    # In solution replace the _votes values with probabilities\n    for col_pre in HBA_names:\n        solution.loc[:, col_pre + \"_vote\"] = Xy_valid_meta[col_pre + \"_prob\"]\n    # Create a 'submission' template\n    submission = solution.copy()\n\nif False:  # if statement in place of the for loop\n##for scan_this in [60,60, 100,100, 200,200, 350,350, 500,500]:\n    # min leaf [3,5,7,9,11,13,16,19, 23, 28, 35]\n    # max feats [0.1,0.1,0.2,0.2,0.35,0.35,0.5,0.5,0.7,0.7]\n    # max samples [0.35,0.35,0.5,0.5,0.7,0.7,0.8,0.8,0.9,0.9]\n    # n_est.s [40,40, 60,60, 100,100, 200,200, 350,350, 500,500]\n    rfmodscan = RandomForestClassifier(\n                    n_estimators=200,  #300,\n                    min_samples_leaf=3,    #5  odd is better than even\n                    max_features=0.2,\n                    max_samples=0.85, #0.9,\n                    oob_score=True, class_weight=\"balanced_subsample\",\n                    n_jobs=-1, verbose=0, random_state=None).fit(X,y)\n    print(\"   _ _ _ scan value =\",str(scan_this),\":\")\n    print(\"RF train score = {:.1f}%\".format(100*rfmodscan.score(X, y)))\n    # Make and add the predictions to Xy_valid_meta\n    Xy_valid_meta[\"pred_id\"] = rfmodscan.predict(Xy_valid_wLRfeats.drop(\n                                    columns=[\"clust_id\"]))\n    \n    # Use the proba values to weight the cluster centers as prediced probs\n    probas = rfmodscan.predict_proba(Xy_valid_wLRfeats.drop(columns=[\"clust_id\"]))\n    pred_probs = probas @ clust_centers  # matrix multiply\n    # Put the probs into the submission column by column\n    for iprob in range(HBA_number):\n        submission[HBA_votes[iprob]] = pred_probs[: , iprob]\n    print(\"Predicted probabilites are proba-weighted cluster centers.\")\n    \n    # Evaluate the KL divergence\n    this_kl = kld_score(solution, submission)\n    print(\"KL from predicted centers: {:.4f}\\n\".format(this_kl))","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:27.301695Z","iopub.execute_input":"2024-03-19T16:55:27.303180Z","iopub.status.idle":"2024-03-19T16:55:27.324102Z","shell.execute_reply.started":"2024-03-19T16:55:27.303133Z","shell.execute_reply":"2024-03-19T16:55:27.322495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>","metadata":{}},{"cell_type":"markdown","source":"## Create and Output the Test Predictions","metadata":{"papermill":{"duration":0.031892,"end_time":"2024-02-01T20:26:05.168059","exception":false,"start_time":"2024-02-01T20:26:05.136167","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Assemble the submission with probabilites assigned by cluster id\n\n# Generate the test features:\nXy_test_feats = assemble_features(test_meta, traintest=\"test\",\n                            smooth_width=SMOOTH_WIDTH)\n\n# Add LR probas \n# Xy_test_feats has no clust_id, don't need to drop it\nXy_test_wLRfeats = Xy_test_feats.copy()\nXlr = Xy_test_feats.drop(columns=Xy_test_feats.columns[-10:])\n#\nif USE_LR1:\n    Xlr1 = Xlr.iloc[:, 0:83+1]\n    lrprobas = lrmodel1.predict_proba(Xlr1)\n    for iadd in range(NUM_CLUSTS):\n        Xy_test_wLRfeats[\"lrMid\"+str(iadd)] = (lrprobas[: , iadd] +\n                            LR_BLUR*(np.random.rand(len(lrprobas))-0.5))\nif USE_LR2:\n    Xlr2 = Xlr.iloc[:, 84:179+1]\n    lrprobas = lrmodel2.predict_proba(Xlr2)\n    for iadd in range(NUM_CLUSTS):\n        Xy_test_wLRfeats[\"lrFFT\"+str(iadd)] = (lrprobas[: , iadd] +\n                            LR_BLUR*(np.random.rand(len(lrprobas))-0.5))\n\n# Start submission with a dataframe of the eeg_id column from test_meta\ntest_submit = test_meta[[\"eeg_id\"]].copy()\n# Add \"_vote\" columns, initially all equal.\nfor new_col in HBA_votes:\n    test_submit[new_col] = 1/HBA_number\n\nif USE_TAMED:\n    # Use the RF model to assign cluster id predictions\n    pred_ids = rfmodel.predict(Xy_test_wLRfeats)\n    # Put the best (validation) tamed cluster centers/probs into the submission\n    for iprob in range(HBA_number):\n        this_col_probs = best_centers[: , iprob]\n        test_submit[HBA_votes[iprob]] = this_col_probs[pred_ids]    \nelse:\n    # Use the proba values to weight the cluster centers as prediced probs\n    probas = rfmodel.predict_proba(Xy_test_wLRfeats)\n    pred_probs = probas @ clust_centers  # matrix multiply\n    # Put the probs into the submission column by column\n    for iprob in range(HBA_number):\n        test_submit[HBA_votes[iprob]] = pred_probs[: , iprob]\n\nprint(test_submit)\n\n# Output the file\ntest_submit.to_csv(\"submission.csv\", header=True, \n                        index=False, na_rep='', float_format='%.6f')\n# Look at the file\n##!more submission.csv","metadata":{"papermill":{"duration":0.060404,"end_time":"2024-02-01T20:26:05.259614","exception":false,"start_time":"2024-02-01T20:26:05.199210","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-19T16:55:27.326602Z","iopub.execute_input":"2024-03-19T16:55:27.327417Z","iopub.status.idle":"2024-03-19T16:55:27.560042Z","shell.execute_reply.started":"2024-03-19T16:55:27.327368Z","shell.execute_reply":"2024-03-19T16:55:27.558991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<HR>\n### The End","metadata":{"papermill":{"duration":0.030063,"end_time":"2024-02-01T20:26:05.320357","exception":false,"start_time":"2024-02-01T20:26:05.290294","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Version history\n#       LB /\"CV\",Train\n#  v6 2.42 /3.00  Worst case: assign clusters randomly to the samples.\n# v17     Processing scheme for the middle 4s of the spectrograms for a sub id:\n#           average the two middle 2s spectra, divide by a reference spectrum,\n#           calculate the mean and median, and divide by the mean.\n#           smooth by n with a rolling window.\n#         Created one large plot with overplotting of ~450 x 4 of these spectra; \n#         color coding by cluster (for the 6 basic clusters) shows some variation.\n# v18     Using the spectrogram processing, create features from the spectrograms\n#         and use them to fit cluster id with a Random Forest model.\n# v20 1.74 /---- Training overfitting gives 0.37 on X sized 354, clust=8, smooth=7.\n# v21 1.72 /---- Add mean/median for each of the \"LL\",\"RL\",\"LP\",\"RP\" spectra.\n# v23 1.56 /1.58 clusts=6, smooth=7, leaf_nodes=4*clusts, features=0.5, samples=0.7\n# v24     Do the log10() before the smoothing and downsample of spectra.\n#     1.64 /1.74,1.03\n# v26     \"Tamed\" the cluster center loc.s/prob.s for best validation KL --\n#         linear mix between center values and all 1/6, set by tamefac = 0.5\n#     1.01 /1.18(tamed),1.14(not tamed)\n# v27     Scan for the best tamed fraction (same for all clusters) - expect ~ v26\n#     0.98 /1.18,0.88  <-- all tamed, clusts=6 width=7\n# v28 0.99 /1.19,0.79  <-- all tamed, clusts=*8* width=7\n# v29 0.95 /1.17,0.84  <-- all tamed, clusts=6 width=*5*\n# v30 0.99 /1.16,0.82  <-- all tamed, clusts=6 width=*3*\n# v31 0.98 /1.16,0.68  <-- all tamed, clusts=*10* width=*5*\n# v32 ---- /1.17,0.83  <-- all tamed, clusts=6 width=*5*   Added the lowest freqs too.\n# v33 0.96 /1.14,0.82  <-- indiv tamed, clusts=6 width=*5* \n# v34 0.95 /1.18,0.79  <-- indiv tamed, clusts=6 width=*5*  * Added LR features, blur=0.7*\n# v35 0.96 /1.21,0.71  <-- indiv tamed, clusts=*9* width=*5*   with LR features, blur=0.7*\n# v36 0.94 /1.13,1.05  <--   \" \", clusts=*9* width=*5* , LR features(0.7), x5 samples\n# v37 0.95 /1.14,1.06  <--   \" \", clusts=*9* width=*5* , LR-no means,l1,0.42, x5 samples\n# v38     Adjust the RForest hyper parameters roughly, got e.g. min leaf = 3, ...\n#     0.95 /1.04,0.29  <--   \" \", clusts=*9* width=*5* , LR-no means,l1,0.42, x5 samples\n# v39 0.91 /1.00,0.65  <-- clusts=9 width=5, LR,l1,0.20, min_leaf=9, 16k samples\n# v40 0.93 /1.02,0.68  <-- clusts=*6* width=5, LR,l1,0.20, min_leaf=9, 16k samples\n# v41 0.93 /1.00,0.65  <-- clusts=*10* width=5, LR,l1,0.20, min_leaf=9, 16k samples\n# v42 0.94 /0.99,0.65  <-- clusts=9 width=5, No LR feats, min_leaf=9, 16k samples\n# v43     Working on adding more features from spectrogram...\n#         Using the ratio of the middle 2s spectra: KL from tamed centers: 1.3643  :(\n# v44     Try using ONLY the center 8s ratio'ed to ave of 48,80,112 s before/after spectra.\n#     1.08 /1.15,0.72  Worse than the 1.07 /1.38 using average probs.\n# v46 0.91 /0.97,0.60  Like v39 except LR_BLUR=0.1, all features included.\n#                      No ratio feature in the top 25 feats.\n# v47     Split the spectro colored plots into 4: one for each \"LL\",\"RL\",\"LP\",\"RP\":\n#         they all look about the same, perhaps not surprisingly...\n#         Split the plots into 9 by the clusters - shows spectral differences better.\n#         Use 3/4 of the eeg_id=0 plus some non-zero ones for training,\n#         other 1/4 for validation plus some non-zero ones.\n#         The Hybrid clusters look different esp the LRDA-Other. [Changed in v52]\n#     0.89 /0.93,0.60  Used all defined training and validations. (~1 hour+ run.)\n# v48     A few little changes. Put the v47 output feature dataframes into a dataset.\n#         Read those files in to avoid the feature extraction.\n#     0.92 /0.93,0.61  ~ Same as v47 except no lr2 (ratio) features (not used anyway.)\n# v49     More cleanup and tried different cluster sizes, not much difference, use 7.\n#     0.90 /0.92,0.61  ~ Same as v47 with 7 clusters.\n# v50 0.94 /0.90,0.58  ~ Same as v47 with 8 clust.s, LR1_C=0.0005(basically off)\n# v52 0.92 /0.92,0.60  ~ v47 with 9 clusts(from all eeg_id<200), LR1_C=0.1\n# v53 0.92 /0.89,0.37  ~ v47 with 9 clusts(eeg<200), LR1,LR2=0.1, RF: 300,5,0.2,0.9\n# v54 0.95 /0.84,0.36  ~ v47 with 9 clusts(eeg<200),  No LRs,     RF: 300,5,0.2,0.9\n# v55 0.90 /0.90,0.37  ~ v47 with 9 clusts(eeg<200), LRCs=1.0,*L2*, RF: 300,5,0.2,0.9\n# v56     Add code to look at the fft of the spectral amplitude, \"Flare-Rate spectra\"\n#         toward adding features.  Split the LRs onto the halfs of the normal feat.s. :(\n#     0.93 /0.90,0.37  ~ v47 9 clusts, LRs split,LRCs=1.0,*L2*, RF: 300,5,0.2,0.9\n# v57,58,59  Continue adding the ffts and plotting them and spectra together.\n#     60,61  \"\n# v62     Added fft features to the middle 8 second features and the means/medians.\n#     -.-- /0.87,0.35  Run (5300s) using all the train, valid feat.s; the csv files are v62.\n# v63 0.?? /0.87,0.35  Same as v62 but the features are taken from the dataset.\n# v64 0.90 /0.88,0.35  Same as v62 with 6 clusters.\n# v65 0.90 /0.85,0.34  Same as v62, LRs L1 reg with Cs = 0.05 and 0.03.\n# v66 -.-- /0.78,0.34  Some fft feature adjustments and some more train/valid samples.\n# v67 0.96 /0.78,0.34  Use v66 features files. Yeesh, validation is better but LB worse.\n# v68 0.96 /0.82,0.46  For variety, change RF params to: 300, 7, 0.30, 0.80 .\n# v69                  In the words of Lord Peter Wimsey: \n#                                  \"Oh Bunter! I've been such an ass!\"\n#                      The \"taming\" step requires knowing the true solution, and if the\n#                      validation and test samples were similar then using the validation\n#                      best_centers for test might be a reasonable thing to do.\n#                      But I should really use the proba values as weights to sum the\n#                      cluster centers. Implemented that with a simple matrix multiply :) \n#     0.72 /0.78,0.61  Nice!  Interestingly, the train KL increased but the valid went down some.\n# v70 0.72 /0.76,0.61  v69 without LRs  \n# v71 0.72 /0.73,0.53  v69 without LRs and RF params back to: 300, 5, 0.20, 0.90\n# v72 0.71 /0.76,0.55  v69 (with LRs) RF params: 300, 5, 0.20, 0.90\n# v73 0.71 /0.73,0.46  v69 (with LRs) RF params: 200, 3, 0.20, 0.85\n# v74 0.?? /----,0.44  v73 trained on combined data\n ","metadata":{"execution":{"iopub.status.busy":"2024-03-19T16:55:27.562049Z","iopub.execute_input":"2024-03-19T16:55:27.562505Z","iopub.status.idle":"2024-03-19T16:55:27.573637Z","shell.execute_reply.started":"2024-03-19T16:55:27.562461Z","shell.execute_reply":"2024-03-19T16:55:27.572550Z"},"trusted":true},"execution_count":null,"outputs":[]}]}