{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6896fb40",
   "metadata": {},
   "source": [
    "# 실습 1 · 시계열 — 서울 지하철 일별 승차 예측과 역 군집\n",
    "\n",
    "**한 줄 목표** 시간순 분할과 윈도우를 직접 만들어, 딥러닝 시계열 모델이 계절 naive 를 *언제* 이기는지 숫자로 확인한다.\n",
    "\n",
    "템플릿 구조: ① 데이터 로드 → ② 분할 → ③ 베이스라인 → ④ 모델 → ⑤ 평가(seed 반복) → ⑥ 군집·시각화 → ⑦ 보고표 → ⑧ 자기 데이터로 바꾸기.\n",
    "자기 데이터에 쓰려면 **①만 바꾸면 된다** (날짜 열 + 값 열 하나).\n",
    "\n",
    "강의 페이지: https://hufs-ai-lecture.pages.dev/practice/01_timeseries.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6961a48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ① 설치와 데이터 로드 — 서울 지하철 일별 총승차 (서울 열린데이터광장 OA-12914 를 일 단위로 합산한 파일)\n",
    "!pip -q install holidays\n",
    "import pandas as pd, numpy as np, torch, torch.nn as nn, time, holidays\n",
    "import matplotlib.pyplot as plt\n",
    "URL='https://hufs-ai-lecture.pages.dev/practice/data/subway_daily_total.csv'\n",
    "df=pd.read_csv(URL,parse_dates=['date']).sort_values('date').reset_index(drop=True)\n",
    "y=df.board.values/1e6                       # 백만 명 단위\n",
    "kr=holidays.KR(years=range(2015,2027))\n",
    "hol=np.array([d in kr for d in df.date.dt.date]).astype(float)\n",
    "dow=df.date.dt.dayofweek.values\n",
    "print(df.shape, df.date.min().date(), df.date.max().date())\n",
    "df.plot(x='date',y='board',figsize=(11,3),legend=False); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c3e8c5d",
   "metadata": {},
   "source": [
    "## ② 분할 — 시간순으로만, 3개 fold\n",
    "평가 구간을 6개월씩 세 번 두고, 각 구간 **이전** 데이터로만 학습한다(rolling-origin). 무작위 분할을 쓰면 미래가 학습에 새어 들어온다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb4044ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "H=28; L=56                                  # 예측 지평 28일, 입력 창 56일\n",
    "folds=[('2025-03-01','2025-08-31'),('2025-09-01','2026-02-28'),('2026-03-01','2026-08-31')]\n",
    "idx=lambda d:int(np.searchsorted(df.date.values,np.datetime64(d)))\n",
    "def cal_feats(t):                            # 요일 원핫 7 + 공휴일 1\n",
    "    f=np.zeros((len(t),8)); f[np.arange(len(t)),dow[t]]=1; f[:,7]=hol[t]; return f\n",
    "def make_windows(start,end):\n",
    "    X,Y,C,O=[],[],[],[]\n",
    "    for o in range(max(start,L),end-H+1):\n",
    "        X.append(y[o-L:o]); Y.append(y[o:o+H]); C.append(cal_feats(np.arange(o,o+H)).ravel()); O.append(o)\n",
    "    return np.array(X),np.array(Y),np.array(C),np.array(O)\n",
    "def test_set(s,e):\n",
    "    Ote=np.arange(s,e-H+1,H)                  # 28일 간격 비중첩 origin\n",
    "    return (np.array([y[o-L:o] for o in Ote]),np.array([y[o:o+H] for o in Ote]),\n",
    "            np.array([cal_feats(np.arange(o,o+H)).ravel() for o in Ote]),Ote)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "345db970",
   "metadata": {},
   "source": [
    "## ③ 베이스라인 — 계절 naive\n",
    "7일 전 값을 그대로 반복한다. MASE 는 이 베이스라인의 학습구간 MAE 로 나눈 값이라 **1 보다 작아야 의미가 있다**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8616e694",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mase(yt,yp,scale): return np.mean(np.abs(yt-yp))/scale\n",
    "rows=[]\n",
    "for fi,(a,b) in enumerate(folds):\n",
    "    s,e=idx(a),idx(b)+1; Xte,Yte,Cte,Ote=test_set(s,e)\n",
    "    scale=np.mean(np.abs(y[7:s]-y[:s-7]))\n",
    "    p_sn=np.array([np.tile(y[o-7:o],H//7) for o in Ote])\n",
    "    rows.append(dict(fold=fi+1,model='SeasonalNaive7',seed=None,MAE=np.mean(np.abs(Yte-p_sn)),MASE=mase(Yte,p_sn,scale)))\n",
    "pd.DataFrame(rows).round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16d9ce6f",
   "metadata": {},
   "source": [
    "## ④ 모델 — lag-MLP · GRU · DLinear\n",
    "셋 다 입력은 같다: 과거 56일 값 + 앞으로 28일의 요일·공휴일. 출력은 28일 벡터(direct multi-output). 손실은 L1, 조기종료는 학습구간 끝 10%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c77af4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "dev='cuda' if torch.cuda.is_available() else 'cpu'\n",
    "class MLP(nn.Module):\n",
    "    def __init__(s,i=L+H*8,h=256): super().__init__(); s.n=nn.Sequential(nn.Linear(i,h),nn.ReLU(),nn.Dropout(.1),nn.Linear(h,h),nn.ReLU(),nn.Linear(h,H))\n",
    "    def forward(s,x,c,cin=None): return s.n(torch.cat([x,c],1))\n",
    "class GRUm(nn.Module):\n",
    "    def __init__(s,h=64): super().__init__(); s.g=nn.GRU(1+8,h,num_layers=2,batch_first=True,dropout=.1); s.o=nn.Linear(h+H*8,H)\n",
    "    def forward(s,x,c,cin): _,hN=s.g(torch.cat([x.unsqueeze(-1),cin],-1)); return s.o(torch.cat([hN[-1],c],1))\n",
    "class DLinear(nn.Module):                     # 추세(이동평균)와 잔차를 각각 선형으로\n",
    "    def __init__(s,k=7): super().__init__(); s.k=k; s.lt=nn.Linear(L,H); s.ls=nn.Linear(L,H); s.lc=nn.Linear(H*8,H)\n",
    "    def forward(s,x,c,cin=None):\n",
    "        pad=torch.cat([x[:,:1].repeat(1,s.k//2),x,x[:,-1:].repeat(1,s.k//2)],1); tr=pad.unfold(1,s.k,1).mean(-1)\n",
    "        return s.lt(tr)+s.ls(x-tr)+s.lc(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e27840fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_predict(name,Xtr,Ytr,Ctr,Otr,Xte,Cte,Ote,seed,use_cal=True):\n",
    "    torch.manual_seed(seed); np.random.seed(seed)\n",
    "    mu,sd=Xtr.mean(),Xtr.std(); T=lambda a:torch.tensor(a,dtype=torch.float32,device=dev)\n",
    "    xt,yt,ct=T((Xtr-mu)/sd),T((Ytr-mu)/sd),T(Ctr if use_cal else 0*Ctr); xe,ce=T((Xte-mu)/sd),T(Cte if use_cal else 0*Cte)\n",
    "    cin_tr=T(np.stack([cal_feats(np.arange(o-L,o)) for o in Otr])*use_cal); cin_te=T(np.stack([cal_feats(np.arange(o-L,o)) for o in Ote])*use_cal)\n",
    "    m={'MLP':MLP,'GRU':GRUm,'DLinear':DLinear}[name]().to(dev); opt=torch.optim.AdamW(m.parameters(),lr=2e-3,weight_decay=1e-4)\n",
    "    n=len(xt); nv=max(int(n*.1),H); tr,va=np.arange(n-nv),np.arange(n-nv,n); best=(1e9,None); wait=0\n",
    "    for ep in range(300):\n",
    "        m.train(); perm=np.random.permutation(tr)\n",
    "        for i in range(0,len(perm),128):\n",
    "            b=perm[i:i+128]; loss=nn.functional.l1_loss(m(xt[b],ct[b],cin_tr[b]),yt[b]); opt.zero_grad(); loss.backward(); opt.step()\n",
    "        m.eval()\n",
    "        with torch.no_grad(): vl=nn.functional.l1_loss(m(xt[va],ct[va],cin_tr[va]),yt[va]).item()\n",
    "        if vl<best[0]-1e-4: best=(vl,{k:v.clone() for k,v in m.state_dict().items()}); wait=0\n",
    "        elif (wait:=wait+1)>=20: break\n",
    "    m.load_state_dict(best[1]); m.eval()\n",
    "    with torch.no_grad(): return m(xe,ce,cin_te).cpu().numpy()*sd+mu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "813f8060",
   "metadata": {},
   "source": [
    "## ⑤ 평가 — seed 3회 × 3 fold (Colab T4 약 5분)\n",
    "관찰 포인트: (1) 세 모델이 계절 naive(MASE≈0.9)를 이기는가. (2) 세 모델 사이 차이가 seed 표준편차보다 큰가. (3) `use_cal=False` 로 바꾸면 얼마나 잃는가 — 이기는 이유가 순환 구조인지 공휴일 특징인지 여기서 갈린다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "984f3f67",
   "metadata": {},
   "outputs": [],
   "source": [
    "seeds=[0,1,2]\n",
    "for fi,(a,b) in enumerate(folds):\n",
    "    s,e=idx(a),idx(b)+1; Xtr,Ytr,Ctr,Otr=make_windows(0,s); Xte,Yte,Cte,Ote=test_set(s,e); scale=np.mean(np.abs(y[7:s]-y[:s-7]))\n",
    "    for nm in ['MLP','GRU','DLinear']:\n",
    "        for use_cal in [True,False]:\n",
    "            for sd_ in seeds:\n",
    "                p=fit_predict(nm,Xtr,Ytr,Ctr,Otr,Xte,Cte,Ote,sd_,use_cal)\n",
    "                rows.append(dict(fold=fi+1,model=nm+('' if use_cal else '(-cal)'),seed=sd_,MAE=np.mean(np.abs(Yte-p)),MASE=mase(Yte,p,scale)))\n",
    "    print('fold',fi+1,'done')\n",
    "R=pd.DataFrame(rows)\n",
    "tbl=R.groupby('model').agg(MAE=('MAE','mean'),MASE=('MASE','mean'),MASE_sd=('MASE','std')).round(3).sort_values('MASE'); tbl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f1a108a",
   "metadata": {},
   "source": [
    "## ⑥ 군집·시각화 — 역별 이용 패턴 → UMAP → K-means\n",
    "역마다 2019년 요일별 승차 비중(7) + 규모 + 승/하차 비 + 코로나 하락(2020/2019) + 회복(2024/2019) 벡터를 만들고 군집한다. **군집은 원 벡터에서, UMAP 은 그림에만** 쓴다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e107b56d",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip -q install umap-learn\n",
    "import umap\n",
    "from sklearn.preprocessing import StandardScaler; from sklearn.cluster import KMeans; from sklearn.metrics import silhouette_score, adjusted_rand_score\n",
    "st=pd.read_parquet('https://hufs-ai-lecture.pages.dev/practice/data/subway_station_daily.parquet')\n",
    "st['key']=st.line+'·'+st.station; st['dow']=st.date.dt.dayofweek; st['year']=st.date.dt.year\n",
    "keep=st.groupby('key').date.nunique().pipe(lambda c:c[c>=4000].index); d=st[st.key.isin(keep)]; y19=d[d.year==2019]\n",
    "piv=y19.pivot_table(index='key',columns='dow',values='board',aggfunc='mean'); feat=piv.div(piv.sum(1),axis=0); feat.columns=[f'dow{i}' for i in range(7)]\n",
    "g=y19.groupby('key'); feat['board_alight']=g.board.sum()/g.alight.sum(); feat['log_size']=np.log(g.board.mean())\n",
    "ym=lambda yr:d[d.year==yr].groupby('key').board.mean(); feat['covid_drop']=ym(2020)/ym(2019); feat['recover_24']=ym(2024)/ym(2019); feat=feat.dropna()\n",
    "X=StandardScaler().fit_transform(feat.values)\n",
    "print({k:round(silhouette_score(X,KMeans(k,n_init=10,random_state=0).fit_predict(X)),3) for k in range(2,7)})\n",
    "K=3; km=KMeans(K,n_init=10,random_state=0).fit(X); feat['cluster']=km.labels_\n",
    "print('seed 안정성 ARI',[round(adjusted_rand_score(km.labels_,KMeans(K,n_init=10,random_state=s).fit_predict(X)),3) for s in range(1,6)])\n",
    "red=umap.UMAP(n_neighbors=15,min_dist=0.1,random_state=0).fit_transform(X)\n",
    "plt.figure(figsize=(5,4)); plt.scatter(red[:,0],red[:,1],c=feat.cluster,s=8,cmap='viridis'); plt.title('역 이용 패턴 UMAP (색=군집)'); plt.show()\n",
    "feat.groupby('cluster').agg(n=('log_size','size'),weekend=('dow5',lambda v:(v+feat.loc[v.index,'dow6']).mean()),covid_drop=('covid_drop','mean'),recover_24=('recover_24','mean')).round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d1f516a",
   "metadata": {},
   "source": [
    "## ⑦ 보고표 — 논문 본문 표 형식으로\n",
    "평균±표준편차, seed 수, 분할 방식, 베이스라인을 한 표에 담는다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d890df7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "rep=R.groupby('model').agg(MAE=('MAE','mean'),MAE_sd=('MAE','std'),MASE=('MASE','mean'),MASE_sd=('MASE','std'),n=('MASE','size')).reset_index()\n",
    "rep['MAE (백만 명)']=rep.apply(lambda r:f\"{r.MAE:.3f} ± {0 if np.isnan(r.MAE_sd) else r.MAE_sd:.3f}\",axis=1)\n",
    "rep['MASE']=rep.apply(lambda r:f\"{r.MASE:.3f} ± {0 if np.isnan(r.MASE_sd) else r.MASE_sd:.3f}\",axis=1)\n",
    "print(rep[['model','MAE (백만 명)','MASE','n']].to_markdown(index=False))\n",
    "print('\\n분할: rolling-origin 3 fold(6개월), 지평 28일, 입력 56일, 손실 L1, 조기종료, seed',seeds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "299015d2",
   "metadata": {},
   "source": [
    "## ⑧ 자기 데이터로 바꾸기\n",
    "①에서 `df` 를 `date`(날짜), `board`(값) 두 열만 있는 데이터프레임으로 바꾸면 나머지 셀은 그대로 돈다. 확인할 것 세 가지:\n",
    "1. 결측 날짜가 없는가 (`pd.date_range` 로 채우기)\n",
    "2. 계절 주기가 7일이 맞는가 — 월별 데이터면 `H`, `L`, naive 의 7 을 12 로\n",
    "3. 평가 구간(`folds`)이 실제 관심 시점(개입 전후 등)과 겹치지 않는가"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}