Contents My Feed Description 2020-03-15T22:23:00+09:00 Tosh Morry https://tm23forest.com/contents TOEFL iBTとTOEICの換算式が違いすぎる話 https://tm23forest.com/contents/toefl-ibt-toeic-mis-exchange 2020-03-15T22:23:00+09:00 2020-03-15T22:23:00+09:00 このページは単なる文句というか意見をさらしておくだけのものにつき予めご了承ください。タグもOpinionになっているしあくまで意見だということで。

はじめに

少し前の話。これは大学・大学院の入学試験を受けたことのある人ならば良くわかる話だと思うが、このあたりになってくると英語はしばしばTOEICあるいはTOEFLで置き換えられる。このときにあり得るパターンとしては次の3つだろう。

  • TOEICのスコアを英語の代わりに提出する
  • TOEFL iBTのスコアを英語の代わりに提出する
  • TOEICまたはTOEFL iBTのスコアを英語の代わりに提出する

TOEFLにも何種類かあると思うが現状もっともあると感じているのはiBT、というかこのページではiBTの話をしたいのでそのようにしていうる。1個目と2個目は問題ない、問題は3番目。

3番目の時には、何かオフィシャルになってるっぽい方法で両者の得点を換算しますっていう記述が目につく。受験者によって提出してくるスコアの種類が違うのだから何らかの方法で換算して比較することはごく自然なことだろう。もう言いたいことわかると思うが、問題なのは換算することではなくてその換算の方法。具体的には両者のスコアレンジの対応表。

TOEFL iBTスコア過小評価されすぎ問題

いちいち具体的なスコアレンジの対応表を上げることはしないが、複数大学受験するなどしてTOEICもTOEFL iBTも受験したことのある人なら痛感していることだろう。明らかにTOEFL iBTのスコアが過小評価され過ぎている。これは少し前といっても5年以上前の話なので今はマシになっていると思いたいが少なくとも5年前は完全にそうだった。

入学試験の要項から普通のウェブサイトまでいたるところで出回っている換算表を見るたびに本当にこの表をさらしている人間は両方のテストを受けたことがあるんだろうかと毎回エアプレイを疑っていた。というかエアプレイが確定しているので疑ってはいなかった。ほとんどの表は何かしらの情報源をもとにその表を作っているだろけど一番大元というのが存在するわけで、それがオフィシャルであろうと信じられなかった。その理由がいくつかある。

iBTはマジでスコアが取れない

まずはこれ。iBTにはリーディング、リスニング、スピーキング、ライティングの4項目がある。長いので以後R、L、S、Wとよぶ。それぞれの配点は30点で満点は120点となる。仮にTOEIC800点くらいの人間がいたとする。TOEIC800点といえば世間的にはまぁそこそこ英語出来るとみなされておかしくない。それぐらいならばiBTで60点、つまり半分くらいは対策しなくても取れると思うだろう。

絶対に取れない。対策しても厳しいかもしれない。

絶対に取れない理由を簡潔に説明すると、まずTOEICはインプットの能力しか問われないので対策無しでiBTに臨んだ場合アプトプットの能力であるS、Wの点数はゼロとなる。ゼロまではいかんと思われた方がいるかもしれないがそんなことは無い。ゼロだ。対策をして臨んだ場合はどうだろうか?TOEIC800点だとS、Wそれぞれ1桁点数がいいとこだろう。10点超えたら素晴らしい。よく800点の能力で対策したな、と感じる。強い言葉を使っているがこれは受けた人に初めて受けたときの感覚を思い出してわかってほしい。

スコアが取れない理由はSとWだけではない。RとLの点数もTOEICが出来るからといってその地力で臨んでいい点数が望めるとは限らない。まず単純にインプットの量が多い。TOEIC10分余して終わるからインプット早いと思っているかもしれないが、そのレベルでもiBTのR全文余裕をもって読むことは簡単ではないだろう。実際には回答しながら読み進めることにはなるのでそう単純なものではないが。あとTOEFLの問題はTOEICのよくわからんインヴォイスとかチラシのように浅いロジックではないから複数回読みも頻繁に発生する。次に問題の回答形式が4択以外もあるし選択肢そのものの文章が長いこともある。

とはいったもののRに関してはまだ何とかなる。Lに関しては1回に聞き取る量、つまり時間がハンパなく長いのと当然のように1回しか読まれないので本当に意味を理解して聞き進めないと得点には結びつかない。断片的に聞こえたどうこうは全く役に立たない。暗記要素をなくすためなのかメモ用のペーパーと鉛筆が渡されるが非ネイティブにとっては聞きながらメモを取るのが難しい、場合によってはメモ取らない全記憶型フルミーニング理解戦法のほうがまだ点数は取れる。Lに関してはSとWの設問部分にも使われるから聞き取れないとそもそも何もできない問題も存在する。

TOEICはスコアを上げられる

見出しがアホみたいな並びだがこれは本当のことである。TOEIC800点を先ほどから例として挙げているが、なぜこの辺りを取り上げているかというとこの辺で換算表がバグを起こし始める印象があるから。一口にTOEIC800点といっても実は2種類存在していて、

  • テスト特化対策によって地力より高い点数を出して800点
  • そこそこのインプットの地力のある人が何となく受けて800点

があり得る。まぁどんなテストでも対策は重要だ。この2種類はいつでも存在する。TOEICの場合はこの幅がすんごい広いのである。インプットの速さ、という1点のみの地力をつけて各パートごとの対策をやって余裕を持って時間内に回答を埋めることが出来れば必ず800点は超える。語彙のレンジもiBTほど広くない。リスニングも暗記要素が無いから本当にTOEFLに比べて楽に点数を上げられる。インプットさえあれば何も怖くない世界だ。上げて800の人間がiBT60点不可能なことは誰もが納得するだろう。

後半の地力に関してはインプットの、と書いたのがミソ。この800という点数はネイティブレベルには取れない点数だ。アウトプット込みのネイティブレベルならばTOEICで900を下回ることはほぼあり得ない。従ってここで言う地力というのはまぁ受験英語の延長で頑張って英語勉強したんだなっていうレベルのもので、その辺の試験英語の長文は全文読んでも間に合うし英語のウェブページやリファレンス、英語論文などはスラスラ読める、ネイティブがしゃべっていることもなんとなく聞き取って意味が取れるよっていうくらい。こういうレベルにいる人はテスト特化対策を施して何となくじゃなくちゃんとやれば900点は超える。なんなら各2ミスで975とかもあり得るだろう。上げて800が出来るのならば上げて900後半もあり得るわけだ。とにかくTOEICは点数を伸ばしやすい。

まとめ

TOEFL iBTは頑張っても点数を伸ばしにくくてTOEICは頑張ったら簡単に点数が伸びる。測りたいものが頑張ったかどうかではなく英語の能力なのであるならば頑張ってすぐ伸びてしまうTOEICの謎の課題評価をやめて、英語の地力が確実に反映されるTOEFL iBTの上の方のレンジの評価をもっと引き上げて換算しましょうということがいいたい。そもそもSとWが入っているから比べているものが違うにしてもTOEFLの方が英語をより高いレベルで使いこなさないと同じRとLでも点数は出ないんだ。強い言葉を使っているけどもこれは似たようなシチュエーションで受けたことがある人にわかってほしい。

]]>
Q-tree: 領域四分木の詳解とPythonによる実装 https://tm23forest.com/contents/qtree-algorithm-python-impl-etc 2020-02-09T23:13:00+09:00 2020-02-09T23:13:00+09:00 特にドロー系の画像編集ソフトなどでは、指定された2次元領域に関係しているオブジェクトを素早く特定して処理を行うために空間インデックスといわれる類のデータ構造をしばしば用いている。その内の一つである領域四分木というデータ構造とアルゴリズムを見ていき、Pythonで素朴に実装してみる。具体的なチンパンコード付き。

Q-treeアルゴリズムのイメージ
Q-treeアルゴリズムのイメージ

データ構造とアルゴリズム

四分木とは、次のような問題を解くことのできるアルゴリズムの1つであり、場合にもよるがQ-treeだとかQuad-treeとも呼ばれる。

入力として次を与える。

  • N個のx,yの組
  • xmin, xmax, ymin, ymax 全てのx,yの組はxmin<x<xmaxかつymin<y<ymaxを満たしているものとする。
  • M個のx0, x1, y0, y1の組 が与えられたとき、それぞれのx0, x1, y0, y1に対して、N個のx,yの組の内x0<x<x1かつy0<y<y1をみたすものを全て出力せよ。N, Mは整数、それ以外は実数とする。

脳死アルゴリズムによってこの問題を解くと計算ステップのオーダはO(MN)である。二分木といった木構造を用いると計算ステップのオーダはO((N+M)log N)となる。

二分木とは

四分木にいくまえにその簡単な場合である二分木についてみていく。二分木とは、大小関係を持ったデータをある意味でソートしながら記憶し、ある点に最も近いデータ、あるいはある区間内にあるデータの列挙といったクエリに対して効率良く答えることが出来る。ここでいう大小関係とは、数直線上の点であればそのまま数値の大小関係を用いればよいし、文字列であれば辞書に出てくる順番の前後関係を大小関係と見なせばよい。

具体的にどうするかというと、次の図のようなデータ点の1つ1つが頂点に対応する木構造を構築する。

二分木の適当な例
二分木の一例

一番上から順番に自分との大小関係を比較していき、自分のいる側に移動していって最終的にそれ以上エッジが無くなったらその先に自分を追加すればよい。全ての頂点は高々2つのエッジを持つという特徴がある。クエリに対しては直接的に二分探索を使って応答することが出来る。

四分木とは

二分木は大小関係を持ったデータに対して適用できるとあったが平面上の点というデータに対しては使えるだろうか?とにかく大小関係を与えればよい、という点では例えば「原点との距離」を用いれば可能だろうが先に上げた例のようなクエリに対してはふさわしくないだろう。あるいはx座標のみ参照して大小関係を与えるという方法もあるがやや微妙。そこで四分木というデータ構造がある。やることは簡単で、xの大小関係とyの大小関係のそれぞれ2つによってエッジを生成していく。そのほかのやり方は二分木の場合と同じ。

領域四分木

単純な四分木でも先の問題に答えることは出来そうだが、まだやや問題がある。二分木とその素朴な拡張による四分木のどちらに対しても言えることとして、生成される木構造がデータの入力順に依存しているということがある。二分木の場合の極端な例として、値の小さいものから順番にデータを木構造に追加していくと右側エッジの先にしかノードが生えず、結果として単なるリストと等価なものが出来上がってしまう。これではlogNにはならないので二分木はデータ追加の過程でルートノードを付けかえる木構造のバランシング処理を行うことが普通だ。四分木の場合にも同じ問題は起こり、コーディングの面とバランス判定基準の調整の面の両方でバランシング処理はより面倒になる。

四分木の場合には、木構造をデータ入力の順番に依存させない領域四分木という方法がある(領域二分木というものが無いわけではないが聞いたことが無い、なぜか?)。領域四分木では全データのバウンディングボックスを再帰的に4分割していくことによって木構造を構成していく。分割のイメージはページ上部のものを参考にされたい。

さっきまでのデータ構造とは異なり、それぞれの頂点は領域を表しデータは木構造の「葉」としてグラフの末端に格納される。もし例の問題のように各座標の上下限がわからなくても一度データを舐めればオーダO(N)でバウンディングボックスの算出は可能である。明らかな特徴として、領域四分木ではデータを入力する順番を変えても得られる木構造は同じになる。

PythonによるQ-treeのコーディング

最初に言っておかなければならないのは、これは実験であり車輪の再開発である。どんな言語でもあると思うがPythonにもpyqtreeというQ-treeの実装を与えるパッケージはあるのですぐに何かをやりたい場合はそれを使うのが賢明だろう。

Pythonプログラム

import numpy as np

QmaxD = 50 # maximum depth
QmaxN = 5  # maximum items per boxel
Qtree = {}
Qbbox = [-1, 1, -1, 1] #xmin, xmax, ymin, ymax

def Qreset(bbox, QmaxN=5, QmaxD=50):
  Qtree = {}
  Qbbox = list(bbox)

def Qinsert(point):
  qtree, qbbox = Qtree, list(Qbbox)

  for _ in range(QmaxD):
    key, x0, y0 = 0, (qbbox[0]+qbbox[1])/2, (qbbox[2]+qbbox[3])/2
    if point[0] < x0:
      qbbox[1] = x0
      key += 1
    else:
      qbbox[0] = x0
    if point[1] < y0:
      qbbox[3] = y0
      key += 2
    else:
      qbbox[2] = y0

    if not key in qtree:
      qtree[key] = []
      break
    if type(qtree[key]) is list:
      break

    qtree = qtree[key]

  if len(qtree[key]) < QmaxN:
    qtree[key].append(point)
  else:
    subtree = {}
    x0, y0 = (qbbox[0]+qbbox[1])/2, (qbbox[2]+qbbox[3])/2
    for point in qtree[key]+[point]:
      subkey = 0
      if point[0] < x0: subkey += 1
      if point[1] < y0: subkey += 2
      if subkey in subtree:
        subtree[subkey].append(point)
      else:
        subtree[subkey] = [point]
    qtree[key] = subtree

def Qquery(bbox):
  result = []
  Qstack = [(Qtree, list(Qbbox), 0)]
  pc = 0
  while len(Qstack) > 0:
    pc += 1
    qtree, qbbox, key = Qstack[-1]
    if type(qtree) is list:
      for point in qtree:
        if bbox[0] <= point[0] and point[0] <= bbox[1] and bbox[2] <= point[1] and point[1] <= bbox[3]:
          result.append(point)
      Qstack.pop()
      continue

    if key > 3:
      Qstack.pop()
      continue

    x0, y0, qbbox = (qbbox[0]+qbbox[1])/2, (qbbox[2]+qbbox[3])/2, list(qbbox)
    if key % 2 == 1:
      qbbox[1] = x0
    else:
      qbbox[0] = x0
    if key // 2 == 1:
      qbbox[3] = y0
    else:
      qbbox[2] = y0

    Qstack[-1] = (Qstack[-1][0], Qstack[-1][1], key+1)
    if key in qtree:
      if qbbox[0] <= bbox[1] and bbox[0] <= qbbox[1] and qbbox[2] <= bbox[3] and bbox[2] <= qbbox[3]:
        Qstack.append((qtree[key], qbbox, 0))
  return result

はい。結構長くなってしまっている。読むのは何かしらのトレーニングにはなるだろうが苦痛なのでやめたほうがいいだろう。これをmyqtree.pyとしていろいろビジュアライズしていこう。

プログラムの解説

解説といってもコードを追うようなことはしないが、Pythonのdictを使いまくっている。パフォーマンスを追求することは目的ではないのが体感クエリに関しては相当早く応答する(dictの特性からどんなにPythonの中で変なことしてても速いのは当たり前だが)。ちなみに予想どおりdictをアロケートしまくるインサートは遅い。

領域四分木では一つの頂点(領域)にぶら下がる葉(データ)の数に上限を設け、その上限を超えたときにその領域を分割する。このプログラムではその上限をQmaxN=5としている。ここで1つだけ問題があり、まったく同一の点がQmaxN個以上入力されると、実装にもよるが「ボクセル内のデータ数がQmaxNを超えてしまう」か「無限に領域を分割するループに入る」のいずれかが起こってしまう。後者の対策としてQmaxDとして深さの上限を設けているが50だと全然ヒットしないので放置してある。ついでにこの実装だとおそらく前者が起こるので気にしなくていいはず。

それにしてもこの手の車輪をC/C++で再開発していたころに比べて大分やりやすくなったもんだな。

疑似データによる実行結果

バウンディングボックスは[-1,+1]^2の範囲としてその中に一様乱数を用いて2000個のランダム点を散らす。3つのクエリに対する応答を可視化した結果次のようになる。これに用いたコードは次のようになる。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as ptch
from myqtree import Qreset, Qinsert, Qquery

points = []
for i in range(2000):
  point = np.random.rand(2)*2-1
  points.append(point)
  Qinsert(point)

bboxes = [
  [0.25, 0.4, 0.25, 0.4],
  [-0.5, -0.3, 0.2, 0.9],
  [-0.9, -0.4, -0.5, -0.4]
]

fig, ax = plt.subplots(figsize=(5,5))
ax.set_xlim([-1,1])
ax.set_ylim([-1,1])
ax.set_aspect('equal')

P = np.array(points).T
ax.scatter(P[0], P[1], s=3)

for bbox in bboxes:
  result = Qquery(bbox)
  if len(result) > 0:
    R = np.array(result).T
    ax.scatter(R[0], R[1], s=6)
    ax.add_patch(ptch.Rectangle((bbox[0], bbox[2]), bbox[1]-bbox[0], bbox[3]-bbox[2], fill=False, lw=1))
plt.show()
plt.close(fig)

これに対する結果は以下。

Q-treeアルゴリズムの実行結果
Q-treeアルゴリズムの実行結果

プロットした点の大きさの関係ではみ出ているように見えるやつもあるがクエリで与えられたボックスの中にあるものをもれなく列挙出来ている。

計算ステップの計測実験

ここまでであれば高々2000点だし脳死アルゴリズムとの違いは大してわからない。そこで、クエリ数M=3(先の例で使ったもの)を固定してデータ数Nを大きくしていったときに、計算量の目安としてループカウンタがどこまで大きくなるかを計測した。ついでに予想される計算量の曲線もプロットしていく。素朴に考えるとインサートにNlogN、クエリには3logNのオーダで曲線は伸びていくはずである。

Q-treeの計算ステップ数のグラフ
ステップ数のプロット

これが結果。insertのループカウンタについてはおおむね予想通りの回数が出ている。ところがqueryについては大幅に大きなステップ数(の増加オーダ)がでています。なぜだろうと考えると、このqueryのループカウンタは四分木のエッジを辿るところにも1ループ使っているので本来ならば2^Nのオーダが出るはずの深さ優先のイテレーションです。深さ優先の利点は記憶域のオーダが高々木構造の深さ倍で済むところで、さらに今回の場合はエッジの先にある領域が完全にバウンディングボックスからはみ出ている場合にはそのエッジを切り落とすことができます(枝刈り)。この切り落とした回数をcutとしてプロットしてみると枝刈りによって相当な恩恵を受けていることがわかります。実際、2^NのオーダはプロットしていませんがNlogNよりも急峻なはずですからqueryの計算ステップは2^Nと比べて大幅に削減されていることがわかります。

まとめ

以上でQ-treeのデータ構造とアルゴリズム、さらにPythonのチンパンコードによる軽い実験をおわる。最後のところはM=O(N)としてみるとN^2と(N+M)logN=NlogNの比較になるので領域四分木が速いだろうことはわかるだろう。ついでにチンパンコードから得た所感はクエリに対してはべらぼうに早いのでM=大だともっと強くなるだろうね。対数目盛でわかりにくいがN=5,000,000までPython3.8.1+Windows10の環境で動いているのでPythonはメモリ1G喰っても文句言わずに動いて偉いですな。

]]>
Pythonで地理院地図を画像としてファイルに保存 https://tm23forest.com/contents/python-gsimaps2file-image 2020-02-07T23:29:00+09:00 2020-02-07T23:29:00+09:00 地理院地図はウェブページから普通に閲覧は出来るが、ウェブマップサービスとしてタイル座標からそれぞれの画像にアクセスすることもできる。これを使ってPythonから画像を取得して1枚の画像にするプログラムを置いておく。

Pythonコード

ワールドファイルは生成しないのでただの画像としてしか使えない。QGISからxyzタイル使って保存したほうが画面見ながら範囲指定して保存できるしワールドファイル作れるし再投影できるしこっちはオワコン。

#!python3
# -*- coding:utf-8 -*-

import sys
import math
import csv
import os.path
import requests
from io import BytesIO
from PIL import Image

def tile2latlon(x, y, z):
    lon = (x / 2.0**z) * 360 - 180
    mapy = (y / 2.0**z) * 2 * pi - pi
    lat = 2 * math.atan(e ** (- mapy)) * 180 / math.pi - 90
    return [lon,lat]

def latlon2tile(lon, lat, z):
    x = int((lon / 180 + 1) * 2**z / 2)
    y = int(((-math.log(math.tan((45 + lat / 2) * math.pi / 180)) + math.pi) * 2**z / (2 * math.pi)))
    return [y,x]

def execute(output, type, zoom, latitude, longitude, hsize, vsize):
  session = requests.session()
  headers = {'User-Agent' : 'Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:55.0) Gecko/20100101 Firefox/55.0'}
  basename, ext = os.path.splitext(output)
  y0, x0 = latlon2tile(longitude, latitude, zoom)
  width = 2*hsize + 1
  height = 2*vsize + 1
  image = Image.new("RGB", (256*width,256*height), (255,255,255))
  for i in range(width):
    for j in range(height):
      x = x0 - hsize + i
      y = y0 - vsize + j
      print('downloading tile at '+str([x,y]))
      response = session.get('http://cyberjapandata.gsi.go.jp/xyz/%s/%d/%d/%d%s' % (type, zoom, x, y, ext), headers=headers)
      subimage = Image.open(BytesIO(response.content))
      image.paste(subimage, (256*i,256*j))
  image.save(output)

def main(argv):
  if len(argv) < 8:
    print('usage: python %s [output] [type] [z] [latitude] [longitude] [hsize] [vsize]' % argv[0])
  else:
    try:
      output = argv[1]
      type = argv[2]
      zoom = int(argv[3])
      latitude = float(argv[4])
      longitude = float(argv[5])
      hsize = int(argv[6])
      vsize = int(argv[7])
    except ValueError:
      print('argument error')
      quit()
    execute(output, type, zoom, latitude, longitude, hsize, vsize)

if __name__ == '__main__':
    main(sys.argv)
]]>
Pythonで国土地理院のDEM5A標高タイルからGeoTiffを作成 https://tm23forest.com/contents/python-gdal-cyberjapandata-dem5tile 2020-02-07T23:07:00+09:00 2020-02-07T23:07:00+09:00 国土地理院が提供する標高タイルをダウンロードしてGeoTiff形式にする簡単なPythonプログラムをメモ用に置いておく。

標高タイルとは

https://maps.gsi.go.jp/development/demtile.html

ここで解説されている、国土地理院が提供するウェブメルカトルタイルのそれぞれの点の標高値が取得できるヤツ。基盤地図情報などとまとめて提供される数値標高モデルとの関係(点の密度はどっちが高いとかね、)はまた別ページで考察したいと思うが、アカウント登録も不要だしはっきり言ってコッチの方が使いやすい。

Pythonコード

基本的にはタイルの標高値でピクセルを埋めるが優先度はDEM5A>DEM5Bとなっている。DEM5Aはレーザ測量で街区や緩斜エリアに多く、DEM5Bは写真測量で山間部に多い。解説はしない。

#!python3
# -*- coding: utf-8 -*-

import os
import sys
import codecs
import csv
import requests
import math
import numpy as np
import osgeo.gdal,osgeo.osr

def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2**zoom
  xtile = int((lon_deg + 180.0)/360.0*n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1/math.cos(lat_rad)))/math.pi)/2.0*n)
  return xtile,ytile

def num2deg(xtile, ytile, zoom):
  n = 2**zoom
  lon_deg = xtile/n*360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi*(1-2*ytile/n)))
  lat_deg = math.degrees(lat_rad)
  return lat_deg,lon_deg

def main(argv):
  try:
    output_filename = argv[1]
    lat = float(argv[2])
    lon = float(argv[3])
    hsize = int(argv[4])
    vsize = int(argv[5])
  except IndexError:
    print('usage: python %s [output_filename] [latitude] [longitude] [hsize] [vsize]' % (argv[0]))
    quit()

  nodata = -1

  z = 15
  x,y = deg2num(lat, lon, z)

  raster = np.zeros((vsize*256, hsize*256), dtype=np.float64)

  headers = {'User-Agent' : 'Mozilla/5.0 (Windows NT 10.0; WOW64; rv:54.0) Gecko/20100101 Firefox/54.0'}

  for i in range(vsize):
    for j in range(hsize):
      tx = x - hsize//2 + j
      ty = y - vsize//2 + i

      url5a = 'http://cyberjapandata.gsi.go.jp/xyz/dem5a/%d/%d/%d.txt' % (z,tx,ty)
      url5b = 'http://cyberjapandata.gsi.go.jp/xyz/dem5b/%d/%d/%d.txt' % (z,tx,ty)
      print(i,j,url5a)
      print(i,j,url5b)
      response1 = requests.get(url5a, headers=headers)
      response2 = requests.get(url5b, headers=headers)

      dem5a = list(csv.reader(response1.text.splitlines()))
      dem5b = list(csv.reader(response2.text.splitlines()))

      if len(dem5a) == 256:
        for p in range(256):
          for q in range(256):
            if dem5a[p][q] == 'e' : raster[i*256 + p][j*256 + q] = nodata
            else: raster[i*256 + p][j*256 + q] = float(dem5a[p][q])
      else:
        print('error dem5a',i,j)
        for p in range(256):
          for q in range(256):
            raster[i*256 + p][j*256 + q] = nodata

      if len(dem5b) == 256:
        for p in range(256):
          for q in range(256):
            if raster[i*256 + p][j*256 + q] == nodata:
              if not dem5b[p][q] == 'e': raster[i*256 + p][j*256 + q] = float(dem5b[p][q])
      else: print('error dem5b',i,j)

  lat0,lon0 = num2deg(x-hsize//2,y-vsize//2,z)
  lat1,lon1 = num2deg(x+hsize//2+1,y+vsize//2+1,z)

  print(lat0,lon0)
  print(lat1,lon1)

  trans = [lon0, (lon1-lon0)/(hsize*256), 0, lat0, 0, (lat1-lat0)/(vsize*256)]

  srs = osgeo.osr.SpatialReference()
  srs.ImportFromEPSG(4326)

  driver = osgeo.gdal.GetDriverByName('GTiff')
  output = driver.Create(output_filename, hsize*256, vsize*256, 1, osgeo.gdal.GDT_Float32)
  output.GetRasterBand(1).WriteArray(raster)
  output.GetRasterBand(1).SetNoDataValue(nodata)
  output.SetGeoTransform(trans)
  output.SetProjection(srs.ExportToWkt())
  output.FlushCache()

if __name__ == '__main__':
  main(sys.argv)

免責事項とライセンス

プログラムは各自の責任のもとで利用されるものとし、何が起こっても一切の責任は負えません。実験や趣味に用いる分には自由に参考にしてもらって構いませんが、何かしらの業務やプロダクト生成の過程でこのプログラムを利用する場合はContact Informationから著作者に通知してください。

]]>
PythonとGDAL/OGRによる地理座標の投影座標への変換 https://tm23forest.com/contents/python-gdal-ogr-coordinatetransformation 2020-02-07T23:00:00+09:00 2020-02-07T23:00:00+09:00 GDALとはGeospatial Data Abstraction Libraryの略であり、もともとはC/C++のライブラリであるがPythonやJavaなどで使えるバインディングも用意されている。これを用いて地理座標を投影座標に変換する方法を書いておく。当然逆も可能。

準備

  • Python 3.8.1
  • GDAL 3.0.1

とう状況のもとで行った。Pythonが使えることは前提としてGDALについてはpipかもしくはWindowsであればUnofficial Windows Binaryからダウンロードしてインストールできる。

EPSGコードを使って変換

from osgeo import ogr, osr, gdal

src_srs, dst_srs = osr.SpatialReference(), osr.SpatialReference()
src_srs.ImportFromEPSG(4612) # from EPSG6668
dst_srs.ImportFromEPSG(6674) # to   EPSG6674
trans = osr.CoordinateTransformation(src_srs, dst_srs)

print(trans.TransformPoint(35, 135))
# (-110481.75032585638, -91280.6397650397, 0.0)

EPSG6668はJGD2011の地理座標、EPSG6674はJGD2011の平面直角座標に投影した座標の6系。これで北緯35度、東経135度という地理座標系を投影座標系に変換すると-110481.75032585638, -91280.6397650397という結果になる。

Shapefileから座標系を取得して変換

ファイルの指定する座標系に合わせたいという場合がある。

from osgeo import ogr, osr, gdal

shapefile_name = 'yourshapefile.shp'

shp = ogr.GetDriverByName('ESRI Shapefile').Open(shapefile_name, 0) # '0' means 'read'
src_srs, dst_srs = shp.GetLayer().GetSpatialRef(), osr.SpatialReference()
dst_srs.ImportFromEPSG(6674)
trans = osr.CoordinateTransformation(src_srs, dst_srs)

print(trans.TransformPoint(35, 135))
# (-110481.75032585638, -91280.6397650397, 0.0)

Shapefileといっても実際に座標系の情報が含まれているのは同名ファイルの拡張子.prjであるので、.prjを直接読み込んでImportFromWktで読み込むことも可能だがShapefileありきであればこの方がスマートそうである。少しだけ問題があって、TransformPointの引数の順番が緯度、経度の順番なのかその逆なのかがEPSGコードだけでは定まらない様子である。よって、場合によってはEPSG6668に合わせた順番である緯度、経度の順で引数を指定すると読み込んでいるShapefileによってはinfが返される。その出力側の座標系の範囲外の座標が渡されるとこうなるようだ。

この例では入力側の座標系をShapefileから読んでいるが出力側も同じようにしてShapefileから読み込める。ちなみにz座標はデフォルトで0になっているが指定することもできる。

複数の点を変換

TransformPointsを使うとできる。ただ、こうなってくるとShapefileごと座標系を変換して一時ファイルに書き出すとか、あるいはogr2ogrであらかじめ所望の座標系に変えておくかすればよいと思われる。引数はおそらくタプルの配列だろうけどよくわかっていないので放置しておく。気が向いたら(いずれ使うことがあれば)メモがてら載せたいと思う。

まとめ

ウェブメルカトル座標のバウンディングボックスの4点を変換したり、とか少数の点を変換する際には使えるのではないでしょうか。

]]>
LaTeX/amsmathによる数式環境を見本付きで https://tm23forest.com/contents/latex-amsmath-guide-with-svg-outputexample 2020-02-05T22:39:00+09:00 2020-02-05T22:39:00+09:00 誰よりも正しくLaTeX/amsmathの数式環境を使え!

行中数式

行中数式あるいはインライン数式とは文章の中に現れる数式のこと。使う場合には$ ... $\( ... \)を使う。

例1. $\lambda,x$をそれぞれ行列$A$の固有値,固有ベクトルとする.
例2. \(\lambda,x\)をそれぞれ行列\(A\)の固有値,固有ベクトルとする.

インライン数式の出力見本

1行のディスプレイ数式

ディスプレイ数式はテキストとは別に1行全体使って表示される数式でequation環境か\[ ... \]を用いる。\[ ... \]は数式番号が着かない。

二次方程式の解の公式は
\begin{equation}
x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}
\end{equation}

ディスプレイ数式の出力見本

ディスプレイ数式と直前のテキストの間に空行を入れると行頭1文字空いてしまうため注意する(段落が変わる場合を除く)。行頭が空く以外にもテキストとの間隔がおかしくなってしまう。ソースファイルの見た目上そこに空行を挟みたい場合は%によるコメントアウトを使う(%のみの行を挿入する)。

複数行数式を作るために\[ \]を連続使用しないようにする。数式の行間がおかしくなるので複数行数式にはgather環境などを使う。

以下は使っていはいけない。

  • \begin{math}\end{math}環境
  • \begin{displaymath}\end{displaymath}環境
  • \begin{eqnarray}\end{eqnarray}環境
  • $$ ... $$によるディスプレイ数式 これらはLaTeXの提供する数式環境が残っていて使えるだけであって、amsmathパッケージが意図するコマンドではないので混在させないようにする。数式内や数式間のスペースの広さや数式番号の配置がおかしくなるので使わないに越したことは無い。

複数行のディスプレイ数式

複数行数式にはそれ用のコマンドがあるので\[ \]などを連続して使用しないようにする。複数行数式で改ページしてもいい場合には\allowdisplaybreaks[1]をプリアンブルに書く。個別に指定したい場合はそれぞれの数式の行末に\displaybreak[1]を書く。ここに[1]は、必ずしも改ページしなくてもいいが必要あればやってもいいということを表す。度合いは[0][4]で指定できる。

複数の数式を独立に並べる:gather

ただ単に数式を中央揃えにして並べる、equation環境を連続使用したような配置になる。式番号も1行ずつそれぞれに着く。等号を縦一列にそろえなくても不親切でないような場合に使うとよいだろう。

\begin{gather}
(a+b)^2=a^2+2ab+b^2\\
(a+b+c)^2=a^2+b^2+c^2+2(ab+bc+ca)
\end{gather}

gather環境の出力見本

長い数式を複数行に分割:multline

項の数が多くて数式が横に長くなり過ぎる場合の数式の改行を左、真ん中、右という風にいい感じに配置してくれる。どういう風に?見ればわかる。

\begin{multline}
(x+y+z)^3=x^3+xy^2+xz^2+2x^2y+2xyz+2zx^2\\
+yx^2+y^3+yz^2+2xy^2+2y^2z+2xyz\\
+x^2z+y^2z+z^3+2xyz+2yz^2+2z^2x
\end{multline}

multlineの出力見本

複数の数式を整列させる:align, alignat, flalign

こちらは&(アンド)を指定した場所で複数行の数式を縦に整列させる。

次の連立方程式を解け.
\begin{align}
  5x+5y&=50 & x+y&=10\\
10x-2y &=-8 & 5x-y&=-4
\end{align}

alignの出力見本

次の連立方程式を解け.
\begin{alignat}{2}
  5x+5y&=50 &\qquad x+y&=10\\
10x-2y &=-8 &     5x-y&=-4
\end{alignat}

alignedの出力見本

次の連立方程式を解け.
\begin{flalign}
  5x+5y&=50 & x+y&=10\\
10x-2y &=-8 & 5x-y&=-4
\end{flalign}

flalignの出力見本

alignは縦にそろえて中央よりに配置。alignatは横の間隔を自分で指定する。この例では\quadを指定しているがこれを書かないと密着する。alignatの引数は数式の「カタマリ」の数を指定する(この例では2つのカタマリがある)。flalignはfull lengthの意味で、横幅をめいっぱい使う。

等号などの整列を保ったまま短いテキストを挿入するには\intertext{ ... }を使用する。

入れ子にする:gathered, multlined, aligned, alignedat

領域$D$を次のように定める.
\begin{equation}
D=\left\{
x,y,z\;\left|\;
  \begin{gathered}
    x^2+y^2\le 1\\
    z\ge 0,\;x+y+z\le 1
  \end{gathered}
\right.
\right\}
\end{equation}

gatheredの出力見本

\begin{equation}
\begin{multlined}
(x+y+z)^3
=x^3+xy^2+xz^2+2x^2y+2xyz+2zx^2\\
 +yx^2+y^3+yz^2+2xy^2+2y^2z+2xyz\\
 +x^2z+y^2z+z^3+2xyz+2yz^2+2z^2x
\end{multlined}
\end{equation}

multlinedの出力見本

前の例と比べて数式番号の位置が上下方向の中心なのに加えて真ん中の行の左右の配置が異なる。ちなみにmultline環境の各行の配置は\shoveleft{ ... }(...を左へずらす)、\shoveright{ ... }(...を右へずらす)で制御できる。

multlined環境と後で現れる rcases環境の使用には\usepackage{mathtools}が必要です。エラーが出た場合はtlmgr install mathtoolsしましょう。

\begin{equation}
\delta_{ij}=\left\{
  \begin{aligned}
  &1 &  (i&=j)\\
  &0 &  (i&\neq j)
  \end{aligned}
\right.
,\quad
\delta_{ij}=
  \begin{cases}
  1 &  (i=j)\\
  0 &  (i\neq j)
  \end{cases}
\end{equation}

alignedとcasesの出力見本

cases環境はこのような場合分けをalignedよりも短く描けますが内容は同じです。alignedatの使い方はきわめて上記に同じです。ここで上げた環境は、その部分の横幅が数式のサイズそのままになるので最初の例のD=の部分のように前後に別な数式を配置できます。splitとはやや挙動がことなります。

cases環境では数式がインライン数式モードになる(行中数式で用いられるサイズが小型のもの)。ディスプレイ数式モードにしたければ\displaystyleを明示的に指定するか、もしくはdcases環境を用いる。中括弧を右側にしたい場合にrcases環境もある。これらは\usepackage{mathtools}が必要である。

数式番号の制御

よく知られた基本的なこととしてアスタリスク*を環境名の末尾に付与した以下の環境では数式番号が割り振られない。文書の種類にもよるが基本的に使わなくてよいだろう。アスタリスクは無い前提でいく。

  • equation*
  • align*
  • multline*
  • alignat*
  • flalign*

複数行に対して1つの数式番号:split

次の連立微分方程式を解け.
\begin{equation}
  \begin{split}
  y_1' &= y_1+2y_2\\
  y_2' &= 2y_1+y_2
  \end{split}
\end{equation}

splitの出力見本

&は1つだけ使える。お気づきの通りこれはalignedでも実現できるうえに、現状equation内のsplitの外に別の数式があると横幅がはみ出すことがある。従って入れ子にするなら-ed系の環境を使ったほうが良いだろう。

部分的に数式番号を付けない, 変える:notag, tag

\begin{align}
\dot x(t)&=Ax(t)+Bu(t)\tag{A}\\
y(y)&=Cx(t)+Du(t)\notag
\end{align}

notagとtagの出力見本

数式番号の形式:theequation, numberwithin, subequations

自動的に割り振られる「tag」の形式を制御できる(LaTeXにおけるtagとlabelの使い分けが難しい)。これらは何かしらのテンプレートを使用していればそのテンプレートが指定する形式になっていることもあるのでその場合はいじらないこと。特にsubequationsの形式などを自分で指定したい場合。

\renewcommand{\theequation}{\thesection.\arabic{equation}}
\numberwithin{equation}{section}

\begin{document}
\section{状態空間モデル}
たぶんコレ.
  \begin{align}
  \dot{\bm{x}}(t)&=A\bm{x}(t)+Bu(t)\\
    y(y)&=C\bm{x}(t)+Du(t)
  \end{align}
より具体的に書くと,
  \begin{subequations}
  \begin{align}
    \dot x_1(t)&=x_2(t)\\
    \dot x_2(t)&=x_3(t)\\
    \dot x_3(t)&=-kx_1(t)-Dx_2(t)+b_3u(t)
  \end{align}
  \end{subequations}
\end{document}

subequationsの出力見本

見ての通りsubequationsで入れ子にした数式にはアルファベットによるtagが振られます。ちなみに、この例からは\numberwithinの効果はわかりません。意味は「equation」というカウンタを「section」ごとにリセットせよ、です。これが無いと次のセクションに行ったときに最初の数式が(2.4)となりますが、これによって(2.1)とすることが出来ます。\renewcommandによる\theequationの上書きは見た通りで、他のフォーマットにしたければ好きな区切り文字などが使えます。subsectionsubsubsectionも当然両方のコマンドで使えます。対応するのはthesubsectionthesubsubsectionとなります。subequations内の数式のフォーマットを変えたければ\renewcommand{\theequation}{\theparentequation\roman{equation}}subequations内で使うと「親の数式番号+ローマ数字の子カウンタ番号」とできます。

数式番号の参照:label, eqref

\begin{align}
\dot{\bm{x}}(t)&=A\bm{x}(t)+Bu(t)\label{ss1}\\
y(y)&=C\bm{x}(t)+Du(t)\label{ss2}
\end{align}
\eqref{ss1}, \eqref{ss2}を合わせて状態方程式という.
  \begin{subequations}\label{ss3}
    \begin{align}
    \dot x_1(t)&=x_2(t)\\
    \dot x_2(t)&=x_3(t)\\
    \dot x_3(t)&=-kx_1(t)-Dx_2(t)+b_3u(t)\tag{A}\label{di}
    \end{align}
  \end{subequations}
\eqref{ss3}の\eqref{di}は運動方程式である.

eqrefの使い方

普通の\refと違ってカッコが勝手につくという違いがある。

参照文献

]]>
TeXstudioのインストールと初期設定 https://tm23forest.com/contents/texstudio-install-and-configure 2020-01-28T11:02:00+09:00 2020-01-28T11:02:00+09:00

TeX/LaTeXの環境は通常コンソールからしか使えないため、例えばTeXLiveなどを入れただけではPDFにして開くまでにいくつかのコマンドを実行する必要があるため、逐一それらを実行しなければいけないのは面倒である。多くのディストリビューションにはTeXworksというエディタが同梱されているが、多分機能的にその上位互換でそんなにデメリットも見当たらないTeXstudioをインストールして速攻で使える用にする設定をメモしておく。

TeXstudioのインストール

https://www.texstudio.org/

からDownload Nowしてポチポチするだけ。

TeXstudioのホームページのスクリーンショット

TeXstudioのインストーラ最初の画面

TeXstudioのインストーラ2番目の画面

TeXstudioのインストーラ3番目の画面

TeXstudioのインストーラ4番目の画面

あっさり終わった。ちなみにアップデート時には一度アンインストールしなくてもそのまま新しいバージョンのインストーラを起動すれば上書きしてくれるし前のバージョンの設定もそのまま引き継がれる。

TeXstudioの初期設定【必須】

「オプション」→「TeXstudioの設定」→「コマンド」を開き、「latex」のところを「platex」か「uplatex」に変更する。

オプション、TeXstudioの設定、コマンド

最近のバージョンのTeXstudioはPCのLaTeX環境をチェックしてくれる模様。TeXLiveなどでインストールして各所コマンドにパスが通っている場合にはデフォルトでこんな感じになっている。<unknown>になっているところはインストールされていないのでtlmgrなどで追加すれば使えるようになるだろう。この画面で1点だけ変更が必要な点は日本語の文書をコンパイルするときには使用するLaTeXコマンドを「platex」もしくは「uplatex」に変更すること。「uplatex」の方はUnicode対応のバージョン。利用するテンプレートクラスなどによっては指定がされているので適宜正しい方を使ようにする。

次に左のペインから「ビルド」を選択し「ビルド&表示」を「DVI->PDFチェーン」、既定のコンパイラを「LaTeX」にする。

オプション、TeXstudioの設定、ビルド

これでplatex(uplatex)でコンパイルしてからdvipdfmxでpdfにして表示という動作をF5キーを押すだけでやってくれる。もちろんそれ以外のコンパイル方法を既定にしたい場合はここから再度選択する。そこまで頻繁に変更項目ではないだろう。

以上で最低限必要な設定は終わり。次のように日本語の文書がコンパイルできる。

TeXstudioのテスト画面

TeXstudioの初期設定【オプション】

これだけの状態だとバックスラッシュを打っただけでコマンドの補間が始まって、かつこのコマンドを補間する動作が結構重かったりして使いにくいのでその辺りの設定をいじっていく。

TeXstudioの設定、一般、Meiryo UIにした

UIフォントはデフォルトで「MS UI Gothic」だったので「Meiryo UI」に変更した。ちなみにMeiryo UIにすると次回にこの「TeXstudioの設定」画面を開いたときにウィンドウの初期サイズが画面からはみ出る。

エディタのフォント、字下げ、行番号、インラインチェック

フォントは「源ノ角ゴシック」にした。字下げは維持する、タブはスペースに置換。タブを2文字にできないのが残念。行番号はエラーメッセージが出たときのために表示、インラインチェックの文法と構文はうるさいのでチェックを外す。

コマンドの補間を入力時に開始しない、ツールチップヘルプとプレビューはオフ

「LaTeXコマンド入力時に自動的に補間を開始」は動作が遅くなるのでチェックを外した。VisualStudioCodeなど他の近代的なエディタであれば同等の機能を特に重くならずにできているので何かしら実装上の問題があるのだろう。もし補間したい場合はコマンドを入力している途中でCtrl+Spaceを入力すればよい。「ツールチップヘルプ」と「ツールチッププレビュー」も編集中に重くなってしまうので今のところは外した。

]]>
QGIS3.10 LTS長期サポート版をインストールした https://tm23forest.com/contents/qgis-3.10-lts-installed 2020-01-24T23:10:00+09:00 2020-01-24T23:10:00+09:00

QGIS(https://www.qgis.org/)の3.10をインストールしたことをメモしておく。結論から言うと重複インストールすることなく勝手に前バージョンをアンインストールしてくれるのでHappyだった。

QGIS 3.10.1をインストール

前バージョンが残っていると次のような流れになった。ちなみにCoruñaはn以降が文字化けしている。なんでそんな文字使うねん。

QGISインストーラ1番目の画面

QGISインストーラ2番目の画面

QGISインストーラ3番目の画面

QGISインストーラ4番目の画面

QGISインストーラ5番目の画面

QGISインストーラ6番目の画面

QGISインストーラ7番目の画面

QGISインストーラ8番目の画面

QGISインストーラ9番目の画面

QGISインストーラ10番目の画面

QGISインストーラ11番目の画面

QGISインストーラ12番目の画面

QGISインストーラ13番目の画面

最後に前バージョンをアンインストールしてくれるようだ。ありがたい。

QGIS 3.10.2のインストール

気づいたらLTS(Long Term Support)版の3.10.2がリリースされていたのでインストールした。このときも前バージョンをアンインストールすることなくインストーラを起動したところ次のようになった。前回は無かったがWindows Defender SmartScreenが発動してしまった。

Windows Defender SmartScreenが起動した画面

Windows Defender SmartScreenが起動した画面(詳細情報を押した)

QGIS3.10.2インストーラの既存インストールの確認

QGIS3.10.2インストーラ1番目の画面

QGIS3.10.2インストーラ2番目の画面

QGIS3.10.2インストーラ3番目の画面

QGIS3.10.2インストーラ4番目の画面

QGIS3.10.2インストーラ5番目の画面

QGIS3.10.2インストーラ6番目の画面

QGIS3.10.2インストーラ7番目の画面

QGIS3.10.2インストーラ8番目の画面

QGIS3.10.2インストーラ9番目の画面

QGIS3.10.2インストーラ10番目の画面

QGIS3.10.2インストーラ11番目の画面

QGIS3.10.2インストーラ12番目の画面

QGIS3.10.2インストーラ13番目の画面

QGIS3.10.2インストーラ14番目の画面

こんどは先に前バージョンがアンインストールされた。

]]>
グラフラプラシアンの解釈と接続行列 https://tm23forest.com/contents/graph-laplacian-incidence-matrix 2020-01-24T12:34:00+09:00 2020-01-24T12:34:00+09:00 グラフラプラシアン(Graph Laplacian)がなぜラプラシアンと呼ばれるのかを、普通の関数に対するラプラシアンオペレータと比較して直観的に説明する。ついでに普通の関数のときにラプラシアンとまとめて紹介されるナブラについて、それに対応する接続行列について紹介する。

グラフラプラシアンってカタカナで書くとめちゃくちゃ読みにくいな。Graph Laplacianなら全然読めるけど両方カタカナで空白つけないのが良くない。こういうのGraphが漢字なら日本語にしたとき丁度よいのだが。

Graph Laplacianの定義

無向グラフ\(G=(V,E)\)\(|V|=n\)に対して次の\(n\times n\)行列\(L(G)\)をGraph Laplacianといい、グラフ\(G\)が何を指すか明らかなときは単に\(L\)とかく。もしくは\(\mathcal{L}\)のときもある。

\[ \begin{equation} L(G)_{ij}=\left\{\begin{matrix} \deg(v_i) & (i=j)\\ -1 & (i\neq j\text{ and }\{v_i,v_j\}\in E)\\ 0 & (i\neq j\text{ and }\{v_i,v_j\}\notin E) \end{matrix}\right. \end{equation}\]

\(L(G)_{ij}\)はGraph Laplacianの\((i, j)\)要素を示し、\(v_i\in V\)\(i\)番目の頂点を示す。\(\deg\)は頂点の次数。念のため次数とはその頂点につながっている辺の数のこと。

グラフの次数行列\(D\)と隣接行列\(A\)を用いるとグラフラプラシアンは次のようにかける。

\[ \begin{equation} L=D-A \end{equation}\]

いずれも\(G\)によって決まる行列だが明らかなので省略する。次数行列(degree matrix)とは対角成分に頂点の次数を持つような対角行列(非対角成分がゼロ)、隣接行列(adjacency matrix)とは\(\{0,1\}\)のみを要素にもつ行列で\(i,j\)間のエッジが存在すれば\(1\)、そうでなければ\(0\)を取る。

\[ \begin{gather} D_{ij}=\left\{\begin{matrix} \deg(v_i) & (i=j)\\ 0 & (i\neq j) \end{matrix}\right.\\ A_{ij}=\left\{\begin{matrix} 1 & (\{v_i,v_j\}\in E)\\ 0 & (\{v_i,v_j\}\notin E) \end{matrix}\right.\\ \end{gather}\]

Graphラプラシアンの例

無向グラフの例

こういうグラフの場合を考えると、それぞれ次のようになる。

\[ \begin{gather} D=\begin{bmatrix} 2 & 0 & 0 & 0 & 0 & 0\\ 0 & 4 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 0 & 0 & 0\\ 0 & 0 & 0 & 3 & 0 & 0\\ 0 & 0 & 0 & 0 & 3 & 0\\ 0 & 0 & 0 & 0 & 0 & 2 \end{bmatrix},\quad A=\begin{bmatrix} 0 & 1 & 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 1 & 1 & 0\\ 1 & 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 1 & 1\\ 0 & 1 & 0 & 1 & 0 & 1\\ 0 & 0 & 0 & 1 & 1 & 0 \end{bmatrix}\\ L=\begin{bmatrix} 2 & -1 & -1 & 0 & 0 & 0\\ -1 & 4 & -1 & -1 & -1 & 0\\ -1 & -1 & 2 & 0 & 0 & 0\\ 0 & -1 & 0 & 3 & -1 & -1\\ 0 & -1 & 0 & -1 & 3 & -1\\ 0 & 0 & 0 & -1 & -1 & 2 \end{bmatrix} \end{gather}\]

ここまでは普通の話。(ここからがマグマなんです)

Laplacianといわれる理由

名前からして実関数に対するラプラシアンオペレータと関係とおかしいだろうということで実際関係していることを直観的に示していく。簡単な場合を使って2変数関数\(f(x,y)\)に対するラプラシアンを考える。

\[ \begin{equation} \Delta\equiv\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2} \end{equation}\]

とおくと、これは関数に作用して関数値を返す演算子となり

\[ \begin{equation} \Delta f(x,y)= \frac{\partial^2 f(x,y)}{\partial x^2}+\frac{\partial^2 f(x,y)}{\partial y^2} \end{equation}\]

となる。こっちのラプラシアンをGraph Laplacianと結びつけるために、この値を離散化して数値計算しようと試みてみる。もっとも単純な離散化のやり方は次だろう。

\[ \begin{align} \frac{\partial^2 f(x,y)}{\partial x^2} &=\frac{1}{h}\left( \frac{f(x+h,y)-f(x,y)}{h}-\frac{f(x,y)-f(x-h,y)}{h} \right)\\ &=\frac{f(x+h,y)-2(x,y)+f(x-h,y)}{h^2} \end{align}\]

小さい\(h\)をとって素朴に前進差分と後退差分を使って勾配を計算し、さらにその差から2階微分を求めている。\(y\)についても同じことをやると次がすぐにわかる。

\[ \begin{multline} \frac{\partial^2 f(x,y)}{\partial x^2} +\frac{\partial^2 f(x,y)}{\partial y^2}\\ =\frac{f(x+h,y)+f(x-h,y)+f(x,y+h)+f(x,y-h)-4f(x,y)}{h^2} \end{multline}\]

もう言いたいこと分かったと思うがこれは実平面を拡張したような次のグラフに対応している。

Graph Laplacianのイメージ

直観的に符号が定義と異なるが、このグラフから向きという概念をなくすということを考えればいい。このグラフはもともとの実平面の性質から数字の大小という意味で「向き」をもっているが、無向グラフの世界ではそのような概念は無いのでラプラシアンをグラフに作用させることを考えるにあたって「自分引く相手の総和」を用いることはまぁアナロジーとして許せる範囲だろう。

ついでにグラフ上の調和関数を考える時はゼロかゼロでないかが重要なので符号が逆転していても構わない。ということでGraph Laplacianは普通のラプラシアンのグラフに対する拡張っぽくなっていることが気持ち的には分かった。もっとうまい説明は無数にある気がしているが次にいく。

グラフに対するLaplacianと調和関数

ここまでの話を踏まえてグラフの頂点を定義域とする関数\(f:V\to\mathbb{R}\)を考えると、実関数に対するラプラシアンから得た「自分引く相手の総和」というアナロジーから次の値をおもいつく。

\[ \begin{equation} \Delta f(v)=\sum_{(v,e)\in E}(f(v)-f(e))\qquad(v\in V) \end{equation}\]

書き換えると

\[ \begin{equation} \Delta f(v_i)=\sum_{j=1}^n \mathbb{I}(\{v_i,v_j\}\in E)(f(v_i)-f(v_j)) \end{equation}\]

ここで、\(\mathbb{I}\)は指示関数(indicator function)で引数内が成り立つなら1、それ以外の場合にはゼロとなる便利なヤツ。実関数にたいするLaplacianがある種の演算子であったことを考えるとグラフに対するLaplacianも\(f\)という関数を介することなく定義されるべきだろう。そこで2番目のラプラシアンの定義式より、各\(f(v_j)\)の係数を抜き出しておいたものを\(L\)とおくとこれはGraph Laplacianの定義そのもの。これは式をいじるまでもなく簡単に確かめられる。あとは

\[ \begin{equation} Lf(v_i)\equiv\sum_{j=1}^n L_{ij}f(v_j),\quad i=1,2,\dots,n \end{equation}\]

とでも定義してやれば同じように使える。注意しておくとこれは行列の計算でもなんでもなくてただ\(Lf\)をこう定義しただけ。ようは\(L\)\(f\)によらずに決定されていてこの\(L\)を用いることで\(i=1,\dots,n\)に対する\(f(v_i)\)のラプラシアンが計算できていればいい。

関数\(f:V\to\mathbb{R}\)についても調和関数かどうかがこれによって定義できる。\(f:V\to\mathbb{R}\)がグラフ\(G\)上で調和であるとは\(i=1,2,\dots,n\)について

\[ \begin{equation} Lf(v_i)=0 \end{equation}\]

を満たすこと。まぁそらそうだよな。いくらでも他の書き方は存在します。

勾配に対応する接続行列

実関数の場合には\(\Delta=\nabla\cdot\nabla\)というような式があったがこれに対応するものが有向グラフの場合にはある。もう一度いうと有向グラフの場合にはある。無向グラフには向きという概念が無いから勾配という概念に対応するものも無くてまぁ不思議ではないね。

接続行列の定義

無向グラフ\(G=(V,E)\)\(|V|=n,|E|=m\)に対する接続行列(incidence matrix)\(B\)はそれぞれの要素が次のようになる\(n\times m\)行列

\[ \begin{equation} B_{ij}=\left\{\begin{matrix} 1 & (v_i\in e_j)\\ 0 & (v_i\notin e_j) \end{matrix}\right. \end{equation}\]

\(v_i, e_j\)はそれぞれ\(i\)番目の頂点と\(j\)番目の辺のこと。有向グラフの場合には次のようにやるのが自然だろう。

\[ \begin{equation} B_{ij}=\left\{\begin{matrix} 1 & (e_j=(v, v_i)\text{ for some }v)\\ -1 & (e_j=(v_i, v)\text{ for some }v)\\ 0 & (\text{otherwise}) \end{matrix}\right. \end{equation}\]

少しだけ説明すると、有向グラフでは辺\(e=\{v, w\}\)において\(v\)\(w\)の順番に意味があるので集合ではなくベクトルのような記法\(e=(v,w)\)を用いることが多い。集合の場合には要素の順番は考慮されないがベクトルの場合は考慮するので。ただ、完全にベクトルというわけではなく有向グラフの場合にも\(v\in e\)なんかの使い方はゆるす。

接続行列の例

例をみたほうが早いだろう。Graphラプラシアンのときと同じグラフを用いる。...

]]>
衝突積分の対称関係式 https://tm23forest.com/contents/boltzmann-collision-integral-symmetric 2020-01-16T20:20:00+09:00 2020-01-16T20:20:00+09:00 ボルツマン方程式の衝突項まわりで現れる衝突積分の対象関係式を導出する。

衝突項

分子運動論だとかに現れるボルツマン(Boltzmann)方程式の右辺のことを衝突項という。(重力などの全ての粒子に一様に働くような)外力が無い場合のボルツマン方程式は

\[ \begin{equation} \frac{\partial f}{\partial t}+\xi_i\frac{\partial f}{\partial X_i}=\Delta f_{\mathrm{collision}} \end{equation}\]

みたいな感じで右辺が衝突項になっている。この書き方は総和の規約が使われているので同じ添え字のでてくる左辺第2項は和をとるものとする。2体衝突だとか3体衝突だとかの違いで衝突項の取り方にもいろいろあるらしいがとにかく分子同士の衝突を表す項のことを衝突項といっている。考えるのはこの衝突項について。

衝突積分作用素

ここでは2体衝突のみを考慮した次のタイプの衝突項を考える。

\[ \begin{equation} J(f,f)=\frac{\sigma^2}{2m} \int_{\mathbb{R}^3\times\mathbb{S}^2} [f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}')-f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})] \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_* \end{equation}\]

ここに、\(f(\xi)=f(\boldsymbol{X},\boldsymbol{\xi},t)\)は分子の位置・速度に関する密度関数、\(\boldsymbol{X},\boldsymbol{\xi},t\)はそれぞれ位置、速度、時刻で位置と時刻は出てこないので省略してかかれている。太字は\(\mathbb{R}^3\)である。また、

  • \(\sigma\)は分子の直径、\(m\)は分子の質量
  • \(\boldsymbol{e}\in\mathbb{S}^2\)は単位球上の点であり\(|e|=1\)\(d\Omega(e)\)\(e\)方向の立体角素(面積素の3次元角バージョン)

さらに

\[ \begin{align} \boldsymbol{\xi}'&\equiv\boldsymbol{\xi}'(\boldsymbol{\xi},\boldsymbol{\xi}_*, \boldsymbol{e}) =\boldsymbol{\xi}+[(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot \boldsymbol{e}]\boldsymbol{e}\\ \boldsymbol{\xi}'_*&\equiv\boldsymbol{\xi}'(\boldsymbol{\xi},\boldsymbol{\xi}_*, \boldsymbol{e}) =\boldsymbol{\xi}_*-[(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot \boldsymbol{e}]\boldsymbol{e} \end{align}\]

とする。\(\boldsymbol{\xi}'\)\(\boldsymbol{\xi}'_*\)の2つはいずれも\(\boldsymbol{\xi},\boldsymbol{\xi}_*,\boldsymbol{e}\)によってきまる。\(J(f,f)\)の式に直接入れ込むこともできるが式が長くなるので。ついでに\(J(f,f)\)自体はいぜんとして\(\boldsymbol{X},\boldsymbol{\xi},t\)についての関数であるので注意しておく。この衝突項を用いた(重力などの全ての粒子に一様に働くような)外力が無い場合のボルツマン方程式はすごく丁寧にかくと次のようにかける。

\[ \begin{equation} \frac{\partial f(\boldsymbol{X},\xi,t)}{\partial t} +\sum_{i=1}^3\xi_i\frac{\partial f(\boldsymbol{X},\xi,t)}{\partial X_i} =J[f,f](\boldsymbol{X},\xi,t) \end{equation}\]

外力無しの場合なので他で見るかもしれない形とくらべて左辺の第3項が無い。こう見ると確かに\(J\)\(f\)に作用する積分作用素。ちなみに\(J(f)\)でなく\(J(f,f)\)になっている理由は

\[ \begin{equation} J(f,g)\equiv \frac{1}{4m}\int_{\mathbb{R}^3\times\mathbb{S}^2} [f(\boldsymbol{\xi}')g(\boldsymbol{\xi}'_*) +f(\boldsymbol{\xi}'_*)g(\boldsymbol{\xi}') -f(\boldsymbol{\xi})g(\boldsymbol{\xi}_*) -f(\boldsymbol{\xi}_*)g(\boldsymbol{\xi})] \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_* \end{equation}\]

という一般化された衝突積分を考えると都合がいい場合がまぁ別な機会にあるから。\(f=g\)としてみると確かに最初にあげた衝突積分(collision integral)に一致していることがわかる。

対称関係式

任意の関数\(\phi=\phi(\boldsymbol{\xi})\)についてのつぎの関係を対称関係式という。正確にはここでいうところの衝突積分作用素についての対称関係式だが例えば線形化された衝突作用素(linearized collision integral)なんかでも同じことが出来る。

\[ \begin{multline} \int_{\mathbb{R}^3}\phi(\boldsymbol{\xi})J(f,f)d\boldsymbol{\xi}\\ =\frac{\sigma^2}{8m}\int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*) -\phi(\boldsymbol{\xi}')-\phi(\boldsymbol{\xi}'_*)) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi} \end{multline}\]

これを示す。まずは左辺を書き下す。

\[ \begin{align} =\frac{\sigma^2}{2m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} \phi(\boldsymbol{\xi}) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi} \end{align}\]

この式は\(\boldsymbol{\xi},\boldsymbol{\xi}_*\)の両方で積分しているので、被積分関数の\(\boldsymbol{\xi}\)\(\boldsymbol{\xi}_*\)を形式的に入れ替えたものも同じ値になるはずである。ところが\(\phi(\boldsymbol{\xi})\)に関する部分以外はいずれも\(\boldsymbol{\xi},\boldsymbol{\xi}_*\)の対称式の形をしている。これに注目するとこの式は次のように変形できる。

\[ \begin{align} =\frac{\sigma^2}{2m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} \frac{\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)}{2} (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi} \end{align}\]

つまり、同じ項を2で割って分割し、片方の変数を形式的に入れ替えてからもう一度足し合わせてくくったものである。この操作をもう一度行う。

\[ \begin{align} =&\frac{\sigma^2}{2m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} \frac{\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)+\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)}{4} (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi}\\ =&\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*) +\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi} \end{align}\]

結構近づいてきたことがわかる。ここで項を前半と後半に分割する。

\[ \begin{align} =&\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*) +\textcolor{red}{\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)}) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi} \end{align}\]

前半は完成しているのでこの赤で示した部分に注目して\(\boldsymbol{\xi},\boldsymbol{\xi}_*\)の積分を\(\boldsymbol{\xi}',\boldsymbol{\xi}'_*\)の積分に書き換えるという操作を行う。そのためにはまずヤコビアンを計算するのだが、実はこれは6変数の積分でめんどうなので結果をしめすと

\[ \begin{gather} \left|\frac{\partial(\boldsymbol{\xi}',\boldsymbol{\xi}'_*)}{\partial(\boldsymbol{\xi},\boldsymbol{\xi}_*)}\right| =1\\ \therefore\quad d\boldsymbol{\xi}_*d\boldsymbol{\xi} =d\boldsymbol{\xi}'_*d\boldsymbol{\xi}' \end{gather}\]

ちゃんとやるなら物理的な対称性から一般性を失わずに\(\boldsymbol{e}\)を特定の成分だけ1のベクトルにするとか、極座標つかって真面目にやるとかすればいいんじゃなかろうか。これについては物理的な意味を考えた方がわかると思うが図が必要なのでやめる。次に

\[ \begin{align} \boldsymbol{\xi}'&\equiv\boldsymbol{\xi}'(\boldsymbol{\xi},\boldsymbol{\xi}_*, \boldsymbol{e}) =\boldsymbol{\xi}+[(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot \boldsymbol{e}]\boldsymbol{e}\\ \boldsymbol{\xi}'_*&\equiv\boldsymbol{\xi}'(\boldsymbol{\xi},\boldsymbol{\xi}_*, \boldsymbol{e}) =\boldsymbol{\xi}_*-[(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot \boldsymbol{e}]\boldsymbol{e} \end{align}\]

という式の差を取って両辺に\(\boldsymbol{e}\)かけると

\[ \begin{align} \boldsymbol{\xi}'-\boldsymbol{\xi}'_* &=\boldsymbol{\xi}-\boldsymbol{\xi}_* +2[(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot \boldsymbol{e}]\boldsymbol{e}\\ (\boldsymbol{\xi}'-\boldsymbol{\xi}'_*)\cdot\boldsymbol{e} &=-(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot\boldsymbol{e} +2(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot \boldsymbol{e}\\ (\boldsymbol{\xi}'-\boldsymbol{\xi}'_*)\cdot\boldsymbol{e} &=(\boldsymbol{\xi}_*-\boldsymbol{\xi})\cdot\boldsymbol{e}\\ \end{align}\]

\(\boldsymbol{e}\cdot\boldsymbol{e}=1\)をつかえるのでこうなる。地味にアスタリスクがついてるほうとの符号の付き方が入れ替わっていることに注目。これを踏まえると

\[ \begin{align} \boldsymbol{\xi}& =\boldsymbol{\xi}'+[(\boldsymbol{\xi}'_*-\boldsymbol{\xi}')\cdot \boldsymbol{e}]\boldsymbol{e}\\ \boldsymbol{\xi}_*& =\boldsymbol{\xi}'_*-[(\boldsymbol{\xi}'_*-\boldsymbol{\xi}')\cdot \boldsymbol{e}]\boldsymbol{e} \end{align}\]

がすぐにわかる。移項しただけ。これで赤字のところを変数変換する準備はできた。

\[ \begin{align} &\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\textcolor{red}{\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)}) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}_*-\boldsymbol{\xi})\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi}\\ &=\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\textcolor{red}{\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)}) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}'-\boldsymbol{\xi}'_*)\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}'_*d\boldsymbol{\xi}' \end{align}\]

といっても書き下さなければ最後の方がかわっただけでそんなに変数変換をやった感じはない。ここでも対称性を利用して全てのダッシュ付きとダッシュ無しを入れ替える。すると

\[ \begin{align} &=\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\textcolor{red}{\phi(\boldsymbol{\xi})+\phi(\boldsymbol{\xi}_*)}) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}'-\boldsymbol{\xi}'_*)\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}'_*d\boldsymbol{\xi}'\\ &=\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\textcolor{red}{\phi(\boldsymbol{\xi}')+\phi(\boldsymbol{\xi}'_*)}) (f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi}) -f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}')) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}-\boldsymbol{\xi}_*)\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi}\\ &=\frac{\sigma^2}{8m} \int_{\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{S}^2} (\textcolor{red}{-\phi(\boldsymbol{\xi}')-\phi(\boldsymbol{\xi}'_*)}) (f(\boldsymbol{\xi}'_*)f(\boldsymbol{\xi}') -f(\boldsymbol{\xi}_*)f(\boldsymbol{\xi})) \left|\boldsymbol{e}\cdot(\boldsymbol{\xi}-\boldsymbol{\xi}_*)\right| d\Omega(\boldsymbol{e})d\boldsymbol{\xi}_*d\boldsymbol{\xi} \end{align}\]

ほんとに形式的に逆にしただけね。一応\(\boldsymbol{\xi},\boldsymbol{\xi}_*\)の中身もいれかわっていることに注意。すべて入れ替えたうえでもの変数の使い方と変わっていなくて唯一変わったのは赤字の部分だけ。これをもとの部分にいれると最初の対称関係式を得る。

]]>