{"cells":[{"metadata":{"_uuid":"d6a3b5c2b7374ad33a71bda5854e7bd43ded657a"},"cell_type":"markdown","source":"During the competition I got curious about possible ways how to reconstruct representative light curves of hard classes. I was lucky to somehow reconstruct those. Features based on it likely helped me to jump from the 335th place on public LB place to the 182th place on private LB, which seems to be one of the largest jumps in this competition.\n\nIt was very exciting for me to see how “signatures of supernovae” emerge from noisy data! 😊 I’ve even created a kernel to share some of routines I tried out. However, I made it private almost immediately. Because when gratefully reading about all the great solutions by top-scoring teams I just felt humble and didn’t want to distract anyone from all those fantastic solutions. After half a year I’d like to briefly share a couple of things I did, just in case one day it might be of any practical use to someone in the astroinformatics community. \n\nThe main idea was to match and then to robustly aggregate light curves of the same class. For instance, after some rescaling and shifting we match two instances of the class 62 like this:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from IPython.display import Image\nImage(\"../input/plasstic-intro/first_iter_plasstic.png\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And after aggregation we have an updated representative light curve for the class:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"Image(\"../input/plasstic-intro/agg_1.png\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Repeating the match-aggregate, 80 iterations later we have:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"Image(\"../input/plasstic-intro/later_iter_plasstic.png\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"Image(\"../input/plasstic-intro/agg_later.png\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This way I match-aggregated light curves which likely had maximum. Next, I match-aggregated light curves without apparent maximum, where I had to be slightly more careful with rescaling:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"Image(\"../input/plasstic-intro/partial_match.png\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the end, I've got representative light curves for different classes, and for the first time I saw slight differences in the hard classes. E.g., we dealt with the class 62 in the previous figures. In the figure below we can compare representative light curves of two different classes with the mentioned 62-class light curve:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"Image(\"../input/plasstic-intro/comparison.png\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, I computed some features based on the representative light curves. The main idea was to try to match a light curve with every class prototype, looking for the best rescaling (both time and flux scales). The matching was done using either squared error or cosine error (so, either forcing points to be close to each other, or forcing segments to have similar directions). As CPMP advised to write efficient code, I wrote an efficient implementation in Cython.\n\nIf interested, you can find the implementation with an example of usage below in this kernel.\n\nP.S. I forgot to mention that before starting the match-aggregate thing, I combined passbands into a single light curve. The details can be also found below as well.\n\n\nI'd like to thank the organizers and the community for fantastic opportunity to learn more about the exciting project and to see all the cool techniques used for the problem!"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nplt.rcParams.update({'figure.max_open_warning': 0, 'axes.titlesize':22, 'axes.labelsize': 22, 'lines.linewidth' : 2, 'lines.markersize' : 7})\nimport os\nimport pickle\nfrom scipy.stats import sigmaclip\nfrom sklearn.metrics import mean_squared_error\nfrom multiprocessing import Pool, cpu_count\nfrom functools import partial\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"### Preprocessing"},{"metadata":{"_uuid":"1c9b4632ce8c9f75a4d9a869e9047c3562f9d722"},"cell_type":"markdown","source":"Reading data, using class 62 for illustration. Doing initial data aggregation of observations from the same night."},{"metadata":{"trusted":true,"_uuid":"8c7fe9c16962919dd2e78ba93b59172f15be2933"},"cell_type":"code","source":"df = pd.read_csv('../input/PLAsTiCC-2018/training_set.csv')\ndf_meta = pd.read_csv('../input/PLAsTiCC-2018/training_set_metadata.csv')\nclass_to_process = 62\ndf_meta = df_meta[df_meta['target']==class_to_process]\n\n# dealing with two measurements during the same night based on CPMP's kernel\ndf['mjd'] += 0.3\n\n# time scale correction based on https://physics.stackexchange.com/questions/172474/redshift-of-supernova-light-curve, https://physics.stackexchange.com/questions/156618/tired-light-red-shift-hypothesis\ndf = df.merge(df_meta[['object_id', 'hostgal_specz']], on='object_id')\ndf['mjd'] /= df['hostgal_specz'] + 1\n\n# aggregating measurements from the same night\ndf['mjd_int'] = df['mjd'].map(lambda x: np.modf(x)[1])\ndf.sort_values(by=['object_id', 'passband', 'mjd_int'], inplace=True)\ndf[['object_id_lag', 'passband_lag', 'mjd_int_lag', 'flux_lag', 'flux_err_lag']] = df[\n    ['object_id', 'passband', 'mjd_int', 'flux', 'flux_err']].shift(1)\ndf[['object_id_lead', 'passband_lead', 'mjd_int_lead', 'flux_lead', 'flux_err_lead']] = df[\n    ['object_id', 'passband', 'mjd_int', 'flux', 'flux_err']].shift(-1)\nprevious_flux_the_same_night = (df['object_id'] == df['object_id_lag']) & (df['passband'] == df['passband_lag']) & (\n        df['mjd_int'] == df['mjd_int_lag'])\nnext_flux_the_same_night = (df['object_id'] == df['object_id_lead']) & (df['passband'] == df['passband_lead']) & (\n        df['mjd_int'] == df['mjd_int_lead'])\ndf.loc[previous_flux_the_same_night, 'flux'] = (df.loc[previous_flux_the_same_night, 'flux'] + df.loc[\n    previous_flux_the_same_night, 'flux_lag']) / 2\ndf.loc[previous_flux_the_same_night, 'flux_err'] = (df.loc[previous_flux_the_same_night, 'flux_err'] ** 2 + df.loc[\n    previous_flux_the_same_night, 'flux_err_lag'] ** 2).map(np.sqrt)\ndf = df[~next_flux_the_same_night]\ndf.drop(\n    ['object_id_lag', 'passband_lag', 'mjd_int_lag', 'flux_lag', 'flux_err_lag', 'object_id_lead', 'passband_lead',\n     'mjd_int_lead', 'flux_lead', 'flux_err_lead'], axis=1, inplace=True)\n\nfirst_timestamp_per_object = df.groupby('object_id')['mjd'].agg(lambda x: np.modf(np.min(x))[1])\nfirst_timestamp_per_objects = df['object_id'].map(lambda x: first_timestamp_per_object[x])\n\n# timeorder stand for order of the observation night\ndf['timeorder'] = df['mjd_int'] - first_timestamp_per_objects\n\ndf = df[['object_id', 'passband', 'timeorder', 'flux', 'flux_err']]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f48407ee6ad181bf4e84ed0e6a51f7aaa157e991"},"cell_type":"markdown","source":"Let's visualize initial time series data."},{"metadata":{"trusted":true,"_uuid":"0772b8b94d566ee2e54710bbb26e64ebb8a705c5"},"cell_type":"code","source":"def plot_object(object_id, class_id='', with_errors=True):\n    object_df = df[df['object_id'] == object_id]\n    min_timeorder = df.loc[(df['object_id'] == object_id) & ~np.isnan(df['flux']), 'timeorder'].min()\n    max_timeorder = df.loc[(df['object_id'] == object_id) & ~np.isnan(df['flux']), 'timeorder'].max()\n    for passband in range(6):\n        plt.figure(figsize=(30,5))\n        object_df_passband = object_df[object_df['passband'] == passband]\n        if with_errors:\n            plt.errorbar(object_df_passband['timeorder'], object_df_passband['flux'], yerr=object_df_passband['flux_err'])\n        else:            \n            plt.plot(object_df_passband['timeorder'], object_df_passband['flux'])\n        plt.scatter(object_df_passband['timeorder'], object_df_passband['flux'])\n        plt.xlim(min_timeorder, max_timeorder)\n        plt.title(f'object id: {object_id}, passband: {passband}')\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n\nid_to_check = df['object_id'].unique()[33]        \nplot_object(id_to_check)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c8b44b73c181bd12e14ae2460d4d715c71582035"},"cell_type":"markdown","source":"Considering slower flux changes in light curves for the classes with single maximum,  I first aggregated and slightly smoothed the time series. "},{"metadata":{"trusted":true,"_uuid":"44ac42c9dc8bc68cf10a94aef2e83638ecffb6a0"},"cell_type":"code","source":"def moving_mean_smoothed_aggregations(df, column_names=['flux', 'flux_err'], agg_window_width=10,\n                                      exp_decaying_tail_width=7, exp_decaying_const=0.3):\n    # imputing anchor timestamps for future aggregations\n    existing_ts_entries = set(map(tuple,\n                                  df.loc[df['timeorder'] % agg_window_width == 0, ['object_id', 'passband',\n                                                                                   'timeorder']].astype(\n                                      np.int).drop_duplicates().values.tolist()))\n\n    max_timeorder = df['timeorder'].max()\n\n    triples_to_add = np.array([[object_id, passband, timeorder, np.nan, np.nan]\n                               for object_id, passband in df[['object_id', 'passband']].drop_duplicates().values\n                               for timeorder in np.arange(0, int(max_timeorder) + 1, agg_window_width)\n                               if (object_id, passband, timeorder) not in existing_ts_entries])\n    if len(triples_to_add) > 0:\n        df_to_append = pd.DataFrame(triples_to_add, columns=df.columns)\n        df = df.append(df_to_append, ignore_index=True)\n        df.sort_values(by=['object_id', 'passband', 'timeorder'], inplace=True)\n\n    # lag/lead pandas-efficiency-friendly computation\n    list_of_shifts = []\n    values_counts = [(~df[column_name].isnull()).map(int) for column_name in column_names]\n    for shift_step in range(1, agg_window_width + exp_decaying_tail_width + 1):\n        shifted_values = df.shift(-1 * shift_step)\n        shifted_values.loc[(shifted_values['object_id'] != df['object_id']) | (\n                shifted_values['passband'] != df['passband']) |\n                           (shifted_values['timeorder'] - df[\n                               'timeorder'] >= agg_window_width + exp_decaying_tail_width), column_names] = np.nan\n        # if exp_decaying_const != 0:\n        decaying_weights = np.exp(-exp_decaying_const *\n                                  (shifted_values['timeorder'] - df['timeorder']).map(\n                                      lambda x: 0 if x < agg_window_width\n                                      else x - agg_window_width + 1))\n        for column_name in column_names:\n            shifted_values[column_name] *= decaying_weights\n        list_of_shifts.append(shifted_values[column_names])\n\n        for i, column_name in enumerate(column_names):\n            values_counts[i] += (~shifted_values[column_name].isnull()).map(int) * decaying_weights.fillna(1)\n\n    for shift_step in range(-exp_decaying_tail_width, 0):\n        shifted_values = df.shift(-shift_step)\n        shifted_values.loc[(shifted_values['object_id'] != df['object_id']) | (\n                shifted_values['passband'] != df['passband']) | (\n                                   shifted_values['timeorder'] - df[\n                               'timeorder'] < -exp_decaying_tail_width), column_names] = np.nan\n        decaying_weights = np.exp(exp_decaying_const * (shifted_values['timeorder'] - df['timeorder']))\n        for column_name in column_names:\n            shifted_values[column_name] *= decaying_weights\n        list_of_shifts.append(shifted_values[column_names])\n        for i, column_name in enumerate(column_names):\n            values_counts[i] += (~shifted_values[column_name].isnull()).map(int) * decaying_weights.fillna(1)\n\n    df[column_names] = df[column_names].fillna(0)\n    for shifted_series in list_of_shifts:\n        for column_name in column_names:\n            df[column_name] += shifted_series[column_name].fillna(0)\n    for i, column_name in enumerate(column_names):\n        df.loc[df[column_name] == 0, column_name] = np.nan\n        df[column_name] /= values_counts[i]\n\n    df = df[df['timeorder'] % agg_window_width == 0]\n    df['timeorder'] /= agg_window_width\n    return df\n\ndf = moving_mean_smoothed_aggregations(df, agg_window_width=7, exp_decaying_tail_width=3)\nplot_object(id_to_check)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ae81984c5f6676b7adf84dfa75ebf9f55a8d6b23"},"cell_type":"markdown","source":"Next, the most reliable passbands were chosen using heuristics based on overall uncertainty of observations in a passband and relative strength of the signal in the passband."},{"metadata":{"trusted":true,"_uuid":"7634274e9200040cc6354b754f9f83eeddf4d5c2"},"cell_type":"code","source":"flux_err_sum_per_passband = df.groupby(['object_id', 'passband'])['flux_err'].agg(np.nansum)\ndf['flux_abs'] = df['flux'].map(abs)\nflux_sum_per_passband = df.groupby(['object_id', 'passband'])['flux_abs'].agg(np.nansum)\n\nflux_err_sum_per_passband = flux_err_sum_per_passband.reset_index()\nflux_sum_per_passband = flux_sum_per_passband.reset_index()\nflux_sum_per_passband['flux_err_rel_sum'] = flux_err_sum_per_passband['flux_err']/flux_sum_per_passband['flux_abs']\n\nmax_per_object_passband = df.groupby(['object_id', 'passband'])['flux_abs'].max().reset_index()\nmax_per_object_passband.columns = ['object_id', 'passband', 'flux_abs_max']\nmax_per_object_passband['flux_object_abs_max'] = max_per_object_passband['object_id'].map(\n    max_per_object_passband.groupby('object_id')['flux_abs_max'].max())\nmax_per_object_passband['passband_to_total_max_ration'] = max_per_object_passband['flux_abs_max']/max_per_object_passband['flux_object_abs_max']\nmax_per_object_passband['object_id_passband'] = max_per_object_passband['object_id'].map(str) + '_' + max_per_object_passband['passband'].map(str)\nmax_per_object_passband_map = max_per_object_passband.set_index('object_id_passband')['passband_to_total_max_ration']\nflux_sum_per_passband['object_id_passband'] = flux_sum_per_passband['object_id'].map(str) + '_' + flux_sum_per_passband['passband'].map(str)\nflux_sum_per_passband_map = flux_sum_per_passband.set_index('object_id_passband')['flux_err_rel_sum']\n\naccepted_flux_err_rel_sum_per_object = flux_sum_per_passband.groupby(['object_id'])['flux_err_rel_sum'].agg(lambda x: sorted(x)[2])\naccepted_passband_to_total_max_ration_per_object = max_per_object_passband.groupby(['object_id'])['passband_to_total_max_ration'].agg(lambda x: sorted(x,key=lambda i: -i)[3])\n\nflux_sum_per_passband = flux_sum_per_passband[\n    (flux_sum_per_passband['flux_err_rel_sum'] <= flux_sum_per_passband['object_id'].map(accepted_flux_err_rel_sum_per_object)) & (\n    (flux_sum_per_passband['object_id'].map(str) + '_' + flux_sum_per_passband['passband'].map(str)).map(max_per_object_passband_map) >= flux_sum_per_passband['object_id'].map(accepted_passband_to_total_max_ration_per_object))]\n\nreliable_passbands_per_object = flux_sum_per_passband.groupby('object_id')['passband'].agg(set)\nobject_passband_max_flux = max_per_object_passband[['flux_abs_max', 'object_id', 'passband']]\nobject_passband_max_flux['object_id_passband'] = object_passband_max_flux['object_id'].map(str) + '_' + object_passband_max_flux['passband'].map(str)\nobject_passband_max_flux['flux_abs_max_sqrt'] = object_passband_max_flux['flux_abs_max'].map(np.sqrt)\nobject_passband_max_flux = object_passband_max_flux.set_index('object_id_passband')['flux_abs_max_sqrt']\n\nreliable_passbands = df['object_id'].map(reliable_passbands_per_object)\nreliable_passbands_without_current_passband = reliable_passbands - df['passband'].map(lambda x: {x})\nreliable_passbands_entries = reliable_passbands_without_current_passband.map(len) != reliable_passbands.map(len)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"da763d0def48975679925d6f855f02f89489c219"},"cell_type":"markdown","source":"Next, the passband scales were normalized so that the maximum flux is equal to one."},{"metadata":{"trusted":true,"_uuid":"f321dfd6f109c8651a7c6753417cfa904b51a6ec"},"cell_type":"code","source":"max_per_object_passband = df.groupby(['object_id', 'passband'])['flux_abs'].max().reset_index()\nmax_per_object_passband.columns = ['object_id', 'passband', 'flux_abs_max']\nmax_per_object_passband['object_id_passband'] = max_per_object_passband['object_id'].map(str) + '_' + max_per_object_passband['passband'].map(str)\nmax_per_object_passband_map = max_per_object_passband[['object_id_passband', 'flux_abs_max']].set_index('object_id_passband')['flux_abs_max']\ndf_object_id_passband = df['object_id'].map(str) + '_' + df['passband'].map(str)\nmax_flux_per_object_passband = df_object_id_passband.map(max_per_object_passband_map)\ndf['flux'] /= max_flux_per_object_passband","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e7614b67cbd1640584e8269c3ff0ba357931aec9"},"cell_type":"markdown","source":"Reliable passbands were aggregated."},{"metadata":{"trusted":true,"_uuid":"ff475662b4977fb1ec65b75115bcc851825ec078"},"cell_type":"code","source":"object_passband = df['object_id'].map(str) + '_' + df['passband'].map(str)\ndf['passband_max_flux_sqrt'] = object_passband.map(max_per_object_passband_map)/(1 + object_passband.map(flux_sum_per_passband_map))\ndf['flux_weight'] = list(zip(df['flux'], df['passband_max_flux_sqrt']))\n\ndef weighted_mean(timestamp_values):\n    non_nans = list(filter(lambda x: not np.isnan(x[0]), timestamp_values))\n    if len(non_nans) == 0:\n        return np.nan\n    values = list(map(lambda x: x[0], non_nans))\n    if len(values) > 2:\n        mean = np.mean(values)\n        std = np.std(values)\n        lb = mean - std\n        ub = mean + std\n        non_nans = list(filter(lambda x: x[0] >= lb and x[0] <= ub, non_nans))\n    elif len(values) == 2:\n        val_range = max(values) - min(values)\n        if val_range > 0.5:\n            return np.nan\n    elif len(values) == 1:\n        return values[0]\n    \n    if len(non_nans) == 0:\n        return np.nan\n    \n    values = list(map(lambda x: x[0], non_nans))\n    weights = list(map(lambda x: x[1], non_nans))\n    result = np.average(values, weights=weights)\n    return result\n    \ndf = df.groupby(['object_id', 'timeorder'], as_index=False)['flux_weight'].agg(weighted_mean)\ndf.columns = ['object_id', 'timeorder', 'flux']\n\n# omitting tailing missing values\nmax_timeorder_per_object = df[~df['flux'].isnull()].groupby(['object_id'])['timeorder'].max()\nmin_timeorder_per_object = df[~df['flux'].isnull()].groupby(['object_id'])['timeorder'].min()\ndf = df[(df['timeorder'] <= df['object_id'].map(max_timeorder_per_object)) & \n        (df['timeorder'] >= df['object_id'].map(min_timeorder_per_object))]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"969247d7c298d7d473b218989fef16da91e7222f"},"cell_type":"markdown","source":"Visualizing the aggregated light curve for the object."},{"metadata":{"trusted":true,"_uuid":"e0e5ece0639827dd4b1da0487adcd93b978a8481"},"cell_type":"code","source":"def plot_agg_light_curve(object_id=None, object_df=None, class_id='', save=False):\n    assert object_id is not None or object_df is not None\n    if object_df is None:\n        object_df = df[df['object_id']==object_id]\n    else:\n        object_id = object_df['object_id'].iloc[0]\n    plt.figure(figsize=(30,5))\n    plt.plot(object_df['timeorder'], object_df['flux'])\n    plt.scatter(object_df['timeorder'], object_df['flux'])\n    plt.title(f'object id: {object_id}', fontsize=30)\n    if save:\n        if not os.path.exists(f'../input/class_{class_to_proceed}'):\n            os.makedirs(f'../input/class_{class_to_proceed}')\n        plt.savefig(f'../input/class_{class_to_proceed}/{object_id}.png')\n    \nplot_agg_light_curve(id_to_check)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f2ca3ac07ce0702e98de510daee4b2f4612b6fab"},"cell_type":"markdown","source":"Hopefully, the presented aggregation is not completely off, the aggregation definitely can be improved."},{"metadata":{"_uuid":"d918fab075c588d90e535a81a7eb18cde63f037c"},"cell_type":"markdown","source":"### Reconstruction of a representative class light curve from instances with apparent maximum"},{"metadata":{"_uuid":"f5f6b3bbee1947a47ad2107c7ad9fb0845a1f452"},"cell_type":"markdown","source":"For this part we need light curves which likely have the pick. I believe it should not be too hard to come up with reasonable automatic selection of such curves. Given the small amount of training data I handpicked some. Bringing the handpicked object ID's to the kernel."},{"metadata":{"trusted":true,"_uuid":"849b2f1a144d5b4b59c6f7eeaa4edd10e994d6ec"},"cell_type":"code","source":"# object_ids_with_full_pick = [int(name.split('.')[0]) for name in os.listdir(f'../input/class_{class_to_process}_init_base')]\nwith open('../input/plasticcwip/object_ids_with_full_pick_class_62.pkl', 'rb') as f:\n    object_ids_with_full_pick = pickle.load(f)\nobject_ids_class_train = set(df_meta.loc[(df_meta['target'] == class_to_process), 'object_id'].values)\nobject_ids_with_full_pick_train = [object_id for object_id in object_ids_with_full_pick if object_id in object_ids_class_train]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ef63c104872ff117c24570c3e910b4f72a7e4ed1"},"cell_type":"markdown","source":"Next we match light curve instances maximum-to-maximum and average the observations."},{"metadata":{"trusted":true,"_uuid":"c98deeb19a81c41ce6d02cadc1b8bfea204c2439"},"cell_type":"code","source":"def rob_mean(timestamp_values):\n        values = list(filter(lambda x: not np.isnan(x), timestamp_values))\n        if len(values) == 1:\n            return values[0]\n        if len(values) == 0:\n            return np.nan\n        if len(values) > 2:\n            mean = np.mean(values)\n            std = np.std(values)\n            lb = mean - std\n            ub = mean + std\n            values = list(filter(lambda x: x >= lb and x <= ub, values))\n        elif len(values) == 2:\n            val_range = max(values) - min(values)\n            if val_range > 0.5:\n                return np.nan\n        if len(values) == 0:\n            return np.nan\n        result = np.mean(values)\n        return result\n    \ndef shift_array_and_vector(ts_1_input, ts_collection, shift_0, original_ts_collection_region=False):\n    # shift = 'ts_collection' - 'ts_1'\n    if not original_ts_collection_region:\n        if shift_0 < 0:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((ts_1_input,\n                                                 np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n                ts_collection = np.concatenate((np.full((-shift_0, ts_collection.shape[1]), np.nan),\n                                                ts_collection,\n                                                np.full((max(0, len(ts_1_input) - ts_collection.shape[0] + shift_0),\n                                                         ts_collection.shape[1]), np.nan)), axis=0)\n            else:\n                ts_1_shifted_0 = np.concatenate((ts_1_input, np.full(-shift_0 - len(ts_1_input) + ts_collection.shape[0], np.nan)))\n                ts_collection = np.concatenate((np.full((-shift_0, ts_collection.shape[1]), np.nan), ts_collection), axis=0)\n        else:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan), ts_1_input))\n                ts_collection = np.concatenate((ts_collection,\n                                                np.full((shift_0 + len(ts_1_input) - ts_collection.shape[0], ts_collection.shape[1]), np.nan)), axis=0)\n            else:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan),\n                                                 ts_1_input,\n                                                 np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n                ts_collection = np.concatenate((ts_collection, np.full((max(0, shift_0 + len(ts_1_input) - ts_collection.shape[0]),\n                                                                        ts_collection.shape[1]), np.nan)), axis=0)\n    else:\n        if shift_0 < 0:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((ts_1_input[-shift_0: len(ts_1_input) - max(0, shift_0 + len(ts_1_input) - ts_collection.shape[0])],\n                                                np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n            else:\n                ts_1_shifted_0 = np.concatenate((ts_1_input[-shift_0:],\n                                                 np.full(-shift_0 - len(ts_1_input) + ts_collection.shape[0], np.nan)))\n        else:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan),\n                                                 ts_1_input[:-shift_0 - len(ts_1_input) + ts_collection.shape[0]]))\n            else:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan),\n                                                 ts_1_input[:len(ts_1_input) - max(0, shift_0 + len(ts_1_input) - ts_collection.shape[0])],\n                                                 np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n\n    return ts_1_shifted_0, ts_collection\n\ndef match_full_picks(ts_1_input, ts_collection, debug=False, iteration=''):\n    \n    ts_2_input = np.apply_along_axis(rob_mean, 1, ts_collection)\n\n    ts_1 = np.nan_to_num(ts_1_input)\n    ts_2 = np.nan_to_num(ts_2_input)\n    idx_max_ts_1 = np.argmax(ts_1)\n    idx_max_ts_2 = np.argmax(ts_2)\n        \n    shift_0 = idx_max_ts_2 - idx_max_ts_1\n    if shift_0 < 0:\n        if len(ts_1) > len(ts_2):\n            ts_1_shifted_0 = np.concatenate((ts_1, np.zeros(max(0, -shift_0 - (len(ts_1) - len(ts_2))))))\n            ts_2_shifted_0 = np.concatenate((np.zeros(-shift_0), ts_2, np.zeros(max(0, len(ts_1) - len(ts_2) + shift_0))))\n        else:            \n            ts_1_shifted_0 = np.concatenate((ts_1, np.zeros(-shift_0 - len(ts_1) + len(ts_2))))\n            ts_2_shifted_0 = np.concatenate((np.zeros(-shift_0), ts_2))\n    else:\n        if len(ts_1) > len(ts_2):\n            ts_1_shifted_0 = np.concatenate((np.zeros(shift_0), ts_1))\n            ts_2_shifted_0 = np.concatenate((ts_2, np.zeros(shift_0 + len(ts_1) - len(ts_2))))\n        else:            \n            ts_1_shifted_0 = np.concatenate((np.zeros(shift_0), ts_1, np.zeros(max(0, -shift_0 - len(ts_1) + len(ts_2)))))\n            ts_2_shifted_0 = np.concatenate((ts_2, np.zeros(max(0, shift_0 + len(ts_1) - len(ts_2)))))\n        \n    if debug:\n        plot_args = list(range(len(ts_1_shifted_0)))\n        plt.figure(figsize=(20,4))\n        plt.plot(plot_args, ts_2_shifted_0)\n        plt.plot(plot_args, ts_1_shifted_0)\n        plt.scatter(plot_args, ts_2_shifted_0)\n        plt.scatter(plot_args, ts_1_shifted_0)\n        plt.title(f'Class representative so far + next curve')\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.legend(['class representative', 'light curve'])\n        \n    min_msa = mean_squared_error(ts_1_shifted_0, ts_2_shifted_0)\n    min_shift_added = 0\n    ts_1_shifted, ts_2_shifted = ts_1_shifted_0, ts_2_shifted_0\n    for shift_added in range(1, 6):\n        ts_1_shifted = np.concatenate((np.zeros(shift_added), ts_1_shifted))\n        ts_2_shifted = np.concatenate((ts_2_shifted, np.zeros(shift_added)))\n        msa = mean_squared_error(ts_1_shifted, ts_2_shifted)\n        if msa < min_msa:\n            min_msa = msa\n            min_shift_added = shift_added\n            \n    ts_1_shifted, ts_2_shifted = ts_1_shifted_0, ts_2_shifted_0  \n    for shift_added in range(-5, 0):\n        ts_1_shifted = np.concatenate((ts_1_shifted, np.zeros(-shift_added)))\n        ts_2_shifted = np.concatenate((np.zeros(-shift_added), ts_2_shifted))\n        msa = mean_squared_error(ts_1_shifted, ts_2_shifted)\n        if msa < min_msa:\n            min_msa = msa\n            min_shift_added = shift_added\n            \n    shift_0 += min_shift_added  \n    \n    ts_1_shifted_0, ts_collection = shift_array_and_vector(ts_1_input, ts_collection, shift_0, original_ts_collection_region=True)\n    \n    ts_collection = np.concatenate((ts_collection, ts_1_shifted_0.reshape(-1, 1)), axis=1)\n    \n    if debug:\n        ts_final = np.apply_along_axis(rob_mean, 1, ts_collection)\n\n        plt.figure(figsize=(20,4))\n        plt.plot(range(len(ts_final)), ts_final)\n        plt.scatter(range(len(ts_final)), ts_final)        \n        plt.title(f'Class {class_to_process} representative after {iteration} iterations')\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.figure(figsize=(20,4))\n        plt.plot(range(len(ts_final)), np.ones(len(ts_final)))\n    return ts_1_shifted_0, ts_collection","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"53b48088857eabbf0f561d4825a9a0e2c48b71b8"},"cell_type":"markdown","source":"Let's visualize the process of reconstructing the class representative. *I found it beautiful, seeing how trimmed mean uncovered structure from the noise. I just hope the uncovered ligth curve is meaningful.*"},{"metadata":{"trusted":true,"_uuid":"482975bca2ec8a1dca7aa0cf1516e74bde0f53df"},"cell_type":"code","source":"def rob_mean(timestamp_values):\n        values = list(filter(lambda x: not np.isnan(x), timestamp_values))\n        if len(values) == 1:\n            return values[0]\n        if len(values) == 0:\n            return np.nan\n        if len(values) > 2:\n            mean = np.mean(values)\n            std = np.std(values)\n            lb = mean - std\n            ub = mean + std\n            values = list(filter(lambda x: x >= lb and x <= ub, values))\n        elif len(values) == 2:\n            val_range = max(values) - min(values)\n            if val_range > 0.5:\n                return np.nan\n        if len(values) == 0:\n            return np.nan\n        result = np.mean(values)\n        return result\n    \ndef shift_array_and_vector(ts_1_input, ts_collection, shift_0, original_ts_collection_region=False):\n    # shift = 'ts_collection' - 'ts_1'\n    if not original_ts_collection_region:\n        if shift_0 < 0:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((ts_1_input,\n                                                 np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n                ts_collection = np.concatenate((np.full((-shift_0, ts_collection.shape[1]), np.nan),\n                                                ts_collection,\n                                                np.full((max(0, len(ts_1_input) - ts_collection.shape[0] + shift_0),\n                                                         ts_collection.shape[1]), np.nan)), axis=0)\n            else:\n                ts_1_shifted_0 = np.concatenate((ts_1_input, np.full(-shift_0 - len(ts_1_input) + ts_collection.shape[0], np.nan)))\n                ts_collection = np.concatenate((np.full((-shift_0, ts_collection.shape[1]), np.nan), ts_collection), axis=0)\n        else:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan), ts_1_input))\n                ts_collection = np.concatenate((ts_collection,\n                                                np.full((shift_0 + len(ts_1_input) - ts_collection.shape[0], ts_collection.shape[1]), np.nan)), axis=0)\n            else:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan),\n                                                 ts_1_input,\n                                                 np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n                ts_collection = np.concatenate((ts_collection, np.full((max(0, shift_0 + len(ts_1_input) - ts_collection.shape[0]),\n                                                                        ts_collection.shape[1]), np.nan)), axis=0)\n    else:\n        if shift_0 < 0:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((ts_1_input[-shift_0: len(ts_1_input) - max(0, shift_0 + len(ts_1_input) - ts_collection.shape[0])],\n                                                np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n            else:\n                ts_1_shifted_0 = np.concatenate((ts_1_input[-shift_0:],\n                                                 np.full(-shift_0 - len(ts_1_input) + ts_collection.shape[0], np.nan)))\n        else:\n            if len(ts_1_input) > ts_collection.shape[0]:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan),\n                                                 ts_1_input[:-shift_0 - len(ts_1_input) + ts_collection.shape[0]]))\n            else:\n                ts_1_shifted_0 = np.concatenate((np.full(shift_0, np.nan),\n                                                 ts_1_input[:len(ts_1_input) - max(0, shift_0 + len(ts_1_input) - ts_collection.shape[0])],\n                                                 np.full(max(0, -shift_0 - len(ts_1_input) + ts_collection.shape[0]), np.nan)))\n\n    return ts_1_shifted_0, ts_collection\n\ndef match_full_picks(ts_1_input, ts_collection, debug=False, iteration=''):\n    \n    ts_2_input = np.apply_along_axis(rob_mean, 1, ts_collection)\n\n    ts_1 = np.nan_to_num(ts_1_input)\n    ts_2 = np.nan_to_num(ts_2_input)\n    idx_max_ts_1 = np.argmax(ts_1)\n    idx_max_ts_2 = np.argmax(ts_2)\n        \n    shift_0 = idx_max_ts_2 - idx_max_ts_1\n    if shift_0 < 0:\n        if len(ts_1) > len(ts_2):\n            ts_1_shifted_0 = np.concatenate((ts_1, np.zeros(max(0, -shift_0 - (len(ts_1) - len(ts_2))))))\n            ts_2_shifted_0 = np.concatenate((np.zeros(-shift_0), ts_2, np.zeros(max(0, len(ts_1) - len(ts_2) + shift_0))))\n        else:            \n            ts_1_shifted_0 = np.concatenate((ts_1, np.zeros(-shift_0 - len(ts_1) + len(ts_2))))\n            ts_2_shifted_0 = np.concatenate((np.zeros(-shift_0), ts_2))\n    else:\n        if len(ts_1) > len(ts_2):\n            ts_1_shifted_0 = np.concatenate((np.zeros(shift_0), ts_1))\n            ts_2_shifted_0 = np.concatenate((ts_2, np.zeros(shift_0 + len(ts_1) - len(ts_2))))\n        else:            \n            ts_1_shifted_0 = np.concatenate((np.zeros(shift_0), ts_1, np.zeros(max(0, -shift_0 - len(ts_1) + len(ts_2)))))\n            ts_2_shifted_0 = np.concatenate((ts_2, np.zeros(max(0, shift_0 + len(ts_1) - len(ts_2)))))\n        \n    if debug:\n        plot_args = list(range(len(ts_1_shifted_0)))\n        plt.figure(figsize=(20,4))\n        plt.plot(plot_args, ts_2_shifted_0)\n        plt.plot(plot_args, ts_1_shifted_0)\n        plt.scatter(plot_args, ts_2_shifted_0)\n        plt.scatter(plot_args, ts_1_shifted_0)\n        plt.title(f'Class representative so far + next curve')\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.legend(['class representative', 'light curve'])\n        \n    min_msa = mean_squared_error(ts_1_shifted_0, ts_2_shifted_0)\n    min_shift_added = 0\n    ts_1_shifted, ts_2_shifted = ts_1_shifted_0, ts_2_shifted_0\n    for shift_added in range(1, 6):\n        ts_1_shifted = np.concatenate((np.zeros(shift_added), ts_1_shifted))\n        ts_2_shifted = np.concatenate((ts_2_shifted, np.zeros(shift_added)))\n        msa = mean_squared_error(ts_1_shifted, ts_2_shifted)\n        if msa < min_msa:\n            min_msa = msa\n            min_shift_added = shift_added\n            \n    ts_1_shifted, ts_2_shifted = ts_1_shifted_0, ts_2_shifted_0  \n    for shift_added in range(-5, 0):\n        ts_1_shifted = np.concatenate((ts_1_shifted, np.zeros(-shift_added)))\n        ts_2_shifted = np.concatenate((np.zeros(-shift_added), ts_2_shifted))\n        msa = mean_squared_error(ts_1_shifted, ts_2_shifted)\n        if msa < min_msa:\n            min_msa = msa\n            min_shift_added = shift_added\n            \n    shift_0 += min_shift_added  \n    \n    ts_1_shifted_0, ts_collection = shift_array_and_vector(ts_1_input, ts_collection, shift_0, original_ts_collection_region=True)\n    \n    ts_collection = np.concatenate((ts_collection, ts_1_shifted_0.reshape(-1, 1)), axis=1)\n    \n    if debug:\n        ts_final = np.apply_along_axis(rob_mean, 1, ts_collection)\n\n        plt.figure(figsize=(20,4))\n        plt.plot(range(len(ts_final)), ts_final)\n        plt.scatter(range(len(ts_final)), ts_final)        \n        plt.title(f'Class {class_to_process} representative after {iteration} iterations')\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.figure(figsize=(20,4))\n        plt.plot(range(len(ts_final)), np.ones(len(ts_final)))\n    return ts_1_shifted_0, ts_collection\n\n# the matching\nbase_prototype = df.loc[df['object_id']==object_ids_with_full_pick_train[0], 'flux'].values.reshape(-1, 1)\n\nfor i, object_id in enumerate(object_ids_with_full_pick_train[1:]):\n    _, base_prototype = match_full_picks(df.loc[df['object_id']==object_id, 'flux'].values, base_prototype, debug=i%40==0, iteration=i+1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"42f76cb60ccbb71a8139001408b24c6466da9bf2"},"cell_type":"code","source":"base_prototype = base_prototype[37:67,:]\nfinal_base_prototype_ts = np.apply_along_axis(rob_mean, 1, base_prototype)\nplt.figure(figsize=(20,4))\nplt.plot(range(len(final_base_prototype_ts)), final_base_prototype_ts)\nplt.scatter(range(len(final_base_prototype_ts)), final_base_prototype_ts)\nplt.xlabel('timestamp order')\nplt.ylabel('flux')\nplt.title(f'Class {class_to_process} representative light curve')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2b9638078ee17aa302748e2fd7f803b0aac771ce"},"cell_type":"markdown","source":"Next, we use light curves without a certain maximum, rescale and shift it appropriately, then use for the class representative reconstruction. \n\nTo find the best match (timestamp shift position + flux rescaling), the following is minimized: MSE against representative-curve divided by mean of squares of the light curve. The intuition: if the representative curve is positioned in the worst possible way having only zeros below the considered light curve instance, then we have maximum value of MSE, which is equal to the mean of light curve fluxes squared. We start from there and would like to improve as much as possible. Direct MSE optimization would result in preferring smaller flux scales."},{"metadata":{"trusted":true,"_uuid":"d25fce0ec53cce95b851702f5cef0ef62fb4d62e"},"cell_type":"code","source":"object_ids_with_full_pick_train = set(object_ids_with_full_pick_train)\nobject_ids_without_full_pick_train = [object_id for object_id in object_ids_class_train if object_id not in object_ids_with_full_pick_train]\n\ndef squared_error_sum(x, y):\n    return np.sum([(x[i] - y[i])**2 for i in range(len(x))])\n\ndef match_partials(ts_1_input, ts_collection, debug=False, iteration=''):\n    \n    ts_2_input = np.apply_along_axis(rob_mean, 1, ts_collection)    \n    ts_2 = np.nan_to_num(ts_2_input)\n    \n    squared_error_min = float('inf')\n    scale_factor_opt = -1\n    ts_1_pos_under_ts_2_start_opt = -1\n    \n    for scale_factor in np.arange(0.1,1.01,0.1):\n        ts_1 = scale_factor*ts_1_input\n            \n        sum_squares_current = np.nansum(ts_1[:-4]**2)\n        sum_squares_total = np.nansum(ts_1**2)\n        \n        # shifting a ts_1_input, abstracting ts_2 as a solid, not moving part  \n        for ts_1_pos_under_ts_2_start in range(len(ts_1) - 5, -(len(ts_2) - 5), -1):\n            ts_1_intersection_part = ts_1[max(0, ts_1_pos_under_ts_2_start):\n                                          min(len(ts_1), \n                                              len(ts_2) + ts_1_pos_under_ts_2_start)]\n            ts_2_intersection_part = ts_2[max(0, -ts_1_pos_under_ts_2_start): \n                                          min(len(ts_2), \n                                             len(ts_1) - ts_1_pos_under_ts_2_start)]         \n                            \n            mask = ~np.isnan(ts_1_intersection_part)\n            \n            squared_error_intersection = squared_error_sum(ts_1_intersection_part[mask], ts_2_intersection_part[mask])\n            \n            # managing non-intersection part of ts_1\n            if ts_1_pos_under_ts_2_start > 0:\n                val = ts_1[ts_1_pos_under_ts_2_start]**2 if not np.isnan(ts_1[ts_1_pos_under_ts_2_start]) else 0\n                sum_squares_current -= val\n                \n            if len(ts_1) - ts_1_pos_under_ts_2_start > len(ts_2):\n                val = ts_1[len(ts_2) + ts_1_pos_under_ts_2_start]**2 if not np.isnan(ts_1[len(ts_2) + ts_1_pos_under_ts_2_start]) else 0\n                sum_squares_current += val\n                \n            squared_error_total = (squared_error_intersection + sum_squares_current)/sum_squares_total\n            if squared_error_total < squared_error_min:\n                squared_error_min = squared_error_total\n                scale_factor_opt = scale_factor\n                ts_1_pos_under_ts_2_start_opt = ts_1_pos_under_ts_2_start\n                \n    ts_1_input *= scale_factor_opt\n    \n    ts_1_shifted_0, ts_collection = shift_array_and_vector(ts_1_input, ts_collection, -ts_1_pos_under_ts_2_start_opt, \n                                                           original_ts_collection_region=True)\n    t2_shifted = np.apply_along_axis(rob_mean, 1, ts_collection)\n                \n    if debug:\n        plt.figure(figsize=(20,4))\n        plt.plot(range(len(t2_shifted)), t2_shifted)\n        plt.plot(range(len(ts_1_shifted_0)), ts_1_shifted_0)\n        plt.scatter(range(len(t2_shifted)), t2_shifted)\n        plt.scatter(range(len(ts_1_shifted_0)), ts_1_shifted_0)\n        plt.title(f'Class representative so far + next curve')\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.legend(['class representative', 'light curve'])\n    \n    ts_collection = np.concatenate((ts_collection, ts_1_shifted_0.reshape(-1, 1)), axis=1)\n    \n    if debug:\n        ts_final = np.apply_along_axis(rob_mean, 1, ts_collection)\n        plt.figure(figsize=(20,4))\n        plt.plot(range(len(ts_final)), ts_final)\n        plt.scatter(range(len(ts_final)), ts_final)  \n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.title(f'Class {class_to_process} representative after {iteration} iterations')\n    return ts_collection\n\nfor i, object_id in enumerate(object_ids_without_full_pick_train):\n    base_prototype = match_partials(df.loc[df['object_id']==object_id, 'flux'].values, base_prototype, debug=i%70==0, iteration=i)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"efdd400d65b74ae34d1f5e70554ab4d4244b4667"},"cell_type":"code","source":"final_prototype_ts = np.apply_along_axis(rob_mean, 1, base_prototype)\nplt.figure(figsize=(20,4))\nplt.plot(range(len(final_prototype_ts)), final_prototype_ts)\nplt.scatter(range(len(final_prototype_ts)), final_prototype_ts)\nplt.xlabel('timestamp order')\nplt.ylabel('flux')\nplt.title(f'Final class {class_to_process} representative light curve')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f60e94170a5a70b50c48256e5a8f9434fe38010d"},"cell_type":"markdown","source":"Analogously, representatives for all other hard classes are computed. "},{"metadata":{"trusted":true,"_uuid":"1262ed61bdf6bae120679f857d9f7695301d6b81"},"cell_type":"code","source":"fold_idx = 0\nclasses = [42, 52, 62, 67, 90, 6, 15, 95, 64]\nclass_prototypes = dict()\nfor class_name in classes:\n    with open(f'../input/plasticcwip/final_prototype_ts_class_{class_name}_fold_{fold_idx}.pkl', 'rb') as f:\n        plt.figure(figsize=(20,4))\n        final_prototype_ts = pickle.load(f)\n        final_prototype_ts, base_prototype = match_full_picks(final_prototype_ts, base_prototype)\n        class_prototypes[class_name] = np.nan_to_num(final_prototype_ts)\n        plt.plot(range(len(final_prototype_ts)), final_prototype_ts)\n        plt.scatter(range(len(final_prototype_ts)), final_prototype_ts)\n        plt.xlabel('timestamp order')\n        plt.ylabel('flux')\n        plt.title(f'Final {class_name} class representative light curve')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b4f4332983e6c1e34c2ffcf5e1f19e42b4951d83"},"cell_type":"markdown","source":"### Faster implementation"},{"metadata":{"_uuid":"8dfc8ddf009682de1a3534cde27ae5e0489b8ecf"},"cell_type":"markdown","source":"Just in case someone might experiment with a similar approach in the future, I'd like share a faster implementation in cython. Hopefully, it can be useful as a starting point. In the provided functions cosine-distance-based matching is also implemented. \n\nHaving a small length of representative class curves, the computational complexity of algorithm is in O(n) of light curve length."},{"metadata":{"trusted":true,"_uuid":"41505bf5d3ab99f502a76d256ca77289c8524e6c"},"cell_type":"code","source":"%load_ext Cython","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1c2dd5480bb1cf5eb8270e4a7a71e3b059b35e9d"},"cell_type":"code","source":"%%cython\n\ncimport numpy as cnp\nfrom numpy cimport ndarray\nctypedef unsigned char uint8\nimport numpy as np\nimport matplotlib.pyplot as plt\n\ncnp.import_array()\ncdef double[:,:] return_empty_2d(int dim1, int dim2):\n    cdef cnp.npy_intp* dims = [dim1, dim2]\n    return cnp.PyArray_SimpleNew(2, dims, cnp.NPY_DOUBLE)\n\ncdef double[:] return_empty_1d_double(int dim):\n    cdef cnp.npy_intp dims = dim\n    return cnp.PyArray_SimpleNew(1, &dims, cnp.NPY_DOUBLE)\n\ncdef int[:] return_empty_1d_int(int dim):\n    cdef cnp.npy_intp dims = dim\n    return cnp.PyArray_SimpleNew(1, &dims, cnp.NPY_INT)\n\ncdef extern from \"math.h\":\n    double hypot(double x, double y)\n\ncdef extern from \"math.h\":\n    double INFINITY\n\ncdef extern from \"math.h\":\n    bint isnan(double x)\n\ncdef extern from \"math.h\":\n    double fabs(double x)\n\ncdef int INT_NAN = -999\ncdef int BIG_M_SEQ_LEN = 10000 # no sequence can be longer\ncdef double NOISE_UB = 0.04\ncdef double NOISE_LB = -0.01\ncdef int INITIAL_INTERSECTION_WIDTH = 5\n\n#cdef double weight_direction_error = 70\n\n\ncdef double[:] squared_error_sum(double[:] pred, double[:] true):\n    cdef int N = pred.shape[0]\n    cdef double sum_err = 0\n    cdef double sum_sq_true = 0\n    cdef int i\n    for i in range(N):\n        if not isnan(pred[i]): # prototypes have no missing values, therefore there is no need to check y\n            sum_err += (pred[i] - true[i])**2\n            sum_sq_true += true[i]**2\n    cdef double[:] result = return_empty_1d_double(2)\n    result[0] = sum_err\n    result[1] = sum_sq_true\n    return result\n\ncdef double[:,:] get_segments(double[:,:] ts_points, NUMBER_OF_POINTS_TS_1=-1):\n    cdef int N\n    if NUMBER_OF_POINTS_TS_1 == -1:\n        N = ts_points.shape[0] - 1\n    else:\n        N = NUMBER_OF_POINTS_TS_1 - 1\n\n    if N < 0:\n        return return_empty_2d(0, 2)\n    cdef double[:,:] segments = return_empty_2d(N, 2)\n    cdef int i\n    for i in range(N):\n        segments[i, 0] = ts_points[i + 1, 0] - ts_points[i, 0]\n        segments[i, 1] = ts_points[i + 1, 1] - ts_points[i, 1]\n    return segments\n\n\ncdef double cos_dist(double x1, double y1, double x2, double y2):\n    cdef double vector_product = x1 * x2 + y1 * y2\n    cdef double cos = vector_product / (hypot(x1, y1) * hypot(x2, y2))\n    return 1 - cos\n\ncdef double cos_dist_ts_to_be_masked(double[:,:] segments_ts_1_intersection, double[:] ts_2_intersection_part, double[:] ts_1_intersection_part):\n    if len(segments_ts_1_intersection) == 0:\n        return 0\n    cdef int number_of_points = len(segments_ts_1_intersection) + 1\n    cdef double[:,:] ts_2_points = return_empty_2d(number_of_points, 2)\n    cdef int current_point_idx = 0\n    cdef int i\n    cdef int N = ts_1_intersection_part.shape[0]\n    for i in range(N):\n        if not isnan(ts_1_intersection_part[i]):\n            ts_2_points[current_point_idx,0] = i\n            ts_2_points[current_point_idx,1] = ts_2_intersection_part[i]\n            current_point_idx += 1\n    cdef double[:,:] ts_2_segments = get_segments(ts_2_points)\n    cdef double cos_dist_sum = 0\n    for i in range(len(segments_ts_1_intersection)):\n        cos_dist_sum += cos_dist(segments_ts_1_intersection[i,0], segments_ts_1_intersection[i,1], ts_2_segments[i,0], ts_2_segments[i,1])\n    return cos_dist_sum\n\n\ndef estimate_scale(double[:] ts_1_input, double[:] ts_2, int debug=1, class_name='', time_scale=''):\n\n    cdef double error_sq_min = INFINITY\n    cdef double error_cos_min = INFINITY\n    cdef double scale_factor_opt_sq = -1.0\n    cdef double scale_factor_opt_cos = -1.0\n    cdef double[:] scale_factor_array = np.arange(0.1, 1.001, 0.1)\n    cdef int LEN_TS_1 = len(ts_1_input)\n    cdef int LEN_TS_2 = len(ts_2)\n    cdef double[:] ts_1 = return_empty_1d_double(LEN_TS_1)\n    cdef double[:] ts_1_points_orig = return_empty_1d_double(LEN_TS_1)\n    cdef double[:,:] ts_1_points = return_empty_2d(LEN_TS_1, 2)\n    cdef double[:,:] ts_1_segments, ts_1_segment_values, segments_ts_1_before_intersection, segments_ts_1_after_intersection, segments_ts_1_non_intersection\n    cdef double[:] segment_base_cos_dists = return_empty_1d_double(LEN_TS_1)\n    cdef double[:] sq_errors, segments_cos_dists_before_intersection, segments_cos_dists_after_intersection, segments_cos_dists_non_intersection, cos_dist_before_intersection, cos_dist_after_intersection\n    cdef double segment, sum_squares_outside_intersection_current, sum_squares_total, sum_segment_cos_dists_total, sum_segment_cos_dists_outside_intersection_current, scale_factor, val, cos_dist_intersection, squared_error_intersection\n    cdef int[:] timestamp_segment_start_2_segment_idx = np.ones(LEN_TS_1, dtype=np.int32) * INT_NAN\n    cdef int[:] timestamp_segment_end_2_segment_idx = np.ones(LEN_TS_1, dtype=np.int32) * INT_NAN\n    cdef int[:] timestamp_2_start_segment_idx_higher = return_empty_1d_int(LEN_TS_1)\n    cdef int[:] timestamp_2_end_segment_idx_lower = return_empty_1d_int(LEN_TS_1)\n    cdef int[:] timestamp_2_point_idx_higher = return_empty_1d_int(LEN_TS_1)\n    cdef int i, segment_idx, start_point_idx, scale_factor_idx, ts_1_pos_under_ts_2_start, current_min, current_max, number_of_segments_outside_intersection_current, initial_intersection_point_idx, ts_1_start_idx, ts_1_end_idx, NUMBER_OF_POINTS_TS_1\n\n    # init\n    NUMBER_OF_POINTS_TS_1 = 0\n    for i in range(LEN_TS_1):\n        if not isnan(ts_1_input[i]):\n            ts_1_points[NUMBER_OF_POINTS_TS_1, 0] = i\n            ts_1_points_orig[NUMBER_OF_POINTS_TS_1] = ts_1_input[i]\n            NUMBER_OF_POINTS_TS_1 += 1\n\n    for segment_idx in range(NUMBER_OF_POINTS_TS_1 - 1):\n        start_point_idx = <int>ts_1_points[segment_idx, 0]\n        timestamp_segment_start_2_segment_idx[start_point_idx] = segment_idx\n    for i in range(1, NUMBER_OF_POINTS_TS_1):\n        end_point_idx = <int>ts_1_points[i, 0]\n        timestamp_segment_end_2_segment_idx[end_point_idx] = i - 1\n\n    current_min = BIG_M_SEQ_LEN\n    timestamp_2_start_segment_idx_higher[-1] = current_min\n    # the wrongly filled \"tail\" of timestamp_2_start_segment_idx_higher/ \"head\" of timestamp_2_end_segment_idx_lower don't bother us, it would result into an empty slice selection\n    for i in range(LEN_TS_1 - 2, -1, -1):\n        if timestamp_segment_start_2_segment_idx[i] != INT_NAN:\n            current_min = timestamp_segment_start_2_segment_idx[i]\n        timestamp_2_start_segment_idx_higher[i] = current_min\n    current_max = -1\n    timestamp_2_end_segment_idx_lower[0] = current_max\n    for i in range(LEN_TS_1):\n        if timestamp_segment_end_2_segment_idx[i] != INT_NAN:\n            current_max = timestamp_segment_end_2_segment_idx[i]\n        timestamp_2_end_segment_idx_lower[i] = current_max\n    current_min = NUMBER_OF_POINTS_TS_1\n    timestamp_2_point_idx_higher[-1] = current_min\n    for i in range(LEN_TS_1 - 2, -1, -1):\n        if not isnan(ts_1_input[i]):\n            current_min -= 1\n        timestamp_2_point_idx_higher[i] = current_min\n\n    # main loop checking different scales\n    for scale_factor_idx in range(10):\n        scale_factor = scale_factor_array[scale_factor_idx]\n        for i in range(LEN_TS_1):\n            ts_1[i] = scale_factor * ts_1_input[i]\n            if ts_1[i] < NOISE_UB and ts_1[i] > NOISE_LB:\n                ts_1[i] = 0\n\n        for i in range(NUMBER_OF_POINTS_TS_1):\n            ts_1_points[i, 1] = scale_factor * ts_1_points_orig[i]\n\n        ts_1_segments = get_segments(ts_1_points, NUMBER_OF_POINTS_TS_1)\n\n        # a helping variable to keep track of values of interest outside of intersection\n        sum_segment_cos_dists_total = 0.0\n        number_of_segments_outside_intersection_current = ts_1_segments.shape[0]\n        for i in range(ts_1_segments.shape[0]):\n            segment_base_cos_dists[i] = cos_dist(ts_1_segments[i, 0], ts_1_segments[i, 1], ts_1_segments[i, 0], 0.0)\n            sum_segment_cos_dists_total += segment_base_cos_dists[i]\n\n        sum_squares_total = 0.0\n        for i in range(NUMBER_OF_POINTS_TS_1):\n            sum_squares_total += ts_1_points[i, 1]**2\n        if sum_squares_total == 0:\n            continue\n\n        initial_intersection_point_idx = NUMBER_OF_POINTS_TS_1 - 1\n\n        sum_squares_outside_intersection_current = sum_squares_total\n        while initial_intersection_point_idx >= 0 and ts_1_points[initial_intersection_point_idx, 0] >= LEN_TS_1 - (INITIAL_INTERSECTION_WIDTH - 1):\n            sum_squares_outside_intersection_current -= ts_1_points[initial_intersection_point_idx, 1]**2\n            initial_intersection_point_idx -= 1\n\n        sum_segment_cos_dists_outside_intersection_current = sum_segment_cos_dists_total\n        for i in range(LEN_TS_1 - (INITIAL_INTERSECTION_WIDTH - 1), LEN_TS_1):\n            if timestamp_segment_start_2_segment_idx[i] != INT_NAN:\n                segment_idx = timestamp_segment_start_2_segment_idx[i]\n                sum_segment_cos_dists_outside_intersection_current -= segment_base_cos_dists[segment_idx]\n                number_of_segments_outside_intersection_current -= 1\n\n        # shifting ts_1_input against ts_2\n        for ts_1_pos_under_ts_2_start in range(LEN_TS_1 - INITIAL_INTERSECTION_WIDTH, -(LEN_TS_2 - INITIAL_INTERSECTION_WIDTH), -1):\n            ts_1_start_idx, ts_1_end_idx = max(0, ts_1_pos_under_ts_2_start), min(LEN_TS_1,\n                                                                                  LEN_TS_2 + ts_1_pos_under_ts_2_start)\n            ts_1_intersection_part = ts_1[ts_1_start_idx: ts_1_end_idx]\n            ts_2_intersection_part = ts_2[max(0, -ts_1_pos_under_ts_2_start):\n                                          min(LEN_TS_2,\n                                              LEN_TS_1 - ts_1_pos_under_ts_2_start)]\n\n            sq_errors = squared_error_sum(ts_1_intersection_part, ts_2_intersection_part)\n            squared_error_intersection = sq_errors[0]\n            squared_baseline_intersection = sq_errors[1]\n            segments_ts_1_intersection = ts_1_segments[timestamp_2_start_segment_idx_higher[ts_1_start_idx]:\n                                                       timestamp_2_end_segment_idx_lower[ts_1_end_idx - 1] + 1]\n\n            segment_cos_dists_intersection = cos_dist_ts_to_be_masked(segments_ts_1_intersection, ts_2_intersection_part, ts_1_intersection_part)\n            # managing non-intersection part of ts_1\n            if ts_1_pos_under_ts_2_start > 0:\n                if not isnan(ts_1[ts_1_pos_under_ts_2_start]):\n                    val = ts_1[ts_1_pos_under_ts_2_start] ** 2\n                    sum_squares_outside_intersection_current -= val\n                    if timestamp_segment_start_2_segment_idx[ts_1_pos_under_ts_2_start] != INT_NAN:\n                        segment_idx = timestamp_segment_start_2_segment_idx[ts_1_pos_under_ts_2_start]\n                        sum_segment_cos_dists_outside_intersection_current -= segment_base_cos_dists[segment_idx]\n                        number_of_segments_outside_intersection_current -= 1\n            if LEN_TS_1 - ts_1_pos_under_ts_2_start > LEN_TS_2:\n                if not isnan(ts_1[LEN_TS_2 + ts_1_pos_under_ts_2_start]):\n                    val = ts_1[LEN_TS_2 + ts_1_pos_under_ts_2_start] ** 2\n                    sum_squares_outside_intersection_current += val\n                    if timestamp_segment_end_2_segment_idx[LEN_TS_2 + ts_1_pos_under_ts_2_start] != INT_NAN:\n                        segment_idx = timestamp_segment_end_2_segment_idx[LEN_TS_2 + ts_1_pos_under_ts_2_start]\n                        sum_segment_cos_dists_outside_intersection_current += segment_base_cos_dists[segment_idx]\n                        number_of_segments_outside_intersection_current += 1\n\n            squared_error_total = squared_error_intersection + sum_squares_outside_intersection_current\n            squared_error_unreduced_portion = squared_error_total / sum_squares_total\n\n            segment_cos_dists_unreduced_portion = (segment_cos_dists_intersection + sum_segment_cos_dists_outside_intersection_current) / sum_segment_cos_dists_total\n\n            if squared_error_unreduced_portion < error_sq_min:\n                error_sq_min = squared_error_unreduced_portion\n                scale_factor_opt_sq = scale_factor\n                squared_error_unreduced_portion_opt_sq, segment_cos_dists_unreduced_portion_opt_sq = squared_error_unreduced_portion, segment_cos_dists_unreduced_portion\n\n            if segment_cos_dists_unreduced_portion < error_cos_min:\n                error_cos_min = segment_cos_dists_unreduced_portion\n                scale_factor_opt_cos = scale_factor\n                squared_error_unreduced_portion_opt_cos, segment_cos_dists_unreduced_portion_opt_cos = squared_error_unreduced_portion, segment_cos_dists_unreduced_portion\n\n    return scale_factor_opt_cos, error_cos_min, scale_factor_opt_sq, error_sq_min, segment_cos_dists_unreduced_portion_opt_sq, squared_error_unreduced_portion_opt_cos","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9cff153b922bbebc069e45e76a7946af9fb3bdbf"},"cell_type":"markdown","source":"Examples of possible features computed using the function above:"},{"metadata":{"trusted":true,"_uuid":"8712331342186cd3fe5979a039dd87c6e3cdf3c6"},"cell_type":"code","source":"def compute_scale_based_on_class_prototype(class_name, df_grouped, time_scale, class_prototypes, object_ids):\n    class_ts_2 = class_prototypes[class_name]\n    mapping = df_grouped.agg(lambda x: estimate_scale(x.values, class_ts_2))\n    added_columns = [f'x_sc_{class_name}_t_sc_{time_scale:.2f}_cos', f'err_{class_name}_t_sc_{time_scale:.2f}_cos',\n                     f'x_sc_{class_name}_t_sc_{time_scale:.2f}_sq',\n                     f'err_{class_name}_t_sc_{time_scale:.2f}_sq',\n                     f'err_{class_name}_t_sc_{time_scale:.2f}_cos_from_sq',\n                     f'err_{class_name}_t_sc_{time_scale:.2f}_sq_from_cos']\n    df_added_columns = pd.DataFrame(np.column_stack(object_ids.map(lambda x: mapping[x])).T,\n                                    columns=added_columns, index=object_ids.index)\n    return df_added_columns\n\ndf_grouped = df.groupby('object_id')['flux']\n\ntime_scales = [1.0, 0.99] # for illustration\nfor time_scale in time_scales:\n    compute_scale_based_on_class_prototype_current = partial(compute_scale_based_on_class_prototype,\n                                                             time_scale=time_scale,\n                                                             df_grouped=df_grouped,\n                                                             class_prototypes=class_prototypes,\n                                                             object_ids=df_meta['object_id'])\n\n    with Pool(cpu_count() - 1) as p:\n        results_list = p.map(compute_scale_based_on_class_prototype_current, classes)\n    all_new_columns = pd.concat(results_list, axis=1)\n    df_meta = pd.concat((df_meta, all_new_columns), axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5087860650b6e5c1370e3f9cfed3aff863eec10c"},"cell_type":"code","source":"################################################\n# Counting matches of cos/sq best similarities\n################################################\nfor class_name in classes:\n    df_meta[f'cos_sq_matches_{class_name}'] = (df_meta[f'err_{class_name}_t_sc_{time_scales[0]:.2f}_cos']==df_meta[f'err_{class_name}_t_sc_{time_scales[0]:.2f}_cos_from_sq']).map(int)\n    for time_scale in time_scales[1:]:\n        df_meta[f'cos_sq_matches_{class_name}'] += (df_meta[f'err_{class_name}_t_sc_{time_scale:.2f}_cos']==df_meta[f'err_{class_name}_t_sc_{time_scale:.2f}_cos_from_sq']).map(int)\n\n################################################\n# Time scale estimation\n################################################\ndef gather_time_scale_info(suffix):\n    gathered_time_scale_values = dict()\n    for time_scale in time_scales:\n        classes_iterator = iter(classes)\n        class_name = next(classes_iterator)\n        gathered_time_scale_values[time_scale] = np.expand_dims(np.row_stack(df_meta[[f'err_{class_name}_t_sc_{time_scale:.2f}_{suffix}',\n                                                                                      f'x_sc_{class_name}_t_sc_{time_scale:.2f}_{suffix}']].values), axis=1)\n        for class_name in classes_iterator:\n            helping_array = np.expand_dims(np.row_stack(df_meta[[f'err_{class_name}_t_sc_{time_scale:.2f}_{suffix}',\n                                                                 f'x_sc_{class_name}_t_sc_{time_scale:.2f}_{suffix}']].values),\n                                           axis=1)\n            gathered_time_scale_values[time_scale] = np.concatenate(\n                (gathered_time_scale_values[time_scale], helping_array), axis=1)\n    return gathered_time_scale_values, np.concatenate(\n        list(map(lambda x: np.expand_dims(x, 0), gathered_time_scale_values.values())))\n\n# gathered_time_scale_values_cos, axes: time_scale, object, class, [err, x_scale]\ngathered_time_scale_values_cos_dict, gathered_time_scale_values_cos = gather_time_scale_info('cos')\ngathered_time_scale_values_sq_dict, gathered_time_scale_values_sq = gather_time_scale_info('sq')\n\ndef get_trimmed_means_per_time_scale(gathered_time_scale_values, number_of_best_scored_classes_to_take=4):\n    gathered_time_scale_values_sorted_per_classes = np.sort(gathered_time_scale_values, axis=2)\n    trimmed_means_per_time_scale = np.apply_along_axis(np.mean, axis=2,\n                                                       arr=gathered_time_scale_values_sorted_per_classes[:, :, :number_of_best_scored_classes_to_take, 0])\n    return trimmed_means_per_time_scale\n\ndef normalize_scores(object_trimmed_means_per_time_scale):\n    max_score = np.max(object_trimmed_means_per_time_scale)\n    object_trimmed_means_per_time_scale = object_trimmed_means_per_time_scale / max_score\n    return object_trimmed_means_per_time_scale\n\ntrimmed_means_per_time_scale_sq = get_trimmed_means_per_time_scale(gathered_time_scale_values_sq)\ntrimmed_means_per_time_scale_sq_normed = np.apply_along_axis(normalize_scores, axis=0,\n                                                             arr=trimmed_means_per_time_scale_sq)\n\ndef neighb_smooth(trimmed_means_per_time_scale_normed):\n    max_idx = trimmed_means_per_time_scale_normed.shape[0] - 1\n\n    def neighb_mean(i):\n        if i > 0 and i < max_idx:\n            return (0.3 * trimmed_means_per_time_scale_normed[i - 1]\n                    + trimmed_means_per_time_scale_normed[i]\n                    + 0.3 * trimmed_means_per_time_scale_normed[i + 1]) / 1.6\n        elif i == 0:\n            return (trimmed_means_per_time_scale_normed[0] + 0.3 * trimmed_means_per_time_scale_normed[1]) / 1.3\n        else:\n            return (trimmed_means_per_time_scale_normed[-1] + 0.3 * trimmed_means_per_time_scale_normed[-2]) / 1.3\n\n    trimmed_means_per_time_scale_normed_smoothed = np.array(\n        list(map(neighb_mean, np.arange(trimmed_means_per_time_scale_normed.shape[0]))))\n    return trimmed_means_per_time_scale_normed_smoothed\n\ntrimmed_means_per_time_scale_sq_normed_smoothed = np.apply_along_axis(neighb_smooth, axis=0,\n                                                                      arr=trimmed_means_per_time_scale_sq_normed)\ntrimmed_means_per_time_scale_sq_smoothed = np.apply_along_axis(neighb_smooth, axis=0, arr=trimmed_means_per_time_scale_sq)\ndf_meta['estimated_t_sc_idx_orig'] = np.argmin(trimmed_means_per_time_scale_sq_smoothed, axis=0)\ndf_meta['estimated_t_sc'] = df_meta['estimated_t_sc_idx_orig'].map(lambda i: time_scales[i])\n################################################\n# red shift estimation\n################################################\ndf_meta['z_approx'] = df_meta['estimated_t_sc']*(df_meta['hostgal_photoz'] + 1) - 1\n\n################################################\n# Amplitude scale estimation\n################################################\n# time_scale, object, class, [err, x_scale]\ndef get_best10matches_amplitude_scales(gathered_time_scale_values):\n    gathered_err = np.swapaxes(gathered_time_scale_values[:, :, :, 0], 1, 2).reshape(-1,gathered_time_scale_values.shape[1])\n    gathered_scale = np.swapaxes(gathered_time_scale_values[:, :, :, 1], 1, 2).reshape(-1,gathered_time_scale_values.shape[1])\n    gathered_val_pos = np.argsort(gathered_err, axis=0)\n    # https://stackoverflow.com/questions/6155649/sort-a-numpy-array-by-another-array-along-a-particular-axis\n    second_scale_indices = list(range(gathered_scale.shape[1]))\n    gathered_scale_sorted_by_err = gathered_scale[gathered_val_pos, second_scale_indices]\n    gathered_scale_best10 = gathered_scale_sorted_by_err[:10, :]\n    return gathered_scale_best10\n\ngathered_scale_best10_cos = get_best10matches_amplitude_scales(gathered_time_scale_values_cos)\ngathered_scale_best10_sq = get_best10matches_amplitude_scales(gathered_time_scale_values_sq)\ngathered_scale_best20 = np.concatenate((gathered_scale_best10_sq, gathered_scale_best10_cos))\ndf_meta['estimated_x_sc'] = np.median(gathered_scale_best20, axis=0)\n\n################################################\n# Tournament scores\n################################################\ndef err_to_rank_points(err_array):\n    results = np.zeros_like(err_array)\n    if np.all(np.isnan(err_array)):\n        return results\n    pos_sorted = np.argsort(err_array)\n    for i, points in enumerate([4, 3, 2, 1]):\n        results[pos_sorted[i]] += points\n    return results\n\ndef compute_tournament_points(gathered_time_scale_values, suffix):\n    global df_meta\n    gathered_time_scale_err = gathered_time_scale_values[:, :, :, 0] # axes: time_scale, object, class, [err, x_scale]\n    gathered_time_scale_tournament_per_scale_points = np.apply_along_axis(err_to_rank_points, axis=2, arr=gathered_time_scale_err)\n    gathered_time_scale_tournament_total_points = np.sum(gathered_time_scale_tournament_per_scale_points, axis=0)\n\n    columns_to_add = [f'tournament_points_{class_name}_{suffix}' for class_name in classes]\n    df_tournament_points = pd.DataFrame(gathered_time_scale_tournament_total_points,\n                                       columns=columns_to_add, index=df_meta.index)\n    df_meta = pd.concat((df_meta, df_tournament_points), axis=1)\n\ncompute_tournament_points(gathered_time_scale_values_cos, 'cos')\ncompute_tournament_points(gathered_time_scale_values_sq, 'sq')\n\n################################################\n# std of amplitude matching scores\n################################################\ndef std_matching_scores(gathered_time_scale_values, suffix):\n    global df_meta\n    gathered_time_scale_err = gathered_time_scale_values[:, :, :, 1]  # axes: time_scale, object, class\n    gathered_time_scale_err_sorted_per_time_scale = np.sort(gathered_time_scale_err, axis=0)\n    gathered_time_scale_err_2best = gathered_time_scale_err_sorted_per_time_scale[:5, :, :]\n    gathered_time_scale_err_2best_mean = np.std(gathered_time_scale_err_2best, axis=0)\n    # gathered_time_scale_err_2best_sum_min = np.min(gathered_time_scale_err_2best_sum, axis=1)\n    gathered_time_scale_err_2best_sum_normed = gathered_time_scale_err_2best_mean  # gathered_time_scale_err_2best_sum/ (np.ones_like(gathered_time_scale_err_2best_sum) *\n    #    gathered_time_scale_err_2best_sum_min.reshape(-1, 1))\n\n    columns_to_add = [f'best2amp_std_{class_name}_{suffix}' for class_name in classes]\n    df_gathered_time_scale_err_2best_sum = pd.DataFrame(gathered_time_scale_err_2best_sum_normed,\n                                                        columns=columns_to_add, index=df_meta.index)\n    df_meta = pd.concat((df_meta, df_gathered_time_scale_err_2best_sum), axis=1)\n\nstd_matching_scores(gathered_time_scale_values_cos, 'cos')\nstd_matching_scores(gathered_time_scale_values_sq, 'sq')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e3392849431e745b7e358e9eac966dec96ba9c57"},"cell_type":"markdown","source":"I'd be grateful for any potential feedback/corrections/suggestions, comparison to similar results obtained by others. Thank you in advance."},{"metadata":{"trusted":true,"_uuid":"ea0a134d8e2b9a12c69a5ac7a68099017801e124"},"cell_type":"markdown","source":""}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}