「サイゼリヤで1000円あれば最大何kcal摂れるのか」をSMTソルバー(Z3)で解いてみた。

前書き

• サイゼリヤ1000円ガチャをつくってみた(Heroku + Flask + LINEbot)
• 「サイゼリヤで1000円あれば最大何kcal摂れるのか」を量子アニーリング計算(Wildqat)で解いてみた。

完全に二番煎じですが、古典コンピューターが好きなので、個人的に古典コンピューター最強のなんだかよく分からないけどよく分からないものをよく分からないうちに解いてくれるソフト、z3を使ってサイゼリア問題を解いてみました。

問題

サイゼリヤのメニューを重複無しで合計1000円以下になるように選んだときに、最大の総カロリーになるようなメニューの組み合わせを求めよ。

サイゼリヤのメニューは https://github.com/marushosummers/Saizeriya_1000yen こちらを使わせて使わせて頂きました。メニューは100種類ぐらいみたいで、カロリーは整数値で、プロコンとは違ってメニューは現実の食べ物なのでカロリーは一品あたり多くても1000kcalぐらいみたいです。

解き方

ナップザック問題で、カロリーのレンジも価格のレンジも狭いので、動的計画法を使えば効率よく解けますが、動的計画法とか言われても難しいので、何も考えずにSMTソルバーに解いてもらいます。

SMTソルバーとは

Satisfiability modulo theories (SMT) とは、ブール式の充足可能性問題（SAT）のブールのところを整数やら実数やらリストやらいろいろなものを取れるようにして、命題論理から述語論理に一般化したもののようです（厳密な定義はよく知りませんすみません）。SMTソルバーとはそういう問題を解いてくれるソルバーです。例えば $x^2+x-1=0$ などと入力してやると、$x=-1.6180339887498948...$ と解いてくれます。もっと変数が多くても充足可能な解を1つ見つけてくれます。解がないときは解がないと言ってくれます。特定の制約の下で特定の式を最大にするような問題も解いてくれます。ブール式に限ればかなり変数が多くても解いてくれますし、整数や実数を含む場合も結構解いてくれますが、変数が多いと現実的な時間で解けないこともあります。

本質的に指数的に時間がかかるものを解かせているので万能というわけではありませんが、とりあえずサイゼリヤのメニューのように100個ぐらいの組み合わせとか、それぐらいの規模であれば問題なく解いてくれます。普通にナイーブな探索だと現実的な時間で終わらないと思うので、ドメイン特化の枝狩りをやらずとも勝手に枝狩りしてくれるイメージでしょうか。

実装

SMTソルバーとして、おそらく現時点で最高性能の Z3 を使います。

フロントエンドしては、Haskellの sbv を使います。z3にはPythonのバインディングとかもあるみたいですが、個人的にHaskellのほうが面白いのでHaskellを使います。先ほどの例だと、replに適当に入力すると、

ghci> sat $ \x -> x^2 + x - 1 .== (0 :: SReal)
Satisfiable. Model:
  s0 = root(1, x^2+x = 1) = -1.6180339887498948... :: Real

こうやって解いてくれます。

ではサイゼリアを解くコードを実装していきます。まずはメニューのデータを定義します。メニューのDBファイルからcsvに変換したものを読み込ませるためにカラムに合わせて型を定義してあります。実際使うのはidと名前と値段とカロリーです。

data Menu = Menu
    { menuId :: Int
    , menuName :: String
    , menuCategory :: String
    , menuType :: String
    , menuPrice :: Int
    , menuCalorie :: Int
    , menuSalt :: Double
    } deriving (Generic, Show)

instance FromRecord Menu

次に、最適化させる式を定義します。メニューのリストから、各メニューについて購入するかどうかのブール変数を作って、その変数の割り当てが条件（今回の問題だと合計金額が1000円以下かどうか）を満たすかチェックして、総カロリーを最大化するような条件を書きます。

maxCalorie :: [Menu] -> Symbolic ()
maxCalorie ms = do
    ss <- sBools (map (show . menuId) ms)

Symbolic ()というのは、ソルバーに食わせられる型の一つであり、モナドになっていて、ここで生成した変数や制約式を満たすものを解いてもらえます。sBoolsは、[String] -> Symbolic [SBool] の型を持つ関数で、引数に与えられた名前のリストでブール変数をそれぞれ作ります。要するにここでは、メニューのリストのメニューIDを名前とするブール変数をそれぞれのメニューに対して作っています。

    let totalPrice :: SInt32 =
            sum [ ite s (fromIntegral $ menuPrice menu) 0 | (menu, s) <- zip ms ss]

総額を計算する式です。ite というのはSBool -> a -> a -> aな関数で、要するにSBool版のifです。Trueが割り当てられているメニューだけ価格を合計しています。fromIntegralが挟まっているのは、IntからSInt32に型を合わせる必要があるためです。

    let totalCalorie :: SInt32 =
            sum [ ite s (fromIntegral $ menuCalorie menu) 0 | (menu, s) <- zip ms ss]

総カロリーの計算です。金額の部分とほとんど同じです。

最後に制約と最大化対象を書きます。

    constrain $ totalPrice .<= 1000
    maximize "total-calorie" totalCalorie

constrainはSBoolを引数にとって、それが満たされる必要があるという制約を追加します。maximizeは与えられた式を最大化するという条件を追加します。

問題の方が書き終わったので、ドライバーの部分を書いていきます。

main :: IO ()
    ...
    LexicographicResult result <- optimize Lexicographic $ maxCalorie menus

optimizeに先ほど定義したmaxCalorieを渡すと、充足結果を返してくれます。optimizeの第一引数は最大化する条件が複数ある場合にどういう最大値を求めるのかを指定するもので、Lexicographicは辞書順最大を指定するオプションですが、今回は「総カロリー」1つしかないので、何でも良いです。

この結果をそのままprintしてもよいのですが、何を選んだのかが把握しづらいので、わかりやすい結果を表示するために値を取り出します。

    let ss = [ menuId m
             | m <- menus
             , let Just s = getModelValue (show $ menuId m) result
             , s ]
    let selectedMenus =
             [ m | s <- ss, let Just m = find (\m -> menuId m == s) menus]

getModelValue という関数で特定の名前の変数の割り当てを取得できます。解が求まっているならNothingではないのでJustでバインドしています。割り当てがTrueのものだけフィルターして、選んだメニューのIDのリストを作っています。それから、そのIDをルックアップしてMenuのリストにしています。

    let totalPrice = sum $ map menuPrice selectedMenus
    let totalCalorie = sum $ map menuCalorie selectedMenus
    printf "* Total: ￥%d %dkcal\n" (length bests + 1) totalPrice totalCalorie
    forM_ selectedMenus $ \menu -> do
            printf "    * %s ￥%d %dkcal\n" (menuName menu) (menuPrice menu) (menuCalorie menu)

総額と総カロリー、そして選んだメニューを表示して、完成です。

$ stack run
...
* Total: ￥992 1940kcal
    * ポテトのグリル ￥199 366kcal
    * アーリオ・オーリオ(Wサイズ) ￥574 1120kcal
    * ラージライス ￥219 454kcal

出ました！

量子アニーリングで求められた結果と同じになりました。こちらはSATソルバーが厳密解を返しているはずなので、おそらくたぶん100％正しい答えです。

2位以降

せっかくなので2位以降も求めてみます。
どうやるかはいろいろあるかもしれませんが、とりあえず「1位の組み合わせを除くものの中での最大」を求めれば、それはおそらく2位ということになるはずです。3位以降も同様に、「1～N-1位までの組み合わせ」を除いたものから最大を求めれば求まりそうです。

これをナイーブに書くとなんとなく式が大きくなりそうな気がしたので、「選んだものを包含する」ものを排除する、ということにしました。問題の性質上、包含するものはより金額が大きくなるはずだし、それがより望ましい解なのであれば、カロリーが増えてより上位に来ているはずなので、これは正しいはずです。

ということで、maxCalorie関数に、これまで選んだ組み合わせを引数に加えます。

maxCalorie :: [Menu] -> [[Int]] -> Symbolic ()
maxCalorie ms doNotUses = do
   ...

そして、その組み合わせを包含するような割り当てはダメだという制約条件を追加します。

    forM_ doNotUses $ \doNotUse ->
        constrain $ bnot $ bAnd [ s | (m, s) <- zip ms ss, menuId m `elem` doNotUse ]

これで、引数に1～N-1位までの組み合わせを与えてやることで、N位組み合わせが求まるはずです。

実験結果

さて、これでサイゼリヤで1000円あれば最高何kcal採れるのかランキングの厳密解が求められます。どうなるでしょうか！結果発表行ってみましょう！

第1位

イメージ図：

• ポテトのグリル 199円 366kcal
• アーリオ・オーリオ(Wサイズ) 574円 1120kcal
• ラージライス 219円 454kcal
• 合計 992円 1940kcal

というわけで、第1位は、先駆者様と同じになりました。炭水化物の力強さを感じさせます。

第2位

イメージ図：

• フォッカチオ 119円 214kcal
• セットプチフォッカ 79円 107kcal
• アーリオ・オーリオ(Wサイズ) 574円 1120kcal
• ラージライス 219円 454kcal
• 合計 991円 1895kcal

先駆者様の第2位よりも、わずかに1kcal高い組み合わせが見つかりました。アーリオ・オリオ3人前よりもフォッカチオに追いフォッカするのが若干コスパが高かったようです。

第3位

イメージ図：

• フォッカチオ 119円 214kcal
• アーリオ・オーリオ 299円 560kcal
• アーリオ・オーリオ(Wサイズ) 574円 1120kcal
• 合計 992円 1894kcal

先駆者様の2位と同じものになります。なんか全体的に白い料理ばかりですね。

第4位

せっかくなので4位以降も。

イメージ図：

• フォッカチオ 119円 214kcal
• アーリオ・オーリオ(Wサイズ) 574円 1120kcal
• トッピング半熟卵 69円 90kcal
• ラージライス 219円 454kcal
• 981円 1878kcal

卵という救世主。それでいてカロリーも申し分ない。

第5位

イメージ図：

• アーリオ・オーリオ(Wサイズ) 574円 1120kcal
• ライス 169円 303kcal
• ラージライス 219円 454kcal
• 合計 962円 1877kcal

先駆者様の3位と同じ、炭水化物オンリーでフィニッシュです。

低カロリーメニュー

見ているだけで太りそうな記事になってしまったので、逆に「サイゼリアで最低1000円は食べて、最もカロリーが低い組み合わせ」も探してみることにしました。SMTソルバー的には符号を逆にするだけで求められます。制約の部分を次のように書き換えるだけです。

    constrain $ totalPrice .>= 1000
    minimize "total-calorie" totalCalorie

では、ランキング行ってみます！

第1位

イメージ図：

• わかめサラダ 299円 92kcal
• セロリのピクルス(季節限定) 199円 52kcal
• 熟成ミラノサラミ(Wサイズ) 598円 190kcal
• 合計 1096円 334kcal

なんとなくサラダ3連打にされるかと思ったら、思いのほかまともそうな組み合わせが1位になりました。痩せたいならサラミを食え。

第2位

イメージ図：

• セロリのピクルス(季節限定) 199円 52kcal
• 熟成ミラノサラミ 299円 95kcal
• 熟成ミラノサラミ(Wサイズ) 598円 190kcal
• 合計 1096円 337kcal

こんなにサラミを食べていいんですか？いいんです。

第3位

イメージ図：

• 小エビのサラダ 349円 115kcal
• わかめサラダ 299円 92kcal
• キャベツとアンチョビのソテー 199円 80kcal
• セロリのピクルス(季節限定) 199円 52kcal
• 合計 1046円 339kcal

ここにきて人気No.1の小エビのサラダが登場。こんなにヘルシーだったなんて。これからはサイゼリヤに行ったらこれを頼みます。

まとめ

SMTソルバーはちょっとした問題を解くのにも便利だと思うのですが、あんまり多くの人が使っているイメージがなかったので、こういう記事を書いてみました。環境の準備がめんどいという方（自分含む）のために、Visual Studio Code の Remote - Containers でサクッと使えるようなレポジトリを作ったので、良ければお試しください。

cloneしてVSCode Remoteで開くだけでsbvが使える環境になっているはずです。

コード

今回のコード全体を載せておきます。

{-# LANGUAGE ScopedTypeVariables, DeriveGeneric #-}

import Control.Monad
import Data.SBV
import Data.Csv
import Data.List
import Data.Maybe
import qualified Data.ByteString.Lazy as B
import qualified Data.Vector as V
import GHC.Generics (Generic)
import Text.Printf

-- CREATE TABLE menu(id integer, name text, category text, type text, price integer, calorie integer, salt real);
data Menu = Menu
    { menuId :: Int
    , menuName :: String
    , menuCategory :: String
    , menuType :: String
    , menuPrice :: Int
    , menuCalorie :: Int
    , menuSalt :: Double
    } deriving (Generic, Show)

instance FromRecord Menu

main :: IO ()
main = do
    bs <- B.readFile "saizeriya.csv"
    let Right (menus :: [Menu]) = fmap V.toList $ decode NoHeader bs
    menus <- return $ filter (\m -> menuCalorie m > 0) menus

    let go bests
            | length bests >= 5 = return ()
            | otherwise = do
                LexicographicResult result <-
                    optimize Lexicographic $ maxCalorie menus bests

                let ss = [ menuId m
                         | m <- menus, let Just s = getModelValue (show $ menuId m) result, s ]
                let selectedMenus = [ m | s <- ss, let Just m = find (\m -> menuId m == s) menus]
                let totalPrice = sum $ map menuPrice selectedMenus
                let totalCalorie = sum $ map menuCalorie selectedMenus

                putStrLn ""
                printf "* Rank %d ￥%d %dkcal\n" (length bests + 1) totalPrice totalCalorie
                forM_ selectedMenus $ \menu -> do
                    printf "    * %s ￥%d %dkcal\n" (menuName menu) (menuPrice menu) (menuCalorie menu)

                go $ ss : bests
    go []

maxCalorie :: [Menu] -> [[Int]] -> Symbolic ()
maxCalorie ms doNotUses = do
    ss <- sBools (map (show . menuId) ms)
    let totalPrice :: SInt32 =
            sum [ ite s (fromIntegral $ menuPrice menu) 0 | (menu, s) <- zip ms ss]
    let totalCalorie :: SInt32 =
            sum [ ite s (fromIntegral $ menuCalorie menu) 0 | (menu, s) <- zip ms ss]
    forM_ doNotUses $ \doNotUse ->
        constrain $ bnot $ bAnd [ s | (m, s) <- zip ms ss, menuId m `elem` doNotUse ]
    constrain $ totalPrice .<= 1000
    maximize "total-calorie" totalCalorie

おまけ

サイゼって行ったことないので行ってみたいんですが、とりあえず万札握りしめて行けばいいですか？

— 加藤公一（はむかず） (@hamukazu) May 17, 2019

サイゼに万札握りしめて行ったらどんな恐ろしいことになってしまうのかを計算してみてあげましょう！

Z3的には、制約式の1000と書いてある部分を10000に変えるだけ。ですが、買えるものの組み合わせが爆発的に増えるからなのか、全然計算が終わらなくなってしまいました。進捗もよくわからないし、何分待てばいいのか…。

とりあえず、何かを最適化したいのであれば、その変数に対して2分探索を行えば何とかなります。
それに最悪厳密解が求まらずとも、何らかの近似解は得られます。

何らかの定数$x$に対して、制約式を、$ 総額≦10000 && 総カロリー≧x $ とすれば、$x$よりカロリーが高い組み合わせが存在するか？という判別式になります。そうすれば、この$x$を二分探索すれば、そこそこ高速に答えが求まるはずです。

まずはソルバーに渡す式を修正します。

maxCalorieBS :: [Menu] -> Int -> Symbolic ()
maxCalorieBS ms x = do
    ss <- sBools (map (show . menuId) ms)
    let totalPrice :: SInt32 =
            sum [ ite s (fromIntegral $ menuPrice menu) 0 | (menu, s) <- zip ms ss]
    let totalCalorie :: SInt32 =
            sum [ ite s (fromIntegral $ menuCalorie menu) 0 | (menu, s) <- zip ms ss]
    constrain $ totalPrice .<= 10000
    constrain $ totalCalorie .>= fromIntegral x

最後の最大値を求めるという条件の部分を。

    maximize "total-calorie" totalCalorie

$x$より大きいというものに書き換えています。

    constrain $ totalCalorie .>= fromIntegral x

そしてこれを二分探索で回します。

search :: [Menu] -> IO ()
search menus = do
    let go lo hi
            | hi - lo > 1 = do
                let mi = (lo + hi) `div` 2
                res <- sat z3 $ maxCalorieBS menus mi
                if modelExists res
                    then do
                        printModel res -- いい感じに解を表示する部分を関数にまとめた
                        go mi hi
                    else go lo mi
            | otherwise = do
                return ()
    go 0 100000

これで、じっくり動作を見守ってあげると、答えが出てきます。

もはや何も考えずにコンピューターに解いてもらうという趣旨からは外れてきている気もしてきたり、
SMTソルバーの内部でこれぐらいやってくれてるんじゃないのとかいう疑問もあり。

まあ何はともあれ結果発表です！

10000円コース、第1位

イメージ図：

• ポテトのグリル 199円 366kcal
• フォッカチオ 119円 214kcal
• プチフォッカ 139円 214kcal
• パンチェッタのピザ 399円 646kcal
• 野菜ときのこのピザ 399円 610kcal
• アラビアータ 399円 591kcal
• アーリオ・オーリオ 299円 560kcal
• キャベツのペペロンチーノ 399円 686kcal
• タラコソースシシリー風 399円 605kcal
• パルマ風スパゲッティ 399円 700kcal
• カルボナーラ 499円 865kcal
• アラビアータ(Wサイズ) 770円 1182kcal
• アーリオ・オーリオ(Wサイズ) 574円 1120kcal
• キャベツのペペロンチーノ(Wサイズ) 770円 1372kcal
• タラコソースシシリー風(Wサイズ) 770円 1210kcal
• パルマ風スパゲッティ(Wサイズ) 770円 1400kcal
• カルボナーラ(Wサイズ) 976円 1730kcal
• ミラノ風ドリア 299円 542kcal
• 半熟卵のミラノ風ドリア 368円 632kcal
• セットプチフォッカ付きミラノ風ドリア 378円 649kcal
• ライス 169円 303kcal
• ラージライス 219円 454kcal
• スモールライス 119円 151kcal
• シナモンフォッカチオ 169円 246kcal
• 合計 10000円 17048kcal

定番のポテト、フォッカチオ、トリプルライスに加えて、いろんな味のパスタが総計15人前！待望のピザ、看板メニューのミラノ風ドリアも存分に楽しめる。さらには、今回初登場のデザートフォッカ！これはもう大満足のコースになりました。これぞサイゼリヤ、いや、むしろイタリアそのものであると言い切っても過言ではないかもしれない。しかも1万円ジャスト。これが万札を握りしめてサイゼリヤにった時の、最高のメニューです。一度はこういう豪遊でイタリアの洪水におぼれてみたい。そんな気がするメニューでした。