{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy  as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport gc\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nfrom sklearn import preprocessing\nimport os\n\nfrom lightgbm import LGBMRegressor\nfrom sklearn.base import BaseEstimator, RegressorMixin, clone\nfrom sklearn.metrics import mean_squared_log_error\n\n#将图表嵌入到notebook中\n%matplotlib inline ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n`feature include:` \n* building_id, meter,site_id, primary_use, square_feet（log1p）,year_built, floor_count,<br>\n* air_temperature, cloud_coverage, dew_temperature, precip_depth_1_hr, sea_level_pressure,wind_direction, wind_speed,<br>   \n* tm_day_of_week, tm_hour_of_day\n\n`find_bad_zeros` identifies rows with zero-readings that have the following characteristics:\n\n* all electrical readings with zero values.\n* 48+ hour runs of steam and hotwater zero-readings *except* for those which are entirely contained within what we identify as typical \"core summer months\".\n* 48+ hour runs of chilledwater zero-readings *except* for those which occur simultaneously at the start and end of the year (i.e. \"core winter months\").\n\n`fixed timestamp` align timestamp by EDA\n* site_GMT_offsets = [-5, 0, -9, -6, -8, 0, -6, -6, -5, -7, -8, -6, 0, -7, -6, -6]\n\n`exclude steam data of 1099 building`\n* return X[(X.building_id==1099)&(X.meter==2)].index\n\n`Model separately by meter`\n* total 4 model of all data\n\n`fill Null data of weather with linear interpolate,relatived with following feature：`\n* air_temperature, cloud_coverage, dew_temperature, precip_depth_1_hr, sea_level_pressure,wind_direction, wind_speed"},{"metadata":{"trusted":true},"cell_type":"code","source":"def input_file(file):\n    path = f\"/kaggle/input/ashrae-energy-prediction/{file}\"\n    if not os.path.exists(path): return path + \".gz\"\n    return path\n## Function to reduce the DF size\ndef compress_dataframe(df):\n    result = df.copy()\n    for col in result.columns:\n        col_data = result[col]\n        dn = col_data.dtype.name\n        if dn == \"object\":\n            result[col] = pd.to_numeric(col_data.astype(\"category\").cat.codes, downcast=\"integer\")\n        elif dn == \"bool\":\n            result[col] = col_data.astype(\"int8\")\n        elif dn.startswith(\"int\") or (col_data.round() == col_data).all():\n            result[col] = pd.to_numeric(col_data, downcast=\"integer\")\n        else:\n            result[col] = pd.to_numeric(col_data, downcast='float')\n    return result","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"这个函数难到极限的python代码，只是理解了实现的功能，没有理解是如何实现的"},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_is_bad_zero(Xy_subset, min_interval=48, summer_start=3000, summer_end=7500):\n    \"\"\"Helper routine for 'find_bad_zeros'.\n    \n    This operates upon a single dataframe produced by 'groupby'. We expect an \n    additional column 'meter_id' which is a duplicate of 'meter' because groupby \n    eliminates the original one.\"\"\"\n    \n    meter = Xy_subset.meter_id.iloc[0]\n    is_zero = Xy_subset.meter_reading == 0\n    #print(is_zero)\n    \n    if meter == 0:\n        # Electrical meters should never be zero. Keep all zero-readings in this table so that\n        # they will all be dropped in the train set.\n        return is_zero\n    \n    transitions = (is_zero != is_zero.shift(1))\n    \n    #print(transitions)\n    all_sequence_ids = transitions.cumsum()\n    \n    \n    #print(all_sequence_ids)\n    ids = all_sequence_ids[is_zero].rename(\"ids\")\n    #print(ids)\n    \n    \n    if meter in [2, 3]:\n        # It's normal for steam and hotwater to be turned off during the summer\n        keep = set(ids[(Xy_subset.timestamp < summer_start) |\n                       (Xy_subset.timestamp > summer_end)].unique())\n        is_bad = ids.isin(keep) & (ids.map(ids.value_counts()) >= min_interval)\n    elif meter == 1:\n        time_ids = ids.to_frame().join(Xy_subset.timestamp).set_index(\"timestamp\").ids\n        is_bad = ids.map(ids.value_counts()) >= min_interval\n\n        # Cold water may be turned off during the winter\n        jan_id = time_ids.get(0, False)\n        dec_id = time_ids.get(8283, False)\n        if (jan_id and dec_id and jan_id == time_ids.get(500, False) and\n                dec_id == time_ids.get(8783, False)):\n            is_bad = is_bad & (~(ids.isin(set([jan_id, dec_id]))))\n    else:\n        raise Exception(f\"Unexpected meter type: {meter}\")\n\n    result = is_zero.copy()\n    result.update(is_bad)\n    return result\n\ndef find_bad_zeros(X,y):\n    \"\"\"Return an Index object containing only the rows which should be deleted.\"\"\"\n    Xy=X.assign(meter_reading=y,meter_id=X.meter)\n    is_bad_zero=Xy.groupby([\"building_id\",\"meter\"]).apply(make_is_bad_zero)\n    return is_bad_zero[is_bad_zero].index.droplevel([0,1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_train():\n    df=pd.read_csv('/kaggle/input/ashrae-energy-prediction/train.csv',parse_dates=[\"timestamp\"])\n    df.timestamp=(df.timestamp-pd.to_datetime(\"2016-01-01\")).dt.total_seconds()//3600\n    return compress_dataframe(df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_test():\n    df=pd.read_csv('/kaggle/input/ashrae-energy-prediction/test.csv',parse_dates=[\"timestamp\"])\n    df.timestamp=(df.timestamp-pd.to_datetime(\"2016-01-01\")).dt.total_seconds()//3600\n    \n    \n    return compress_dataframe(df).set_index(\"row_id\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_building_metadata():\n    return compress_dataframe(pd.read_csv('/kaggle/input/ashrae-energy-prediction/building_metadata.csv').fillna(-1).set_index(\"building_id\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#从EDA 时间偏移分析出来的时间偏移量\nsite_GMT_offsets = [-5, 0, -9, -6, -8, 0, -6, -6, -5, -7, -8, -6, 0, -7, -6, -6]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#这里取消了 add_na_indicators\ndef read_weather_train(fix_timestamps=True,interpolate_na=True):\n    df=pd.read_csv('/kaggle/input/ashrae-energy-prediction/weather_train.csv',parse_dates=[\"timestamp\"])\n    df[\"timestamp\"]=(df.timestamp-pd.to_datetime(\"2016-01-01\")).dt.total_seconds()//3600\n    \n    if fix_timestamps:\n        GMT_offset_map={site:offset for site,offset in enumerate(site_GMT_offsets)}\n        df.timestamp=df.timestamp+df.site_id.map(GMT_offset_map)\n    if interpolate_na:\n        site_dfs=[]\n        \n        for site_id in df.site_id.unique():\n            #训练集2016年有366天\n            site_df=df[df.site_id==site_id].set_index(\"timestamp\").reindex(range(8784))  \n            site_df.site_id=site_id\n            \n            for col in [c for c in site_df.columns if c!=\"site_id\"]:\n                \n                #经过EDA测试，天气的相关缺失值使用线性填补\n                site_df[col]=site_df[col].interpolate(limit_direction='both',method='linear')\n            site_dfs.append(site_df)\n        df=pd.concat(site_dfs).reset_index()\n    return compress_dataframe(df).set_index([\"site_id\",\"timestamp\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#这里取消了 add_na_indicators\ndef read_weather_test(fix_timestamps=True,interpolate_na=True):\n    df=pd.read_csv('/kaggle/input/ashrae-energy-prediction/weather_test.csv',parse_dates=[\"timestamp\"])\n    df[\"timestamp\"]=(df.timestamp-pd.to_datetime(\"2016-01-01\")).dt.total_seconds()//3600\n    \n    if fix_timestamps:\n        GMT_offset_map={site:offset for site,offset in enumerate(site_GMT_offsets)}\n        df.timestamp=df.timestamp+df.site_id.map(GMT_offset_map)\n    \n    if interpolate_na:\n        site_dfs=[]\n        \n        for site_id in df.site_id.unique():\n            #测试集2017年有365天\n            site_df=df[df.site_id==site_id].set_index(\"timestamp\").reindex(range(8760))  \n            site_df.site_id=site_id\n            \n            for col in [c for c in site_df.columns if c!=\"site_id\"]:\n                \n                #经过EDA测试，天气的相关缺失值使用线性填补\n                site_df[col]=site_df[col].interpolate(limit_direction='both',method='linear')\n            site_dfs.append(site_df)\n        df=pd.concat(site_dfs).reset_index()\n    return compress_dataframe(df).set_index([\"site_id\",\"timestamp\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def combined_train_data(fix_timestamps=True,interpolate_na=True,add_na_indicators=True):\n    Xy=compress_dataframe(read_train()\n                        .join(read_building_metadata(),on=\"building_id\")\n                        .join( read_weather_train(fix_timestamps,interpolate_na),on=[\"site_id\",\"timestamp\"])\n                        .fillna(-1))\n    \n    return Xy.drop(columns=[\"meter_reading\"]), Xy.meter_reading","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def _add_time_features(X):\n    return X.assign(tm_day_of_week=(X.timestamp//24)%7,tm_hour_of_day=(X.timestamp%24))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def combined_test_data(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    X = compress_dataframe(read_test().join(read_building_metadata(), on=\"building_id\").join(\n        read_weather_test(fix_timestamps, interpolate_na),\n        on=[\"site_id\", \"timestamp\"]).fillna(-1))\n    return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#根据EDA显示去掉电力中的  building_id  0-104 是site_id==0\ndef find_bad_sitezero(X):\n    return X[(X.timestamp<3378)&(X.site_id==0)&(X.meter==0)].index","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#经过EDA 发现1099号建筑的蒸汽数据有问题，需要去掉\ndef find_bad_building1099(X,y):\n    return X[(X.building_id==1099)&(X.meter==2)].index","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 移除经过EDA发现的异常数据\n#包括四种能源数据的季节性清洗以及大块0数据以及异常高的数据\ndef find_bad_rows(X,y):\n    return find_bad_zeros(X,y).union(find_bad_sitezero(X)).union(find_bad_building1099(X,y))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class CatSplitRegressor(BaseEstimator,RegressorMixin):\n    def __init__(self,model,col):\n        self.model=model\n        self.col=col\n    \n    def fit(self,X,y):\n        self.fitted={}\n        importances=[]\n        \n        for val in X[self.col].unique():\n            X1=X[X[self.col]==val].drop(columns=[self.col])\n            self.fitted[val]=clone(self.model).fit(X1,y.reindex_like(X1))\n            importances.append(self.fitted[val].feature_importances_)\n            del X1\n        \n        fi=np.average(importances,axis=0)\n        col_index=list(X.columns).index(self.col)\n        self.feature_importances_=[*fi[:col_index],0,*fi[col_index:]]\n        return self\n    \n    def predict(self,X):\n        result=np.zeros(len(X))\n        for val in X[self.col].unique():\n            ix=np.nonzero((X[self.col]==val).to_numpy())\n            predictions=self.fitted[val].predict(X.iloc[ix].drop(columns=[self.col]))\n            result[ix]=predictions\n        return result","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"categorical_columns=[\"building_id\",\"meter\",\"site_id\",\"primary_use\",\"tm_day_of_week\",\"tm_hour_of_day\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class LGBMWrapper(BaseEstimator,RegressorMixin):\n    def __init__(self,categorical_feature=None,**params):\n        self.model=LGBMRegressor(**params)\n        self.categorical_feature=categorical_feature\n    def fit(self,X,y):\n        with warnings.catch_warnings():\n            cats = None if self.categorical_feature is None else list(\n                X.columns.intersection(self.categorical_feature)\n            )\n            warnings.filterwarnings(\"ignore\",\"categorical_feature in Dataset is overriden\".lower())\n            self.model.fit(X,y,**({} if cats is None else {\"categorical_feature\":cats}))\n            self.feature_importances_=self.model.feature_importances_\n            \n            return self\n    \n    def predict(self,X):\n        return self.model.predict(X)\n    \n    def get_params(self,deep=True):\n        return {**self.model.get_params(deep),\"categorical_feature\":self.categorical_feature}\n    \n    def set_params(self,**params):\n        ctf=params.pop(\"categorical_feature\",None)\n        if ctf is not None: self.categorical_feature=ctf\n        self.model.set_params(params)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"X,y=combined_train_data()\nbad_rows=find_bad_rows(X,y)\npd.Series(bad_rows.sort_values()).to_csv(\"rows_to_drop.csv\",header=False,index=False)\n\nX=X.drop(index=bad_rows)\ny=y.reindex_like(X)\n\n#Additional preprocessing\n\n#增加了 dayofweek 以及hourofday\nX=compress_dataframe(_add_time_features(X))\n\n#square_feet 平滑处理\nX['square_feet']=np.log1p(X['square_feet'])\n\n#对 primary_use 进行转换\n#labelEncoder = preprocessing.LabelEncoder()\n#X['primary_use'] = labelEncoder.fit_transform(X['primary_use'].astype(str))\n\n#原始的时间戳特征在预测的时候并没有什么作用\nX=X.drop(columns=\"timestamp\")\n\ny=np.log1p(y)\n\nmodel=CatSplitRegressor(\n    LGBMWrapper(random_state=0,n_jobs=-1,categorical_feature=categorical_columns),\"meter\"\n)\n\nmodel.fit(X,y)\ndel X,y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X=combined_test_data()\n\n#增加了 dayofweek 以及hourofday\nX=compress_dataframe(_add_time_features(X))\n\n#原始的时间戳特征在预测的时候并没有什么作用\nX=X.drop(columns=\"timestamp\")\n\n#对 primary_use 进行转换\n#labelEncoder = preprocessing.LabelEncoder()\n#X['primary_use'] = labelEncoder.fit_transform(X['primary_use'].astype(str))\n\n#square_feet 平滑处理\nX['square_feet']=np.log1p(X['square_feet'])\n\npredictions=pd.DataFrame({\n    \"row_id\":X.index,\n    \"meter_reading\":np.clip(np.expm1(model.predict(X)),0,None)\n    \n})\n\ndel X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions.to_csv(\"submission.csv\",index=False,float_format=\"%.4f\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.4"}},"nbformat":4,"nbformat_minor":4}