Advent of Code 2022 を完走した

毎年12月に開催されている Advent of Code に、2019年から参加している。

過去記事:

2022年のAdvent of Codeにも挑戦していて、年が明けてしまったが先日ようやく25日すべての問題に解答して 50 個のスターを集めることができた。

2022年の問題もどれもとても面白かった。

データ構造を扱う問題、ある程度のアルゴリズム知識が求められる問題、ちょっとした数学的な考え方が必要な問題、ひたすら実装が大変な問題、視覚化が面白い問題、などが揃っているのは毎年だが、特に2022年はバランス良くバラエティに富むものが出題されていて取り組み甲斐があったように思う。

工夫のしどころも多くあって、自分では上手く解いたつもりでも Reddit で他の人の解答を見ると目から鱗のものがたくさん発見できたり。

基本的にはヒント無しで自力で正解まで辿り着いたが、day19のpart2だけはどうにもならなくて他の人の解答を見たりしてようやく、という感じになってしまった。悔しい。

全部Rustで解いていて、あとでPythonでも解こうとしている。

github.com

とても楽しかったし勉強になったので、もうちょっと同じ問題に取り組んだ人達でわちゃわちゃと「ここが難しかった」「これはこんな手法があって」「こんな書き方もあって」など議論できると嬉しい…。 それはRedditで英語で頑張れば良いのだろうけど、ともかくもう少し日本語圏の開発者の間でも流行って仲間が増えてくれるといいな、と思っている。

という想いもあり、2020のときにもやったけど 日本語の解説本を書いてみている。

github.com

まだ半分もいってないけど… 少しでも多くのプログラマの人たちに届くといいな

2023パズル をRustで解いてみる

tkihiraさんの問題が面白そうだったので挑戦してみた。

既に解説記事が出ているので解答はこちらをどうぞ。

nmi.jp

結局自分は自力では解けなくて 他の人の解法や上記の解説記事を読んでようやくできた、のだけど… 自分なりに理解して改めてRustで実装してみた。

RPN(逆ポーランド記法)の backtracking

まずは型定義。

#[derive(Clone, Copy, PartialEq, Eq)]
enum Op {
    Add,
    Sub,
    Mul,
    Div,
}

#[derive(Clone)]
enum ExpressionElement {
    Operand(i32),
    Operator(Op),
}

この ExpressionElement を stack に詰んでいって、末尾が Operator だったらさらに Operand を 2つ pop() して、演算結果を push() していくことで計算が実現できる。

常に Operand の方が多い状態でないと Operator は挿入できないので、各々の出現回数をカウントしながら条件を満たしているときだけ操作するようにする。

様々な実装が試せるように Trait を定義。

trait Rpn {
    fn traverse(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
    ) {
        if i == j + 1 && self.get(i).is_none() {
            if self.evaluate(expr) {
                results.push(expr.clone());
            }
            return;
        }
        if let Some(&n) = self.get(i) {
            self.backtrack(expr, (i + 1, j), results, ExpressionElement::Operand(n));
        }
        if i > j + 1 {
            for op in &[Op::Add, Op::Sub, Op::Mul, Op::Div] {
                self.backtrack(expr, (i, j + 1), results, ExpressionElement::Operator(*op));
            }
        }
    }
    fn evaluate(&self, expr: &[ExpressionElement]) -> bool;
    fn get(&self, i: usize) -> Option<&i32>;
    fn backtrack(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
        e: ExpressionElement,
    ) {
        expr.push(e);
        self.traverse(expr, (i, j), results);
        expr.pop();
    }
}

self.get(i)i番目に使う数字を呼び出す。すべて使い切っていて Operatorも適切な回数使われていたら self.evaluate(expr) で計算結果を確認し、条件に合致していた場合のみ その exprresults に格納する。

exprpush / pop して再帰的に self.traverse() を呼ぶ backtracking の部分は別の実装で上書きできるようにあえて独立のmethodで定義。

何も難しいことを考えずにこれを満たす実装を作ると、以下のようになる。

struct Fraction(i32, i32);

struct DefaultSearcher {
    nums: Vec<i32>,
    target: i32,
}

impl DefaultSearcher {
    pub fn new(nums: Vec<i32>, target: i32) -> Self {
        Self { nums, target }
    }
}

impl Rpn for DefaultSearcher {
    fn evaluate(&self, expr: &[ExpressionElement]) -> bool {
        let mut stack = Vec::new();
        for e in expr {
            match e {
                ExpressionElement::Operand(n) => stack.push(Fraction(*n, 1)),
                ExpressionElement::Operator(op) => {
                    if let (Some(n0), Some(n1)) = (stack.pop(), stack.pop()) {
                        if let Some(n) = op.apply(&n1, &n0) {
                            stack.push(n);
                        } else {
                            return false;
                        }
                    }
                }
            }
        }
        stack
            .last()
            .map(|n| n.1 * self.target == n.0)
            .unwrap_or(false)
    }
    fn get(&self, i: usize) -> Option<&i32> {
        self.nums.get(i)
    }
}

impl Op {
    fn apply(&self, lhs: &Fraction, rhs: &Fraction) -> Option<Fraction> {
        match self {
            Self::Add => Some(Fraction(lhs.0 * rhs.1 + rhs.0 * lhs.1, lhs.1 * rhs.1)),
            Self::Sub => Some(Fraction(lhs.0 * rhs.1 - rhs.0 * lhs.1, lhs.1 * rhs.1)),
            Self::Mul => Some(Fraction(lhs.0 * rhs.0, lhs.1 * rhs.1)),
            Self::Div if rhs.0 != 0 => Some(Fraction(lhs.0 * rhs.1, lhs.1 * rhs.0)),
            _ => None,
        }
    }
}

解説記事では浮動小数点でも十分ということでやっていたけど、自分的には整数値だけで処理したいので Fraction(i32, i32) を定義して分数で計算するようにした。約分のことはとりあえず考えなくてよくて、最終的な値の計算結果の確認は 分子が分母の 2023 倍になっているか否か、で判定できる。 各Opの適用結果を Op::apply で計算。ゼロ除算の可能性があるので返り値は Option<Fraction> に。

self.evaluate() 内で stackを用意し、exprを順に見ながら Operand だったら Fraction に変換してpushし、 Operator だったら 2つを pop して演算結果をまた格納。理論上 正しい形式のRPN式であれば最終的にはstackには1要素だけ残るはずなので、その値を確認すれば良い。

これでとりあえず、(98の間の割り算を無視して)計算結果が 2023 になる式の組み合わせを列挙できる。

let nums = (1..=10).rev().collect();
let mut results = Vec::new();
DefaultSearcher::new(nums, 2023).traverse(&mut Vec::new(), (0, 0), &mut results);

手元の環境で実行すると 3370 件の結果が約2分弱で列挙された。

...

10-(9-(((8*(7*6))+(5-4))*(3*(2*1))))
10-(9-(((8*(7*6))+(5-4))*(3*(2/1))))
10-(9-(((8*(7*6))+(5-4))*(3+(2+1))))
10-(9-(8*(((7*(6+(5/4)))*(3+2))-1)))
10-(9-(8*((7*((6+(5/4))*(3+2)))-1)))
Completed 3370 results in 111.933449708s

探索の高速化

計算結果を保持

前述の実装だと、式が完成するたびにstack用意して計算結果を確認して… というのが非常に無駄になる。解説記事によると 1,274,544,128 通り生成されるようなので、ここはもっと計算量を削りたい。

途中までの計算結果は変わらないものなので、それを保持しながら探索した方が効率的になる。 内部に別のstackを持って、 Operator を受け取った時点で演算処理結果を格納しておくようにしておく。

前述の self.evaluate() でやっていた計算を self.backtrack() の中で逐次やっておく、という感じ。

pub struct FastSearcher {
    nums: Vec<i32>,
    target: i32,
    stack: Vec<Fraction>,
}

impl FastSearcher {
    pub fn new(nums: Vec<i32>, target: i32) -> Self {
        let stack = Vec::with_capacity(nums.len());
        Self {
            nums,
            target,
            stack,
        }
    }
}

impl Rpn for FastSearcher {
    fn evaluate(&self, _: &[ExpressionElement]) -> bool {
        self.stack[0].1 * self.target == self.stack[0].0
    }
    fn get(&self, i: usize) -> Option<&i32> {
        self.nums.get(i)
    }
    fn backtrack(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
        e: ExpressionElement,
    ) {
        match e {
            ExpressionElement::Operand(n) => {
                self.stack.push(Fraction(n, 1));
                expr.push(e);
                self.traverse(expr, (i, j), results);
                expr.pop();
                self.stack.pop();
            }
            ExpressionElement::Operator(op) => {
                if let (Some(n0), Some(n1)) = (self.stack.pop(), self.stack.pop()) {
                    if let Some(n) = op.apply(&n1, &n0) {
                        self.stack.push(n);
                        expr.push(e);
                        self.traverse(expr, (i, j), results);
                        expr.pop();
                        self.stack.pop();
                    }
                    self.stack.push(n1);
                    self.stack.push(n0);
                }
            }
        }
    }
}

これによって self.evaluate() に到達した時点で計算結果は出ているので、その値だけを確認すれば良い、ということになる。また、途中でゼロ除算が生じた場合にその先の探索をしなくなるので無駄を省く効果もありそうだ。

これで、前述の 3370 件の列挙が 約30秒に短縮された。約4倍の高速化。

探索結果のキャッシュ

また、この内部stackは同じ状態を経過することも多い。例えば

  • 10*(9-8)+...
  • 10/(9-8)+...

8 までの計算結果はどちらも変わらない。どちらか一方で全探索して1件も見つからなかったなら、もう片方でも1件も見つからないことが分かる。

ので、探索して条件に合致する式が見つかったか否かを返すようにし、それが false だった場合には同様の状態からの探索をスキップするように変更。

use std::collections::HashSet;

type Key = (Vec<Fraction>, (usize, usize));

pub struct FastSearcher {
    nums: Vec<i32>,
    target: i32,
    stack: Vec<Fraction>,
    seen: HashSet<Key>,
}

impl FastSearcher {
    pub fn new(nums: Vec<i32>, target: i32) -> Self {
        let stack = Vec::with_capacity(nums.len());
        Self {
            nums,
            target,
            stack,
            seen: HashSet::new(),
        }
    }
}

impl Rpn for FastSearcher {
    fn traverse(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
    ) -> bool {
        if i == j + 1 && self.get(i).is_none() {
            let found = self.evaluate(expr);
            if found {
                results.push(expr.clone());
            }
            return found;
        }
        let key = (self.stack.clone(), (i, j));
        if self.seen.contains(&key) {
            return false;
        }
        let mut ret = false;
        if let Some(&n) = self.get(i) {
            ret |= self.backtrack(expr, (i + 1, j), results, ExpressionElement::Operand(n));
        }
        if i > j + 1 {
            for op in &Op::ALL {
                ret |= self.backtrack(expr, (i, j + 1), results, ExpressionElement::Operator(*op));
            }
        }
        if !ret {
            self.seen.insert(key);
        }
        ret
    }
}

これで、約12秒くらいに短縮された。さらに2倍以上速くなった。

実際にはこの key に使っている self.stack は約分されていない計算結果なので、同じ状態を経由しているはずなのに 値が異なっているという扱いになってしまうケースがある。 gcd を使って既約分数を格納するようにしてみる。

impl FastSearcher {
    fn gcd(a: i32, b: i32) -> i32 {
        if b == 0 {
            a
        } else {
            Self::gcd(b, a % b)
        }
    }
}

impl Rpn for FastSearcher {
    ...

    fn backtrack(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
        e: ExpressionElement,
    ) -> bool {
        let mut ret = false;
        match e {
            ExpressionElement::Operand(n) => {...}
            ExpressionElement::Operator(op) => {
                if let (Some(n0), Some(n1)) = (self.stack.pop(), self.stack.pop()) {
                    if let Some(n) = op.apply(&n1, &n0) {
                        let gcd = Self::gcd(n.0, n.1);
                        self.stack.push(Fraction(n.0 / gcd, n.1 / gcd));

                        ...
                    }
                    self.stack.push(n1);
                    self.stack.push(n0);
                }
            }
        }
        ret
    }
}

これで、3370 件の列挙が 約7.5秒程度まで短縮された。 このキャッシュ手法は探索目標の値などによって効率が変わってくるので一概には言えないかもしれないが、少なくともこの例においては数倍の高速化が実現できた。

割り算を考慮した条件つき生成

ここまでは計算結果が 2023 になるものを全列挙することだけを考えていたが、件のパズル問題では「98の間には既に÷が入っています」という制約がある。

前述した 3370 件から条件に合致するものを抽出しても良いが、そもそも探索する時点で条件に合致しない式になるものを除外していくほうが効率が良い。

これは自分ではまったく思いつかなくて解説記事を読んで目から鱗だったのだけど、 8の数字が出た後、数字と演算子の数が一致していた場合」の演算子を割り算に固定 することによって実現できる。

解説記事を参考に、専用のSearcherを実装してみる。

struct Div8Searcher {
    nums: Vec<i32>,
    target: i32,
    stack: Vec<Fraction>,
    eight_depth: i32,
}

impl Div8Searcher {
    pub fn new(nums: Vec<i32>, target: i32) -> Self {
        let stack = Vec::with_capacity(nums.len());
        Self {
            nums,
            target,
            stack,
            eight_depth: -1,
        }
    }
}

impl Rpn for Div8Searcher {
    fn traverse(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
    ) -> bool {

        ...

        if i > j + 1 {
            for op in &Op::ALL {
                if self.eight_depth == 0 && op != &Op::Div {
                    continue;
                }
                ret |= self.backtrack(expr, (i, j + 1), results, ExpressionElement::Operator(*op));
            }
        }
        ret
    }
    fn backtrack(
        &mut self,
        expr: &mut Vec<ExpressionElement>,
        (i, j): (usize, usize),
        results: &mut Vec<Vec<ExpressionElement>>,
        e: ExpressionElement,
    ) -> bool {
        let mut ret = false;
        match e {
            ExpressionElement::Operand(n) => {
                let orig = self.eight_depth;
                if n == 8 {
                    self.eight_depth = 0;
                } else if self.eight_depth >= 0 {
                    self.eight_depth += 1;
                }
                self.stack.push(Fraction(n, 1));
                expr.push(e);
                ret |= self.traverse(expr, (i, j), results);
                expr.pop();
                self.stack.pop();
                self.eight_depth = orig;
            }
            ExpressionElement::Operator(op) => {
                if let (Some(n0), Some(n1)) = (self.stack.pop(), self.stack.pop()) {
                    if let Some(n) = op.apply(&n1, &n0) {
                        self.eight_depth -= 1;
                        self.stack.push(n);
                        expr.push(e);
                        ret |= self.traverse(expr, (i, j), results);
                        expr.pop();
                        self.stack.pop();
                        self.eight_depth += 1;
                    }
                    self.stack.push(n1);
                    self.stack.push(n0);
                }
            }
        }
        ret
    }
}

eight_depth の値を保持し、8 が出現したら 0 にセットして OperandOperator かによって値を増減。値が 0 のときには次に push する OpeatorOp::Div だけに制限される。

前述したような探索結果のキャッシュも使いたいところだが、新たな状態としてこの self.eight_depth もキーに含める必要があり、その割にはそれほど効果も無さそうなのでここでは省いた。

これによって 「98の間には必ず/が入る」式だけが探索されるようになる。

条件に合致する 530 件が、約4秒強で列挙されるようになった。全列挙してから抽出するよりも明らかに速いことが分かる。

...

(10*(9/((8-7)/((6*5)/(4/3)))))-(2/1)
(10*(9/((8-7)/(6*((5/4)*3)))))-(2*1)
(10*(9/((8-7)/(6*((5/4)*3)))))-(2/1)
(10*(9/((8-7)/(6*(5/(4/3))))))-(2*1)
(10*(9/((8-7)/(6*(5/(4/3))))))-(2/1)
Completed 530 results in 4.187013041s

中置記法、正規化

ここまでで探索して収集された resultsVec<Vec<ExpressionElement>> であり、これを中置記法に書き換える操作もまたstackを使って実現していくことになる。

impl Display for Op {
    fn fmt(&self, f: &mut Formatter<'_>) -> Result {
        f.write_char(match self {
            Self::Add => '+',
            Self::Sub => '-',
            Self::Mul => '*',
            Self::Div => '/',
        })
    }
}

fn solve(searcher: &mut impl Rpn) {
    let mut results = Vec::new();
    searcher.traverse(&mut Vec::new(), (0, 0), &mut results);
    for result in results {
        let mut stack = Vec::new();
        for e in &result {
            match e {
                ExpressionElement::Operand(n) => stack.push((n.to_string(), None)),
                ExpressionElement::Operator(op) => {
                    if let (Some((mut s0, o0)), Some((mut s1, o1))) = (stack.pop(), stack.pop()) {
                        if o1.is_some() {
                            s1 = format!("({s1})");
                        }
                        if o0.is_some() {
                            s0 = format!("({s0})");
                        }
                        stack.push((format!("{s1}{op}{s0}"), Some(*op)));
                    }
                }
            }
        }
        println!("{}", stack[0].0);
    }
}

連結していくための String と、どの Operator で得られたものか(もしくは Operand か)、を保持しながら処理していく。簡単にすべて括弧をつけてしまうのなら上述のように書ける。

不要な括弧を除去して出力する場合は、左operandの場合と右operandの場合で処理が変わる(引き算・割り算は右operandに注意する必要がある)。matches! で判定してやってみる。

        ...

        ExpressionElement::Operator(op) => {
            if let (Some((mut s0, o0)), Some((mut s1, o1))) = (stack.pop(), stack.pop()) {
                if matches!((op, o1), (Op::Mul | Op::Div, Some(Op::Add | Op::Sub))) {
                    s1 = format!("({s1})");
                }
                if matches!(
                    (op, o0),
                    (Op::Sub | Op::Mul, Some(Op::Add | Op::Sub)) | (Op::Div, Some(_))
                ) {
                    s0 = format!("({s0})");
                }
                stack.push((format!("{s1}{op}{s0}"), Some(*op)));
            }
        }

これで、得られた結果を HashSet などで重複排除することで、解説記事と同様の 81 件を約4秒強で列挙することができた。

((10+(9/8+7)*6*5)*4-3)*2-1
(10*9/((8-7)/(6*5))/(4/3)-2)*1
(10*9/((8-7)/(6*5))/(4/3)-2)/1
(10*9/((8-7)/(6*5))/4*3-2)*1
(10*9/((8-7)/(6*5))/4*3-2)/1
(10*9/((8-7)/(6*5)*4)*3-2)*1
(10*9/((8-7)/(6*5)*4)*3-2)/1
(10*9/((8-7)/(6*5)*4/3)-2)*1
(10*9/((8-7)/(6*5)*4/3)-2)/1
(10*9/((8-7)/(6*5/(4/3)))-2)*1
(10*9/((8-7)/(6*5/(4/3)))-2)/1
(10*9/((8-7)/(6*5/4))*3-2)*1
(10*9/((8-7)/(6*5/4))*3-2)/1
(10*9/((8-7)/(6*5/4)/3)-2)*1
(10*9/((8-7)/(6*5/4)/3)-2)/1
(10*9/((8-7)/(6*5/4*3))-2)*1
(10*9/((8-7)/(6*5/4*3))-2)/1
(10*9/((8-7)/6)*5/(4/3)-2)*1
(10*9/((8-7)/6)*5/(4/3)-2)/1
(10*9/((8-7)/6)*5/4*3-2)*1
(10*9/((8-7)/6)*5/4*3-2)/1
(10*9/((8-7)/6/(5/(4/3)))-2)*1
(10*9/((8-7)/6/(5/(4/3)))-2)/1
(10*9/((8-7)/6/(5/4))*3-2)*1
(10*9/((8-7)/6/(5/4))*3-2)/1
(10*9/((8-7)/6/(5/4)/3)-2)*1
(10*9/((8-7)/6/(5/4)/3)-2)/1
(10*9/((8-7)/6/(5/4*3))-2)*1
(10*9/((8-7)/6/(5/4*3))-2)/1
(10*9/((8-7)/6/5)/(4/3)-2)*1
(10*9/((8-7)/6/5)/(4/3)-2)/1
(10*9/((8-7)/6/5)/4*3-2)*1
(10*9/((8-7)/6/5)/4*3-2)/1
(10*9/((8-7)/6/5*4)*3-2)*1
(10*9/((8-7)/6/5*4)*3-2)/1
(10*9/((8-7)/6/5*4/3)-2)*1
(10*9/((8-7)/6/5*4/3)-2)/1
(10*9/(8-7)*6*5/(4/3)-2)*1
(10*9/(8-7)*6*5/(4/3)-2)/1
(10*9/(8-7)*6*5/4*3-2)*1
(10*9/(8-7)*6*5/4*3-2)/1
10*9/((8-7)/(6*5))/(4/3)-2*1
10*9/((8-7)/(6*5))/(4/3)-2/1
10*9/((8-7)/(6*5))/4*3-2*1
10*9/((8-7)/(6*5))/4*3-2/1
10*9/((8-7)/(6*5)*4)*3-2*1
10*9/((8-7)/(6*5)*4)*3-2/1
10*9/((8-7)/(6*5)*4/3)-2*1
10*9/((8-7)/(6*5)*4/3)-2/1
10*9/((8-7)/(6*5/(4/3)))-2*1
10*9/((8-7)/(6*5/(4/3)))-2/1
10*9/((8-7)/(6*5/4))*3-2*1
10*9/((8-7)/(6*5/4))*3-2/1
10*9/((8-7)/(6*5/4)/3)-2*1
10*9/((8-7)/(6*5/4)/3)-2/1
10*9/((8-7)/(6*5/4*3))-2*1
10*9/((8-7)/(6*5/4*3))-2/1
10*9/((8-7)/6)*5/(4/3)-2*1
10*9/((8-7)/6)*5/(4/3)-2/1
10*9/((8-7)/6)*5/4*3-2*1
10*9/((8-7)/6)*5/4*3-2/1
10*9/((8-7)/6/(5/(4/3)))-2*1
10*9/((8-7)/6/(5/(4/3)))-2/1
10*9/((8-7)/6/(5/4))*3-2*1
10*9/((8-7)/6/(5/4))*3-2/1
10*9/((8-7)/6/(5/4)/3)-2*1
10*9/((8-7)/6/(5/4)/3)-2/1
10*9/((8-7)/6/(5/4*3))-2*1
10*9/((8-7)/6/(5/4*3))-2/1
10*9/((8-7)/6/5)/(4/3)-2*1
10*9/((8-7)/6/5)/(4/3)-2/1
10*9/((8-7)/6/5)/4*3-2*1
10*9/((8-7)/6/5)/4*3-2/1
10*9/((8-7)/6/5*4)*3-2*1
10*9/((8-7)/6/5*4)*3-2/1
10*9/((8-7)/6/5*4/3)-2*1
10*9/((8-7)/6/5*4/3)-2/1
10*9/(8-7)*6*5/(4/3)-2*1
10*9/(8-7)*6*5/(4/3)-2/1
10*9/(8-7)*6*5/4*3-2*1
10*9/(8-7)*6*5/4*3-2/1
Completed 81 results in 4.229359166s

CLIアプリケーション化

ここまで幾つかの実装ができたので、折角なので引数でパラメータを変えながら全列挙できる CLI アプリケーションを作成してみた。

$ ./puzzle2023 --help
Solvers for https://twitter.com/tkihira/status/1609313732034965506

Usage: puzzle2023 [OPTIONS]

Options:
  -t, --target <TARGET>      [default: 2023]
  -m, --max-num <MAX_NUM>    [default: 10]
  -r, --rev <REV>            [default: true] [possible values: true, false]
  -s, --searcher <SEARCHER>  [default: fast] [possible values: default, fast, div8]
  -n, --normalize
  -v, --verbose
  -h, --help                 Print help information
  -V, --version              Print version information

全パターン列挙

$ ./puzzle2023 -v

...

10-(9-(((8*(7*6))+(5-4))*(3*(2*1))))
10-(9-(((8*(7*6))+(5-4))*(3*(2/1))))
10-(9-(((8*(7*6))+(5-4))*(3+(2+1))))
10-(9-(8*(((7*(6+(5/4)))*(3+2))-1)))
10-(9-(8*((7*((6+(5/4))*(3+2)))-1)))
Completed 3370 results in 7.338064208s

9と8の間の割り算に限定

$ ./puzzle2023 -v --searcher div8

...

(10*(9/((8-7)/(6*((5/4)*3)))))-(2/1)
(10*(9/((8-7)/(6*(5/(4/3))))))-(2*1)
(10*(9/((8-7)/(6*(5/(4/3))))))-(2/1)
Completed 530 results in 4.219031208s

正規化

$ ./puzzle2023 -v --searcher div8 --normalize

...

10*9/(8-7)*6*5/(4/3)-2/1
10*9/(8-7)*6*5/4*3-2*1
10*9/(8-7)*6*5/4*3-2/1
Completed 81 results in 4.239637041s

計算結果の確認

$ ./puzzle2023 --searcher div8 --normalize | python3 -c 'for s in open(0): print(f"{eval(s):.1f}")' | sort | uniq -c
  81 2023.0

使うのを 1 から 8 の昇順にして 2024 を作るようにしてみる

$ ./puzzle2023 -m8 -rfalse -t2024 -nv
(((1+2)*3*4+5)*6+7)*8
(1*2*(3+4*5*6)+7)*8
(1+(2+3*4)*(5+6+7))*8
(1+(2+3-(4-5))*6*7)*8
(1+(2+3-4+5)*6*7)*8
(1+2*(3+4)*(5+6+7))*8
(1+2*(3+4+5+6)*7)*8
(1+2*3+4)*(5*6-7)*8
(1+2/(3/((4+5)*6))*7)*8
(1+2/(3/((4+5)*6)/7))*8
(1+2/(3/((4+5)*6*7)))*8
(1+2/(3/(4+5))*6*7)*8
(1+2/(3/(4+5)/(6*7)))*8
(1+2/(3/(4+5)/6)*7)*8
(1+2/(3/(4+5)/6/7))*8
(1+2/3*(4+5)*6*7)*8
(1-(2-3*4))*(5*6-7)*8
(1-2*(3*4-5*6)*7)*8
(1-2*3*(4-5)*6*7)*8
(1-2*3/((4-5)/(6*7)))*8
(1-2*3/((4-5)/6)*7)*8
(1-2*3/((4-5)/6/7))*8
(1-2*3/(4-5)*6*7)*8
(1-2+3*4)*(5*6-7)*8
1*(2*(3+4*5*6)+7)*8
Completed 25 results in 131.605458ms

「8の前を割り算に限定」にしているのも、コマンドラインオプションで違う値で指定できるようにしたかったけど、複雑になりすぎる感じがしたのでそれは断念…

改善?

tkihiraさんの最適化ではstackを push / pop するのではなく固定長配列を用意して見るべきindexだけをズラしていく、といった工夫をしていた。真似して実装してみたが本実装においてはそこはあまり効果なさそうだった。

あとはこうして再帰的に探索した結果を一度 results という Vec に格納してしまっているのが実は非効率なのでは、と思ったが… Iterator として使いたかったが そうするためには探索の途中の状態を保持しておいて next() が呼ばれるたびに続きを探索していく、という処理にする必要があり、そのあたりの実装がうまくできなかった。何か良い方法あるのだろうか…

Repository

github.com

Rubyでバイナリデータに対するrindex検索の挙動でハマったので調べたことメモ

自分の手元の環境でこんなことが起きた。

$ ruby -v
ruby 3.1.2p20 (2022-04-12 revision 4491bb740a) [arm64-darwin21]
$ irb
irb(main):001:0> "\x01\x80\x00\x00".index("\x01")
=> 0
irb(main):002:0> "\x01\x80\x00\x00".rindex("\x01")
=> 1

\x010 番目にしかないのだから、 .index でも .rindex でも 0 が返ってくるはずではないの??

先に結論

バイナリデータを扱うときには必ずEncodingを ASCII-8BIT に指定しておくこと!

きっかけ

roo というgem (記事時点で 2.9.0)を使って、Excelファイルを開こうとした。

require 'roo'

p Roo::Excelx.new("hoge.xlsx")

これは問題ないが、 #initialize の引数は file_or_strem ということでファイルを open したものを渡しても良いはず、と

require 'roo'

File.open("hoge.xlsx") do |f|
  p Roo::Excelx.new(f)
end

ということをすると以下のような謎のエラーを吐いて落ちる。

/Users/sugyan/.rbenv/versions/3.1.2/lib/ruby/gems/3.1.0/gems/rubyzip-2.3.2/lib/zip/entry.rb:365:in `check_c_dir_entry_static_header_length': undefined method `bytesize' for nil:NilClass (NoMethodError)

      return if buf.bytesize == ::Zip::CDIR_ENTRY_STATIC_HEADER_LENGTH
                   ^^^^^^^^^

で困っていたので他の人に聞いてみたところ「自分のところではそんな問題は起きていない」と言われ。 試しにDockerでRuby環境作って実行してみると、確かに問題なく動く…何故? と調べはじめた。

String#rindex の謎挙動

自分の環境と問題なく動く環境とで比較しながら見てみると、roo が使っている rubyzip (記事時点で 2.3.2) の中でファイル内容を読み取った結果の値が違っていた。

module Zip
  class CentralDirectory
    include Enumerable

    END_OF_CDS             = 0x06054b50

    ...

    def get_e_o_c_d(buf) #:nodoc:
      sig_index = buf.rindex([END_OF_CDS].pack('V'))
      ...

buf の中から 0x06054b50 をpackしたもの つまり "PK\x05\x06" のデータを末尾から探すために .rindex() を呼んでいるのだが、手元の環境と 問題なく動く環境ではその値が 1 ズレている。。

何故?と思って色々試しているうちに冒頭の例のようなものに辿り着いた。

もう少し深く追う

少なくとも \x80 などを含まないASCIIだけのものであればおかしな挙動にはならない。

irb(main):001:0> "\x01\x7F\x00\x00".rindex("\x01")
=> 0

また、これが増えるとズレもどんどん大きくなる。

irb(main):001:0> "\x01\x80".rindex("\x01")
=> 0
irb(main):002:0> "\x01\x80\x80".rindex("\x01")
=> 1
irb(main):003:0> "\x01\x80\x80\x80".rindex("\x01")
=> 2
irb(main):004:0> "\x01\x80\x80\x80\x80".rindex("\x01")
=> 3

逆から走査する際に異常な文字が検出された際に読み飛ばし、その結果の走査数を文字列長から引いたものを返したことにより値がおかしくなる、というのが予想される。

ソースを読んでみる。 rstring 関連はこのへんのようだ。

#ifdef HAVE_MEMRCHR によって実装が分かれている。 memrchr というのがglibcに入っているもので、自分のmacOS環境では使えなくて Docker内のLinux環境などでは使える、これによって環境によって挙動が変わったりする、のかもしれない。

で、使えない環境の方での実装は。 対象が見つかるまで while loop の中で rb_enc_prev_char で前の文字を走査している、という感じのようだ。encodingの影響を受けそう…?

    while (s) {
        if (memcmp(s, t, slen) == 0) {
            return pos;
        }
        if (pos == 0) break;
        pos--;
        s = rb_enc_prev_char(sbeg, s, e, enc);
    }

Encodingと実行環境

Encodingが関係しているようだったので色々調べてみた。

Rubyではバイナリデータ(バイト列)も String で扱う。

バイナリの取扱い

Ruby の String は、文字の列を扱うためだけでなく、バイトの列を扱うためにも使われます。しかし、Ruby M17N には直接にバイナリを表すエンコーディングは存在しません。このため、バイナリを String で扱う際には、ASCII 互換オクテット列を意味する ASCII-8BIT を用います。これにより、ASCII 互換であるこの String は 7bit クリーンな文字列と比較・結合が可能となります。

https://docs.ruby-lang.org/ja/latest/doc/spec=2fm17n.html

ということで、バイナリの場合は本来は ASCII-8BIT のEncodingを持つ文字列として検索しなければならないのに、手元の環境では UTF-8 のままで検索しようとしていた、というところに「使い方の問題」があったようだ。

irb(main):001:0> "\x01\x80\x00\x00".encoding
=> #<Encoding:UTF-8>

Encodingが ASCII-8BIT になっていれば、rindex でも正しい値を得ることができそうだ。 全然知らなかったが String#bASCII-8BIT の複製を得ることもできるらしい。

https://docs.ruby-lang.org/ja/latest/method/String/i/b.html

irb(main):001:0> "\x01\x80\x00\x00".rindex("\x01")
=> 1
irb(main):002:0> "\x01\x80\x00\x00".force_encoding(Encoding::ASCII_8BIT).rindex("\x01")
=> 0
irb(main):003:0> "\x01\x80\x00\x00".b.rindex("\x01")
=> 0

なるほど〜。

ではこの Encoding は何で決まるのか?というのも上記リンクに書いてある。

リテラルエンコーディング

文字列リテラル正規表現リテラルそしてシンボルリテラルから生成されるオブジェクトのエンコーディングスクリプトエンコーディングになります。

またスクリプトエンコーディングが US-ASCII である場合、7bit クリーンではないバックスラッシュ記法で表記されたリテラルエンコーディングは ASCII-8BIT になります。

さらに Unicode エスケープ (\uXXXX) を含む場合、リテラルエンコーディングUTF-8 になります。

複雑。。。

とにかく通常はスクリプトエンコーディングがまず大事のようだ。

スクリプトエンコーディング

スクリプトエンコーディングとは Ruby スクリプトを書くのに使われているエンコーディングです。スクリプトエンコーディングは マジックコメントを用いて指定します。スクリプトエンコーディングには ASCII 互換エンコーディングを用いることができます。 ASCII 非互換のエンコーディングや、ダミーエンコーディングは用いることができません。

現在のスクリプトエンコーディング__ENCODING__ により取得することができます。

さらに、magic commentによってこのスクリプトエンコーディングを決定でき、それが無い場合の挙動も書かれている。

マジックコメントが指定されなかった場合、コマンド引数 -K, RUBYOPT およびファイルの shebang からスクリプトエンコーディングは以下のように決定されます。左が優先です。

magic comment(最優先) > -K > RUBYOPTの-K > shebang

上のどれもが指定されていない場合、通常のスクリプトなら UTF-8、-e や stdin から実行されたものなら locale がスクリプトエンコーディングになります。 -K オプションが複数指定されていた場合は、後のものが優先されます。

通常のスクリプト-e かどうか、でも変わったりするのね…。そして特に指定ない場合は最終的にはlocaleが使われる。 ので、 LC_CTYPE などでも挙動が変わり得るようだ。

$ ruby -e 's="\x01\x80\x00\x00"; p __ENCODING__, s.encoding, s.rindex("\x01")'
#<Encoding:UTF-8>
#<Encoding:UTF-8>
1
$ ruby -Kn -e 's="\x01\x80\x00\x00"; p __ENCODING__, s.encoding, s.rindex("\x01")'
#<Encoding:ASCII-8BIT>
#<Encoding:ASCII-8BIT>
0
$ RUBYOPT="-Kn" ruby -e 's="\x01\x80\x00\x00"; p __ENCODING__, s.encoding, s.rindex("\x01")'
#<Encoding:ASCII-8BIT>
#<Encoding:ASCII-8BIT>
0
$ LC_CTYPE=C ruby -e 's="\x01\x80\x00\x00"; p __ENCODING__, s.encoding, s.rindex("\x01")'
#<Encoding:US-ASCII>
#<Encoding:ASCII-8BIT>
0

という具合に、何もしないと UTF-8 文字列として扱ってしまうが、オプションや環境変数によって リテラルの Encoding を変えることもできる。

こうやって正しく ASCII-8BIT のEncodingを持つ文字列として検索すれば、 String#rindex の値がズレるということはなさそうだ。

つまり再現条件は

手元の環境で String#rindex の値がズレたのは2つの要因があって

  1. memrchr が使えない環境でビルドしたRubyで実行していて
  2. ASCII-8BIT でないEncodingの文字列に対して rindex をかけていた

ということになる。2つ揃っていないと起きないものと思われる。

Rooの問題

前述の、Rooでエラーが起きる件について考える。

そもそもEncodingが不明で渡ってくるものに対して安易に rindex を使うものではない、ということで rubyzip 側に落ち度がありそうではある。が 現在の master branch では リリースされている 2.3.2 とは大きく変わっていて、もう同じことは起きないのかもしれない。詳しくは追っていないので分からないが、今回のと関連するissueがあった。

既にCloseされているが、今回調べた Zip::CentralDirectory は直接使用するものではなく Zip::File を使ってください、ということのようだ。

つまりこの場合、 Roo 側の使い方が悪い、ということになりそう。 で、Rooの方も調べてみると ちゃんとそれに対応しようとしていると思われる修正があった。

これが取り込まれれば手元の環境でも問題なく動くようになるかな…?

もしくは現時点でも、使う側が渡す引数のEncodingを明示的に指定してやることで一応回避できそうではある。

require 'roo'

File.open("hoge.xlsx", "r:ASCII-8BIT") do |f|
  p Roo::Excelx.new(f)
end

Rubyのバグではないの?

しかし Encodingを正しく指定していなかったために起きているとしても、 String#indexString#rindex も検索した対象の開始indexを返すものであるはずなのだから、それがズレるのはおかしいのではないのか…?という気はする。

さらに memrchr が使えるか否かのビルド環境によってだけで結果が変わる(ちゃんと確かめてないけど おそらくそう…)、というのも気持ち悪い。

かなり昔から現在の実装になっているようだし 今からでは挙動を変えづらそうではある。 仕様といってしまえばそれまでだけど…。 このあたりはRuby開発陣の方々の見解も聞いてみたいところではあります。

3.2

ちなみに Ruby 3.2 からは String#byteindexString#byterindex が追加されるそうです。今回のようなバイナリデータに対する検索にピッタリ使えそうですね。

spherical linear interpolation(slerp)によるlatent spaceでのnoise補間

memo.sugyan.com

の記事を書いてから、先行事例の調査が足りていなかったなと反省。 Latent Seed の Gaussian noise 間での morphing はあんまりやっている人いないんじゃないかな、と書いたけど、検索してみると普通に居た。

Stable Diffusion が公開されるよりも前の話だった…。

そしてこの morphing を作るためのコードも gist で公開されている

読んでみると、2 つの noise の間の値を得る方法として slerp という関数が使われている。

        for i, t in enumerate(np.linspace(0, 1, num_steps)):
            init = slerp(float(t), init1, init2)

https://gist.github.com/karpathy/00103b0037c5aaea32fe1da1af553355#file-stablediffusionwalk-py-L179-L180

この定義は以下のようになっている。

def slerp(t, v0, v1, DOT_THRESHOLD=0.9995):
    """ helper function to spherically interpolate two arrays v1 v2 """

    if not isinstance(v0, np.ndarray):
        inputs_are_torch = True
        input_device = v0.device
        v0 = v0.cpu().numpy()
        v1 = v1.cpu().numpy()

    dot = np.sum(v0 * v1 / (np.linalg.norm(v0) * np.linalg.norm(v1)))
    if np.abs(dot) > DOT_THRESHOLD:
        v2 = (1 - t) * v0 + t * v1
    else:
        theta_0 = np.arccos(dot)
        sin_theta_0 = np.sin(theta_0)
        theta_t = theta_0 * t
        sin_theta_t = np.sin(theta_t)
        s0 = np.sin(theta_0 - theta_t) / sin_theta_0
        s1 = sin_theta_t / sin_theta_0
        v2 = s0 * v0 + s1 * v1

    if inputs_are_torch:
        v2 = torch.from_numpy(v2).to(input_device)

    return v2

https://gist.github.com/karpathy/00103b0037c5aaea32fe1da1af553355#file-stablediffusionwalk-py-L101-L125

不勉強で知らなかったが、これは "spherical linear interpolation"、日本語では「球面線形補間」と呼ばれるのかな、という手法で、3DCGの世界などではクォータニオン(Quaternion)の補間としてよく使われていたりする、ようだ。

en.wikipedia.org

書かれている通り、元々はクォータニオンの補間の文脈で紹介されていたが、次元数に関係なく適用することができるということらしい。 空間上での原点から2つの点へのベクトルを考え、その2つがなす角  \Omega をドット積とノルムを使って求めることができる。 始点から終点まで、 0 から  \Omega へと角度を線形に変化させながら2点を結ぶ円弧上を移動させていく、という感じだろうか。

 \displaystyle Slerp(p_0,p_1;t) = \frac{\sin{[(1-t)\Omega]}}{\sin{\Omega}}p_0 + \frac{\sin{[t\Omega]}}{\sin{\Omega}}p_1

画像生成モデルでの latent space における補間としては 2016年くらいには提案されていたようだ。 そういえば GAN で遊んでいたときにちょっとそういう話題を聞いたことがあったような気もする…。

arxiv.org

しかしこれによって Gaussian noise として分散を保ったまま変化させられる、ことになるのか…? 数学的なことはちょっとよく分からない。 まぁ球面上を角度だけ変えて動いていると考えればノルムが変化しないから大丈夫なのかな??くらいの感覚でしかない…。

実際にこれを使って補間していったときにどのような変化になるか、前回の記事で書いていた sqrt を使うものと上述の slerp を使うもので ランダムな2点を繋ぐ軌跡を見てみた。2次元、3次元 空間上でそれぞれプロットすると以下のような違いがあった。

間隔が異なるだけでなく、明らかに軌道も異なるところを通るようになるようだ。

実装

引用したコードでも良いが、自分でもPyTorch前提で Slerp を書いてみる。 NumPy への変換はしなくても計算できそうだった。

def myslerp(
    t: float, v0: torch.Tensor, v1: torch.Tensor, DOT_THRESHOLD: float = 0.9995
) -> torch.Tensor:
    u0 = v0 / v0.norm()
    u1 = v1 / v1.norm()
    dot = (u0 * u1).sum()
    if dot.abs() > DOT_THRESHOLD:
        return (1 - t) * v0 + t * v1
    omega = dot.acos()
    return (((1.0 - t) * omega).sin() * v0 + (t * omega).sin() * v1) / omega.sin()

前述の数式の定義通りに実装するだけではあるが、やはり特別な場合に注意する必要はある。2点の原点からの方向がまったく同じか、もしくは正反対の方向のとき、理論上では単位ベクトルのドット積は 1-1 になる。 が、float値の多次元配列で計算してみると 多少の誤差が生じることがある。

for _ in range(10):
    v0 = np.random.normal([1, 4, 64, 64])
    v1 = 1.0 * v0
    print((v0 * v1 / (np.linalg.norm(v0) * np.linalg.norm(v1))).sum())
0.9999999999999997
1.0
1.0
0.9999999999999998
0.9999999999999999
1.0
1.0
0.9999999999999998
1.0
1.0000000000000002
for _ in range(10):
    v0 = np.random.normal([1, 4, 64, 64])
    v1 = -1.0 * v0
    print((v0 * v1 / (np.linalg.norm(v0) * np.linalg.norm(v1))).sum())
-1.0000000000000004
-1.0
-0.9999999999999998
-1.0
-1.0
-1.0
-0.9999999999999999
-1.0
-1.0000000000000002
-1.0

slerp では このドット積の値から arccos() で角度を求める必要があるが、例えば np.arccos() は 引数が -1.01.0 の間に収まっていない場合に nan を返すことになってしまう。 また、後でこの角度の sin() で割る処理があるので、ドット積が 1.0 で返せていたとしてもその後の計算で今度は inf が出てきたりする。 ので、こういうケースでは球面ではなく普通に線形で繋いでしまった方が良いだろうね、ということで DOT_THRESHOLD というものが使われて処理を分岐させているようだ。


※追記

指摘いただいたが、まったく同じ向きのときはともかく 正反対の方向の場合は線形補間が正しいとは限らない。 例えば2次元において (1, 0)(-1, 0) を補間する場合は原点を通らず (0, 1) もしくは (0, -1) を通る円軌道になった方が良かったり 3次元の球でいうと北極と南極を繋ぐのは赤道上を通る球面軌道であって欲しい、など。 これは2点のどちらかを適度にブレさせれば算出できる場合もあるかもしれないし、上記のように明らかに通るべき点があればまず両点からそこにまず補間してやるなど、色々な考え方がありそう。これもまた用途や場合によると思われる。 今回の Stable Diffusion で使っている乱数においては正反対向きはまず有り得ないものと思ってよさそうだ。


この小数点誤差がどれくらいになるのかは、要素数だったり 計算の順番だったり numpy で計算するか torch で計算するか などによっても変わってくるようだ。 どれくらいの値になったら線形と見做すか、の閾値(上述では 0.9995)は場合によると思われる。 少なくとも Stable Diffusion で使う [4, 64, 64] の shape で torch.randn() で生成される乱数同士の場合は、通常はドット積は 0.1 にすら届くことがまずないので、極端な話ではそれくらい低くても問題なさそうではある。 morphing の2点間で同一の noise を使う可能性がない、という前提であれば分岐なくしても問題ないくらい。

比較

で、実際に slerp を使うように変更すると morphing はどう変化するのか。 以前に載せた anime girl での morphing の slerp 版を作ってみた。

並べて比較

ちょっと変化のタイミングが変わったかな…? という程度で 大きな差は感じられない。

考察

結局のところ、Stable Diffusion においては Gaussian noise である限り何らかの画像を生成できるようになっているし、指定の2点の間をどう遷移しようと違いはないのかもしれない。 ただ slerp は線形に角度変化していくという観点でも、その遷移の移動量としてはより安定した変化が期待できるかな、という気はする。 むしろ前回の記事の手法が適当に思い付きでやった割にはそこそこ近いものが出来ていたのすごかったのでは…? というくらい。

また、今回は Gaussian noise 間でということを考えていたが slerp の手法自体は別にどこでも使えるものだとは思うので知っておいて損はないはず。また StyleGAN とかで遊ぶことがあったらそこでも使ってみても良さそう。

あと、prompt の embedding 結果に対する morphing でも slerp を使うことはできるはず、だが どうなんだろう。実験してみていないけど、A から B を遷移するのに まったく無関係な C を通過することになってしまったり 余計な変化が増えてしまうだけのような気もするので、こちらは線形に最短距離で遷移して刻み幅だけ調整するくらいが良いのではないかな、と思っている。実際にはやはりどちらでも大して変わらないかもしれないし prompt 次第かもしれない。試行錯誤するときの選択肢として持っておくと良さそうではある。

Stable Diffusionでmorphing

ということで少し触って遊んでみたのでメモ。

Stable Diffusion をザックリ理解

先月公開された Stable Diffusion。

stability.ai

高精度で美しい画像を出力できる高性能なモデルながら、Google Colab などでも手軽に動かせるし、 Apple silicon でもそれなりに動かせる、というのが魅力だ。

中身については 以下の記事の "How does Stable Diffusion work?" 以降のところが分かりやすい。

huggingface.co

図をそのまま引用させていただくと

という仕組みになっていて、受け取る入力は "User Prompt" と "Latent Seed" のみ。 前者が「どのような画像を生成するか」を決めて、後者で「どんなバリエーションでその画像を生成するか」を決めるような感じ。

User Prompt は [77, 768] の空間にエンコードされて、これを使って [4, 64, 64] の Gaussian noise を scheduler によって繰り返し denoise していくことで目標の画像のための "latent image representations" を生成していく。最後はこれを VAE(Variational Auto-Encoder) で拡大していくことで最終的な [512, 512, 3] の画像を得る。 この途中の scheduler による sampling がアルゴリズムによって結果の質が変わってきたりするし、回数が少なすぎると denoise が足りなくて汚い出力になったりする、というわけだ。

ともかく、重要なのはこの 2 つの入力だけで出力が決まるということ、そして prompt の入力は [77, 768] に embedding されたものが使われる、ということ。 prompt の文字列を工夫していくのも良いが、そこから embedding されて渡すものを直接指定してしまっても良いわけだ。 また、 Latent Seed の noise の方も少しずつ変えていくことで少しだけ違う感じの出力を得たりすることができそうだ。

自分は以前に StyeleGAN で latent space を変化させて生成画像の morphing などをやっていたので、それと同じようなことをやってみることにした。

Prompt 間での interpolation, morphing

まずは 2 つの異なる prompt から生成される画像の間を補間して繋いでみる。

(当然ながら、京都と東京の町並みを補間するからといって愛知や静岡の風景が生成されたりはしない。)

scripts/txt2img.py の中にある、prompt から画像を生成する部分のメインはここにある。

github.com

    c = model.get_learned_conditioning(prompts)
    shape = [opt.C, opt.H // opt.f, opt.W // opt.f]
    samples_ddim, _ = sampler.sample(
        S=opt.ddim_steps,
        conditioning=c,
        batch_size=opt.n_samples,
        shape=shape,
        verbose=False,
        unconditional_guidance_scale=opt.scale,
        unconditional_conditioning=uc,
        eta=opt.ddim_eta,
        x_T=start_code
    )
    x_samples_ddim = model.decode_first_stage(samples_ddim)

model.get_learned_conditioning() で与えられた prompt 文字列から [N, 77, 768]Tensor に embedding されたものが得られる。 これを sampler (ここではデフォルトで使われる DDIMSampler を使用している) に与えて ddim_steps 回数の sampling を実行して [N, 4, 64, 64] の denoise された結果が得られる。これを model.decode_first_stage() に与えることで最終的な画像に使われる値が得られるようだ。

sampler.sample() には色々なパラメータがあるが、とにかく重要なのは conditioning (c) と x_T (start_code) だけ。これを変化させることで生成画像をコントロールしていく。

start_code の方はここでは固定した値を使うようにすることで、「何を描くか」だけを徐々に変化させていく様子を作れる。一度だけ乱数を生成してそれを繰り返し使うようにすると良い。

で、 c の方は 2 つの異なる prompt からそれぞれ embedding された値を取り出して、線形に変化させていく。

指定した cstart_code から生成画像だけを得るような関数を書いておくとやりやすい。 model のロード方法などについては割愛。

from contextlib import nullcontext

import numpy as np
import torch
from PIL import Image

from ldm.models.diffusion.ddim import DDIMSampler
from ldm.models.diffusion.ddpm import LatentDiffusion


def get_device() -> torch.device:
    ...


def load_model() -> LatentDiffusion:
    ...


model = load_model(...)


def generate(
    c: torch.Tensor, start_code: torch.Tensor, ddim_steps: int = 50
) -> Image:
    batch_size = 1
    device = get_device()
    precision_scope = torch.autocast if device.type == "cuda" else nullcontext
    with torch.no_grad():
        with precision_scope("cuda"):
            with model.ema_scope():
                uc = model.get_learned_conditioning(batch_size * [""])
                shape = [4, 64, 64]
                samples_ddim, _ = DDIMSampler(model).sample(
                    S=ddim_steps,
                    conditioning=c,
                    batch_size=batch_size,
                    shape=shape,
                    verbose=False,
                    unconditional_guidance_scale=7.5,
                    unconditional_conditioning=uc,
                    eta=0.0,
                    x_T=start_code,
                )

                x_samples_ddim = model.decode_first_stage(samples_ddim)
                x_samples_ddim = torch.clamp(
                    (x_samples_ddim + 1.0) / 2.0, min=0.0, max=1.0
                )
                image = (
                    255.0 * x_samples_ddim.cpu().permute(0, 2, 3, 1).numpy()[0]
                ).astype(np.uint8)
                return Image.fromarray(image)

ちなみに、 sampler.sample() された結果の方を線形に繋いで変化させていくという手法もあるのだけど、これはもはやどんな画像が生成されるかほぼ決定された後の値なので、 morphing しても単なる画像合成のような感じにしかならなくて面白くはない。

指定した cstart_code で画像を生成する準備ができたら、あとはその入力を作っていくだけ。

def morph_prompts(prompts: Tuple[str, str], steps: int) -> None:
    start_code = torch.randn([1, 4, 64, 64], device=get_device())
    c0 = model.get_learned_conditioning(prompts[0])
    c1 = model.get_learned_conditioning(prompts[1])
    for i in range(steps + 1):
        x = i / steps
        c = c0 * (1.0 - x) + c1 * x
        img = generate(c, start_code)
        img.save(f"morphing_{i:03d}.png")

これらを繋げてアニメーションさせれば、2 つの異なる prompt 間の morphing が出来上がる。

ただ、似ているものならまだあまり違和感ないが あまりに異なる 2 つを morphing させようとすると、急激に変化してしまって面白くない。

embedding された空間がどんなものかは未知だが、ともかく A と B の 2 点間には必ず「denoise された結果 A になるもの」と「denoise された結果 B になるもの」が分断される地点がどこかに存在してしまう。 それは A B の中心かもしれないし、少しズレたところかもしれないが、そのあたりで急激な変化が起こり得る。 ので、中点に近い位置は出来るだけ細かい step で刻んだ方がよりシームレスな morphing になりやすいように感じた。 ので、単純な線形に繋ぐのではなく双曲線関数で刻み幅を微妙に変えながら作ってみることにした。

    a = np.arccosh(5.0)
    for i in range(steps + 1):
        t = i / steps
        x = sinh(a * (t * 2.0 - 1.0)) / sinh(a) / 2.0 + 0.5
        c = c0 * (1.0 - x) + c1 * x
        ...

それでもやっぱり急激な変化は捉えられないことが多々あるけれども…。

Seed 間での interpolation, morphing

今度は、同一の prompt で異なる Latent Seed を使用した 2つの画像間での morphing。

prompt の方は固定して、 torch.randn() で生成していた Gaussian noise の方を徐々に変えていく。生成する画像の「お題」は一緒だが、違うバリエーションのものになっていく、という morphing。

prompt のときと同じように変化させていけば良いだけ、と思ったが そうはいかない。実際やってみると中間点あたりはボヤけた画像になってしまうようだ。 最初 何故だろう…?と思ったが どうやらこの noise は "Gaussian noise" であることが重要で、標準正規分布として  {\mu} = 0, {\sigma}^2 = 1 になっていなければならない、ということらしい。 単純に v0 * (1.0 - x) + v1 * x のように単純な線形結合で変化させていくと、中心に近づくにつれてその標準偏差は小さくなってしまう。

それを防ぐために、足し合わせる前にそれぞれの倍率の sqrt をとるようにすると、合成された noise は標準偏差を保持したまま遷移することができそうだ。

すると今度は 0 付近と 1 付近で急激な変化が起こりやすそうなので、 prompt morphing のときのように刻み幅を調節する。

def morph_noises(prompt: str, steps: int) -> None:
    c = model.get_learned_conditioning([prompt])
    n0 = torch.randn([1, 4, 64, 64], device=get_device())
    n1 = torch.randn([1, 4, 64, 64], device=get_device())
    for i in range(steps + 1):
        t = i / steps
        x = 2.0 * t**2 if t < 0.5 else 1.0 - 2.0 * (1.0 - t) ** 2
        start_code = n0 * math.sqrt(1.0 - x) + n1 * math.sqrt(x)
        img = generate(c, start_code)
        img.save(f"morphing_{i:03d}.png")

これで、同じお題(prompt)に対して複数のバリエーションで描かれたものを連続的に変化させていくことができる。 好みの画像を出力する seed を幾つかピックアップして繋いでみたりするとより好みのものが見つかるかもしれないし、意外とブレンドされたものは好みではないものになるかもしれない。


※追記

memo.sugyan.com


まとめ

以上の2つができれば、その応用として prompt と noise を同時に変化させていったり、交互に変化させていったり、数回変化させた後にまた元の画像に戻ってきたり、といったものも作っていける。

雑に試行錯誤しながら実行できるように Google Colab でscriptを書いていたけど、ちょっと整理して公開する予定(需要あるかどうか分からないけど)。 prompt 間の morphing はやっている人結構いるけど、noise 間のものはまだあんまり見かけないような気はする?

※追記: GitHub - sugyan/stable-diffusion-morphing でとりあえず公開しておきました。多分動くはず?

Slackの家庭内日記をはてなブログに同期していく

夫婦で日記の習慣を始めたのは3年ほど前。簡単な記録でもつけておくと後で色々見返せて良いよね、と。

とはいえなかなか習慣をつけるのは難しいなと思って

  • ブログサービスは最も日付と関連づけて文書が保存されるし良い
  • しかしブログ投稿、となるとちょっと気合いが要る
  • SNS感覚で気軽に投稿できると良い
  • が、下手なSNS誤爆したりすると大変

といったものがあり、色々試行錯誤もした結果 家庭内Slackで日記用channelを作り、1日1投稿ずつ書き込む、という運用にしたのだった。

  • 誤爆の心配は少ない (少なくとも妻は家庭以外でSlackを使わないし)
  • 専用channelで投げていくだけなので他の話題と混ざったりもしない
  • 気軽に書けて編集なども簡単
  • 検索もそれなりにできて過去の記録を見返しやすい

ということで気に入っていた。

ところが90日以上前のものが見られなくなるということで移行を検討した。 べつに課金すれば良いだけではあるのだけど、せっかくならバックアップも兼ねて別のブログに同期しておいたら良いのでは、と。

ので、過去の日記投稿たちをSlackから取得してブログ記事として投稿していく、というものを作った。

なんとなく書きやすいのはPythonかな?と思ってPythonで。ライブラリも充実していて数十行程度でそれぞれのロジックが書けてラクだったので正解だったと思う。

APIでの取得と投稿

SlackのAPIからchannelの履歴を遡り取得していく。 text要素だけ取ると絵文字が :bow: などで表示されてしまうので rich_text_section から復元していってSlackに書かれた通りの文字列でとれるようにした。

同期先のブログははてなブログにした(自分が使っていて多少は勝手が分かるので)。自分と妻のアカウントだけが見られる限定公開にして、そこに AtomPubAPIで投稿していく。 updated を指定すれば日時を指定できるので過去のものも過去の日記として記録されていく。

同期をするタイミングは確定していないけど、複数回実行するし 一度投稿したものも編集する可能性はゼロではないので、既に記事があるものは PUT で更新、まだ無いものは POST で新規投稿する。が 投稿しようとしているものが既に記事あるかどうかを調べるのが難しく、ドキュメントには entry_id がepochと書いてあるけどそうでもない… まぁ投稿日時が一致するようにするから updated に使う日時と同じ epoch の記事が既にあるかどうかを調べれば良いか、ということでそうした。

はてなブログの制限

いざ過去の日記をどんどん記事として投稿していこうとしたら50日ぶんくらいで止まってしまった。 知らなかったけどどうやらはてなブログには「1日100件まで」という上限があるらしい…。 まぁ3年分なら毎日コツコツやっていけば 365 * 3 * 2人 = 2190件 だから約21日で終わるはず。ギリギリ8月中には。。

もっと急ぎで大量の記事を移したかったら インポートできる形式のファイル を作るとかすることになるんだろうか。

定期実行

過去のはボチボチやっていくとして、今後これから書くものをどうするか。 すっかり3年の習慣になっているので今からブログに全面移行するつもりもなくて、日記は今まで通りにSlackで書き続けて、定期的に直近の日記をブログに同期していけば良い、ということにした。

Cloud Functionsにdeployして Cloud Schedulerから定期実行して Pub/Sub経由で自動で定期実行してくれるように設定。これであとは勝手に同期されるようになってくれるはず。

書くときはSlackに書き続け、90日以上前のものは調べたり検索したりできなくなるけど そのときはブログの方で検索すればもっと昔まで遡ることができる、という運用。

まとめ

まだ完全に同期に成功しているわけではないけど、この運用でしばらくやってみる予定。そもそも家庭内の日記がどこまで続くか…という問題もあるけど。まぁ夫婦ともに3年くらい続いているわけだし しばらくは続くんじゃないかな? 子育てで大変な日々も数年後には笑って読み返せるといいな

Rustで将棋棋譜変換ライブラリを作った

というわけで作った記録

将棋棋譜の形式色々

プロ対局、コンピュータ将棋、対戦アプリ、…と様々な将棋対局についてのデータが世の中には存在しているが、その形式は様々だったりする。

KIF

将棋ソフト「柿木将棋」 で使われている形式だそうで、現在でも多くの将棋ソフトなどで使われているようだ。プロ対局の中継サイトなどでもこの形式が使われている。

例:

# ---- Kifu for Windows95 V3.53 棋譜ファイル ----
開始日時:1999/07/15(木) 19:07:12
終了日時:1999/07/15(木) 19:07:17
手合割:平手
先手:先手の対局者名
後手:後手の対局者名
手数----指手---------消費時間-- # この行は、なくてもいい
1 7六歩(77) ( 0:16/00:00:16)
2 3四歩(33) ( 0:00/00:00:00)
3 中断 ( 0:03/ 0:00:19)

初期盤面や指し手など日本語で記述されていて、人間が読むぶんには直感的に分かりやすい。

が、コンピュータで解析するにはなかなか厄介…。 書式の定義については以下のリンクだけが存在している状態のようで、これを頼りに解析していくしかない。しかし「省略できる」「任意のものを追加できる」など なかなか自由な項目が多くて対応するのは大変。

優れている点としては、上記リンクでは言及されていないが「指し手の変化・分岐が表現できる」、という点で、枝分かれしたものも含めてすべて記録として残せるようになっているようだ。自分は柿木将棋などのシリーズ製品を持っていないけど、おそらくこれらがちゃんと分岐を含む棋譜も読み取って再現できるようになっていると思われる。

CSA

これは Computer Shogi Association で定めている形式のようで、KIFと比較すると圧倒的にコンピュータで解析しやすい形になっている。

例:

'----------棋譜ファイルの例"example.csa"-----------------
'バージョン
V2.2
'対局者名
N+NAKAHARA
N-YONENAGA
'棋譜情報
'棋戦名
$EVENT:13th World Computer Shogi Championship
'対局場所
$SITE:KAZUSA ARC
'開始日時
$START_TIME:2003/05/03 10:30:00
'終了日時
$END_TIME:2003/05/03 11:11:05
'持ち時間:25分、切れ負け
$TIME_LIMIT:00:25+00
'戦型:矢倉
$OPENING:YAGURA
'平手の局面
P1-KY-KE-GI-KI-OU-KI-GI-KE-KY
P2 * -HI *  *  *  *  * -KA * 
P3-FU-FU-FU-FU-FU-FU-FU-FU-FU
P4 *  *  *  *  *  *  *  *  * 
P5 *  *  *  *  *  *  *  *  * 
P6 *  *  *  *  *  *  *  *  * 
P7+FU+FU+FU+FU+FU+FU+FU+FU+FU
P8 * +KA *  *  *  *  * +HI * 
P9+KY+KE+GI+KI+OU+KI+GI+KE+KY
'先手番
+
'指し手と消費時間
+2726FU
T12
-3334FU
T6
%CHUDAN
'---------------------------------------------------------

主にASCII文字で表現されるようになっていて 駒種は2文字、場所を示すのも2桁の数字、といった一定の決まりがあるので解析しやすい。 これと同様の表現を用いたCSAサーバプロトコルがあり、世界コンピュータ将棋選手権などではこれが使われているらしい。

USI/SFEN

これは棋譜形式とはちょっと違うかもしれないけど、 将棋GUIソフト「将棋所」 で使われて広く普及しているUSIという通信プロトコルがある。 GUIと思考エンジンの間の通信をするために局面情報や指し手の情報を SFEN(Shogi Forsyth-Edwards Notation)という形式で表現していて、これによって棋譜を再現できる。

例:

position sfen lnsgkgsnl/9/ppppppppp/9/9/9/PPPPPPPPP/1B5R1/LNSGKGSNL w - 1 moves 5a6b 7g7f 3a3b

文字表記

そもそもコンピュータで扱われる前は文字情報だけで記録が残っているはずで、これは指し手をそれぞれ数字と駒種と相対位置や動作で表現している。KIF形式はこれらの表現をカスタマイズした形式といえる。

76歩
52銀右上成

その他

その他にも幾つかの形式があるようで、どれがどのくらい普及しているのかよく分からない…。 おそらく多くはKIFかCSAが使われていると思われる。CSAはむしろコンピュータ将棋の人にしか使われていないか…?

GUIソフトや将棋ソフトの開発者はそれぞれこれらの形式を読み書きできるように対応していくしかない、という状況のようだ…

JSON棋譜フォーマット(JKF, json-kifu-format)

上述した棋譜形式はそれぞれ読みやすさや表現できる情報が違っているわけだが、実際のところ指し手には多くの情報が必要で、単独の指し手の表記だけを取り上げたときに「どの駒が」「どこから動いて」「どうなったか」などの情報をすべて知ることは出来ない。直前の指し手や現在の盤面の情報がないと分からないものが多かったりする。

そこで、json-kifu-formatという形式が提案されて公開されている。

github.com

READMEに書かれている通り、各指し手について必要な情報をすべて含めたJSON形式で表現し、また分岐にも対応している。

例:

{
  "header": {
    "先手": "na2hiro",
    "後手": "うひょ"
  },
  "moves": [
    {},
    {"move":{"from":{"x":7,"y":7},"to":{"x":7,"y":6},"color":0,"piece":"FU"}},
    {"move":{"from":{"x":3,"y":3},"to":{"x":3,"y":4},"color":1,"piece":"FU"}},
    {"move":{"from":{"x":8,"y":8},"to":{"x":2,"y":2},"color":0,"piece":"KA","capture":"KA","promote":false}},
    {"move":{"from":{"x":3,"y":1},"to":{"x":2,"y":2},"color":1,"piece":"GI","capture":"KA","same":true}},
    {"move":{"to":{"x":4,"y":5},"color":0,"piece":"KA"}},

    {"special": "CHUDAN"}
  ]
}

さらにこの形式の提案だけでなく、KIFやCSAからJKF形式に変換するためのparserと、JKFを再生するためのプレイヤーも同梱している。

実装としては それぞれの棋譜形式を一度parseしたのち、初期局面から一手ずつ実際に進めて足りない情報を補っていくことで最終的なJKFを生成しているようだ。

自作parser, converter

上記 json-kifu-parser はJavaScriptで実装されているが、自分は最近はRustで将棋関連のものを作っていてRustでも同様のものが欲しいと思った。

具体的には、自作の詰将棋ソルバ に色々な問題を食わせてみようと思ったときにKIFファイルが多かったり 自分で書くにはCSAで書いたりで色々な入力に対応したいというのがあった。

CSAのparseには csa crateが既に存在していて問題ないが、KIFは自力でやるしかなく、いちおう自前である程度のparseできるところまで作っていたが、折角だから分岐もちゃんと対応してJKFに互換のデータを作れるParserを作ってみよう、という気持ちになった。 ので作ることにした。

github.com

KIFのparseには csa と同様に nom というparser combinators libraryを使って頑張って実装した。optmany0 が多くてなかなか大変だったが、どうにか様々なケースに対応できたと思う。 パフォーマンス的にちょっと問題があるかもしれないけど… 大量の棋譜を食わせて高速に処理をしなければならない、という場面でない限りは大丈夫か。

分岐は一度最後まで読み込んでstackに積んでおいてから後で forks に詰め替えていく必要があってなかなか頭を使った。

CSAのparseはcsa crateに任せてparse結果からJKFに変換する処理だけを書いて作ったが、csa crateは指し手のコメントなどを捨ててしまっていて復元できなかったりもあるので、こちらも自分で実装しなおしても良いのかもしれない。

それぞれparseして結果を JsonKifuFormat に詰めたら、足りない情報を補うために normalize() を呼ぶ。 この中で実際に局面を動かして確かめる必要があるので内部で shogi_core crateを使った。forks再帰的に処理していく。

また、 relative 情報は他の同種駒と比較して相対的なものを調べる必要があり、計算が大変めんどう…なので shogi_official_kifu crateで表記結果から取り出して詰めることでお茶を濁した。本来ならこの棋譜表記を作成するための情報なので本末転倒な感もあるが…。そもそも relative はそういった棋譜表記をするときくらいにしか使わない(棋譜変換などでは必要ない)ので計算もoptionalで良いのかもしれない。

一度 JsonKifuFormat が作れれば、そこからKIF形式やCSA形式に書き出すのは比較的容易にできる。(KIFの分岐はまだ対応できていないが…)。 また shogi_core::Position にも変換できるし shogi_core::ToUsi traitを利用してUSIで書き出すのもやりやすい。

ので、もし対応する棋譜形式を増やそうと思ったらそれぞれ JsonKifuFormat と相互変換するためのparser/converterだけを用意していけば良い、ということになる。 記事冒頭に載せた図の通りで、変換の中継のための形としてとても良いデータ型になっていると思う。

クレート公開

せっかく便利ライブラリとして実装できたので、(需要はまったく無いかもしれないが)crates.ioに公開してみることにした。 できるだけドキュメントをちゃんと書いて、初の cargo publish

docs.rs

Rustで将棋棋譜を扱いたい方が居たら使ってみていただけると嬉しいです。

WASM化してブラウザ上でclient sideだけで棋譜変換できるWebApp、なら多少は使いどころあるかな…?