{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"In this kernel, I’m sharing my Gaussian Process Feature Engineering method with celerite used on my 19th place solution ( https://www.kaggle.com/c/PLAsTiCC-2018/discussion/75167 ) \n\nGaussian Process is used to interpolate flux curve, then extract features from the light curves. Some of Gaussian Process model tuning refinement from my original code was done with reference @CPMP 's solution as the following.\n\nThanks @CPMP.\nCPMP’s Solution detail: https://www.kaggle.com/c/PLAsTiCC-2018/discussion/75050  \nCPMP’s Github repo:  https://github.com/jfpuget/Kaggle_PLAsTiCC  \n\nRemaining issue:  Some of the flux max can not be captured with current implementation of gaussian process curve.\n"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"%load_ext autoreload\n%autoreload 2\n\n# Basic Library\nimport pandas as pd\nimport pandas.io.sql as psql\nimport numpy as np\nimport numpy.random as rd\nimport gc\nimport multiprocessing as mpa\nimport os\nimport sys\nimport pickle\nfrom collections import defaultdict\nfrom glob import glob\nimport math\nfrom datetime import datetime as dt\nfrom pathlib import Path\nimport scipy.stats as st\nimport re\nfrom scipy.stats.stats import pearsonr\nfrom itertools import combinations\n\n# Matplotlib\nimport matplotlib\nfrom matplotlib import font_manager\nimport matplotlib.pyplot as plt\nimport matplotlib.cm as cm\nfrom matplotlib import rc\n\nfrom matplotlib import animation as ani\nfrom IPython.display import Image\n\nfrom joblib import Parallel, delayed\nimport multiprocessing as mp\n\nplt.rcParams[\"patch.force_edgecolor\"] = True\n#rc('text', usetex=True)\nfrom IPython.display import display # Allows the use of display() for DataFrames\nimport seaborn as sns\nsns.set(style=\"whitegrid\", palette=\"muted\", color_codes=True)\nsns.set_style(\"whitegrid\", {'grid.linestyle': '--'})\nred = sns.xkcd_rgb[\"light red\"]\ngreen = sns.xkcd_rgb[\"medium green\"]\nblue = sns.xkcd_rgb[\"denim blue\"]\n\npd.set_option(\"display.max_colwidth\", 100)\npd.set_option(\"display.max_rows\", None)\npd.set_option(\"display.max_columns\", None)\npd.options.display.float_format = '{:,.5f}'.format\n\n%matplotlib inline\n#%config InlineBackend.figure_format='retina'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"502a5baf536cebcd8c042084e723936abb395012"},"cell_type":"code","source":"# for gaussian process\nfrom scipy.optimize import minimize\nimport celerite\nfrom celerite import terms\nfrom celerite.modeling import Model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"65b5e75c6327568bb2784c67d61ecee8d1aa23f0"},"cell_type":"code","source":"def get_data_for_gp(df_pb, offset = 11):\n    # original code of this part of code is CPMP's following notebook\n    # https://github.com/jfpuget/Kaggle_PLAsTiCC/blob/master/code/celerite_003.ipynb\n    x_min = df_pb.mjd.min()\n    x_max = df_pb.mjd.max()\n\n    yerr_mean = df_pb.flux_err.mean()\n    x = df_pb.mjd.values\n    y = df_pb.flux.values\n    yerr = df_pb.flux_err\n\n    x = np.concatenate((np.linspace(x_min-250, x_min -200, offset), x, np.linspace(x_max+200, x_max+250, offset),))\n    y = np.concatenate((np.random.randn(offset) * yerr_mean, y, np.random.randn(offset) * yerr_mean))\n    yerr = np.concatenate((yerr_mean * np.ones(offset), yerr, yerr_mean * np.ones(offset) ))\n    return x, y, yerr\n\n\noptimizer_list = [\"L-BFGS-B\",\n                  \"Nelder-Mead\",\n                  \"Powell\",\n                  \"CG\",\n                  \"BFGS\",\n                  \"Newton-CG\",\n                  \"TNC\",\n                  \"COBYLA\",\n                  \"SLSQP\",\n                  \"dogleg\",\n                  \"trust-ncg\",]\n\ndef neg_log_like(params, y, gp):\n    gp.set_parameter_vector(params)\n    return -gp.log_likelihood(y)\n\ndef grad_neg_log_like(params, y, gp):\n    gp.set_parameter_vector(params)\n    return -gp.grad_log_likelihood(y)[1]\n\ndef build_gp_model(x, y, yerr, n_param = 2):\n    log_sigma = 0\n    log_rho = 0\n    eps = 0.001\n    bounds = dict(log_sigma=(-15, 15), log_rho=(-15, 15))\n    kernel = terms.Matern32Term(log_sigma=log_sigma, \n                                log_rho=log_rho, \n                                eps=eps, \n                                bounds=bounds)\n\n    gp = celerite.GP(kernel, mean=0)\n    gp.compute(x, yerr) \n\n    initial_params = gp.get_parameter_vector()\n    bounds = gp.get_parameter_bounds()\n\n    # depend on a combination of optimizer and dataset, minimize function throw exception,\n    # so trying all type of method.\n    for opt in optimizer_list:\n        try:\n            r = minimize(neg_log_like, \n                         initial_params, \n                         jac=grad_neg_log_like, \n                         method=opt, #\"L-BFGS-B\", \n                         bounds=bounds, \n                         args=(y, gp))\n            return gp\n\n        except Exception as e:\n            pass\n    raise Exception(\"[build_gp_model] can`t optimize\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"aa66eea8822f5223160b98a3509ae14e9ec29a01"},"cell_type":"code","source":"ZIP = False\nif ZIP:\n    train = pd.read_csv('../input/training_set.csv.zip', compression='zip')\nelse:\n    train = pd.read_csv('../input/training_set.csv')\n    \nmeta_train = pd.read_csv('../input/training_set_metadata.csv')\nmeta_train[\"gal\"] = (meta_train['hostgal_specz'] == 0).astype(int)\nclasses = np.sort(meta_train.target.unique())\n########################################\nn_row = 2 # with changing this variable, you can set the number of graphs for each target showing in this notebook.\n########################################\nimport numpy.random as rd\nrd.seed(71)\ndf_list = []\nfor u in np.sort(meta_train.target.unique()):\n    df_list.append(meta_train[meta_train.target==u].sample(frac=1).iloc[:n_row,:]) # with shuffle\ndf_target_short = pd.concat(df_list)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e8c2aaa820501c54942108bba6563b3d24fe7963"},"cell_type":"markdown","source":"# Visualization\n"},{"metadata":{"trusted":true,"_uuid":"30dd64839eb2997033adf220a5f09f234dda4e95"},"cell_type":"code","source":"def draw_gp_result(object_id_list):\n    for object_id in object_id_list:\n        meta = meta_train[meta_train.object_id==object_id]\n        df = train[train.object_id==object_id]\n\n        for pb in range(6):\n            df_pb = df[df.passband == pb]\n\n            flux_err_mean = df_pb.flux_err.mean()\n            flux_err_std = df_pb.flux_err.std()\n\n            # using flux_err within mean±6std\n            df_pb = df_pb[df_pb.flux_err <= flux_err_mean + 6*flux_err_std]\n\n            x, y, yerr = get_data_for_gp(df_pb)\n            gp = build_gp_model(x, y, yerr)\n\n            # graph drowing\n            colors = [\"r\",\"b\", \"g\", \"k\", \"purple\", \"orange\"]\n            n_xx = 800\n            gp_xx = np.linspace(x.min(), x.max(), n_xx)\n            mu, var = gp.predict(y, gp_xx, return_var=True) \n            # N/A interpolation\n            mu = np.interp(gp_xx, gp_xx[~np.isnan(mu)], mu[~np.isnan(mu)], )\n            \n            # Draw Graph\n            ax = plt.subplot(111)\n            df_pb.plot.scatter(\"mjd\", \"flux\", c=colors[pb], figsize=(18,4), ax=ax) #cmap=cm.rainbow,\n            ax.errorbar(df_pb.mjd, df_pb.flux, df_pb.flux_err, lw=1, fmt=\"none\", ecolor=\"k\", zorder=-1, label=\"flux_err\")\n            plt.plot(gp_xx, mu, color=\"gray\")\n\n            plt.fill_between(gp_xx, mu+np.sqrt(var), mu-np.sqrt(var), color=\"orange\", alpha=0.3, edgecolor=\"none\")\n            plt.title(f\"oid:{object_id}, passband:{pb}, class:{meta.target.values[0]}\")\n            plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dbe5bb0cc124bd5285fb372db7955957ee0db925","scrolled":false},"cell_type":"code","source":"for c in classes:\n    print(\"=\"*30, f\" class_{c} \", \"=\"*30)\n    oid_list = df_target_short[df_target_short.target==c].object_id.values\n    draw_gp_result(oid_list)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ac6db6902fadb7e2bfaba142e624b57468f1b52e"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a13583b8ebfb199ac384f89c71f45ac3036c9710"},"cell_type":"markdown","source":"# Create Features"},{"metadata":{"trusted":true,"_uuid":"7f60227a27e7196dea433bd8ce69124b49658244"},"cell_type":"code","source":"\ndef applyParallel(dfGrouped, func):\n    retLst = Parallel(n_jobs=mp.cpu_count())(delayed(func)(group) for name, group in dfGrouped)\n    return pd.concat(retLst)\n\ndef gp_features(df):\n    prefix = \"gp001:\"\n    feature = {}\n    feature[\"object_id\"] = [df.object_id.values[0]]\n\n    mu_interp_list = []\n    for pb in range(6):\n        df_pb = df[df.passband == pb]\n\n        flux_err_mean = df_pb.flux_err.mean()\n        flux_err_std = df_pb.flux_err.std()\n\n        df_pb = df_pb[df_pb.flux_err <= flux_err_mean + 6*flux_err_std]\n        x, y, yerr = get_data_for_gp(df_pb)\n\n        try:\n            gp = build_gp_model(x, y, yerr)\n\n            n_xx = 800\n            dx = (x.max() - x.min()) / n_xx\n            gp_xx = np.linspace(x.min(), x.max(), n_xx)\n            mu, var = gp.predict(y, gp_xx, return_var=True) \n            mu = np.interp(gp_xx, gp_xx[~np.isnan(mu)], mu[~np.isnan(mu)], )\n            mu_interp_list.append(mu)    \n\n            feature[f\"{prefix}mean_var_pb{pb}\"] = [np.mean(var)]\n            feature[f\"{prefix}skew_pb{pb}\"] = [st.skew(mu)]\n            feature[f\"{prefix}range_pb{pb}\"] = [mu.max() - mu.min()]\n\n            mu_slope = pd.Series(np.diff(mu) / dx).rolling(window=30).mean()\n            feature[f\"{prefix}slope_max_pb{pb}\"] = [mu_slope.max()]\n            feature[f\"{prefix}slope_min_pb{pb}\"] = [mu_slope.min()]\n\n        except Exception as e:\n            print(e)\n            mu_interp_list.append(np.zeros_like(gp_xx))\n\n            feature[f\"{prefix}skew_pb{pb}\"] = [np.nan]\n            feature[f\"{prefix}range_pb{pb}\"] = [np.nan]\n            feature[f\"{prefix}slope_max_pb{pb}\"] = [np.nan]\n            feature[f\"{prefix}slope_min_pb{pb}\"] = [np.nan]\n\n    for c1, c2 in combinations(range(6), 2):\n        ratio = mu_interp_list[c1] / mu_interp_list[c2]\n        ratio = ratio[~np.isnan(ratio)]\n        feature[f\"{prefix}ratio_pb{c1}_{c2}\"] = [np.mean(ratio)]\n\n        feature[f\"{prefix}corr_pb{c1}_{c2}\"] = [pearsonr(mu_interp_list[c1], mu_interp_list[c2])[0]]\n\n    return pd.DataFrame(feature)\n\nimport pickle\ndef unpickle(filename):\n    with open(filename, 'rb') as fo:\n        p = pickle.load(fo)\n    return p\n\ndef to_pickle(filename, obj):\n    with open(filename, 'wb') as f:\n        pickle.dump(obj, f, -1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ad4d90a853cf126224097eeca2eb99cf8991b433"},"cell_type":"code","source":"# Test one object\noid = 23822\ndf = train[train.object_id==oid]\ndf_gp_feat = gp_features(df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"125d4686b71e9a5b02668dc5229f14a09ba23acb"},"cell_type":"code","source":"df_gp_feat","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9885d8ea66a960ec80bca7a1bfb3d059a199760e"},"cell_type":"code","source":"%%time\n# all object ids\ndf_gp_feature = applyParallel(train.groupby(\"object_id\"), gp_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"26f217460f85c1561d6a2fc23c2b73d93fd4415c"},"cell_type":"code","source":"df_gp_feature.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8870dc37637a38ca7a719a5e10611c24ac8471e4"},"cell_type":"markdown","source":"# EDA"},{"metadata":{"trusted":true,"_uuid":"3242648407e54ba4805bae7fa6d2743b41c0ed01"},"cell_type":"code","source":"df_gp_feature = df_gp_feature.merge(meta_train, on=\"object_id\", how=\"left\")\ndf_gp_feature_dropna = df_gp_feature.replace(np.inf, np.nan).replace(-np.inf, np.nan).dropna()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"5f86e3af6fb4ca00d3788dcc03e37926830985f6"},"cell_type":"code","source":"upper_limit = 99999999\ndf_gp_feature_dropna[\"gp001:ratio_pb1_3\"] = np.where(np.abs(df_gp_feature_dropna[\"gp001:ratio_pb1_3\"])>upper_limit, \n                                                     np.sign(df_gp_feature_dropna[\"gp001:ratio_pb1_3\"])*upper_limit, \n                                                     df_gp_feature_dropna[\"gp001:ratio_pb1_3\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"837039d609b9ce2bf94f40e5c39fefc95ec33b2b","scrolled":false},"cell_type":"code","source":"for i, c in enumerate(df_gp_feature_dropna.loc[:,df_gp_feature_dropna.columns.str.startswith(\"gp001:\")]):\n    try:\n        print(c)\n        df_gp_feature_dropna[c] = np.where(np.abs(df_gp_feature_dropna[c]) > 99999999, \n                                                     np.sign(df_gp_feature_dropna[c])*99999999, \n                                                     df_gp_feature_dropna[c])\n        \n        plt.figure(figsize=[20, 4])\n        plt.subplot(1, 2, 1)\n        sns.violinplot(x='target', y=c, data=df_gp_feature_dropna)\n        plt.grid()\n\n        plt.subplot(1, 2, 2)\n        sns.distplot(df_gp_feature_dropna[c], kde=False)\n        plt.yscale('log')\n        plt.legend(['train', 'test'])\n        plt.grid()\n        plt.show();\n    except Exception as e:\n        print(e)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"db59ff4f80cbd4c9e30897c84f1598c5b91df52c"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}