JuliaのパッケージGeometricIntegrators.jlを書いて公開してみた

# JuliaのパッケージGeometricIntegrators.jlを書いて公開してみた
[Julia Advent Calendar 2016](http://qiita.com/advent-calendar/2016/julialang)の４日目として書いてみます。
著者は、ふだんはFortranを使って[物性のシミュレーション](http://loto.sourceforge.net/feram/)などをしている、Julia初心者です。
[Julia練習帳](https://gist.github.com/t-nissie/261896dfc388da246ce3)にメモを取りながら練習しています。
アドバイスなどをいただけないかと思い、[Julia Advent Calendar 2016](http://qiita.com/advent-calendar/2016/julialang)に挑戦してみます。

新型コロナウイルスが流行って外出が難しくなった2020年4月にGeometricIntegrators.jlをJulia 1.4に対応させたので、本文章も更新してみました。

## パッケージとは
Juliaは簡単かつ強力に機能を拡張できる[パッケージ・マネージャPkg](https://docs.julialang.org/en/v1/stdlib/Pkg/index.html)を持っています。例えば`using Pkg; Pkg.add("Winston")`などで`$HOME/.julia/`以下にそのパッケージがインストールされます。

## パッケージの作りかた
JuliaではGitHubに自分で書いたパッケージを置くことで簡単に公開することができます。パッケージの作りかたと公開の方法はどこかに書いてありますが、他人のパッケージのマネをするのが簡単です。

## Geometric integratorとは
[Geometric integrator](https://en.wikipedia.org/wiki/Geometric_integrator)とは常微分方程式の解法（数値積分法）の一つです。[調和振動子](https://ja.wikipedia.org/wiki/%E8%AA%BF%E5%92%8C%E6%8C%AF%E5%8B%95%E5%AD%90)や[振り子](https://ja.wikipedia.org/wiki/%E6%8C%AF%E3%82%8A%E5%AD%90)や[惑星の運動](https://ja.wikipedia.org/wiki/%E3%82%B1%E3%83%97%E3%83%A9%E3%83%BC%E3%81%AE%E6%B3%95%E5%89%87)のシミュレーションをするのに便利です。たいてい系のエネルギーが原理的によく保存されます。

## GeometricIntegrators.jl
### 紹介
[GeometricIntegrators.jl](https://github.com/t-nissie/GeometricIntegrators.jl)はJuliaでそのGeometric integrationをするモジュールを提供します。Julia-1.4で使えるはずです。ホームページは https://github.com/t-nissie/GeometricIntegrators.jl です。

### セットアップ方法
依存するパッケージなどはないので

    using Pkg; Pkg.add(PackageSpec(url="https://github.com/t-nissie/GeometricIntegrators.jl", rev="master"));

とするだけでインストールできます。

### 使用例
#### モジュールの中の例
モジュール`src/GeometricIntegrators.jl`の最後のほうに２つの独立な調和振動子の例が書いてあります。Julia-1.4と[Winston](https://github.com/nolta/Winston.jl)をお使いなら、

    $ cd src/
    $ julia GeometricIntegrators.jl

とすると`GeometricIntegrators.eps`というプロットが得られます。

#### １次元調和振動子
`test/harmonic_oscillator.jl`は１次元調和振動子の例です。`harmonic_oscillator.jl`を実行すると`harmonic_oscillator_??.eps`などのいくつかのプロットが得られます。

#### ２体問題
`test/twobody.jl`は２体問題の例です。`twobody.jl`を実行すると、図１のような`twobody_??.eps`のプロットが得られます。

![twobody](https://raw.githubusercontent.com/t-nissie/GeometricIntegrators.jl/master/docs/twobody.jpg "two body problem")

図１ ２体問題。質量M=2とm=1の質点の軌跡。

### テスト
GeometricIntegrators.jlパッケージをクローンして、テストして、消去するだけなら

    using Pkg; Pkg.add(PackageSpec(url="https://github.com/t-nissie/GeometricIntegrators.jl", rev="master")); Pkg.test("GeometricIntegrators"); Pkg.rm("GeometricIntegrators")
とすればOKです。

### 使用方法
座標と運動量の微分を`qdot`と`pdot`として用意してください。両方とも関数の配列です。ハミルトニアンを作って、その正準方程式を考えるのがよいでしょう。あとは初期値とともに例にならって使ってみてください。

## Travis CI で Continuous Integration
[travis-ci.org](https://travis-ci.org)にサインインして[.travis.yml](https://github.com/t-nissie/GeometricIntegrators.jl/blob/master/.travis.yml)というファイルを適宜用意しておけば`git-push`のたびにテストを実行してくれるように設定できます。Cron Jobsとして設定しておけば定期的かつ自動的に最新のJuliaでパッケージが動くかどうかチェックしてもらえます。さらに[coveralls.io](https://coveralls.io)にサインインして[.travis.yml](https://github.com/t-nissie/GeometricIntegrators.jl/blob/master/.travis.yml)の`after_success:`にデータを送るように設定しておけばcode coverageをチェックすることもできます。[travis-ci.org](https://travis-ci.org)と[coveralls.io](https://coveralls.io)のどちらもGitHubのアカウントでサインインできます。[.travis.yml](https://github.com/t-nissie/GeometricIntegrators.jl/blob/master/.travis.yml)などは次のドキュメントを参考に設定しました。

* [GitHubと連携できる継続的インテグレーションツール「Travis CI」入門](https://knowledge.sakura.ad.jp/3754/)
* http://docs.travis-ci.com/user/languages/julia/
* https://github.com/JuliaLang/Example.jl/blob/master/.travis.yml

## パッケージを作っていて困ったこと
### Julia-0.5での変更にハマりました

Julia-0.4.6:

    julia> VERSION
    v"0.4.6"
    
    julia> typeof(x -> x[1])
    Function

Julia-0.5:

    julia> VERSION
    v"0.5.0"
    
    julia> typeof(x -> x[1])
    ##1#2

これはJulia-0.5の新しい仕様みたいです→[Julia言語の0.5の変更点](http://bicycle1885.hatenablog.com/entry/2016/09/23/084524)

こんなかんじで困っています：

    julia> VERSION
    v"0.5.0"
    
    julia> function g(f::Array{Function,1})
               x = [1,2,3]
               map(f_each -> f_each(x), f)
           end
    g (generic function with 1 method)
    
    julia> g([x -> 1x[1], x -> 2x[2], x -> 3x[3]])
    3-element Array{Int64,1}:
      1
      4
      9
    
    julia> g([x -> 1x[1], x -> 2x[2]])
    2-element Array{Int64,1}:
      1
      4
    
    julia> g([x -> 1x[1]])
    ERROR: MethodError: no method matching g(::Array{##13#14,1})
    Closest candidates are:
      g(::Array{Function,1}) at REPL[1]:2

`function g(f::Array{Function,1})`と引数の型を指定しないで、ただ`function g(f)`と書けば動きました。

### Julia-1だかの変更にハマりました
`where` を使って書くことになったみたいです。

    function energy(q::Array{T,2},
                    p::Array{T,2}, m, k, n::Integer) where {T<:AbstractFloat}

あと、Project.tomlを用意することになったみたいです。[パッケージ・マネージャPkg](https://docs.julialang.org/en/v1/stdlib/Pkg/index.html)にリンクのある文章を参考にしてください。

## 他の常微分方程式を解くパッケージ
* [RungeKutta.jl](https://github.com/timothyrenner/RungeKutta.jl) （ただし、Runge-Kutta法はGeometricではない、もしくは、[Symplectic](https://en.wikipedia.org/wiki/Symplectic_integrator)ではないです）更新が止まってしまったみたいです。
* [Sundials.jl](https://github.com/JuliaDiffEq/Sundials.jl)
    * [#MonthOfJulia Day 21: Differential Equations](http://www.exegetic.biz/blog/2015/09/monthofjulia-day-21-differential-equations-2/)
* [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl)
    * [Introducing DifferentialEquations.jl](http://www.stochasticlifestyle.com/introducing-differentialequations-jl/)
* その他[JuliaDiffEq](https://github.com/JuliaDiffEq/)にあるもの