{"cells":[{"metadata":{},"cell_type":"markdown","source":"\n## Generalized mean\nWe use the [weighted generalized mean](https://en.wikipedia.org/wiki/Generalized_mean) to blend our predictions.  \n\n$$\n\\bar{x}\n= \\left( \\sum_{i=1}^n w_i x_i^p \\right)^{1/p}\n$$\n\nWe tuned the weights and $p$ using the optuna library."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import os\nimport optuna\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom functools import partial\nfrom sklearn.metrics import mean_squared_error","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nroot = '../input/ashrae-feather-format-for-fast-loading'\ntest = pd.read_feather(f'{root}/test.feather')\nmeta = pd.read_feather(f'{root}/building_metadata.feather')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"leak = pd.read_feather('../input/ashrae-leak-data-station/leak.feather')\n\nleak.fillna(0, inplace=True)\nleak = leak[(leak.timestamp.dt.year > 2016) & (leak.timestamp.dt.year < 2019)]\nleak.loc[leak.meter_reading < 0, 'meter_reading'] = 0 # remove negative values\nleak = leak[leak.building_id != 245]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Leak Validation for public kernels(not used leak data)"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_list = [   \n    \"20191214-catboost-no-split-1-100-v2/2019-12-14_catboost-no_split_1-100_v2\",\n    \"aggregate-models-v5/weighted_blend_2019-12-14_gbm_split_primary_use\",\n    \"aggregate-models-v5/weighted_blend_2019-12-10_gbm_with_trend\",\n]\n\nfor i,f in enumerate(submission_list):\n    x = pd.read_csv(f'../input/{f}.csv', index_col=0).meter_reading\n    x[x < 0] = 0\n    test[f'pred{i}'] = x\n\ndel  x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"leak = pd.merge(leak, test[['building_id', 'meter', 'timestamp', *[f\"pred{i}\" for i in range(len(submission_list))], 'row_id']], \"left\")\nleak = pd.merge(leak, meta[['building_id', 'site_id']], 'left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(len(submission_list)):\n    sns.distplot(np.log1p(leak[f\"pred{i}\"]))\n    sns.distplot(np.log1p(leak.meter_reading))\n    leak_score = np.sqrt(mean_squared_error(np.log1p(leak[f\"pred{i}\"]), np.log1p(leak.meter_reading)))\n    print(f'score{i}={leak_score}')    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Leak Validation for Blending"},{"metadata":{"trusted":true},"cell_type":"code","source":"# log1p then mean\nlog1p_then_mean = np.mean(np.log1p(leak[[f\"pred{i}\" for i in range(len(submission_list))]].values), axis=1)\nleak_score = np.sqrt(mean_squared_error(log1p_then_mean, np.log1p(leak.meter_reading)))\nprint('log1p then mean score =', leak_score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# mean then log1p\nmean_then_log1p = np.log1p(np.mean(leak[[f\"pred{i}\" for i in range(len(submission_list))]].values, axis=1))\nleak_score = np.sqrt(mean_squared_error(mean_then_log1p, np.log1p(leak.meter_reading)))\nprint('mean then log1p score=', leak_score)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Tune with Optuna"},{"metadata":{"trusted":true},"cell_type":"code","source":"class GeneralizedMeanBlender():\n    \"\"\"Combines multiple predictions using generalized mean\"\"\"\n    def __init__(self, p_range=(-2,2)):\n        \"\"\"\"\"\"\n        self.p_range = p_range\n        self.p = None\n        self.weights = None\n                \n    def _objective(self, trial, X, y):\n                    \n        # create hyperparameters\n        p = trial.suggest_uniform(f\"p\", *self.p_range)\n        weights = [\n            trial.suggest_uniform(f\"w{i}\", 0, 1)\n            for i in range(X.shape[1])\n        ]\n\n        # blend predictions\n        blend_preds, total_weight = 0, 0\n        if p <= 0:\n            for j,w in enumerate(weights):\n                blend_preds += w*np.log1p(X[:,j])\n                total_weight += w\n            blend_preds = np.expm1(blend_preds/total_weight)\n        else:\n            for j,w in enumerate(weights):\n                blend_preds += w*X[:,j]**p\n                total_weight += w\n            blend_preds = (blend_preds/total_weight)**(1/p)\n            \n        # calculate mean squared error\n        return np.sqrt(mean_squared_error(y, blend_preds))\n\n    def fit(self, X, y, n_trials=10): \n        # optimize objective\n        obj = partial(self._objective, X=X, y=y)\n        study = optuna.create_study()\n        study.optimize(obj, n_trials=n_trials)\n        # extract best weights\n        if self.p is None:\n            self.p = [v for k,v in study.best_params.items() if \"p\" in k][0]\n        self.weights = np.array([v for k,v in study.best_params.items() if \"w\" in k])\n        self.weights /= self.weights.sum()\n\n    def transform(self, X): \n        assert self.weights is not None and self.p is not None,\\\n        \"Must call fit method before transform\"\n        if self.p == 0:\n            return np.expm1(np.dot(np.log1p(X), self.weights))\n        else:\n            return np.dot(X**self.p, self.weights)**(1/self.p)\n    \n    def fit_transform(self, X, y, **kwargs):\n        self.fit(X, y, **kwargs)\n        return self.transform(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = np.log1p(leak[[f\"pred{i}\" for i in range(len(submission_list))]].values)\ny = np.log1p(leak[\"meter_reading\"].values)\n\ngmb = GeneralizedMeanBlender()\ngmb.fit(X, y, n_trials=20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.sqrt(mean_squared_error(gmb.transform(X), np.log1p(leak.meter_reading))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submit"},{"metadata":{"trusted":true},"cell_type":"code","source":"# make test predictions\nsample_submission = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/sample_submission.csv\")\nX_test = test[[f\"pred{i}\" for i in range(len(submission_list))]].values\nsample_submission['meter_reading'] = np.expm1(gmb.transform(np.log1p(X_test)))\nsample_submission.loc[sample_submission.meter_reading < 0, 'meter_reading'] = 0\n\n# fill in leak data\nleak = leak[['meter_reading', 'row_id']].set_index('row_id').dropna()\nsample_submission.loc[leak.index, 'meter_reading'] = leak['meter_reading']\n\n# save submission\nsample_submission.to_csv('submission.csv', index=False, float_format='%.4f')\nsample_submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(np.log1p(sample_submission.meter_reading))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}