{"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\nimport lightgbm as lgb\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, is_holiday<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 2 model of all data by month%2==0\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\n\n`is holiday`\n* holidays = \n            [\"2016-01-01\", \"2016-01-18\", \"2016-02-15\", \"2016-05-30\", \"2016-07-04\",\n            \"2016-09-05\", \"2016-10-10\", \"2016-11-11\", \"2016-11-24\", \"2016-12-26\",\n            \"2017-01-02\", \"2017-01-16\", \"2017-02-20\", \"2017-05-29\", \"2017-07-04\",\n            \"2017-09-04\", \"2017-10-09\", \"2017-11-10\", \"2017-11-23\", \"2017-12-25\",\n            \"2018-01-01\", \"2018-01-15\", \"2018-02-19\", \"2018-05-28\", \"2018-07-04\",\n            \"2018-09-03\", \"2018-10-08\", \"2018-11-12\", \"2018-11-22\", \"2018-12-25\",\n            \"2019-01-01\"]"},{"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        \n        if col==\"original_timestamp\":\n            continue\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['original_timestamp']=df.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['original_timestamp']=df.timestamp\n    df.timestamp=(df.timestamp-pd.to_datetime(\"2016-01-01\")).dt.total_seconds()//3600\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,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),tm_month=(X.original_timestamp.dt.month),\n                   is_holiday=((X.timestamp.isin(holidays)).astype(int)))","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":"#美国的重大节假日 每年10个\nholidays = [\"2016-01-01\", \"2016-01-18\", \"2016-02-15\", \"2016-05-30\", \"2016-07-04\",\n            \"2016-09-05\", \"2016-10-10\", \"2016-11-11\", \"2016-11-24\", \"2016-12-26\",\n            \"2017-01-02\", \"2017-01-16\", \"2017-02-20\", \"2017-05-29\", \"2017-07-04\",\n            \"2017-09-04\", \"2017-10-09\", \"2017-11-10\", \"2017-11-23\", \"2017-12-25\",\n            \"2018-01-01\", \"2018-01-15\", \"2018-02-19\", \"2018-05-28\", \"2018-07-04\",\n            \"2018-09-03\", \"2018-10-08\", \"2018-11-12\", \"2018-11-22\", \"2018-12-25\",\n            \"2019-01-01\"]\n    ","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 以及monthofday用来切分数据集\nX=compress_dataframe(_add_time_features(X))\n\n#square_feet 平滑处理\nX['square_feet']=np.log1p(X['square_feet'])\n\n#删掉冗余的特征\ndrop_features=[\"timestamp\",\"sea_level_pressure\",\"wind_direction\",\"wind_speed\",\"original_timestamp\"]\nX.drop(drop_features, axis=1, inplace=True) #0 index,1 columns\n\n#根据奇偶月份将数据集划分成2份\nX_half_1=X[X['tm_month']%2==0]\nX_half_2=X[X['tm_month']%2==1]\n\ny_half_1=np.log1p(X_half_1['meter_reading'])\ny_half_2=np.log1p(X_half_2['meter_reading'])\n\n#month会造成不太好的扰动，因此要去掉\nX_half_1.drop([\"meter_reading\",\"tm_month\"],axis=1,inplace=True)\nX_half_2.drop([\"meter_reading\",\"tm_month\"],axis=1,inplace=True)\n\ncategorical_features = [\"building_id\", \"site_id\", \"meter\", \"primary_use\", \"tm_hour_of_day\", \"tm_day_of_week\",\"is_holiday\"]\nX_half_1.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"d_half_1=lgb.Dataset(X_half_1,label=y_half_1,categorical_feature=categorical_features, free_raw_data=False)\nd_half_2=lgb.Dataset(X_half_2,label=y_half_2,categorical_feature=categorical_features, free_raw_data=False)\n\nwatchlist_1 = [d_half_1, d_half_2]\nwatchlist_2 = [d_half_2, d_half_1]\n\nparams = {\n    \"objective\": \"regression\",\n    \"boosting\": \"gbdt\",\n    \"num_leaves\": 45,\n    \"learning_rate\": 0.05,\n    \"feature_fraction\": 0.85,\n    \"reg_lambda\": 2,\n    \"metric\": \"rmse\"\n}\n\nprint(\"Building model with first half and validating on second half:\")\nmodel_half_1 = lgb.train(params, train_set=d_half_1, num_boost_round=1000, valid_sets=watchlist_1, verbose_eval=200, early_stopping_rounds=200)\n\nprint(\"Building model with second half and validating on first half:\")\nmodel_half_2 = lgb.train(params, train_set=d_half_2, num_boost_round=1000, valid_sets=watchlist_2, verbose_eval=200, early_stopping_rounds=200)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del  X_half_1, X_half_2, y_half_1, y_half_2, d_half_1, d_half_2, watchlist_1, watchlist_2\ngc.collect()","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\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#删掉冗余的特征\ndrop_features=[\"timestamp\",\"sea_level_pressure\",\"wind_direction\",\"wind_speed\",\"original_timestamp\",\"tm_month\"]\nX.drop(drop_features, axis=1, inplace=True) #0 index,1 columns\n\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred = np.expm1(model_half_1.predict(X, num_iteration=model_half_1.best_iteration)) / 2\n\ndel model_half_1\ngc.collect()\n\npred += np.expm1(model_half_2.predict(X, num_iteration=model_half_2.best_iteration)) / 2\n    \ndel model_half_2\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame({\"row_id\": row_ids, \"meter_reading\": np.clip(pred, 0, a_max=None)})\nsubmission.to_csv(\"submission2.csv\", index=False)","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}