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で使えるか分からないので、こちらの人のブログ記事を参考にさせて頂きます。サンキューです!

vincenty 0.1.4

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「1」は計算した距離)でソートする。
    # 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.


Responsive image

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)等工夫すればもっと高速に動きそうな気配です。 自分の環境だとメモリが少なく、(4GB)すぐメモリエラーになりましたけど…。 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をさわってみておこう。

Amazon Cloud9では、gccとvimも使おう。

High performance tensor library in Nim

Nimで競プロをするための基本的な機能

Arraymancer

Tutorial: First steps

nim

Module nimdata

Nimで文字列から他の型への相互変換の仕方

Nim ファイルの読み書き

Nim 文字コードを指定したファイルの読み書き

csvtools

Nimで文字コード変換

【Nim】文字列操作

NimData

nimのJSON操作

Nim Tutorial (Part I)

NimのライブラリArraymancer を触ってみます.

Nimの制御構文(if, switch case)

For Loops & Iterators

regex

times

strutils

Nim for Python Programmers

Nimのオブジェクト指向の整理

Nim】クラス定義マクロを作ってみる

nim two key table with generics

tables

Macros

nimコレクション操作

【Nim】seq型とstring型のリファレンスの罠【追記:v0.19.0にて解決済み】

json

nim_syntax.md

Nim basics

Nim に入門してみて予想通り最高のプログラミング体験ができている件 (1) ~ 変数・制御構文・プロシージャ編

algorithm

Nimのデータ型

Files

vincenty 0.1.4

[Python]緯度経度から2地点間の距離と方位角を計算する

長々とお付き合い頂きありがとうございます。

Thank you for reading this long blog.




See You Next Page!