12. 粒子データ入出力
赤木 ：ここからは、割合本格的に多体計算しようみたいな話にはいろうかと思って。
学生C：今までとどう違うんですか？
赤木 ：今までのプログラムは、N体とかでもあんまり粒子数大きいのじゃなく
て、論文というか研究に使える感じじゃないんだけど、それを、まあ、研究レベルにね。
学生C：なんか大変そうですが、、、
赤木 ：まあ、その辺を、 Crystal で書くことでどれくらい楽にできるかをみ
ていこうみたいなね。
学生C：はあ。
赤木 ：まず、ある程度の粒子数を扱う、とすると、結果をファイルに書いた
り、入力をファイルから読んだり、としたくなるのね。今回その話。
学生C：ずーっと前の章の練習問題になんかなかったでしたっけ？FDPSのサン
プルプログラムのフォーマットとか。
赤木 ：粒子数とか時刻を最初に書いて、あと1行1粒子みたいなのね？
学生C：それです。
赤木 ：あれはあれでいけないわけじゃないんだけど、このデータは何か、ど
うやって作ったものか、というのが記録されてないじゃない？例えば、どうい
うコマンドで作ったか、とか、そもそもこの数字は何か、とか。
学生C：えー、データファイルってそういうものじゃなんですか？数字が並ぶ
わけですよね？
赤木 ：それだとなんだかわからなくなるから、わかるようにしておこうよ、
という話が色々なあるわけ。IT業界だと XML とか JSON とか YAML とか。
天文学だと FITS というのが観測データではよく使われているし、
大規模数値計算では HDF かしらね。
学生C：じゃあそのどれか使えばいいんでは?
赤木 ： それでもいいんだけど、そういうの使うのなんというか結構面倒でプ
ログラムも大変だから、なるべく簡単なのですませたいのね。まあ、他の
人が書いたプログラムに、だと変換プログラム書かないといけないけど、
それはまあしょうがないとして。
学生C：そうすると、どうするんですか？
赤木 ：最近だと楽に使えて、あとできたファイルも扱いやすいのは YAML か
な。作者氏とあと色々な人で昔論文書いてて、標準形式にしようと提案はして
るの。
学生C：でもそんなには使われてないと、、、、
赤木 ：まあそうなんだけど、この際だから。少なくとも論文があるから、こ
れに従ってやってね、といえるのはメリットね。YAML っぽいものはすで
使ってて、コマンドラインオプションを書く形式がそれなんだけど、あそこで
は割合簡単なのしか使ってなかったらから、論文に従って粒子データの例をま
ず出すわね。例えば。番号(名前)、位置、速度、質量をもつ1つの粒子を
以下のように表しましょう、というのが主張で、これに PSDF (Particle
Stream Data Format) っていう名前つけてるの。
--- !Particle
id: 0
r:
- 0.1
- 0.2
- 0.3
v:
- -1
- -2
- -3
m: 1.0
基本的に
key: value
みたいな、名前+コロン、スペース、値 である名前の変数なり物理量の値を指
定、位置とか速度みたいに複数の値の時には、上の例みたいに次の行を "-"
で始めて、それで始まる行が続く限り、配列に値がはいるみたいな。最初の
「---] は、新しく何かの始まりで、そのあとの 「!Particle」は
「タグ」っていう仕掛けになるのね。この辺、YAML のちゃんとした話は
少し複雑なんだけど、とりあえず1つの粒子のデータが
--- !Particle
で始まる、なので、次に 「---」 がでてくるかファイルの終わりがくるかし
たらその粒子は終わりで別にの粒子なりなんか違うものになるわけ。
PSDF の論文では、 粒子の名前、物理量として以下を使おう、としてるわね。
id 名前(整数でもテキストでもよい)
m 質量
t 粒子の時刻
t_max この粒子が有効な最大時刻
r 位置 3次元ベクトル
v 速度 3次元ベクトル
pot 重力ポテンシャル
acc 加速度 3次元ベクトル
jerk 加速度の一階時間導関数 3次元ベクトル
snap 2階 3次元ベクトル
crackle 3階 3次元ベクトル
pop 4階 3次元ベクトル
これに従っておけば、とりあえずファイルの意味は明らか、ということね。
本当は、!Particle のところで、このデータ構造が定義されている URI、ウェ
ブ上の場所ね、を書くとかあるんだけど、論文ではそこまでやってないからと
りあえずこれで。
学生C：なんか一杯ありますね。全部ないといけないんですか？
赤木 ：良い質問ね。もちろん、使うものだけあればいい、というのと、プロ
グラムの作り方として、例えばこの粒子データファイル、みんなスナップショッ
トというからここでもそう呼ぶけど、スナップショットをなんか加工するプロ
グラムがあったとして、それは、自分が扱うデータのことはしか知らないとし
て、他のデータは読んだものをそのまま書き出して欲しいわけ。で、この時に、
何なのか知らないデータがあっても、それを消したり変なものにしないような
仕掛けがほしいのね。
あと、データファイルに、こういう粒子のデータだけでなくて実行したコマン
ドとかオプションとかの記録があったら付け加える、というのも。
学生C：注文が多い気がしますが、、、なんか複雑なので、少し整理させて下
さい。まず、ここで作るのは、あくまでも多体シミュレーションのスナップ
ショットファイルである、ということでいいですね？
赤木 ：いいわよ。でも、それに、粒子データじゃない、実行したコマンドの
名前とかオプションとか、あとそれ以外の色々なプログラムのログ出力とかも
書けるようにしたいわけ。
学生C：それはだから、 Particle じゃない何かが必要だってことですね。そ
うしたら、なんかこんなふうなファイルですかね？
---- !CommandLog
command: mkplummer -n 10 ...
log:
(if any)
---- !Particle
id: 0
....
---- !Particle
id: 1
....
---- !Particle
id: 9
....
赤木 ：そうね。
学生C：なんか、クラスから YAML に変換は色々書くとやってくれるみたいな
ので、とりあえず上の CommandLog と Particle に対応するクラス定義と、
それでなんか書くプログラム作ってみました。こっちがクラス定義
1:#
2:# nacsio0.cr
3:#
4:# Basic YAML-based IO library for
5:# nacs (new art of computational science)
6:# Copyright 2020- Jun Makino
7:
8:require "yaml"
9:require "./vector3.cr"
10:
11:struct Vector3
12: def initialize(ctx : YAML::ParseContext, node : YAML::Nodes::Node)
13: a=Array(Float64).new(ctx, node)
14: @x=a[0]; @y=a[1]; @z=a[2]
15: end
16:end
17:
18:module Nacsio
19: class CommandLog
20: YAML.mapping(
21: command: {type: String, default: "",},
22: log: {type: String, default: "",},
23: )
24: def self.new
25: CommandLog.from_yaml("")
26: end
27: def self.new(logstring : String)
28: c=CommandLog.new
29: c.command = ([PROGRAM_NAME]+ ARGV).join(" ")
30: c.log = logstring
31: c
32: end
33: def to_nacs
34: self.to_yaml.gsub(/---/, "--- !CommandLog")
35: end
36: def add_command
37: @command += "\n"+([PROGRAM_NAME]+ ARGV).join(" ")
38: self
39: end
40:
41: end
42:
43: class Particle
44: def to_nacs
45: self.to_yaml.gsub(/---/, "--- !Particle")
46: end
47: def self.new
48: Particle.from_yaml("")
49: end
50: YAML.mapping(
51: id: {type: Int64, default: 0i64,},
52: time: {type: Float64, key: "t", default: 0.0,},
53: mass: {type: Float64, key: "m",default: 0.0,},
54: pos: {type: Vector3, key: "r",default: [0.0,0.0,0.0].to_v,},
55: vel: {type: Vector3, key: "v",default: [0.0,0.0,0.0].to_v,},
56: )
57: end
58:end
こっちがテストプログラム
1:# nacsiotest0.cr
2:#
3:# Test code for Basic YAML-based IO library for
4:# nacs (new art of computational science)
5:
6:require "./nacsio0.cr"
7:include Nacsio
8:
9:print CommandLog.new("Sample log message").to_nacs
10:(0..2).each{|id|
11: obj=Particle.new
12: obj.id =id.to_i64
13: obj.pos[0]=id*0.1
14: print obj.to_nacs
15:}
で、実行結果が
gravity> crystal nacsiotest0.cr
[2mShowing last frame. Use --error-trace for full trace.[0m
In [4mnacsio0.cr:20:10[0m
[2m 20 | [0m[1mYAML.mapping([0m
[32;1m^------[0m
[33;1mError: undefined method 'mapping' for YAML:Module[0m
です。 nacs という名前は、作者が元々使っていた acs に new の意味で n つけました。
赤木 ：なるほど。コードの中身説明してみて。
学生C：まず、CommandLog クラスですね。その前に Vector3 をなんかしてま
すがこの話はあとで。
YAML.mapping(
command: {type: String, default: "",},
log: {type: String, default: "",},
)
の、YAML.mapping がクラスとYAML 書式の関係つける仕掛けみたいです。最初
の command とか log がメンバー変数の名前、{} の中は type が型を、
default が(あれば) 指定がなかった時のデフォルト値、 key が(あれば)
YAML の中でのキー(なければメンバー変数の名前)、あと nilable (値を nil
にできるかどうか)と、もっと色々なことの指定ができるみたいですがまあ普
通これくらいで。
ちょっと面白かったのは、これで型とデフォルトの値を書いておくと、
initialize を別に定義しなくても from_yaml というクラス関数が new の代
わりになることで。でも new はないので、
def self.new
CommandLog.from_yaml("")
end
でデータがない文字列を解釈して、というのでデフォルトの値がはいった
インスタンスをつくれます。もっとちゃんとしたやり方もありそうですが、動
いたからこれでいいかなと。 new とか書換えられるんですね。
あと property も勝手に入るみたいで、
property :log
とか書かなくてもクラスの外から読み書きできました。あと、 new で文字列
渡すと、 command のほうに実行したコマンド、 log のほうにその文字列を
入れるようにしました。
赤木 ：実行結果の command のところが謎文字列だけど？
学生C：これ多分、 crystal foo.cr とすると実際にはコンパイルされるわけ
で、コンパイルされたものがそれだということなんではないかと、、、
赤木 ：そうね。多分。
学生C：あと、 to_yaml でYAML の文字列にしてくれるんですが、!CommandLog
のとこがつかないので gsub でつける関数が to_nacs です。この辺なんか
もっと YAML 的に正しい作法とかあるんじゃないかと思いますが、、、、
赤木 ：まあとりあえずこれでいいわよ。Particle は？
学生C：YAML.mapping の中身と、あと、つけるタグが CommandLog じゃなくて
Particle だというだけですね。あ、ややこしいのが、 pos, vel で、
これ Vector3 にしてるじゃないですか？最初は Array(Float64) で書いてた
んですが、それだとあとで不便な気がしたのでこっちにしました。
で、コンパイルしたら、なんか
> 197 |
> 198 |
> 199 | Vector3.new(ctx, value_node)
^--
Error: no overload matches 'Vector3.new' with types YAML::ParseContext, YAML::Nodes::Node+
Overloads are:
- Vector3.new(x : Float64 = 0, y : Float64 = 0, z : Float64 = 0)
みたいな、謎なエラーなんですが、要するに Vector3.new(ctx :
YAML::ParseContext, node : YAML::Nodes::Node) という関数がないといって
るっぽくて、これが Array では通ってたのは Array にはそういうものがある
のかな？とレファレンスみたらあったので、それを Array で実行した結果を
Vector3 にする initialize 関数を Vector3 に追加したら上手くいきました。
struct Vector3
def initialize(ctx : YAML::ParseContext, node : YAML::Nodes::Node)
a=Array(Float64).new(ctx, node)
@x=a[0]; @y=a[1]; @z=a[2]
end
end
のとこです。
赤木 ： ctx, node って何？
学生C： そこ調べてないです。動くからいいんじゃないかと。
赤木 ： そう、、、ねえ、、、
Crystal 公式のYAMLドキュメント
も正直なところよくわからないしね、、、
学生C：ちょっとこの辺まだドキュメント書けてない感じしますね。
赤木 ：そうね。これ、新しく作って書くのはこうでいいとして、読んでなん
かして書く、というのをやってみてくれない？つまり、例えば、このファイル
読んで、 id を1増やす、みたいな、粒子に変更加えるもの。
この時に、自分が知らない粒子クラスの YAML データでも、 id のデータがあ
ればそれをいじる、それ以外はさわらない、として欲しいの。
学生C： それは、この YAML 用の Particle クラスは使わない、ということに
なります？
赤木 ：そうね。まずはそれで。
学生C：CommandLog はどうしましょうか？
赤木 ：これは、使って、あったコマンドログの最後に自分を付け加えるでどう？
学生C：考えてみます。
(数日後)
一応できたことにしたいですが、、、、
赤木 ：あら、早いわね。動かすとどんな感じ？
学生C：
gravity> crystal nacsiotest0.cr > test.yml
[2mShowing last frame. Use --error-trace for full trace.[0m
In [4mnacsio0.cr:20:10[0m
[2m 20 | [0m[1mYAML.mapping([0m
[32;1m^------[0m
[33;1mError: undefined method 'mapping' for YAML:Module[0m
gravity> crystal nacsreadwrite3.cr < test.yml
[2mShowing last frame. Use --error-trace for full trace.[0m
In [4mnacsreadwrite3.cr:10:8[0m
[2m 10 | [0m[1mYAML.mapping([0m
[32;1m^------[0m
[33;1mError: undefined method 'mapping' for YAML:Module[0m
と、こんな感じで、 id が 0-2 だったのを1増やして 1-3 になってます。
赤木 ：それっぽいわね。プログラム本体は？
学生C： とりあえずこれです
1:# nacsreadwrite3.cr
2:#
3:# Test code for Basic YAML-based read/write for
4:# nacs (new art of computational science)
5:
6:require "./nacsio.cr"
7:include Nacsio
8:
9:class IDParticle
10: YAML.mapping(
11: id: {type: Int64, default: 0i64,},
12: )
13:end
14:
15:update_commandlog
16:
17:while (sp= CP(IDParticle).read_particle).y != nil
18: p = sp.p
19: p.id += 1
20: sp.p = p
21: sp.print_particle
22:end
ここでは、6行目で nacsio.cr を読み込んで、9-13行で IDParticle という id だけもったク
ラスを定義します。クラスは YAML.mapping だけです。
で、 nacsio.cr のほうでは
update_commandlog
でログ部分を読んで書く、
CP(T).read_particle
でクラス T の粒子を読み込む、 print_particle で CP(T) 型をYAML で書く、
となるわけです。
赤木 ：CP(T) 型って何？
学生C： これは、T のとこがクラスで、そのクラスの変数をメンバーとしても
つクラスを作ります。ライブラリ本体は
1:#
2:# nacsio.cr
3:#
4:# Basic YAML-based IO library for
5:# nacs (new art of computational science)
6:# Copyright 2020- Jun Makino
7:
8:require "yaml"
9:require "./vector3"
10:
11:# struct Vector3
12:# def initialize(ctx : YAML::ParseContext, node : YAML::Nodes::Node)
13:# a=Array(Float64).new(ctx, node)
14:# @x=a[0]; @y=a[1]; @z=a[2]
15:# end
16:# end
17:
18:module Nacsio
19: extend self
20: class CommandLog
21: include YAML::Serializable
22: @[YAML::Field(key: "command")]
23: property command : String
24: @[YAML::Field(key: "log")]
25: property log : String
26:
27: # YAML.mapping(
28: # command: {type: String, default: "",},
29: # log: {type: String, default: "",},
30: # )
31: def self.new
32: CommandLog.from_yaml(%(
33: command:
34: log:
35:))
36: end
37: def self.new(logstring : String, progname = PROGRAM_NAME, options = ARGV)
38: c=CommandLog.new
39: c.command = ([progname]+ options).join(" ")
40: c.log = logstring
41: c
42: end
43: def to_nacs
44: self.to_yaml.gsub(/---/, "--- !CommandLog")
45: end
46: def add_command(progname = PROGRAM_NAME, options = ARGV)
47: @command += "\n"+([progname]+ options).join(" ")
48: self
49: end
50:
51: end
52:
53: class Particle
54: def to_nacs
55: self.to_yaml.gsub(/---/, "--- !Particle")
56: end
57: def self.new
58: Particle.from_yaml(%(
59: id: 0
60: t: 0.0
61: m: 0.0
62: r:
63: x: 0.0
64: y: 0.0
65: z: 0.0
66: v:
67: x: 0.0
68: y: 0.0
69: z: 0.0
70:))
71: end
72: include YAML::Serializable
73: @[YAML::Field(key: "id")]
74: property id : Int64
75: @[YAML::Field(key: "t")]
76: property time : Float64
77: @[YAML::Field(key: "m")]
78: property mass : Float64
79: @[YAML::Field(key: "r")]
80: property pos : Vector3
81: @[YAML::Field(key: "v")]
82: property vel : Vector3
83: # YAML.mapping(
84: # id: {type: Int64, default: 0i64,},
85: # time: {type: Float64, key: "t", default: 0.0,},
86: # mass: {type: Float64, key: "m",default: 0.0,},
87: # pos: {type: Vector3, key: "r",default: [0.0,0.0,0.0].to_v,},
88: # vel: {type: Vector3, key: "v",default: [0.0,0.0,0.0].to_v,},
89: # )
90: end
91:
92: def update_commandlog(progname = PROGRAM_NAME, options = ARGV,
93: of = STDOUT)
94: s=gets("---")
95: s=gets("---") if s == "---"
96: a=s.to_s.split("\n")
97: a.pop if a[a.size-1]=="---"
98: if a[0] != " !CommandLog"
99: raise("Input line #{a[0]} not the start of Commandlog")
100: else
101: a.shift
102: ss = (["---\n"] + a).join("\n")
103: of.print CommandLog.from_yaml(ss).add_command(progname,options).to_nacs
104: end
105: end
106: def read_commandlog(ifile= STDIN)
107: c=CommandLog.new
108: s=ifile.gets("---")
109: s=ifile.gets("---") if s == "---"
110: a=s.to_s.split("\n")
111: a.pop if a[a.size-1]=="---"
112: if a[0] != " !CommandLog"
113: raise("Input line #{a[0]} not the start of Commandlog")
114: else
115: a.shift
116: ss = (["---\n"] + a).join("\n")
117: c= CommandLog.from_yaml(ss)
118: end
119: c
120: end
121:
122: class CP(T)
123: property :p, y
124: def initialize(p : T, y : YAML::Any)
125: @p=p
126: @y=y
127: end
128: def self.read_particle(ifile = STDIN)
129: s=ifile.gets("---")
130: retval=CP(T).new(T.new,YAML::Any.new(nil))
131: if s != nil
132: a=s.to_s.split("\n")
133: a.pop if a[a.size-1]=="---"
134: if a.size > 0
135: if a[0] != " !Particle"
136: raise("Input line #{a[0]} not the start of Particle")
137: else
138: a.shift
139: ystr = (["--- \n"] + a).join("\n")
140: retval=CP(T).new(T.from_yaml(ystr), YAML.parse(ystr))
141: end
142: end
143:
144: end
145: retval
146: end
147:
148: def print_yamlpart(of = STDOUT)
149: of.print @y.to_yaml.gsub(/---/, "--- !Particle")
150: end
151:
152: def print_particle(of = STDOUT)
153: yy=@y.as_h.to_a
154: ycore = YAML.parse(@p.to_yaml).as_h.to_a
155: ycore.each{|core|
156: yy.reject!{|x| x[0]== core[0]}
157: }
158: yy = ycore + yy
159: newstring = YAML.build{|yaml|
160: yaml.mapping{
161: yy.each{|x|
162: yaml.scalar x[0].as_s
163: if x[1].raw.class == Int64
164: yaml.scalar x[1].as_i64
165: elsif x[1].raw.class == Float64
166: yaml.scalar x[1].as_f
167: elsif x[1].raw.class == String
168: yaml.scalar x[1].as_s
169: elsif x[1].raw.class == Array(YAML::Any)
170: xx = x[1].as_a
171: yaml.sequence { xx.each{|v| yaml.scalar v.as_f}}
172: end
173: }
174: }
175: }
176: of.print newstring.gsub(/---/, "--- !Particle")
177: end
178: end
179:
180: def repeat_on_snapshots(p = Particle.new, read_all = false)
181: sp= CP(typeof(p)).read_particle
182: no_time = (sp.y["t"] == nil )
183: while sp.y != nil
184: pp=[sp]
185: time = sp.y["t"].as_f unless read_all
186: while (sp= CP(typeof(p)).read_particle).y != nil &&
187: (read_all || no_time || sp.y["t"].as_f == time)
188: pp.push sp
189: end
190: yield pp
191: end
192: end
193:end
194:
で、CP(T) は
def initialize(p : T, y : YAML::Any)
@p=p
@y=y
end
が new で呼ぶ関数なので、 T型の p と、YAML::Any という謎型の y をメン
バ変数に持つわけです。で、 read_particle はこの型のクラスメソッドで
sp= CP(IDParticle).read_particle
で標準入力、これ変更できたほうがいい気もしますが、から1粒子読んで、
引数で与えた粒子クラスで解釈できるものを sp.p に、それを含めた全部を
Crystal の YAML ライブラリ のデフォルトの型である YAML::Any にいれる、
ということをします。あと、ファイルが終わってるとかの時には sp.y のほう
を nil にします。こっちは Any なので nil にしてもコンパイラが文句をい
かなかったので。
プログラム本体のほうでは、 sp.p が IDParticle型になってるので、 id を
変更したりできるわけです。色々してから、 print_particle をします。
これは、色々ややこしいことしてますが、要するに IDParticle 型の p にあるデー
タはそっちで、ないものは y のほうにあるのをそのまま書く、というもので
す。ただし、 YAML のデータ構造を触るのがあんまり上手くできてなくて、一
旦配列にして、また戻すとかしているので、ちょっと汎用性がなくて、今のところ
しか扱えないです。
赤木 ：まあ、そうわかってればいいかしらね。
学生C：
赤木 ：
学生C：
赤木 ：
学生C：
赤木 ：
学生C：
赤木 ：
学生C：
12.1. 課題
-
上の nacsreadwrite3.cr を参考に、位置座標を2倍にして出力するプログ
ラムを作って動作を確認せよ。
-
さらに、コマンドラインオプションで指定した値で 位置、速度、質量をス
ケールするプログラムに拡張してみよ。
12.2. まとめ
-
YAML 型式で粒子データを読み書きするライブラリを作成した。
-
このライブラリは、プログラムの粒子データクラスにはいってないものは
そのまま残して、おいて、粒子データクラスと合わせてあとで出力することもできる
12.3. 参考資料
module YAML
PSDF論文