よっしーの私的空間

機械学習を中心に興味のあることについて更新します

SHAPを使用した回帰問題の機械学習モデルの局所解釈方法

shapライブラリを使用して、回帰問題を解いた機械学習モデルの大局的解釈を行う。

1.SHAPとは

SHAP(SHapley Additive exPlanations)は、機械学習モデルの局所的な解釈可能性を提供するためのフレームワークです。SHAPは、個々の特徴量が予測にどのように寄与しているかを計算することにより、モデルの解釈性を高めます。

2.データセット

今回はsklearnのdiabetes(糖尿病)データセットを例に実装してみる。データセットの概要は以下の通り。

  • 目的変数:1年後の糖尿病の進行に関する測定値
  • 説明変数:11個の特徴量
特徴量 説明
age 患者の年齢
sex 患者の性別
bmi 患者のBMI
bp 患者の平均血圧
S1 T-Cells (白血球の1種)の活性化に関与する化学物質
S2 インスリンのレセプターの活性に関与する化学物質
S3 グルコース代謝に関与する化学物質
S4 トリグリセリド (脂質の一種)の代謝に関与する化学物質
S5 血清アルブミン (タンパク質の一種)のレベル
S6 銅の輸送に関与する化学物質

※S1~S6: 血清の6つの化学指標。

3.実装

3.1.Lightgbmを用いた簡易的なモデルを作成

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
import lightgbm as lgb
import pandas as pd
import shap

# データのロード
diabetes = load_diabetes()
x_train, x_test, y_train, y_test = train_test_split(diabetes.data, diabetes.target, test_size=0.3)

# LightGBMモデルのトレーニング
model = lgb.LGBMRegressor().fit(x_train, y_train)

3.2.SHAPの実装

# SHAP値の計算
df_x_test = pd.DataFrame(x_test, columns=diabetes.feature_names)
explainer = shap.Explainer(model.predict, df_x_test)
shap_values = explainer(df_x_test)

# 特定のサンプルのSHAP値を表示(shap_values[x]:x番目のデータのSHAP値を表示、max_display=n:n個の特徴量を描画)
shap.plots.waterfall(shap_values[0], max_display=10)

上記の結果、以下のような図が表示されます。

E[f(x)]は、モデルの期待出力値(期待値)を表しています。つまり、全体的な平均予測値です。 SHAP値は、個々の特徴量がこの期待値からどれだけ上下に影響を与えたかを表します。

PDP/ICEを使用した回帰問題の機械学習モデルの大局的解釈方法

pdpboxライブラリを使用して、回帰問題を解いた機械学習モデルの大局的解釈を行う。

1.データセット

今回はsklearnのボストン住宅価格データセットを例に実装してみる。データセットの概要は以下の通り。

  • 目的変数:ボストンの地域別住宅価格
  • 説明変数:以下13個の特徴量1
特徴の名前 説明
CRIM 人口 1 人当たりの犯罪発生数
ZN 25,000 平方フィート以上の住居区画の占める割合=「広い家の割合」
INDUS 小売業以外の商業が占める面積の割合
CHAS チャールズ川に関わるダミー変数 (1: 川の周辺, 0: それ以外)
NOX 一酸化窒素の濃度
RM 住居の平均部屋数
AGE 1940 年より前に建てられた物件の割合 =「古い家の割合」
DIS ボストン市の5 つの雇用施設からの距離 (重み付け済) =「主要施設への距離」
RAD 高速道路へのアクセスのしやすさ
TAX $10,000 ドルあたりの固定資産税率
PTRATIO 町毎の生徒と教師の比率
B 町毎の黒人 (Bk) の比率
LSTAT 低所得者人口の割合

2.実装

2.1.Lightgbmを用いた簡易的なモデルを作成

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
import pandas as pd
import lightgbm as lgb

boston = load_boston()
x_train, x_test, y_train, y_test = train_test_split(boston.data, boston.target, test_size=0.3)

#LGB用のデータに変形
lgb_train = lgb.Dataset(x_train, y_train)
lgb_eval = lgb.Dataset(x_test, y_test)

params = {
    'objective': 'regression',
    'boosting_type': 'gbdt',
    'metric': 'mae'
    }

model = lgb.train(
    params = params, 
    train_set = lgb_train,
    valid_sets= [lgb_eval,lgb_train],
)

2.2.PDP/ICEの実装

# pdpboxはpandas.DataFrameしか受け付けない
test_df = pd.DataFrame(x_test, columns=boston.feature_names)

# 特徴量毎にpdp/iceを描画
for col in train_df.columns:
    pdp_boston = pdp.pdp_isolate(
        model=model, dataset=test_df, model_features=test_df.columns, feature=col
    )
    fig, axes = pdp.pdp_plot(pdp_boston, col, plot_lines=True, frac_to_plot=1.0, plot_pts_dist=True)

特徴量毎に以下のような図が描画される。


上図より、一酸化窒素の濃度が0.7までは住宅価格にあまり影響がなく、0.7以上あたりから住宅価格を下げる要因と機械学習モデルが判断していることが分かる。

2.3.PDP(2変数)の実装

ついでに2変数の相互作用を等高線として図示する方法もまとめておく。

import itertools
columns = test_df.columns
combinations = itertools.combinations(columns, 2)

for combination in combinations:
    interaction = pdp.pdp_interact(model=model, dataset=test_df, model_features=test_df.columns, features=combination)
    fig, axes = pdp.pdp_interact_plot(pdp_interact_out=interaction, feature_names=combination)

特徴量の組み合わせ毎に以下のような図が描画される。


上図より、固定資産税率が350を超えたあたりから住宅価格にあまり影響がなく、犯罪発生率の方が影響していることが分かる。

Pandasでグループ毎に過去データを集約する方法

PythonのPandasを用いて、任意のカラム毎にグルーピングして、過去データを集約する方法についてまとめます。調べても意外と情報が少なくて苦戦しました。見つけても、複数カラムでグルーピングできなかったり、自信のデータを除外できなかったりと。これらの点も対応しているので、良ければ参考にしてください。

例が無いと分かりにくいので、以下のデータについて考えていきます。生徒毎に過去2年分の平均点数を算出する方法をまとめます。

結論ですが、以下で実施可能です。

import pandas as pd
# サンプルデータを作成する。
df = pd.DataFrame(
    [["A君","2020年","60"],
     ["B君","2020年","86"],
     ["C君","2020年","56"],
     ["A君","2021年","90"],
     ["B君","2021年","72"],
     ["C君","2021年","70"],
     ["A君","2022年","95"],
     ["B君","2022年","85"],
     ["C君","2022年","80"]],
    columns=["生徒","試験日","点数"]
)

# 集計する事前準備として、データを生徒毎に試験日で並び替える
df = df.groupby(['生徒']).apply(lambda x: x.sort_values(by=['試験日'], ascending=[True])).reset_index(drop=True)

# 過去2年分データを集計して、新カラム作成
df['過去2年の平均点数'] = pd.Series(
    df.groupby(['生徒'])['点数']\
    .apply(lambda v: v.shift(+1).rolling(2,min_periods=1).mean())
)

display(df.sort_values(by=['試験日','生徒']))

結果は以下の通りになります。想定通りですね。

もう少し解説していこうと思います。まずは並び替え部分について。

続いて、集計部分について。

競馬のレーティングをしてみた(Elo Rating)

前回はGlicko2 Rating Systemを使用してレーティングを行いましたが、今回はElo Ratingを使用してみました。今回は競走馬だけではなくジョッキーやトレイナーについてもレーティングを行っています。Elo Ratingは通常チェス等の1対1の対戦において使用されるレーティングシステムですが、今回はそれを複数人対戦用に拡張したライブラリを使用してみました。英語ですが解説はこちらです。

1.レーティング結果

2015年以降の全レース結果を対象にイロレーティングを適用しました。結果は以下の通りでした。上位5位までを掲載します。6位以下についても気になる人がいるようであれば開示します。
f:id:t-yoshi-book:20220301005110p:plain

2.ソースコード

from multielo import MultiElo
import pickle
import pandas as pd
pd.set_option('display.max_columns', 100)
import numpy as np
import datetime
from tqdm.notebook import tqdm as tqdm

# 事前に取得した競馬レースデータをロードする。競馬データの取得は以下を参考にした。
# http://houdoukyokucho.com/2020/08/31/post-1617/
with open('df.pickle', 'rb') as f:
    df = pickle.load(f)

# ジョッキーの名前から減量記号を削除
df['jockey_name'] = df['jockey_name'].str.replace("▲","").replace("△","").replace("☆","").replace("★","").replace("◇","")
# 中止レース等を除外
df = df[(df["rank"]!="中止") & (df["rank"]!="除外") & (df["rank"]!="取消") & (df["rank"]!="失格")]
# レートを格納する列を追加
df["horse_rate"]=np.nan
df["jockey_rate"]=np.nan
df["trainer_rate"]=np.nan

# レーティング結果を格納するデータフレームを作成
def createRatingTable(df, colName):
    df = df[[colName]].copy()
    df = df.drop_duplicates()
    df["rating"] = 1500
    return df

df_horse_rate = createRatingTable(df, "horse_name")
df_jockey_rate = createRatingTable(df, "jockey_name")
df_trainer_rate = createRatingTable(df, "trainer")

# レーティング処理を高速化するためにPandasをリストに変換
horse_rates = df_horse_rate.values.tolist()
jockey_rates = df_jockey_rate.values.tolist()
trainer_rates = df_trainer_rate.values.tolist()
race_datas = df.values.tolist()

# 繰り返し計算用に辞書定義
list_dict = {
    "horse_name":horse_rates,
    "jockey_name":jockey_rates,
    "trainer":trainer_rates
}
rating_col_dict = {
    "horse_name":"horse_rate",
    "jockey_name":"jockey_rate",
    "trainer":"trainer_rate"
}

# レース単位に繰り返し計算を行う
elo = MultiElo(k_value=32, d_value=400, score_function_base=2)
finished = []

for i,race_data in enumerate(tqdm(race_datas)):
    target_race_id = race_data[12]
    if target_race_id not in finished:
        finished.append(target_race_id)
        single_race = [rd for rd in race_datas if rd[12]==target_race_id]
        ranks = [int(horse[15]) for horse in single_race]
        
        # 馬・ジョッキー・トレイナーのレートを計算
        for rating_target in ["horse_name","jockey_name","trainer"]:
            # 馬・ジョッキー・トレイナーの名前のリストを作成
            name_col_no = df.columns.get_loc(rating_target)
            names = [single_record[name_col_no] for single_record in single_race]
            # 最新レートを格納するリストを取得
            rating_list=list_dict[rating_target]
            # レートを格納する列番号を取得
            rate_col_no = df.columns.get_loc(rating_col_dict[rating_target])
            
            # 将来的に機械学習等を行うときの説明変数とするために、レース前のレートを記録しておく
            players=[]
            for j, name in enumerate(names):
                # 過去レートを取得
                rate = [rate for rate in rating_list if rate[0]==name][0][1]
                players.append(rate)
                # レースデータのリストに登録
                race_datas[i+j][rate_col_no] = rate
            
            # レートを計算する
            new_player_rates  = elo.get_new_ratings(players, ranks)
            
            # 最新のレート情報を記録
            for name, new_player_rate in zip(names,new_player_rates):
                for j,rate in enumerate(rating_list):
                    if rate[0] == name:
                        rating_list[j][1]=new_player_rate
                        break

# 見やすいようにDataFrameに戻す
def createDataFrame(rated_list, original_df):
    columns=original_df.columns
    df = pd.DataFrame(rated_list, columns=columns)
    return df

df = createDataFrame(race_datas, df)
df_horse_rate = createDataFrame(horse_rates, df_horse_rate)
df_jockey_rate = createDataFrame(jockey_rates, df_jockey_rate)
df_trainer_rate = createDataFrame(trainer_rates, df_trainer_rate)

# 結果表示
display(df_horse_rate.sort_values("rating",ascending=False).head(5))
display(df_jockey_rate.sort_values("rating",ascending=False).head(5))
display(df_trainer_rate.sort_values("rating",ascending=False).head(5))

競走馬をレーティングしてみた(Glicko2 Rating System)

競走馬の強さを評価するために競走馬のレーティングをしてみました。レーティングはGlicko2 Rating Systemというのを使用してみました。Glicko2の概要やPythonでの実装方法については以下の記事にまとめたので、良ければ見てください。
book-read-yoshi.hatenablog.com

レーティングは2015年以降のレース結果を基に実施しました。結果は以下の通りです。上位10頭まで表示してます。感覚と一致しますか??
f:id:t-yoshi-book:20220226223440p:plain

レーティング用に作成したPythonコードは以下の通りです。レースデータの取得はこちらを参考にしました。以下ではレースデータの取得方法は割愛してます。

from glicko2 import glicko2
import pickle
import pandas as pd
pd.set_option('display.max_columns', 100)
import numpy as np
from tqdm.notebook import tqdm as tqdm

# レーティングを行う関数
def calcRatings(players,ranks):
    newPlayers=[]
    for i, (target_player,target_rank) in enumerate(zip(players,ranks)):
        new_target_player = copy.deepcopy(target_player)
        ratings=[]
        rds=[]
        outcomes=[]
        for j, (player, rank) in enumerate(zip(players,ranks)):
            if not i==j:
                ratings.append(player.rating)
                rds.append(player.rd)
                if rank>target_rank:
                    outcomes.append(1)
                elif rank<target_rank:
                    outcomes.append(0)
                elif rank==target_rank:
                    outcomes.append(0.5)
        
        new_target_player.update_player(ratings, rds, outcomes)
        newPlayers.append(new_target_player)
    
    return newPlayers

# 事前に取得した競馬レースデータをロードする。競馬データの取得は以下を参考にした。
# http://houdoukyokucho.com/2020/08/31/post-1617/
with open('df_with_race_data.pickle', 'rb') as f:
    df = pickle.load(f)

# 中止レース等を除外
df = df[df["rank"]!="中止"]
df = df[df["rank"]!="除外"]
df = df[df["rank"]!="取消"]
df = df[df["rank"]!="失格"]

# レート等を格納する列を作成
df["horse_rate"]=np.nan
df["horse_rd"]=np.nan
df["horse_vol"]=np.nan

# レーティング結果を格納するデータフレームを作成
def createRatingTable(df, colName):
    df = df[[colName]].copy()
    df = df.drop_duplicates()
    df["rating"] = 1500
    df["rd"] = 350
    df["vol"] = 0.06
    return df

df_horse_rate = createRatingTable(df, "horse_name")

# レーティング処理を高速化するためにPandasをリストに変換
horse_rates = df_horse_rate.values.tolist()   # 各競走馬の最新のレート情報等を登録するリスト
race_datas = df.values.tolist()    # レースデータを蓄積したリスト

# レース単位に繰り返し計算を行う
finished = []
for i,race_data in enumerate(tqdm(race_datas)):
    target_race_id = race_data[12]  # [12]は列番号なのでデータによって異なる
    if target_race_id not in finished:
        finished.append(target_race_id)
        single_race = [race_data for race_datain race_datas if race_data[12]==target_race_id]
        ranks = [int(horse[15]) for horse in single_race]    # [15]は列番号なのでデータによって異なる
        
        horse_names  = [single_record[3] for single_record in single_race]   # [3]は列番号なのでデータによって異なる

        # 将来的に機械学習等を行うときの説明変数とするために、レース前のレートを記録しておく
        horses = []
        for j, horse_name in enumerate(horse_names):
            # 過去レートを取得
            horse_rate,horse_rd,horse_vol = [horse_rate for horse_rate in horse_rates if horse_rate[0]==horse_name][0][1:4]
            player = glicko2.Player(horse_rate,horse_rd,horse_vol)
            horses.append(player)
            # レースデータのリストに登録 ※以下の[25:28]は列番号なので作成データによって異なる
            race_datas[i+j][25:28] = [horse_rate,horse_rd,horse_vol]

        # 競走馬のレートを計算する
        new_horses  = calcRatings(horses, ranks)

        # 各競走馬の最新のレート情報を記録
        for horse_name, new_horse in zip(horse_names,new_horses):
            for j,horse_rate in enumerate(horse_rates):
                if horse_rate[0] == horse_name:
                    horse_rates[j][1:4]=[new_horse.rating, new_horse.rd, new_horse.vol]

# レートを高い順に表示する
df_rate = pd.DataFrame(horse_rates,columns=["horse_name","rating","rd","vol"])
df_rate.sort_values("rating",ascending=False)

多人数対戦におけるレーティング方法(glicko2 rating systemを用いて)

競走馬の強さをレーティングをしたくて、レーティングシステムについて調べたところ、glicko-2 rating systemというのが良さそうでした。本記事ではこれの直感的な理解やPythonによる実装方法についてまとめます。

Glicko-2 rating systemの計算式について知りたい場合は、以下の記事がよくまとまっています。以下の記事に書いてあったのですが、glicko-2 rating systemというのはsplatoon2スプラトゥーン2)でも使用されているようです。
レーティングについてと,グリコ2レーティング(Glicko-2 System)におけるレーティング算出方法 - 機械学習、たまにゲーム

なお、本記事は執筆にあたり、以下を参考にしています。

1.glicko-2 rating systemの直感的な理解

Glicko-2のパラメータは以下3つで構成されていて、対戦毎にこれらのパラメータが更新されます。

  • rating(レーティング):r(初期値:1500)
  • rating deviation(レーティング偏差):RD(初期値:350)
  • rating volatility(レーティング変動率):σ(初期値は適用対象による。論文では0.06と設定。)

ratingは解説するまでもないですが、高ければ高いほど強いことを表し、初期値は1500です。対戦で強い相手に勝つほどratingスコアが高くなります。逆に負けるとスコアが下がります。

rating deviationはratingスコアの振れ幅を表します。ratingは単一の値で表現されますが、これが本当に真の強さを表現できているとは限りません。そこで、スコアの振れ幅を設定して真の強さが収まる範囲を設定します。論文では、ratingが1850、rating deviationが50のとき、真の強さの95%信頼区間は1750~1950と言えると記載されていました。これはつまり、ratingスコアが正規分布に従うと仮定していて、その幅を以下の通り定義しているものと思われます。
 r-1.96RD ≦ 真の強さ ≦ r+1.96RD

最後にrating volatilityについてですが、これは強さが安定しているかを表す指標です。恒常的に同じパフォーマンスを発揮している場合はこの値が低くなり、安定していない場合はこの値が高くなるようです。

2.glicko-2 rating systemの実装

実装に当たり、deepy/glicko2を使用しました。論文執筆者Mark Glickmanの実装(ryankirkman/pyglicko2)も試してみたのですが、稀にエラーが発生するので、今回はdeepy/glicko2を使用しています。
さっそくですが、実装したコードは以下の通りです。

# ライブラリのインポート
from glicko2 import glicko2
import copy

# プレイヤーのパラメータを設定する関数
## 変数rating, rd, volはそれぞれrating、rating deviation、rating volatilityの初期値
## playerオブジェクトを1つ返す
def setPlayer(rating=1500, rd=350, vol=0.06):
    player = glicko2.Player()
    player.rating=rating
    player.rd=rd
    player.vol=vol
    return player

# 複数プレイヤーのrating等パラメータを計算する関数
## 変数playersは複数のplayerオブジェクトを格納したリスト
## 変数ranksは複数のplayerの対戦時における順位を格納したリスト
## 各playerのrating等を更新してplayersリストを返す
def calcRatings(players,ranks):
    newPlayers=[]
    for i, (target_player,target_rank) in enumerate(zip(players,ranks)):
        new_target_player = copy.deepcopy(target_player)
        ratings=[]
        rds=[]
        outcomes=[]
        for j, (player, rank) in enumerate(zip(players,ranks)):
            if not i==j:
                ratings.append(player.rating)
                rds.append(player.rd)
                if rank>target_rank:
                    outcomes.append(1)
                elif rank<target_rank:
                    outcomes.append(0)
                elif rank==target_rank:
                    outcomes.append(0.5)
        
        new_target_player.update_player(ratings, rds, outcomes)
        newPlayers.append(new_target_player)
    
    return newPlayers

使い方は以下の通り。

  1. setPlayer関数でplayerオブジェクトを生成して、それらをリストとして格納する。
  2. 各playerの順位をリストとして格納する。
  3. 上記1,2のリストをcalcRatings関数に代入する。
  4. 各playerのrating, rating deviation, rating volatilityが更新されたリストが返ってくる。

上記の手順を試しに実装すると以下の通り。

player1 = setPlayer(1400,30,0.06)
player2 = setPlayer(1550,100,0.06)
player3 = setPlayer(1700,300,0.06)
player4 = setPlayer(1500,200,0.06)

players=[player1,player2,player3,player4]
ranks=[4,2,1,3]

players=calcRatings(players,ranks)

上記を実行した際の結果は以下の通り各プレイヤーのrating等が更新される。
f:id:t-yoshi-book:20220213224138p:plain

Tensorflow-EfficientnetでGradCAMを実装してみた

前々から気になってたGrad-CAMを実装してみました。実装に当たって以下を参考としてます。
EfficientNet-Keras-GradCam-Visualization/inference_example.ipynb at master · lvisdd/EfficientNet-Keras-GradCam-Visualization · GitHub

上記は訓練済みのEfficientnetモデルをそのまま使っているのですが、この記事ではFinetuneもしてみました。Efficientnetは以下を使用しています。
GitHub - qubvel/efficientnet: Implementation of EfficientNet model. Keras and TensorFlow Keras.

1.個人的に引っかかったポイント

私はよく以下のようにEfficientnetのモデルを作成するのですが、以下の方法だとGrad-CAMで最終レイヤーを取得できず失敗するようです。

from efficientnet.tfkeras import EfficientNetB0

model = Sequential()
model.add(EfficientNetB0(
    include_top=False,
    weights='imagenet',
    input_shape=input_shape))
model.add(GlobalAveragePooling2D())
model.add(Dense(num_classes, activation="softmax"))

上記でモデルを作成して、model.summary()を実行すると以下の通りとなります。efficientnet-b0レイヤーとしてまとまってしまうのが問題なんだと思います

代わりに以下の通りモデルを作成することでうまくいきました。

from efficientnet.tfkeras import EfficientNetB0

eff = EfficientNetB0(
    include_top=False,
    weights='imagenet',
    input_shape=input_shape)
x = tf.keras.layers.GlobalAveragePooling2D()(eff.output)
output = tf.keras.layers.Dense(num_classes, activation='softmax', name='last_output')(x)
model = tf.keras.Model(inputs=eff.inputs, outputs=output, name='model')

model.summary()の実行結果は以下の通りです。長いので途中省略してますが、要は上記でモデルを作成すると、Efficientnet内部のレイヤーが残るようです。

2.Grad-CAMの実装

先述の通り、Grad-CAMの実装に当たってはこちらを参考としていますが、参考サイトではKerasを使用しており、Tensorflow.Kerasを使用している場合、エラーとなるので、多少修正しました。最終版は以下の通りです。tensorflow==2.3.0では問題なく動きました。

# grad cam関数
def grad_cam(input_model, image, cls, layer_name):
    """GradCAM method for visualizing input saliency."""
    y_c = input_model.output[0, cls]
    conv_output = input_model.get_layer(layer_name).output
    g = tf.Graph()
    with g.as_default():
        grads = tf.gradients(y_c, conv_output)[0]
    # Normalize if necessary
    # grads = normalize(grads)
    gradient_function = K.function([input_model.input], [conv_output, grads])

    output, grads_val = gradient_function([image])
    output, grads_val = output[0, :], grads_val[0, :, :, :]

    weights = np.mean(grads_val, axis=(0, 1))
    cam = np.dot(output, weights)

    # Process CAM
    cam = cv2.resize(cam, (image_size, image_size), cv2.INTER_LINEAR)
    cam = np.maximum(cam, 0)
    cam = cam / cam.max()
    return cam

3.CIFAR10で試してみた

EfficientnetのモデルをCIFAR10でFinetuneし、Grad-CAMを適用してみました。結果は以下の通りです。無事、注視範囲の可視化に成功しました。目と耳に反応しているんですかね?

参考:ソースコード全文

データセット:CIFAR10
使用モデル:Efficientnet-v1-b0

#ライブラリインポート
import tensorflow as tf
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, GlobalAveragePooling2D, Flatten

from tensorflow.keras.utils import to_categorical
from tensorflow.keras.preprocessing.image import ImageDataGenerator

from efficientnet.tfkeras import EfficientNetB0
from tensorflow.keras import optimizers
from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping

from sklearn.model_selection import train_test_split
from tensorflow.keras.datasets import cifar10

import datetime

from tensorflow.keras import backend as K
import cv2

# 画像をリサイズ(拡大)する関数を定義
def upscale(image):
    size = len(image)
    data_upscaled = np.zeros((size, RESIZE_TO, RESIZE_TO, 3,))
    for i in range(len(image)):
        data_upscaled[i] = cv2.resize(image[i], dsize=(RESIZE_TO, RESIZE_TO), interpolation=cv2.INTER_CUBIC)
    image = np.array(data_upscaled, dtype='uint8')
    
    return image

# データインポート
(x_train, y_train), (x_test, y_test) = cifar10.load_data()

# 画像リサイズ
IMAGE_SIZE  = 32
RESIZE_TO   = 160
x_train = upscale(x_train)
x_test  = upscale(x_test)

# 変数定義
image_size = RESIZE_TO
input_shape=(image_size,image_size,3)

# 予測するクラスの数
num_classes = 10

# モデル定義
def build_model():
    eff = EfficientNetB0(
        include_top=False,
        weights='imagenet',
        input_shape=input_shape)
    x = tf.keras.layers.GlobalAveragePooling2D()(eff.output)
    output = tf.keras.layers.Dense(num_classes, activation='softmax', name='last_output')(x)
    model = tf.keras.Model(inputs=eff.inputs, outputs=output, name='model')

    model.compile(optimizer=optimizers.Adam(learning_rate=1e-4), loss="categorical_crossentropy", metrics=["accuracy"])
    
    return model

# 訓練用関数の定義
def train_model(X, y, epochs, batch_size, callbacks):
    # 訓練データと評価データの分割
    X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, stratify=y, shuffle=True)
    y_train = to_categorical(y_train)
    y_valid = to_categorical(y_valid)
    
    # Data Augumentation
    datagen = ImageDataGenerator(
        rotation_range=20,
        horizontal_flip=True,
        width_shift_range=0.2,
        zoom_range=0.2,
    )    
        
    #モデル構築
    model = build_model()

    # 学習
    model.fit(
        datagen.flow(X_train, y_train, batch_size=batch_size),
        validation_data=(X_valid, y_valid),
        steps_per_epoch=len(X_train) / batch_size,
        epochs=epochs,
        callbacks=callbacks,
        shuffle=True
    )
    
    return model

# 訓練実施
epochs = 3
batch_size = 32

reduce_lr = ReduceLROnPlateau(
    monitor='val_accuracy',
    factor=0.2,
    patience=1,
    verbose=1,
    min_delta=1e-4,
    min_lr=1e-6,
    mode='max'
)
earlystopping = EarlyStopping(
    monitor='val_accuracy',
    min_delta=1e-4,
    patience=2,
    mode='max',
    verbose=1
)

callbacks = [earlystopping, reduce_lr]

print('開始時間:',datetime.datetime.now())
model = train_model(x_train, y_train, epochs, batch_size, callbacks)
print('終了時間:',datetime.datetime.now())

# 予測
preds = []
X = x_test
pred = model.predict(X)

# 予測結果の確認
df_pred = pd.DataFrame(pred)
pred = np.array(df_pred.idxmax(axis=1))
df_pred = pd.DataFrame(pred)
df_y = pd.DataFrame(y_test)
df_result = pd.concat([df_y, df_pred], axis=1, join_axes=[df_y.index])
df_result.columns = ['y','pred']

from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix
print('Confusion Matrix:')
print(confusion_matrix(df_result['y'],df_result['pred']))
print()
print('Accuracy :{:.4f}'.format(accuracy_score(df_result['y'],df_result['pred'])))
print('Precision:{:.4f}'.format(precision_score(df_result['y'],df_result['pred'],average='macro')))
print('Recall   :{:.4f}'.format(recall_score(df_result['y'],df_result['pred'],average='macro')))
print('F_score  :{:.4f}'.format(f1_score(df_result['y'],df_result['pred'],average='macro')))

# grad cam関数
def grad_cam(input_model, image, cls, layer_name):
    """GradCAM method for visualizing input saliency."""
    y_c = input_model.output[0, cls]
    conv_output = input_model.get_layer(layer_name).output
    g = tf.Graph()
    with g.as_default():
        grads = tf.gradients(y_c, conv_output)[0]
    # Normalize if necessary
    # grads = normalize(grads)
    gradient_function = K.function([input_model.input], [conv_output, grads])

    output, grads_val = gradient_function([image])
    output, grads_val = output[0, :], grads_val[0, :, :, :]

    weights = np.mean(grads_val, axis=(0, 1))
    cam = np.dot(output, weights)

    # Process CAM
    cam = cv2.resize(cam, (image_size, image_size), cv2.INTER_LINEAR)
    cam = np.maximum(cam, 0)
    cam = cam / cam.max()
    return cam

# modelの確認
model.summary()

# cifar10 class labels
labels = {
    0:"airplane",
    1:"automobile",
    2:"bird",
    3:"cat",
    4:"deer",
    5:"dog",
    6:"frog",
    7:"horse",
    8:"ship",
    9:"truck"
}

# grad cam
n = 16 #n番目の画像に対してGradCAM適用する
preprocessed_input = x_test[n:n+1].copy()
predictions = model.predict(preprocessed_input)
cls = np.argmax(predictions)
layer_name='top_activation'

gradcam = grad_cam(model, preprocessed_input, cls, layer_name)

# visualise grad cam
print("TRUE : {}".format(labels[y_test[n,0]]))
print("PRED : {}".format(labels[cls]))
for i in range(num_classes):
    print("  {:<10s} : {:5.2f}%".format(labels[i],predictions[0,i]*100))

plt.figure(figsize=(5, 5))
plt.subplot(121)
plt.title('GradCAM')
plt.imshow(preprocessed_input[0])
plt.imshow(gradcam, cmap='jet', alpha=0.5)

plt.subplot(122)
plt.title('Original')
plt.imshow(preprocessed_input[0])

plt.show()