Asprova コンテストに参加した

年末年始は Asprova Programming Contest をやってました。なんか問題文が軽い気持ちで読むには重すぎる感じなのですが,短くまとめると

  • いくつかの設備を備えた工場があり,そこでいろんな製品を生産しています
  • 各製品はいくつかの工程を経て完成します
  • 各製品の各工程は処理できる設備とできない設備があります
  • いま,どの製品がいつまでにいくつほしいという注文が200件来ました
  • いい感じの稼動スケジュールを求めてください

という感じです。実際にはもっと条件があります。

やったこと(だいたい時系列順)

数字はその時点でのスコアです。ちなみに配布されてるサンプルコードをそのまま提出すると408307056327点です(試してないけどたぶん)。

  • 問題文が何度読んでもよくわからないので,とりあえず配布されてたサンプルコードを読んだり入力を読むコードを自分で書いたりしてみる。
  • サンプルコードを読んで出力を見てみたら,常に最初に見つかった設備に割り付けるようになっており,何も割り付けられない設備がたくさんあることに気付く。ちゃんと分散させたらそれだけでもスコアが上がるはずなので,それを実装してみることにする。
  • 割り付け可能なすべての設備について,すでに割り付けられている一番遅い工程の終了時刻が一番早いところに割り付ける貪欲を書いた(つまりサンプルをすべての設備を見るように変えた)。ついでに注文を納期の早い順にソートしておくことで早く着手したほうがよさそうなものを早い時刻に割り付けるようにした。447088233657
  • 出力を眺めていたら,一部の設備への割り付け時刻がどんどん後ろに遅延していく様子。たぶん単純に後ろに置いていくだけだと前のほうにたくさん隙間ができていそうで,そこに割り付けられる工程がありそう。そこで別の工程の間に入れられそうな隙間があったらそこに詰め込むようにした。段取り時間を計算し直す必要があったりしてこのあたりからだんだん実装がややこしくなる。448737404998
  • 基本的に納期に対して遅れがちになるみたいなので,着手遅延を狙うよりは早めに着手して納期の遅れを減らすほうがスコアはよくなる傾向がありそうだと予想し,なるべく早い時間に工程を詰め込む方針考えることにする。ただし納期遅れが発生していない場合には,最後に納期の範囲内で全体を後ろへずらして着手遅延ボーナスを稼ぎたい。工程が複数ある入力でこれをやるのが面倒だったので,input 1(マスターデータ 1 から作られたやつ)の場合だけ最後にその処理を入れることにした(たぶんスコアへの貢献はかなり小さい)。
  • 今までは開始時刻が最小のところに割り付けていたのだが,段取りのペナルティの係数をちょっといじったものを開始時刻に足した値が最小になるよう割り付けてみた。係数はいろいろ試して一番よかったのを採用する(この評価値,落ち着いて考えてみると質の異なるものを足しているのでどういう意味があるのか疑問だったのだが,最後までこれよりもよいものは見つけられなかったのでしかたなく使い続けた)。449614829656
  • 簡単なビジュアライザを書いた。input 3 では最初の1/3を除いてずっと遅れが発生しているような状態になるらしいことがわかった。この入力だと最初のほうの工程を担っている設備が1--6で他の設備と担当できる工程が disjoint なのだけど,1--6がほとんど休みなしで稼動してこれだけ遅れているので,大幅な改善は難しそうに見える。この時点で上位陣は 4498 に乗っていたので,どこに改良の余地があるのだろうかと不思議に思い始める。ただ,納期の比較的早いものがずっと後回しにされているような印象はあったので,その部分をうまく解決する方法を考えられればよかったのかもしれない。
  • 一つの注文は複数の工程からなっているのだが,注文ごとではなく複数の工程を別々に処理したほうが効率よく割り付けられるのではないか,と思い試してみる。優先度つきキューに処理すべき工程を入れて順に割り付けていくのだが,優先度のつけ方を納期の早い順にすると前よりも少しよくなった。449640963038
  • さらに納期から所要時間の見積りを引いたものを優先度にしてみたり評価関数のパラメータの動かし方を変えてみたりパラメータをひとつ追加したりした(何も割り付けれていない隙間ができることに対してペナルティを課す項を足した)。449730476138
  • 貪欲でこれ以上根本的な改善ができる気がしないので,乱数に頼るかという気持ちになり、ここまでで得られた解から始めて山登りをした。既に割り付けられている注文のいくつかをランダムに除去して割り付け直して,スコアが上がったら採用,を繰り返す。割り付けの制約が複雑なので,制約に違反しないように割り付け直すのがかなり大変で苦労する。assert をたくさん書いた。あと細かいことをいろいろ試して 449800668222
    • ちなみにこのスコアが出たとき提出したコードにはバグがあって手元で動かすとときどき segmentation fault が出るのだけど,このバグ(スコアにはさほど影響するとは思えない)を直すとスコアが1000万点ぐらい落ちるので困惑した。手元では直してもそんなに悪くなっているようには見えず,結局何が起きていたのかわからないままである。もしかしたら,たまたまジャッジサーバで実行したときにいい乱数が出るコードだったのかもしれないが,そんなことあるだろうか?

できるかもしれないと思ったけど結局やらなかったこと

  • 高速化。山登りが2秒では全然収束してなかったのでできたらよかったのだけど。割り付けられそうな隙間を探すときに,線形に探索するのではなく区間木のようなデータ構造を持っておけば効率よく見つけられたかもしれない。線形探索でも平均すると10〜20回ぐらいループを回せば割り付け箇所が見つかっているようだったし,更新のオーバーヘッドも大きくなるはずだし,本当に割り付けられるかどうかは前後の工程とか見ないと判断できないしで,全体として本当に改善するのか確信がなったのでやる気が出なかったところはある。
  • ビームサーチ。使えるタイプの問題なのかなと思ったけど,どういう探索木を考えればいいのかよくわからなかった。そもそもビームサーチを使ったことがないので勝手がわからないというのもある。
  • この問題自体はかなり複雑だけど,適当に問題を単純化したら厳密解が出て,そこから始めてうまく補正したらいい解が作れないかとちらっと思った。そこまでアルゴリズム力がないので詳しく検討せず。

感想

終わってみると,貪欲法で4497が出たあたりから先は問題の特性をうまく捉えて利用するという意味で本質的にいい改善はできなかったように思う。改善できるとしたら納期遅れペナルティを減らすことなんだろうと思ったけどどうすればそれがうまくできるのかわからなかった。

後半でいろいろ細かいことを試せたのはそれなりに勉強になったし悪くはなかったと思う一方で,根拠があるのかないのかわからない思いつきをいろいろ試していただけで意味のある洞察が新しく得られたということはなかったようにも感じる。後者の意味では後半の試行錯誤を振り返ると何をやっていたんだろうなという気もする。

二分探索

二分探索は,値の探索に使うという認識が一般的だと思いますが,実際にはもっと守備範囲が広いです。ということでなるべく一般的な形で書くとどうなるのかを整理した話を書いてみることにしました。why3 で書いて正しいことを検証してみたらけっこうすっきり理解できた気がします。

int 上の述語 f があるとします。f を真にする a と偽にする b がひとつずつわかっていて,a < b であるとしましょう(requires のところ)。そうすると,f x だが f (x+1) ではない x を求めることができる(ensures のところ,result が戻り値),というのが関数の仕様です。

中身はまあ普通の二分探索なのですが,上限 ub は必ず f を満たさないように,下限 lb は必ず f を満たすようにするとわかりやすくて不変条件が簡単に書けます。

ソートされた列 A から値 v のインデックスを求めるよくある二分探索は,f x として A[x] <= v をとった場合になります。

ちなみにコード中では f に特に条件を付けていませんが,真偽が頻繁に入れ替わるような述語に対して使うケースはあまり見た記憶がありません。ある値を境にしてそれより下で真,上で偽(あるいはその逆)になる述語を考えて,その境界の値を求めるのに使うことが多いです。また a, b は自分で与えないといけませんが,たいてい自明な上界と下界があります。

上のコードが正しいことのチェックは,例えば手元の環境では why3 prove -P cvc4 binsearch.mlw を実行すると

binsearch.mlw BinSearch WP_parameter binsearch : Valid (0.03.s)

と出て,正しいことが自動証明できました。int の長さが限られている言語だと足し算がオーバーフローした場合に正しく動かない気がしますが,そこまではチェックしてないようです。

今回は整数上でやりましたが,実数上でも適当なεに対して f x だが f (x+ε) ではない x を求める,のような形で書けるでしょう(たぶん正確にはちょっと違う)。

最短距離と自由豊穣圏

距離空間は豊穣圏であるという Lawvere の論文(http://www.tac.mta.ca/tac/reprints/articles/1/tr1abs.html)がありますが,それに似た感じで重みつきグラフから自由に生成された豊穣圏を考えると hom-object が最短経路になるらしいことに気付いて面白いと思ったので書いてみます.(よく見たら Lawvere の論文にも書いてありましたが,まあいいでしょう.)

R の monoidal structure

R は非負実数の集合に∞を追加した集合とする.これに普通と逆の順序を入れて圏とみなす.実数の普通の足し算はこの圏でのモノイド積になる(0が単位元).また,この圏における直和は inf で与えられる.特に始対象は∞である.モノイド積はすべての直和に対して分配律を満たす.

V-graph と重みつきグラフ

V を可算直和をもつモノイド圏で,分配律  A \otimes \coprod_i B_i \simeq \coprod_i (A \otimes B_i) \left(\coprod_i A_i\right) \otimes B \simeq \coprod_i (A_i \otimes B) が成り立つものとする(Lawvere の論文には詳しい条件が書かれていないけど,後でこれくらい必要になる気がするので最初から仮定しておく).

V-グラフ  \mathcal A は,集合  \mathcal A_0 と,各  u, v \in \mathcal A_0 に対する対象の割り当て  \mathcal A(u, v) \in V からなる.

特に V = R のときを考えると,V-グラフとは非負の重みつき単純有向グラフのことである.逆に非負の重みつき単純有向グラフ G が与えられたとき,辺がないところには重み ∞ の辺があると考えることで G を R-グラフとみなせる.

free V-category と最短距離

V-グラフ \mathcal Aに対して,V-category  \widetilde{\mathcal A}を対象が \mathcal A_0の要素全体で,hom-objectが  \widetilde{\mathcal A}(u, v) = \coprod_{n \ge 0} \mathcal A^n(u, v) となるようなものとする(0は始対象,1はモノイド積の単位元).ただし
 \mathcal A^0(u, v) =
  \begin{cases}
    0 & u \ne v \\ 1 & u = v,
  \end{cases}
  \quad
  \mathcal A^{n+1}(u, v) =
  \coprod_{w \in \mathcal A_0} (\mathcal A^n(w, v)\otimes \mathcal A(u, w))
とする.これは V-category になる(たぶん).Lawvere によると,この構成は V-Cat から V-graph への忘却関手の左随伴になっているそうだ(いかにもなってそうである).

重みつきグラフ G を R グラフとみなして上記の構成を行うと, \widetilde G(u, v)は u から v への最短距離になる( \otimesが足し算で \coprodがinfであることを思い出せばわかる).

Z/nZ の分解周りのメモ

 n = q_1 \dots q_k, q_i = p_i^{e_i} と書いたとき中国剰余定理として知られる同型  \mathbb Z / p^n\mathbb Z \simeq \prod_i \mathbb Z / q_i \mathbb Z の具体的な与え方が毎回思い出せないのでどこかに書いておきたいと思ってメモ。ついでに単数群の生成元も。

同型

 \mathbb Z / n \mathbb Z  \to \prod_i \mathbb Z / q_i \mathbb Z は単純にそれぞれの成分に射影すればよい。逆が少しややこしくて,
 \prod_i \mathbb Z / q_i \mathbb Z \ni (x_1, \dots, x_k) \mapsto \sum_i \left(\frac{n}{q_i} \mod q_i\right)^{-1} \cdot \frac{nx_i}{q_i} \mod n \in \mathbb Z / n \mathbb Z
となる。mod q_i での逆元は互除法を使うか  \phi(q_i) - 1 乗すると求まる。

単数群

p が奇素数のときは  (\mathbb Z / p^n \mathbb Z)^{\times}巡回群で,r を mod p での原始根とすると  (p+1)r で生成される。他にも生成元はたくさんある。

p = 2 のときは n > 2 なら巡回群ではなく, (\mathbb Z / 2^n \mathbb Z)^{\times} \simeq \mathbb Z / 2 \mathbb Z \times \mathbb Z / 2^{n-2} \mathbb Z となる。右から左への同型は  (x \mod 2, y \mod 2^{n-2}) \mapsto (-1)^x \cdot 3^y で与えられる(左辺は乗法群,右辺は加法群なので混乱しないよう注意)。他の同型もたぶんたくさんある。

Graphillion で半順序の個数を数えてみた話

Graphillion で遊んでみた。

やったこと

有限集合上の半順序の個数を数えてみた(cf. OEIS A001035)。これを効率的に求めるのは未解決問題らしい。ちなみに有限集合上の位相の個数を数える問題とだいたい同じ。

素数 n を固定して,Graphillion を使って [1..n] 上の半順序をわりと愚直に列挙した。

実行時間は n=6 で 0.3s, n=7 で 8.6s くらい。n=8 になると三時間ぐらい走らせても終わらなかったのでやめた。ちなみに OEIS を見る限り少なくとも n=18 までは既知らしいので,こういうアプローチで新しい値を見つけるのは無理がありそう。

やりかた

a <= b だったら a から b へ辺を伸ばすことにすると,半順序は有向グラフとみなせる。そこで n 点からなる完全グラフの部分グラフで,半順序になっているようなものを列挙すればよい。Graphillion は基本的な集合演算をサポートしているので,基本的には順序の条件をそのまま書くだけでよい。

ただし,Graphillion は無向グラフを扱うライブラリなのに対して,ここでは有向グラフを扱いたい。そこで,source と target を別の点にして有向グラフをむりやり無向グラフで表すことにした。

一応ちゃんと書くと以下のとおり。有向グラフ  (G, E) に対して無向グラフ  (G+G', E') (a, b') \in E' \iff (a, b) \in E となるように作る。ただし  G' G のコピーで, (-)'\colon G \to G'全単射とする。

実際のコードでは,G の点を正の整数,G' の点を負の整数で表した。あと半順序は反射的反対称的推移的関係だけど,対角線部分はあっても情報が増えないので,非対称的推移的関係(strict order; 日本語だと強順序?)にした。こっちのほうが辺の数が減って少し効率がよくなるようだ。

from graphillion import GraphSet

def posets(n):
    univ = [(i, j) for i in range(1, n+1) for j in range(-n, 0)]
    GraphSet.set_universe(univ)
    x = GraphSet({})

    # asymmetry
    for i in range(1, n+1):
        x = x.excluding((i, -i))

    # transitivity
    for i in range(1, n+1):
        for j in range(1, n+1):
            if i != j:
                for k in range(1, n+1):
                    if j != k:
                        y = x.including([(i, -j), (j, -k)]).excluding((i, -k))
                        x = x - y

    return x

import sys
print(len(posets(int(sys.argv[1]))))

実行してみた。

$ time python3 posets.py 6
130023

real	0m0.328s
user	0m0.293s
sys	0m0.028s

$ time python3 posets.py 7
6129859

real	0m8.648s
user	0m8.251s
sys	0m0.324s

出てきた数字は合ってる。

感想とか

ZDD を使うとコンパクトに表現できるかと思ったが,期待したほどコンパクトにはならなかったようだ。ZDD で効率的に表現できるのは基本的には疎な組合せ集合(正確にはちょっと違うと思う)なのだけど,考えてみれば順序全体ってそれほど疎ではないかもしれない。

強引に無向グラフにしてるところにちょっと無駄がある感じはして,Graphillion を使わないで ZDD を自分で操作すればもうちょっと効率よく列挙できそうには思う。

あとは,順序の表現を単に strict order にしたけど,推移律から推論できる辺を除くなどすればもっと辺の数を節約することはできるかもしれない。あるいは,実は BDD のほうがいいとか,(Z)SDD を使うといいとかはあるかもしれない。

Emacs 上でカーソル位置にある文字の T-Code での入力方法を勝手に表示する

Emacs で日本語を入力するときに T-Code を少し使ってるのですが,ある文字の打ち方を知りたいときにわざわざ 55 って打つのがめんどくさいのでカーソルを合わせたら勝手に表示されるようにしたいとずっと前から定期的に思っていて,今日も思ったので40分くらい調べたらできました。

(defun mylib-show-tc-stroke ()
  (if (and (eq major-mode 'text-mode)
           (not (eq last-command 'skk-insert)))
      (let ((c (char-after)))
        (if (and (integerp c) (< 128 c))
            (tcode-query-stroke (point))))))

(setq mylib-show-tc-stroke-timer
      (run-with-idle-timer 0.2 0.2 'mylib-show-tc-stroke))

関数の中でいろいろ条件つけてるのは適当に書いたのでそんなによくないかも。

あとこれだと勝手にウィンドウが分割されて画面の半分が Help 用のバッファになったりするのであまりスマートでなくて,もうちょっといい感じになったらいいのにと思います。

関数のオーダーの上の順序と足し算

https://twitter.com/yoshihiro503/status/920979507053412353 から始まってなんか議論になっている感じだったのを見て,ちょっと考えてみた話.要するに,関数のオーダー全体の集合は自然な順序が入って join-semilattice になる.あと足し算もできる(実は join と一致する).

以下詳細.

 \mathbb{R}_+は正の実数全体の集合として, f: \mathbb{N} \to \mathbb{R}_+に対してそのオーダーを
 
  O(f) = \{g : \mathbb{N} \to \mathbb{R}_+ \mid \exists n_0, \exists c
  > 0, \forall n \ge n_0, g(n) \le cf(n)\}
とする.オーダー全体の集合を \mathbb{O}と書くことにする.この集合には包含関係で半順序が入る.

 \mathbb{O}は全順序ではない*1が,任意の二元は最小上界をもち,それは O(f) \vee O(g) = O(\max(f, g)) = O(f+g)で与えられる( \max(f, g)は各点ごとのmax).

証明:  O(f), O(g) \subseteq O(\max(f,g))は明らか.もし O(f), O(g) \subseteq O(h)とすると, f \le ch,  g \le ch*2となる cがとれる(この cは共通にとれる).このとき \max(f,g) \le chであるから \max(f,g) \in O(h)すなわち O(\max(f,g)) \subseteq O(h)である. O(\max(f,g)) = O(f+g)であることは \max(f,g) \le f+g \le 2\max(f,g)だから明らか.

オーダーの和を O(f) + O(g) = \{ h \mid \exists h_1 \in O(f), \exists h_2 \in O(g), h = f + g\}と定義すると, O(f) + O(g) = O(f+g)が成り立つ.

証明: 左から右への包含関係は明らか.逆にもし h \in O(f+g)ならば,ある cが存在して h \le c(f+g)となる.このとき h = \frac{f}{f+g}h + \frac{g}{f+g}hと分解してやればそれぞれ \frac{f}{f+g}h \le cf, \frac{g}{f+g}h \le cgとなるので h \in O(f) + O(g)である.

*1:証明は簡単な演習問題

*2:正確には,有限個を除いてすべて点でこの不等式が成り立つ.以下の不等号も同様.