import lightgbm as lgb
import numpy as np
import pandas as pd

n=100000
lam = .01
X = np.floor(np.random.lognormal(size=(n,2))).astype(int)
y = np.maximum(X[:,0],X[:,1])+np.random.poisson(lam=lam, size=n)

train_inds = np.arange(int(n/3))
val_inds =   np.arange(int(n/3), int(2*n/3))
test_inds =   np.arange(int(2*n/3), int(n))

X_test, y_test = X[test_inds,:], y[test_inds]

ds = lgb.Dataset(X,y, categorical_feature=[1])

ds_train = ds.subset(train_inds)
ds_val = ds.subset(val_inds)

params = {'objective':'poisson',
         'metric':'rmse',
         'learning_rate':.1
         }
gbm = lgb.train(params, ds_train, num_boost_round=300, early_stopping_rounds=20, valid_sets=[ds_val, ds_train], 
                verbose_eval=100, categorical_feature=[1])

yhat = gbm.predict(X_test)
print('neg obs:', len(yhat[yhat<0]))

# 字段说明Instant 记录号

Dteday：日期
Season：季节
    1=春天
    2=夏天
    3=秋天
    4=冬天
yr：年份，(0: 2011, 1:2012)
mnth：月份( 1 to 12)
hr：小时 (0 to 23) （只在 hour.csv 有，作业忽略此字段）
holiday：是否是节假日
weekday：星期中的哪天，取值为 0～6
workingday：是否工作日
    1=工作日 （非周末和节假日）
    0=周末
weathersit：天气
1：晴天，多云
2：雾天，阴天
3：小雪，小雨
4：大雨，大雪，大雾
temp：气温摄氏度
atemp：体感温度
hum：湿度
windspeed：风速

y值
-   casual：非注册用户个数
-  registered：注册用户个数
-  cnt：给定日期（天）时间（每小时）总租车人数，响应变量 y

### EXPAND TRAIN DATA ###

X_train['count_so_far'] = X_train.apply(lambda x: np.arange(x.casual), axis=1)
X_train['stop'] = X_train.apply(lambda x: np.append(np.zeros(x.casual-1), 1), axis=1)
X_train = X_train.apply(pd.Series.explode)   # 数据扩充的方式

X_train['count_so_far'] = X_train.count_so_far.astype(int)
X_train['stop'] = X_train.stop.astype(int)

y_train = X_train['stop']
X_train = X_train[columns+['count_so_far']]

print(X_train.shape)
X_train

### EXPAND VALID DATA FOR PREDICTION ###

X_val_surv = df[(df.dteday >= datetime(year=2012, month=1, day=1))&(df.dteday < datetime(year=2012, month=5, day=1))].copy()

X_val_surv['count_so_far'] = X_val_surv.apply(lambda x: np.arange(max_count), axis=1)
X_val_surv['stop'] = X_val_surv.apply(lambda x:np.append(np.zeros(x.casual), np.ones(max_count-x.casual)), axis=1)
X_val_surv = X_val_surv.apply(pd.Series.explode)

X_val_surv['count_so_far'] = X_val_surv.count_so_far.astype(int)
X_val_surv['stop'] = X_val_surv.stop.astype(int)

y_val_surv = X_val_surv['stop']
X_val_surv = X_val_surv[columns+['count_so_far']]

print(X_val_surv.shape)

params = {
         'objective':'poisson',
         'num_leaves':30, 
         'learning_rate': 0.001,
         'feature_fraction': 0.8,
         'bagging_fraction': 0.9,
         'bagging_seed': 33,
         'poisson_max_delta_step': 0.8,
         'metric': 'poisson'
         }

### FIT LGBM WITH POISSON LOSS ### 

trn_data = lgb.Dataset(X_train, label=y_train)
val_data = lgb.Dataset(X_val, label=y_val)

model = lgb.train(params, trn_data, num_boost_round=1000,
                  valid_sets = [trn_data, val_data],
                  verbose_eval=50, early_stopping_rounds=150)

### PREDICT HAZARD FUNCTION ON VALIDATION DATA AND TRANSFORM TO SURVIVAL ###

p_val_hz = model.predict(X_val_surv).reshape(-1,max_count)
p_val = 1-np.exp(-np.cumsum(p_val_hz, axis=1))
X_val_surv.shape
>>> (484000, 12)
p_val.shape
>>> (121, 4000)

### CRPS ON VALIDATION DATA ###
crps(t_val, p_val)
>>> 0.17425293919055515
### CRPS ON VALIDATION DATA WITH BASELINE MODEL ###
crps(t_val, np.repeat(cdf, len(t_val)).reshape(-1,max_count))