数式微分器の作成 と D言語での木走査
モチベーション
無性に微分がしたくなった。
結論
出力が頭悪い。
D言語と振り分け式の外部訪問器は相性がいい。
特に数学的におもしろい話はありません。
はじめに
偏導関数値を求める方法には、数値微分 や (トップダウン/ボトムアップ)型自動微分などがありますが、今回は具体的な値を求めたいわけではないので数式微分を行います。
方法
式木で表される関数を、再帰的に微分することで偏微分関数を表す式木を作ります。
まず、ある程度簡単な関数であれば、三角関数のようなより単純な関数の組み合わせや四則演算で表すことが出来ます。
それらの 関数の組み合わせ に対する微分結果の式木を求めることで微分を行います。
例えば、 という関数の場合、
というふうに偏導関数を求めることが出来ます。以下に幾つかの関数の組み合わせに対応する微分結果を挙げます。
微分前 | 微分結果 |
---|---|
加算: | |
乗算: | |
正弦: | |
余弦: | |
冪乗: | |
対数: |
このような規則に従い式木を走査し組み直すプログラムをD言語で作成しました。
https://github.com/Ryooooooga/DifferentialCulculator
では試しに を微分してみます。
出力は以下の様になります。
((((-1)*sin(x))*sin(x)^{(-1)})+(cos(x)*(sin(x)^{(-1)}*((cos(x)*(-1))*sin(x)^{(-1)}))))
はい、出力の整形を面倒くさがったため括弧地獄になっています。
これを適当に手動で整形すると、
の様になり、 の微分結果である
と一致することがわかります。
まとめ
一応式木を走査して最低限無駄な項を切り落とすようにはしているのですが、出力結果が見れたものではないので、出力の式木に対してより高度な最適化を施す必要がありそうです。
木構造の走査と戦略
デザインパターンっぽい話です。
式木や構文木のような木構造の走査の戦略にはいくつかのパターンがあり、それぞれにメリット、デメリットが存在します。
下に、よく使われる3つのパターンを挙げます。
また、それぞれに例として簡単なC#のソースコードを載せます。
1. 走査メソッド(Interpreter パターン)
最も直感的なパターンです。
式木のノードを表すクラスが走査メソッドを持ち、再帰的に子ノードのメソッドを呼び出すことで木を走査します。
// ノード interface INode { // 走査メソッド double Eval(); } // 加算 class Add: INode { private INode left, right; public Add(INode left, INode right) { this.left = left; this.right = right; } // 走査メソッドの実装 public double Eval() => this.left.Eval() + this.right.Eval(); } // 乗算 class Mul: INode { private INode left, right; public Mul(INode left, INode right) { this.left = left; this.right = right; } // 走査メソッドの実装 public double Eval() => this.left.Eval() * this.right.Eval(); } // 数値 class Num: INode { private double value; public Num(double value) { this.value = value; } // 走査メソッドの実装 public double Eval() => this.value; } class Prog { public static void Main() { // 1*2 + 3*4 var tree = new Add( new Mul(new Num(1), new Num(2)), new Mul(new Num(3), new Num(4)) ); System.Console.WriteLine(tree.Eval()); } }
実行結果
14
上の例では、「加算」、「乗算」、「数値」というノードクラスがそれぞれ Eval
メソッドを実装することで木の評価を行っています。
メリット
わかりやすい
仕組みが単純で動作がわかりやすいです。
引数をとれる、返り値を返せる
基底クラス(インターフェース)の
Eval
メソッドの定義を変更することで自由な個数の引数をとる様にしたり、返り値の型を変更することが出来ます。クラスの内部構造を隠蔽できる
メソッドにより走査するため、クラス内部の情報を外部に公開する必要がありません。
上の例で言えば、
Num
クラスは内部にどのような値を保持しているのかを外部に公開していません。木クラスの追加が容易
上の例に平方根を表すノード
Sqrt
を追加する場合、Sqrt
のEval
メソッドを適切に実装するだけで平方根を含む式木の評価を行えます。
デメリット
処理が分散する
上では比較的単純な例を示しましたが、非正規化非均質抽象構文木のような複雑になりやすい木構造の場合、クラスの数が数十以上になってしまうこともありえます。その場合、それぞれのクラスに処理が分散してしまうため、変更などが行いにくくなります。
また、処理の追加(例: プリティプリント機能 や 型検査機能、最適化機能など)を行うごとに木クラスの定義が膨れあがっていきます。
単純な木構造の走査に向いたパターンと言えます。
2. 外部木訪問器(Visitorパターン)
非均質抽象構文木などの走査によく用いられるパターンです。
ダブルディスパッチにより、訪問した木クラスの型に応じた訪問器のメソッドに処理を振り分けます。
// 外部訪問器インターフェース interface IVisitor<T> { T Visit(Add node); T Visit(Mul node); T Visit(Num node); } // ノード interface INode { // 訪問器の受け入れ T Accept<T>(IVisitor<T> visitor); } // 加算 class Add: INode { public INode Left { get; } public INode Right { get; } public Add(INode left, INode right) { this.Left = left; this.Right = right; } // IVisitor.Visit(Add) を呼び出す public T Accept<T>(IVisitor<T> visitor) => visitor.Visit(this); } // 乗算 class Mul: INode { public INode Left { get; } public INode Right { get; } public Mul(INode left, INode right) { this.Left = left; this.Right = right; } // IVisitor.Visit(Mul) を呼び出す public T Accept<T>(IVisitor<T> visitor) => visitor.Visit(this); } // 数値 class Num: INode { public double Value { get; } public Num(double value) { this.Value = value; } // IVisitor.Visit(Num) を呼び出す public T Accept<T>(IVisitor<T> visitor) => visitor.Visit(this); } // 評価訪問器 class Evaluator: IVisitor<double> { public double Visit(Add node) => node.Left.Accept(this) + node.Right.Accept(this); public double Visit(Mul node) => node.Left.Accept(this) * node.Right.Accept(this); public double Visit(Num node) => node.Value; } class Prog { public static void Main() { // 1*2 + 3*4 var tree = new Add( new Mul(new Num(1), new Num(2)), new Mul(new Num(3), new Num(4)) ); // 評価器 var evaluator = new Evaluator(); System.Console.WriteLine(tree.Accept(evaluator)); } }
実行結果
14
注目して頂きたいのは 各木クラスに定義された Accept
メソッドです。
全く同じ定義の関数がそれぞれのクラスで定義されています。これは一見無駄に見えますが、this
の型によって呼び出すメソッドを解決するため必要になります。
また、式木の評価は訪問器が担当するため、評価に関する処理が木クラスから切り離されていることがわかります。
メリット
木構造と処理を分離できる
処理を訪問器が行うため、木クラスは処理に関する詳細を知る必要がなくなりました。
訪問器が状態を持てる
走査メソッドの場合、走査途中の状態は引数として伝播させる必要がありますが、Visitorパターンの場合は訪問器が処理途中の状態を持つことが出来ます。
機能の追加に強い
機能を追加する際には、新しく処理を行う訪問器を作成すればよいため、既存のコードを変更する必要がありません。
例として、式木を文字列に変換する訪問器
Printer
の定義を示します。
class Printer: IVisitor<string> { public string Visit(Add node) => $"({node.Left.Accept(this)} + {node.Right.Accept(this)})"; public string Visit(Mul node) => $"({node.Left.Accept(this)} * {node.Right.Accept(this)})"; public string Visit(Num node) => $"{node.Value}"; } tree.Accept(new Printer()); // ((1 * 2) + (3 * 4))
デメリット
木の内部情報を訪問器に公開する必要がある
上の例であれば、
Num
クラスの保持している数値に関する情報などの、処理に必要な情報が外部に公開されていることがわかります。引数の型や個数が固定
訪問器によって引数の型や個数を変えることが難しくなります。
更に、上の例の場合 ジェネリクスを用いて返り値の型を変えていますが、C++の場合は (templateと仮想関数の相性が良くないため) 処理の結果を呼び出し元へと伝播させるためには訪問器の内部状態を利用するなどの工夫が必要になります。
処理を行うメソッド名が固定
Interpreterパターンの場合は処理によって適切なメソッド名を付けられますが、上の例の場合は
Visit
に固定されます。
構文木の解析などの比較的複雑な木構造の走査に向いているパターンと言えます。
3. 外部訪問器(型情報などによる振り分け)
メソッドの振り分けを訪問器の内部で行うパターンです。
// ノード interface INode { // Acceptは必要ない } // 加算 class Add: INode { // ... } // 乗算 class Mul: INode { // ... } // 数値 class Num: INode { // ... } // 評価訪問器 class Evaluator { public double Dispatch(INode node) { // C# 6.0 if (node is Add) return this.Eval((Add)node); if (node is Mul) return this.Eval((Mul)node); if (node is Num) return this.Eval((Num)node); // 処理できなかった throw new System.Exception(); } private double Eval(Add node) => this.Dispatch(node.Left) + this.Dispatch(node.Right); private double Eval(Mul node) => this.Dispatch(node.Left) * this.Dispatch(node.Right); private double Eval(Num node) => node.Value; } class Prog { public static void Main() { // 1*2 + 3*4 var tree = new Add( new Mul(new Num(1), new Num(2)), new Mul(new Num(3), new Num(4)) ); // 評価器 var evaluator = new Evaluator(); System.Console.WriteLine(evaluator.Dispatch(tree)); } }
実行結果
14
木クラスの Accept
メソッドがなくなり、代わりに訪問器の Dispatch
メソッドでクラスに応じたメソッドに処理が振り分けられていることがわかります。
mcsが C# 7 に対応していないため if
文を羅列していますが、C# 7 に対応した環境であれば switch
式を用いたパターンマッチングが使える様になるため、よりスマートに記述することが可能になります。
メリット
処理を行うメソッド名を自由に決められる
引数の数や個数が自由
各走査メソッドを
private
に出来る。
デメリット
各訪問器に
Dispatch
メソッドが必要メソッドの実装忘れがコンパイル時に確認できない
Dispatch
メソッド内のノード型に対応した振り分け機構に実装忘れがあっても、コンパイラはそれを検知できません。実行時に例外が発生してようやく実装忘れを確認出来ます。
このパターンの利点は一見少ないように見えますが、今回の数式微分を行うプログラムでは式木の走査にこれを採用しました。
なぜかと言うと、このパターンとD言語の言語機能である mixin / template mixin の相性が非常に良いためです。
それらの言語機能に関する詳細は省きますが、文法的にクリーンなC言語のマクロを想像すればだいたいそんな感じです。
では、D言語での 評価訪問器 の実装例を以下に示します。
// 訪問器ミックスイン template Visitor() { // 振り分けメソッド auto dispatch(string func, Args...)(Node node, Args args) { import std.algorithm: castSwitch; return node.castSwitch!( (Add node) => mixin("this." ~ func ~ "(node, args)"), (Mul node) => mixin("this." ~ func ~ "(node, args)"), (Num node) => mixin("this." ~ func ~ "(node, args)"), )(); } } // 評価訪問器 class Evaluator { mixin Visitor; private double eval(Add node) { return this.dispatch!"eval"(node.left) + this.dispatch!"eval"(node.right); } private double eval(Mul node) { return this.dispatch!"eval"(node.left) * this.dispatch!"eval"(node.right); } private double eval(Num node) { return node.value; } }
この訪問器の振り分けメソッドは、visitor.dispatch!"eval"(node)
のように呼び出します。
上の Visitor
が template mixin と呼ばれる言語機能を使った部分で、mixin Visitor;
と書かれた部分に Visitor
内部の 文 宣言 (この場合は dispatch
メソッド)を展開します。
これによって複数の訪問器を書く場合にも、dispatch
メソッドの定義は一度だけでよく、実装忘れの可能性を低減することができます。
更に、dispatch
のテンプレート引数(上の場合は "eval"
) に適当な文字列を渡せば 振り分け先のメソッドを eval
以外にも変えられる他、node
に続いて 引数を渡すことで 複数の引数を振り分け先のメソッドに渡すことができます。
最後に
D言語はそれなりに良い言語なのでみんなも書けばいいと思います。
C++17 optionalの実装について
C++17で、optionalが標準ライブラリに入ることが決定した。
そして、C++17 optionalはconstexprに対応している。
以前からその実装について気になっていたので調べたついでに少しまとめる。
C++03 Boost.Optionalの実装
C++03時代のBoost.Optionalの実装についても触れておく。
C++03 Boost.Optionalの実装は単純だ。
予めoptionalオブジェクト内にデータを格納するストレージを確保しておき、
適切なタイミングでplacement newをして、初期化する。
また、オブジェクトを破棄する際には、必要があれば明示的デストラクタ呼び出しを行う。
ストレージの確保には、(C++11以降でいうところの)std::aligned_storageのようなものを用いる。
つまり、optionalの内部で領域の動的確保は行われない。
下記のコードはBoost.Optionalの実装のイメージだ(故に正確ではない)。
template <typename T> class optional { // storage typename std::aligned_storage<sizeof(T), alignof(T)>::type storage_; bool init_; // pointer to storage T* ptr() { return static_cast<T*>(static_cast<void*>(&storage_)); } public: // default constructor optional() : init_(false) {} // value constructor optional(const T& val) : init_(true) { // placement new new (&storage_) T(val); } // destructor ~optional() { if (init_) { // explicit destructor call ptr()->~T(); } } };
見て分かる通り、この実装のoptionalクラスはtrivialデストラクタを持たない。
更には引数付きコンストラクタはconstexprコンストラクタの条件を満たさない。
なので、このままではconstexpr対応の条件であるリテラル型クラスの条件を満たさない。
C++17 optionalはどのような実装により、constexpr対応を果たしたのだろうか。
その前に
上のコードでは、Tの型に関わらずoptional<T>型はデストラクタで明示的デストラクタ呼び出しを行っているが、
実際には、Tがtrivialデストラクタをもつ場合は明示的デストラクタ呼び出しを省略できる。
そのため、Tがtrivialデストラクタを持つときには、optional
実装について
さて、ようやくC++17 optionalの実装を確認する。
constexpr対応のoptionalの実装は以下のリポジトリのものを参考にした。
だが、実際のソースコードは一息に読むには長いので、部分部分に分解して見ていく。
まずは、optional<T>型の定義を確認する。
定義は上記のリポジトリ、optional.hppファイルの350行目付近である。
// 350行目付近 template <class T> class optional : private OptionalBase<T> { ... };
optional<T>型はOptionalBase<T>型を継承している。
OptionalBase<T>型はそのすぐ上で定義されている。
// 341行目付近 template <class T> using OptionalBase = typename std::conditional< is_trivially_destructible<T>::value, constexpr_optional_base<typename std::remove_const<T>::type>, optional_base<typename std::remove_const<T>::type> >::type;
std::conditional<B, X, Y>::typeは、B == trueの時 X型になり、B == falseの時 Y型となる。
ようは、条件分岐のようなものだ。
この場合の分岐条件は、
is_trivially_destructible<T>::value
つまり、Tがtrivialデストラクタを持つかどうかだ。
Tがtrivialデストラクタを持つ時、OptionalBase<T>は、constexpr_optional_base<T>となり、
Tが非trivialデストラクタを持つ時、OptionalBase<T>は、optional_base<T>となる。
双方の定義を確認する。
// 296行目付近 template <class T> struct optional_base { bool init_; storage_t<T> storage_; ... ~optional_base() { if (init_) storage_.value_.T::~T(); } }; // 319行目付近 template <class T> struct constexpr_optional_base { bool init_; constexpr_storage_t<T> storage_; ... ~constexpr_optional_base() = default; };
「...」で省略した部分はどちらのクラスもほぼ同じであった。
違いは2点。
1つ目は、storage_ メンバ変数の型、
2つ目は、デストラクタでのメンバの明示的デストラクタ呼び出しの有無だ。
constexpr_optional_baseは、Tがtrivialデストラクタを持つため、明示的デストラクタ呼び出しを省略しているのだろうと予想できる。
では、storage_t と constexpr_storage_t の定義を確認する。
// 266行目付近 template <class T> union storage_t { unsigned char dummy_; T value_; ... ~storage_t(){} }; // 281行目付近 template <class T> union constexpr_storage_t { unsigned char dummy_; T value_; ... ~constexpr_storage_t() = default; };
このoptionalは、ストレージにstd::aligned_storageのような型を使わず、共用体を使って実現しているようだ。
storage_t と constexpr_storage_t の定義の違いは、デストラクタの定義のみである。
ではこの違いは何故必要なのか。
C++11以降では非trivialデストラクタを持つオブジェクトも共用体のメンバにすることが可能になった。
しかし、その時その共用体のデストラクタは暗黙的にdeleteされる。
そのため、非trivialデストラクタを持つT型のメンバ value_ を持つstorage_tには空のデストラクタ定義が必要になる。
さらに、共用体のコンストラクタの初期化子を使って value_ を初期化すれば、T型のコンストラクタが呼び出されるため、
optional<T>のコンストラクタ内での placement new の必要がなくなる。
以上をもって、Tがリテラル型である場合に、constexpr_storage_t<T>、constexpr_optional_base<T>、optional<T>の全てがリテラル型の条件を満たし、コンパイル時のオブジェクト構築が可能となる。
#include <type_traits> #include <string> #include <cassert> template <typename T> union storage_t { char dummy_; T value_; template <typename... Args> constexpr storage_t(Args&&... args) : value_(args...) {} ~storage_t() {} }; template <typename T> union constexpr_storage_t { char dummy_; T value_; template <typename... Args> constexpr constexpr_storage_t(Args&&... args) : value_(args...) {} ~constexpr_storage_t() =default; }; template <typename T> struct optional_base { bool init_; storage_t<T> storage_; template <typename... Args> constexpr optional_base(Args&&... args) : init_(true), storage_(args...) {} ~optional_base() { if (init_) storage_.value_.~T(); } }; template <typename T> struct constexpr_optional_base { bool init_; constexpr_storage_t<T> storage_; template <typename... Args> constexpr constexpr_optional_base(Args&&... args) : init_(true), storage_(args...) {} ~constexpr_optional_base() =default; }; template <typename T> using OptionalBase = std::conditional_t< std::is_trivially_destructible<T>::value, constexpr_optional_base<T>, optional_base<T> >; template <typename T> class optional : OptionalBase<T> { public: constexpr optional(const T& value) : OptionalBase<T>(value) {} constexpr const T& operator*() const { return OptionalBase<T>::storage_.value_; } }; int main() { // constexpr constexpr optional<int> a(42); static_assert(*a == 42); // 非constexpr /*constexpr*/ optional<std::string> b("nyan"); assert(*b == "nyan"); }