Previous ToC Next

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. 課題

  1. 上の nacsreadwrite3.cr を参考に、位置座標を2倍にして出力するプログ ラムを作って動作を確認せよ。
  2. さらに、コマンドラインオプションで指定した値で 位置、速度、質量をス ケールするプログラムに拡張してみよ。

12.2. まとめ

  1. YAML 型式で粒子データを読み書きするライブラリを作成した。
  2. このライブラリは、プログラムの粒子データクラスにはいってないものは そのまま残して、おいて、粒子データクラスと合わせてあとで出力することもできる

12.3. 参考資料

module YAML

PSDF論文
Previous ToC Next