Julia: 東京から世界までの距離

東京から世界までの距離

定数値が何箇所にも出るのとか，plot!() を何回も書くのとかを書き直して，場所も追加変更してみた。

using GMT

# Tokyo
lon  = 139.76
lat  = 35.68
Gdist = gmt("grdmath -Rg -I1 $lon $lat SDIST")

coast(region=:global, land=:lightgreen, ocean=:blue, area=1000,
     frame=(grid=30, title="Distances from Tokyo to the World"),
     shore=:thinnest, proj=(name=:Hammer, center=90),
     figsize=20, portrait=false)
grdcontour!(Gdist, annot=(int=1000, labels=(curved=true, unit="\" km\"",
            font="white")), labels=(line="Z-/Z+",), smooth=8, cont=500,
            pen=((contour=true, pen="thinnest,white,-"),
            (annot=true, pen="thin,white")) )

# Location info for any other cities + label justification
cities = [30.50 50.45; 9.42 0.38; 0 -90; 282.95  -12.1; 151.21 -33.87; 237.67 47.58; 28.20 -25.75; 12.5 41.99];
just_names = ["LT KYIV", "RM LIBREVILLE", "RM SOUTH POLE", "LM LIMA", "LM SYDNEY", "RT SEATTLE", "LM PRETORIA", "RM ROME"];
D = text_record(cities, just_names)

# For each of the cities, plot great circle arc to Rome with gmt psxy
for i in 1:size(cities, 1)
    plot!([lon lat; cities[i, 1] cities[i, 2]], lw=:thickest, lc=:red)
end

# Plot red squares at cities and plot names:
plot!(cities, marker=:square, ms=0.2, fill=:red, markerline=:thinnest)
text!(D, offset=(away=true, shift=(0.15, 0)), font=(12, "Courier-Bold", :red),
      justify=true, no_clip=true)

# Place a yellow star at Tokyo
plot!([lon lat], marker=:star, ms=0.5, fill=:yellow, ml=:thin)

# Sample the distance grid at the cities and use the distance in km for labels
dist = grdtrack(Gdist, D);
text!(round.(dist), offset=(0, -0.5), noclip=true, fill=:white, pen=true, clearance=0.05,
      font=(12,"Helvetica-Bold"), justify=:CT, zvalues="%.0f", show=true)

【数学のひっかけ問題】三角形の面積を求めよ

【数学のひっかけ問題】三角形の面積を求めよ
https://nazesuugaku.com/fake_triangle/
で，「底辺が10cm、高さ6cmの直角三角形の面積を求めなさい。」

に対して，「そのような直角三角形は存在しない」とある。

確かに記事中のように底辺をとると「そのような直角三角形は存在しない」であるが，明らかに別の直角三角形もあるじゃないか？

短い問題文の中に，この三角形を排除する記述はないと思うのだが。

元々は「マイクロソフトの入社試験」だそうなので，原文を見ないとどうとも言えないが。

どうも，ご丁寧に図が描いてあったらしい。これじゃ「そんな三角形は存在しない」わなぁ。

真面目な解答もありました。（すごいなあ）

これ解けますか？（後編：マイクロソフトの三角形ありました）
https://blog.tech-monex.com/entry/2021/08/06/130000

 

簡単なデータ前処理で速度比較 tidyr, BaseR, AWK -- その3

sample.csv は 17035201 行の CSV ファイルである。最初の 7 行は以下のようになっている。

"id","label"
1,"aaa_あああ"
2,"aaa_いいい"
3,"bbb_ううう"
4,"bbb_えええ"
5,"ccc_おおお"
6,"ccc_かかか"

これを入力とし，以下のように変更を加え出力する。

"id",label"
1,"aaa"
1,"あああ"
2,"aaa"
2,"いいい"
3,"bbb"
3,"ううう"
4,"bbb"
4,"えええ"
5,"ccc"
5,"おおお"
6,"ccc"
6,"かかか"

# tidyr

あまり 'tidy' なプログラムでもないが，この 3 行のプログラムでの実行時間は 125.581, 122.979, 121.734 であった。
平均値は 123.431 秒。

しかし，システムが使用する時間が同程度あり，結局経過時間（計算開始から結果が出るまでの時間）は実に388.996, 283.905, 297.495 で，平均値は 323.465 秒，実に 5 分超である。

library(dplyr)
library(tidyr)

system.time({
    sample = read.csv("sample.csv")
    sample_row <- sample %>% separate_rows(label,sep="_")
    write.csv(sample_row, "tidyr.csv", row.names=FALSE)
})

# Base R

最近では tidyr は知っているが，Base R は知らない（書けない）という人も多いようであるが，この 4 行のプログラムでの実行時間は 25.035, 25.480, 25.446 であった。
平均値は 25.320 秒。

経過時間は 24.190 秒である。tidyr によるものより 30 倍速い。

system.time({
    sample = read.csv("sample.csv")
    a = unlist(strsplit(sample$label, "_"))
    sample = data.frame(id=rep(sample$id, each=2), label=a)
    write.csv(sample, "BaseR.csv", row.names=FALSE)
})

# AWK

AWK を知らない人も多いと思うが，この例のような簡単なデータ前処理には十分使える。
この CSV ファイルの構造をフルに利用して，実質 5 行のプログラムを書く。実行時間は 21.380, 22.024, 21.848 であった。
平均値は 21.751 秒。

BaseR によるよりも 16% 速い。tidyr に比べると 5 倍速い。

BEGIN {
    FS = ","
    getline
    print $0
}
{
    split($2, t, "_")
    printf "%s,%s\"\n%s,\"%s\n", $1, t[1], $1, t[2]
}

結果

またまた，AWK が最速という結果になった。

R はデータをすべてメモリ中に保存している。しかも，場合によってはそのコピーを作ったりすることもある。このプログラムが動いているとき R は 2GB ものメモリを消費している。

AWK はデータを「一行読んで，処理して，書き出して」を繰り返すのでメモリは殆ど使わない。

そのあたりがこういう結果になっているのであろう。

特に tidyr は，BaseR に比べてメモリ操作が下手くそ（余計なことをしている）なのであろう。

簡単なデータ前処理で速度比較 tidyr, BaseR, AWK -- その2

sample.csv は 17035201 行の CSV ファイルである。最初の 7 行は以下のようになっている。

"id","label"
1,"aaa_あああ"
2,"aaa_いいい"
3,"bbb_ううう"
4,"bbb_えええ"
5,"ccc_おおお"
6,"ccc_かかか"

これを入力とし，以下のように変更を加え出力する。

"id","new_label1","new_label2"
1,"aaa","あああ"
2,"aaa","いいい"
3,"bbb","ううう"
4,"bbb","えええ"
5,"ccc","おおお"
6,"ccc","かかか"

# tidyr

あまり 'tidy' なプログラムでもないが，この 3 行のプログラムでの実行時間は 111.195, 101.132, 100.276 であった。
平均値は 104.201 秒。

しかし，システムが使用する時間が同程度あり，結局経過時間（計算開始から結果が出るまでの時間）は実に751.059, 686.955, 827.347 で，平均値は 755.120 秒，実に 12 分超である。

library(dplyr)
library(tidyr)

system.time({
    sample = read.csv("sample.csv")
    sample_col <- sample %>% separate(label,into=c("new_label1","new_label2"),sep="_")
    write.csv(sample_col, "tidyr.csv", row.names=FALSE)
})

# Base R

最近では tidyr は知っているが，Base R は知らない（書けない）という人も多いようであるが，この 6 行のプログラムでの実行時間は 23.509, 23.268, 22.892 であった。平均値は 23.223 秒。

経過時間は 24.190 秒である。tidyr によるものより 30 倍速い。

system.time({
    sample = read.csv("sample.csv")
    a = unlist(strsplit(sample$label, "_"))
    sample$new_label1 = a[seq(1, 12, by=2)]
    sample$new_label2 = a[seq(2, 12, by=2)]
    sample = sample[c("id","new_label1","new_label2")]
    write.csv(sample, "BaseR.csv", row.names=FALSE)
})

# AWK

AWK を知らない人も多いと思うが，この例のような簡単なデータ前処理には十分使える。
この CSV ファイルの構造をフルに利用して，実質 5 行のプログラムを書く。実行時間は 20.328, 20.503, 20.126 であった。

平均値は 20.319 秒。

BaseR によるよりも二割方速い。tidyr に比べると 5 倍速い。

BEGIN {
    FS = ","
    getline
    print "\"id\",\"new_label1\",\"new_label2\""
}
{
    split($2, t, "_")
    printf "%s,%s\",\"%s\n", $1, t[1], t[2]
}

結果

AWK が最速という結果になった。

簡単なデータ前処理で速度比較 tidyr, BaseR, AWK

sample.csv は 17035201 行の CSV ファイルである。最初の 7 行は以下のようになっている。

"id","label","label2"
1,"aaa","あああ"
2,"aaa","いいい"
3,"bbb","ううう"
4,"bbb","えええ"
5,"ccc","おおお"
6,"ccc","かかか"

これを入力とし，以下のように変更を加え出力する。

"id","new_label"
1,"aaa_あああ"
2,"aaa_いいい"
3,"bbb_ううう"
4,"bbb_えええ"
5,"ccc_おおお"
6,"ccc_かかか"

# tidyr

あまり 'tidy' なプログラムでもないが，この 3 行のプログラムでの実行時間は 11.633, 11.909, 12.028 であった。
平均値は 11.857 秒。

library(dplyr)
library(tidyr)

system.time({
    sample = read.csv("sample.csv")
    sample_col <- sample %>% unite(new_label,label,label2,remove=TRUE,sep="_")
    write.csv(sample_col, "tidyr.csv", row.names=FALSE)
})

# Base R

最近では tidyr は知っているが，Base R は知らない（書けない）という人も多いようであるが，この 4 行のプログラムでの実行時間は 12.017, 11.987, 12.020 であった。
平均値は 12.008 秒。

system.time({
    sample = read.csv("sample.csv")
    sample$new_label = paste(sample$label, sample$label2, sep="_")
    sample = sample[c("id", "new_label")]
    write.csv(sample, "BaseR.csv", row.names=FALSE)
})

# AWK

AWK を知らない人も多いと思うが，この例のような簡単なデータ前処理には十分使える。
この CSV ファイルの構造をフルに利用して，実質 4 行のプログラムを書く。
これくらい短ければ，使い捨てのプログラムをすぐ書ける。
実行時間は 11.026, 11.087, 10.952 であった。
平均値は 11.0217 秒。

BEGIN {
    getline
    print "\"id\",\"new_label\""
}

{
    sub("\",\"", "_", $0)
    print $0
}

結果

AWK が最速という結果になった。

階乗競争

階乗の限界と計算速度を Python, Julia, R, Octave で比較してみた

追記: 2022/08/18  R での結果を更新。x86-64 版でやったら早々にクラッシュしていた。M1 チップ対応版で Julia と同等以上の速度である。

# Python 3.10.6

from math import factorial
from time import time
s = time()
a = factorial(100000)
print(time() - s)

    0.1885082721710205 sec.

# Julia 1.7.3

@time a = factorial(big(100000));

      0.011688 seconds (302 allocations: 3.076 MiB)

# R 4.2.1 (M1)

library(gmp)
> system.time({a = factorialZ(100000)})
             ユーザ             システム                 経過  
0.009999999999999787 0.000000000000000000 0.009999999999990905 

# Octave 6.2.0

tic
a = factorial(18); # 整数値解を表示するにはこれが限界
toc

    Elapsed time is 0.00777102 seconds.

ということで，Julia と R の圧勝...かな?

Julia も R も，gmp を使っているので，Python も gmpy を使えば同じ速度になるのかも知れない。が，現在の環境では pip install gmpy はできるのだが，いざ import gmpy すると，エラーになるので，確認できない。

ヒルベルト行列--2

ヒルベルト行列

ヒルベルト行列とその逆行列を求め，両者の行列乗算で単位行列になることを，4 種の言語でやってみた。

Octave と Julia はさすが。
R はちょっと誤差が大きい?
Python は R より誤差が大きい?

Octave

hilb
Return the Hilbert matrix of order N.

invhilb
Return the inverse of the Hilbert matrix of order N.

hilb(5)

    ans =
    
       1.0000   0.5000   0.3333   0.2500   0.2000
       0.5000   0.3333   0.2500   0.2000   0.1667
       0.3333   0.2500   0.2000   0.1667   0.1429
       0.2500   0.2000   0.1667   0.1429   0.1250
       0.2000   0.1667   0.1429   0.1250   0.1111

invhilb(5)

    ans =
    
           25     -300     1050    -1400      630
         -300     4800   -18900    26880   -12600
         1050   -18900    79380  -117600    56700
        -1400    26880  -117600   179200   -88200
          630   -12600    56700   -88200    44100

hilb(5) * invhilb(5)

    ans =
    
       1.0000        0        0        0        0
            0   1.0000        0        0        0
            0        0   1.0000  -0.0000        0
            0        0        0   1.0000        0
            0        0        0        0   1.0000
    

Julia

# import Pkg; Pkg.add("Nemo")
using Nemo
M = MatrixSpace(QQ, 5, 5)
h = hilbert(M)
    
    Welcome to Nemo version 0.32.2
    
    Nemo comes with absolutely no warranty whatsoever

    [   1   1//2   1//3   1//4   1//5]
    [1//2   1//3   1//4   1//5   1//6]
    [1//3   1//4   1//5   1//6   1//7]
    [1//4   1//5   1//6   1//7   1//8]
    [1//5   1//6   1//7   1//8   1//9]

inv(h)

    [   25     -300      1050     -1400      630]
    [ -300     4800    -18900     26880   -12600]
    [ 1050   -18900     79380   -117600    56700]
    [-1400    26880   -117600    179200   -88200]
    [  630   -12600     56700    -88200    44100]

h * inv(h)

    [1   0   0   0   0]
    [0   1   0   0   0]
    [0   0   1   0   0]
    [0   0   0   1   0]
    [0   0   0   0   1]

R

# install.packages("matrixcalc")
library(matrixcalc)

h = hilbert.matrix(5)
h

A matrix: 5 × 5 of type dbl
1.0000000 0.5000000 0.3333333 0.2500000 0.2000000
0.5000000 0.3333333 0.2500000 0.2000000 0.1666667
0.3333333 0.2500000 0.2000000 0.1666667 0.1428571
0.2500000 0.2000000 0.1666667 0.1428571 0.1250000
0.2000000 0.1666667 0.1428571 0.1250000 0.1111111

solve(h)

A matrix: 5 × 5 of type dbl
   25   -300    1050   -1400    630
 -300   4800  -18900   26880 -12600
 1050 -18900   79380 -117600  56700
-1400  26880 -117600  179200 -88200
  630 -12600   56700  -88200  44100

h %*% solve(h)

A matrix: 5 × 5 of type dbl
 1.000000e+00 -4.547474e-13 -1.818989e-12  0.000000e+00  0.000000e+00
-1.421085e-14  1.000000e+00  0.000000e+00 -5.456968e-12 -2.728484e-12
 0.000000e+00 -4.547474e-13  1.000000e+00  0.000000e+00 -9.094947e-13
 0.000000e+00  2.273737e-13  0.000000e+00  1.000000e+00 -9.094947e-13
 1.421085e-14 -2.273737e-13  0.000000e+00 -1.818989e-12  1.000000e+00

round(h %*% solve(h), 10)

A matrix: 5 × 5 of type dbl
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1

Python

from scipy.linalg import hilbert, invhilbert
hilbert(5)

    array([[1.        , 0.5       , 0.33333333, 0.25      , 0.2       ],
           [0.5       , 0.33333333, 0.25      , 0.2       , 0.16666667],
           [0.33333333, 0.25      , 0.2       , 0.16666667, 0.14285714],
           [0.25      , 0.2       , 0.16666667, 0.14285714, 0.125     ],
           [0.2       , 0.16666667, 0.14285714, 0.125     , 0.11111111]])

invhilbert(5)

    array([[ 2.500e+01, -3.000e+02,  1.050e+03, -1.400e+03,  6.300e+02],
           [-3.000e+02,  4.800e+03, -1.890e+04,  2.688e+04, -1.260e+04],
           [ 1.050e+03, -1.890e+04,  7.938e+04, -1.176e+05,  5.670e+04],
           [-1.400e+03,  2.688e+04, -1.176e+05,  1.792e+05, -8.820e+04],
           [ 6.300e+02, -1.260e+04,  5.670e+04, -8.820e+04,  4.410e+04]])

hilbert(5) @ invhilbert(5)

    array([[ 1.00000000e+00, -1.39888101e-13,  6.29496455e-13,  2.65876210e-12, -1.32938105e-12],
           [-2.00395256e-14,  1.00000000e+00, -2.34356978e-12,  2.63500333e-12, -1.31750166e-12],
           [ 2.34257058e-14, -1.27453603e-13,  1.00000000e+00, -2.93853830e-12,  5.59774449e-13],
           [ 0.00000000e+00, -2.27373675e-13,  9.09494702e-13,  1.00000000e+00,  0.00000000e+00],
           [-1.80966353e-14,  3.05089287e-13, -3.49720253e-13,  2.36299869e-12,  1.00000000e+00]])

import numpy as np
np.round(hilbert(5) @ invhilbert(5), 10)

    array([[ 1., -0.,  0.,  0., -0.],
           [-0.,  1., -0.,  0., -0.],
           [ 0., -0.,  1., -0.,  0.],
           [ 0., -0.,  0.,  1.,  0.],
           [-0.,  0., -0.,  0.,  1.]])

 

早期リターンではなくなっているぞ

以下の例を上げてきたが，これはもはや「早期リターンではない」

class Student:

  def __init__(self, student_id, name):
    self.student_id = student_id
    self.name = name

  def student_check(self):
    
    if self.student_id != "" or self.name != "":
      if len(self.student_id) > 12:
        flg = False
      else:
        flg = True
    else:
      flg = False
    return flg

student_exe = Student("101122222212","山田花子")
student_exe.student_check()

「if のネストを避ける」というなら，以下のようにすればよい。これで「早期リターン」だ！！

  def student_check(self):
    
    if self.student_id != "" or self.name != "":
      return len(self.student_id) <= 12
    else:
      return False

蛇足だが，flg なんて変数名はダメダメ。flag にしよう。

早期リターンとガード節

早期リターンとガード節

ガード節を使わないプログラムとして挙げられているのが prog1

if value === "hoge"
  result = "hoge"
else
  if value === "fuga"
    result = "fuga"
  else
    if value === "hogefuga"
      result = "hogefuga"
end

return result

ガード節を使うプログラムとして挙げられているのが prog2

if value === "hoge"
  return "hoge"
end

if value === "huga"
  return "fuga"
end

return "hogefuga"

この言語は何かわからないが，やっていることはわかる。

疑問は，「この言語は if-else しかないのか?if-elseif-else はないのか?」ということ

prog2 は，以下のように書くとわかりやすいのではないか? prog3

if value === "hoge"
  return "hoge"
else if value === "huga"
  return "fuga"
else
  return "hogefuga"
end

prog1 では，value が "hoge", "huga", "hogefuga" 以外の場合，未定義の result を返そうとする。言語によっては単に何も返さないと解釈されることもあるが，いずれにせよ，論理エラーである。

prog2, prog3  では，value が "hoge", "huga" 以外の場合，無条件に "hogefuga" を返す。これは論理エラーである。

prog3 はもう一段 else if があるべきだ。 prog4

これは prog1 を else if を使って，更に条件のチェック漏れを補完したプログラムだ。

if value === "hoge"
  return "hoge"
else if value === "huga"
  return "fuga"
else if value === "hogefuga"
  return "hogefuga"
else
  return "ERROR!!"
end

最後に，このバカバカしいプログラムに限ってだが，以下のように書くのがベストだ。prog5

if value === "hoge" || value === "huga" || value === "hogefuga"
  return value
else
  return "ERROR!!"
end

あるプログラミングスタイルを推奨する記事を書くとき，悪い例と良い例を提示するが，提示されたプログラムがクソだと，プログラミングスタイル自体がクソに見える。

早期リターンは必要なら使えばよい。

いずれにしても，プログラムをよく解析して，最良なアルゴリズムを適用することだ。

ヒルベルト行列

Julia で書くと以下のようになる。

一行で書けるが決してこの書き方をすすめるわけではない。普通に2重ループで書くほうがよい。

演算子に // を使ったのは，分数型で答えを出力するためである。

julia> h(n) = [1 // (i + j - 1) for i in 1:n, j in 1:n]

h (generic function with 1 method)

julia> h(5)

5×5 Matrix{Rational{Int64}}:

 1//1  1//2  1//3  1//4  1//5

 1//2  1//3  1//4  1//5  1//6

 1//3  1//4  1//5  1//6  1//7

 1//4  1//5  1//6  1//7  1//8

 1//5  1//6  1//7  1//8  1//9

実数型に変換すると。

julia> float.(h(5))

5×5 Matrix{Float64}:

 1.0       0.5       0.333333  0.25      0.2

 0.5       0.333333  0.25      0.2       0.166667

 0.333333  0.25      0.2       0.166667  0.142857

 0.25      0.2       0.166667  0.142857  0.125

 0.2       0.166667  0.142857  0.125     0.111111

Python だと

>>> a = [[1/(i+j-1) for i in range(1, n+1)] for j in range(1, n+1)]

のようになるが，range(1, n+1) が鬱陶しい。

>>> a = [[1/(i+j+1) for i in range(n)] for j in range(n)]

とすればよい。

しかし，n = 5000 のとき，Julia で // を / にして直接実数型の結果を返す場合に比べて，Python は，ほぼ 24 倍遅い。

追記：

Nemo パッケージに hilbert 関数があった。

 

julia> using Nemo

julia> M = MatrixSpace(QQ, 5, 5)

Matrix Space of 5 rows and 5 columns over Rational Field

 

julia> h = hilbert(M)

[   1   1//2   1//3   1//4   1//5]

[1//2   1//3   1//4   1//5   1//6]

[1//3   1//4   1//5   1//6   1//7]

[1//4   1//5   1//6   1//7   1//8]

[1//5   1//6   1//7   1//8   1//9]

Julia の Plots で日本語

もう，何回も記事にしてきたが，環境が変わるとできなくなったりする。

もし，日本語が使えなくなったら以下を試してみよう。

日本語を使う関数で fontfamily="Helvetica" とか書いてみる。

ただし，正確にそのフォントが使われる保証はない。フォント名をどのように変えても，「なんとかゴシック」が使われるようだ。

でもまあ，トーフが描かれるよりはまし。

julia> using Plots
julia> histogram(randn(1000), title="ヒストグラム", fontfamily="Helvetica")

関数ごとに fontfamily を指定するのが面倒ならば，最初に他のオプションも含めて指定しておく。これ一回で後々全て有効。

gr(     grid=false,
        tick_direction=:out,
        right_margin=15mm,
        fontfamily="IPAMincho",
        foreground_color_legend = nothing)

Julia の小ネタ--046 URL エンコーディングとデコーディング

URL エンコーディングとデコーディング

パーセントエンコーディング (英: percent-encoding)とは、URIにおいて使用できない文字を使う際に行われるエンコード（一種のエスケープ）の名称である。

using URIs

 

エンコーディング

t = escapeuri("エンコーディング")

"%E3%82%A8%E3%83%B3%E3%82%B3%E3%83%BC%E3%83%87%E3%82%A3%E3%83%B3%E3%82%B0"

 

デコーディング

unescapeuri(t)

"エンコーディング"