Published Date : 2019年8月18日6:04
Python学習者のためのNim Part 3
Learning Nim for Python Users Part 3
This blog has an English translation
前回の続きと150時間の壁
Continuation of last blog and 150 hour wall.
前回の続きです。 さっそく、前回ダウンロードした、東京都の緯度経度の情報が入っているCSVファイルを計算しやすいように前処理をしていきます。
This is a continuation of the previous blog. I'm going to preprocess the CSV file that I downloaded last time and contains the latitude and longitude of Tokyo so that it's easy to calculate.
東京都の番地までをキーとして、緯度経度をバリューとしたJSONファイルを作っていきます。
Using the address of Tokyo as a key, we create a JSON file with the value of latitude and longitude.
Python Code [createjson.py]
import pandas as pd
# 文字コードをCP932（Shift-JIS）でCSVファイルを読み込み、データフレームオブジェクトにする。
# Read the CSV file with CP 932 (Shift-JIS) to make it a data frame object.
df = pd.read_csv('13_2018.csv', sep=',', encoding="cp932")
# 都道府県から番地までをstr.catメソッドで一つの文字列として連結させます。
# The str.cat method concatenates the address from the state as one string.
addresses = df["都道府県名"].str.cat(df["市区町村名"]).str.cat(df["大字・丁目名"]).str.cat(df["街区符号・地番"].astype(str)).tolist()
# 緯度経度の座標情報をリストにします。
# Lists latitude-longitude coordinate information.
coordinates = [[df.loc[i, "緯度"], df.loc[i, "経度"]] for i in range(len(addresses))]
# UTF-8にエンコードして、JSONファイルとして書き出す。
# Encode in UTF -8 and export as JSON file.
add_dict = {addresses[i]:[cd[0],cd[1]] for i,cd in enumerate(coordinates)}
with open("tokyoCoordinates.json", "w", encoding="utf-8") as f:
# ensure_ascii=False、これは日本語をJSONファイルに書き出すと、文字コードを出力してしまうのを防ぐためです。
# ensure_ascii=False, which prevents the output of character codes when Japanese is written to a JSON file.
json.dump(add_dict, f, indent=2, ensure_ascii=False)
それではPythonで全ての番地に対して、全ての東京の住所の座標から近い順TOP10を計算して出力するプログラムを書いていきます。
I will write a program in Python that calculates and outputs the top 10 addresses in order of proximity to the coordinates of all addresses in Tokyo for all addresses.
Pythonだと座標から距離を計算できる素晴らしいライブラリがあるので、簡単にできますが、nimで使えるか分からないので、こちらの人のブログ記事を参考にさせて頂きます。サンキューです！
Python has a great library that can calculate distances from coordinates, so it's easy to do, but I don't know if it works with nim, so I'll refer to this person's blog post. Thank you!
[Python]緯度経度から2地点間の距離と方位角を計算する
それから、タイトルを見ても分かる通り、膨大な計算時間になります。ご了承ください。
And as you can see from the title, it takes a lot of calculation time. Thank you for your understanding.
Python Code [vincenty.py]
# Classファイルにするため1％改良を加えましたが、99％は本家のコピペです。
# I made 1% improvements to make it a Class file, but 99% of it was copied and pasted from the original.
from math import radians,atan,sin,cos,sqrt,atan2,degrees,tan,pi
class Vincenty:
def __init__(self):
self.lat1 = 0.0
self.lon1 = 0.0
self.lat2 = 0.0
self.lon2 = 0.0
self.ellipsoid = 1
# 楕円体
# ellipsoid
self.ELLIPSOID_GRS80 = 1 # GRS80
self.ELLIPSOID_WGS84 = 2 # WGS84
# 楕円体ごとの長軸半径と扁平率
# major radius and flatness of each ellipsoid
self.GEODETIC_DATUM = {
self.ELLIPSOID_GRS80: [
6378137.0, # [GRS80]長軸半径
1 / 298.257222101, # [GRS80]扁平率
],
self.ELLIPSOID_WGS84: [
6378137.0, # [WGS84]長軸半径
1 / 298.257223563, # [WGS84]扁平率
],
}
# 反復計算の上限回数
# Maximum Iterations
self.ITERATION_LIMIT = 1000
'''
Vincenty法(逆解法)
Vincenty method (inverse solution)
2地点の座標(緯度経度)から、距離を計算する
Calculate distance from 2 point coordinates (Latitude Longitude)
:param lat1: 始点の緯度 Start Latitude
:param lon1: 始点の経度 Start Longitude
:param lat2: 終点の緯度 End Latitude
:param lon2: 終点の経度 End longitude
:param ellipsoid: 楕円体 ellipsoid
:return: 距離 Distance
'''
def vincenty_inverse(self):
# 差異が無ければ0.0を返す
# Returns 0.0 if there are no differences
if self.lat1 == self.lat2 and self.lon1 == self.lon2:
return 0.0
# 計算時に必要な長軸半径(a)と扁平率(ƒ)を定数から取得し、短軸半径(b)を算出する
# The major axis radius (a) and the flatness rate (ƒ) required for calculation are obtained from the constant, and the minor axis radius (b) is calculated.
# 楕円体が未指定の場合はGRS80の値を用いる
# If no ellipsoid is specified, the value of GRS 80 is used.
a, ƒ = self.GEODETIC_DATUM.get(self.ellipsoid, self.GEODETIC_DATUM.get(self.ELLIPSOID_GRS80))
b = (1 - ƒ) * a
φ1 = radians(self.lat1)
φ2 = radians(self.lat2)
λ1 = radians(self.lon1)
λ2 = radians(self.lon2)
# 更成緯度(補助球上の緯度)
# Further Latitude (Latitude on auxiliary sphere)
U1 = atan((1 - ƒ) * tan(φ1))
U2 = atan((1 - ƒ) * tan(φ2))
sinU1 = sin(U1)
sinU2 = sin(U2)
cosU1 = cos(U1)
cosU2 = cos(U2)
# 2点間の経度差
# Longitude difference between two points
L = λ2 - λ1
# λをLで初期化
# Initialize λ with L
λ = L
# 以下の計算をλが収束するまで反復する
# Repeat the following calculations until λ converges
# 地点によっては収束しないことがあり得るため、反復回数に上限を設ける
# There is an upper limit on the number of iterations because it may not converge at certain points.
for _ in range(self.ITERATION_LIMIT):
sinλ = sin(λ)
cosλ = cos(λ)
sinσ = sqrt((cosU2 * sinλ) ** 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cosλ) ** 2)
cosσ = sinU1 * sinU2 + cosU1 * cosU2 * cosλ
σ = atan2(sinσ, cosσ)
sinα = cosU1 * cosU2 * sinλ / sinσ
cos2α = 1 - sinα ** 2
cos2σm = cosσ - 2 * sinU1 * sinU2 / cos2α
C = ƒ / 16 * cos2α * (4 + ƒ * (4 - 3 * cos2α))
λʹ = λ
λ = L + (1 - C) * ƒ * sinα * (σ + C * sinσ * (cos2σm + C * cosσ * (-1 + 2 * cos2σm ** 2)))
# 偏差が.000000000001以下ならbreak
# Break if deviation is .000000000001 or less
if abs(λ - λʹ) <= 1e-12:
break
else:
# 計算が収束しなかった場合はNoneを返す
# Returns None if the calculation did not converge
return None
# λが所望の精度まで収束したら以下の計算を行う
# Once λ has converged to the desired accuracy, the following calculations are performed
u2 = cos2α * (a ** 2 - b ** 2) / (b ** 2)
A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
Δσ = B * sinσ * (cos2σm + B / 4 * (cosσ * (-1 + 2 * cos2σm ** 2) - B / 6 * cos2σm * (-3 + 4 * sinσ ** 2) * (-3 + 4 * cos2σm ** 2)))
# 2点間の楕円体上の距離
# Distance on ellipsoid between two points
s = b * A * (σ - Δσ)
# 計算した距離を返す
# return calculated distance
return s
続いて、実際に計算をするコードを書いてきます。
Next, I'll write the code to actually do the calculation.
Python Code [calc_distance.py]
from vincenty_class import Vincenty
import time
import json
with open("tokyoCoordinates.json","r",encoding="utf-8") as f:
add_dict = json.load(f)
vincenty = Vincenty()
start1 = time.time()
# 二回ループさせるとう暴挙にでる。
# Loop twice. That's pretty stupid code.
for key, value in add_dict.items():
start2 = time.time()
temp_dict = {}
vincenty.lat1 = value[0]
vincenty.lon1 = value[1]
for key2,value2 in add_dict.items():
if key == key2:
pass
else:
vincenty.lat2 = value2[0]
vincenty.lon2 = value2[1]
temp_dict[key2] = round(vincenty.vincenty_inverse(), 3)
# 辞書のキー（この場合X「１」は計算した距離）でソートする。
# Sort by dictionary key (where X "1" is the calculated distance).
temp_dict = sorted(temp_dict.items(), key=lambda x:x[1])
temp_dict = dict(temp_dict)
f = open("result.txt","a",encoding="utf-8")
for i,(key3, value3) in enumerate(temp_dict.items()):
if i < 10:
if value3 != 0.0:
f.write(f"{key}との距離 TOP{i+1} {key3} : {value3} m.\n")
else:
break
f.write("\n")
f.close()
print(f"Progress in the process. {key} : {round(time.time() - start2)} sec")
print(f"End. {round(time.time() - start1)} sec")
続いて同じコードをNimで書いていきます。
The same code is then written in Nim.
Nim Code [calc_distance.nim]
import json
import os,strutils,tables,times
import Tables
import math
import algorithm
let start = cputime()
type Vincenty = ref object
lat1 : float
lon1 : float
lat2 : float
lon2 : float
ellipsoid : int
ELLIPSOID_GRS80 : int
ELLIPSOID_WGS84 : int
GEODETIC_DATUM : Table[int, seq[float]]
ITERATION_LIMIT : int
proc vincenty_inverse(calculator:Vincenty): float =
# 差異が無ければ0.0を返す
# Returns 0.0 if there are no differences
if calculator.lat1 == calculator.lat2 and calculator.lon1 == calculator.lon2:
return 0.0
# 計算時に必要な長軸半径(a)と扁平率(ƒ)を定数から取得し、短軸半径(b)を算出する
# The major axis radius (a) and the flatness rate (ƒ) required for calculation are obtained from the constant, and the minor axis radius (b) is calculated.
# 楕円体が未指定の場合はGRS80の値を用いる
# If no ellipsoid is specified, the value of GRS 80 is used.
var a = calculator.GEODETIC_DATUM[calculator.ELLIPSOID_GRS80][0]
var ƒ = calculator.GEODETIC_DATUM[calculator.ELLIPSOID_GRS80][1]
var b = (1 - ƒ) * a
var φ1 = degToRad(calculator.lat1)
var φ2 = degToRad(calculator.lat2)
var λ1 = degToRad(calculator.lon1)
var λ2 = degToRad(calculator.lon2)
# 更成緯度(補助球上の緯度)
# Further Latitude (Latitude on auxiliary sphere)
var U1 = arctan((1 - ƒ) * tan(φ1))
var U2 = arctan((1 - ƒ) * tan(φ2))
var sinU1 = sin(U1)
var sinU2 = sin(U2)
var cosU1 = cos(U1)
var cosU2 = cos(U2)
# 2点間の経度差
# Longitude difference between two points
var L = λ2 - λ1
# λをLで初期化
# Initialize λ with L
var λ = L
var sinλ:float
var cosλ:float
var sinσ:float
var cosσ:float
var σ:float
var sinα:float
var cos2α:float
var cos2σm:float
var C:float
var λʹ:float
# 以下の計算をλが収束するまで反復する
# Repeat the following calculations until λ converges
# 地点によっては収束しないことがあり得るため、反復回数に上限を設ける
# There is an upper limit on the number of iterations because it may not converge at certain points.
for count in 0..calculator.ITERATION_LIMIT:
sinλ = sin(λ)
cosλ = cos(λ)
sinσ = sqrt((cosU2 * sinλ) ^ 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cosλ) ^ 2)
cosσ = sinU1 * sinU2 + cosU1 * cosU2 * cosλ
σ = arctan2(sinσ, cosσ)
sinα = cosU1 * cosU2 * sinλ / sinσ
cos2α = 1 - sinα ^ 2
cos2σm = cosσ - 2 * sinU1 * sinU2 / cos2α
C = ƒ / 16 * cos2α * (4 + ƒ * (4 - 3 * cos2α))
λʹ = λ
λ = L + (1 - C) * ƒ * sinα * (σ + C * sinσ * (cos2σm + C * cosσ * (-1 + 2 * cos2σm ^ 2)))
# 偏差が.000000000001以下ならbreak
# Break if deviation is .000000000001 or less
if abs(λ - λʹ) <= 1e-12:
break
# λが所望の精度まで収束したら以下の計算を行う
# Returns None if the calculation did not converge
var u2 = cos2α * (a ^ 2 - b ^ 2) / (b ^ 2)
var A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
var B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
var Δσ = B * sinσ * (cos2σm + B / 4 * (cosσ * (-1 + 2 * cos2σm ^ 2) - B / 6 * cos2σm * (-3 + 4 * sinσ ^ 2) * (-3 + 4 * cos2σm ^ 2)))
# 2点間の楕円体上の距離
# Distance on ellipsoid between two points
var s = b * A * (σ - Δσ)
# 計算した距離を返す
# return calculated distance
return s
var calculator : Vincenty = new Vincenty
# パラメーターの初期化
# parameter initialization
calculator.lat1 = 0.0
calculator.lon1 = 0.0
calculator.lat2 = 0.0
calculator.lon2 = 0.0
calculator.ellipsoid = 1
calculator.ELLIPSOID_GRS80 = 1
calculator.ELLIPSOID_WGS84 = 2
calculator.GEODETIC_DATUM[calculator.ELLIPSOID_GRS80] = @[6378137.0, 1 / 298.257222101]
calculator.GEODETIC_DATUM[calculator.ELLIPSOID_WGS84] = @[6378137.0, 1 / 298.257223563]
calculator.ITERATION_LIMIT = 1000
type TempType = tuple
distance: float
address: string
type CoordinateType = tuple
address: string
distance: float
let jsonObj = parseFile("tokyoCoordinates.json")
let start1 = cputime()
for key1, value1 in jsonObj:
var start2 = cputime()
var temp_dict: seq[CoordinateType]
calculator.lat1 = value1[0].getFloat()
calculator.lon1 = value1[1].getFloat()
for key2, value2 in jsonObj:
calculator.lat2 = value2[0].getFloat()
calculator.lon2 = value2[1].getFloat()
if key1 == key2:
discard
elif round(calculator.vincenty_inverse(), 3) == 0.0:
discard
else:
temp_dict.add((key2,round(calculator.vincenty_inverse(), 3)))
var ranking : seq[TempType]
for i,kv in temp_dict:
ranking.add((kv[1],kv[0]))
ranking.sort
var newkey = key1 & " との距離が近いランキング TOP "
var coordinate : CoordinateType
block addTop10:
var f : File = open("calc_distance.txt" ,fmAppend)
for i,kv in ranking:
if i < 10:
coordinate = (kv[1],kv[0])
f.writeLine(newkey, i+1, " : ", coordinate," m")
else:
defer: f.close()
break addTop10
echo key1, " . ", round(cputime() - start2, 3), " sec"
echo "End. ", round(cputime() - start1, 3), " sec"
Pythonコードとほぼ同じなので、説明を省きますが、 Nimの正攻法でクラスを作ろうとすると、上のようなCの構造体のような感じになります。 import macrosでちゃんとクラスっぽい記述ができますので、良かったら調べてください。 3日目のNimなので大分手こずりました。。。
It's almost the same as Python code, so I won't go into it. If you try to create a class using Nim's straightforward approach, it will look like the C structure shown above. I can write a class-like description with import macros, so please check if you like. It was Nim on the third day, so I had a lot of trouble.
結果がこちらです。 Nimが2倍速いですが、それでも150時間を超える計算です。
Here are the results. Nim is 2 times faster, but still more than 150 hours.
Nim側の文字化けは気にしないでください（笑） テキストファイルは通常に出力されています。
Please do not worry about the phenomenon that the letters on Nim side are not displayed well. The text file is normally output.
取り敢えず、PythonのクラスファイルをNumpyやTensorflowで書き換えて、実際の計算コードも、Numpyのnp.meshgrid()や、np.expand_dims(Array, axis=1) 、np.expand_dims(Array, axis=0)等工夫すればもっと高速に動きそうな気配です。 自分の環境だとメモリが少なく、（４GB）すぐメモリエラーになりましたけど…。 Colaboratoryの出番ですね！
First of all, by rewriting Python class files with Numpy and Tensorflow, and devising Numpy's np.meshgrid(), np.expand_dims(Array, axis=1), np.expand_dims(Array, axis=0), and so on, the actual calculation code seems to work faster. In my environment, there is not enough memory (4 GB) and I got a memory error right away. It's Collaboratory's turn!
しかし、何も考えずFor文を二回回すだけのコードでも、約2倍ほどNimのほうが速くなったのは当然なのでしょうか？ いや、書き方が悪い、全然遅い、もっと速くできる？ （NimやPythonの達人達ならもっと速く動かせるコードを書けるだろうなぁ…）
But is it natural that code that only needs to be turned twice without thinking about it is about twice as fast with Nim? No, the writing is bad, it's too late, can you do it faster? (If you're familiar with Nim or Python, you can probably write code that runs faster...)
参考になったサイト
Useful Sites
最後に色々お世話になったサイトをズラーッと貼り付けておきます。 知恵をくれた先人達に感謝です。
Lastly, I will put up a lot of sites that you helped me a lot. I thank my predecessors for their wisdom.
Python使いで『今後はデータビジネスの現場かも』って人、NimData/Arraymancerをさわってみておこう。
High performance tensor library in Nim
nim two key table with generics
【Nim】seq型とstring型のリファレンスの罠【追記:v0.19.0にて解決済み】
Nim に入門してみて予想通り最高のプログラミング体験ができている件 (1) ~ 変数・制御構文・プロシージャ編
[Python]緯度経度から2地点間の距離と方位角を計算する
長々とお付き合い頂きありがとうございます。
Thank you for reading this long blog.
See You Next Page!