趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

numba.jitの衝撃:5分 ⇒ 3秒

先日お話させて頂きましたnumba.jitの件です.私の手続きが間違ってました. 前回の河床変動の計算時間ですが,タイトルのとおり,超高速化します.100倍です. github更新しておきました.

github.com

numba.jitの使い方

でも,実問題はそんなわけにはいきません.

以下のサイトを参考にしました.

http://yutori-datascience.hatenablog.com/entry/2014/12/10/123157

http://www7a.biglobe.ne.jp/~thor/novel/column/14.html

ポイントは,公式 http://numba.pydata.org/ にもあるように,

@numba.jit(nopython=True)
or
@numba.njit

と書くことです.nopythonオプションは,最適化できないものはエラーとして返すというものです.これがないと最適化できるものはするけど,できないものはしないのでほんとに最適化できているかわからないということです.

特に,後述の高速化できないものがループに含まれるとそのループ自体が最適化されないので,ほとんど早くなりません.

一応並列化もやってみましたが,一次元なので,逆に遅くなりましたが,大規模だと早くなるようです.

@numba.jit(nopython=True, parallel=True)
# 並列化loopをrange(imax)から,
numba.prange(imax)
# に変更

まだ,やってないですがclassの最適化もあるようです.ちょっと面倒です.

https://stackoverflow.com/questions/38682260/how-to-nest-numba-jitclass

numba.jitで使えないもの

文字列型

一切駄目です.例えば,print(str(i))も駄目でした.

lambda式

駄目です.ちょっと不便

pandas

これも駄目です.

dictionary型

駄目です.

通常のリスト型およびその処理

すべて,numpyにしましょう.例えば,通常sumも使えません.np.sumならO.K.です

感想

  • もとのpythonのコードをほとんど手を加えずに,超高速化出来るので便利です.野望が広がります.
  • pythonらしいコードは書けないので,必要なdefだけ最適化すれば良いと思います.それでも,リスト内包表記や配列の演算などは普通に使えると十分です.
  • なので,静的な言語の知識は若干必要かも.我々のように高級言語から入った人間にはそんなに苦労する話じゃないですが.
  • numbaのエラー文がめちゃくちゃわかりにくい.エラー箇所が特定できないので,ちょっとコーディングが大変.

rigid-lid近似

連投になります.

直交格子でLES + 河床変動計算をやりたいと考えており, iricにあるNabiさんのLESモデルについて調べていると,河床変動は出来ないのですが,関連する論文を調べているとすでに河床変動計算も完成しているようです.

この辺りです.

https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2012WR011911 https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1002/wrcr.20303 https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1002/wrcr.20457

水面の処理の話が殆ど書いていないので,よくよく読んでいるとrigid-lid近似としていると書いてありました.

これって,海の計算でよく使う「海面の鉛直変位は認めない」という条件です.要するに,水位は一定ということです.

河床変動って,河床面の摩擦力なので,垂直抗力(重力)が重要なので,このような近似を設けても良いのかなと疑問に思います.

一次元混合粒径掃流砂河床変動計算の更新

いろいろ修正しました.

だいぶ計算も安定してきました.

今度はgithubに上げました.

ちょっとおしゃれな図化にしました.

nbviewerのリンクは以下です.

http://nbviewer.jupyter.org/github/computational-sediment-hyd/1DBedVariationMixed/blob/master/make_graph_case01.ipynb

計算時間は5分くらいだと思います.

github.com

札幌市清田区里塚の航空写真,古地図,地質図を並べてみた

2018/9/24追記

清田区里塚で液状化した場所と住所は?被災前後で画像を比較 | わしろぐtwitter画像の位置の住所がまとめられていたのでポップアップを作ってみました.

Leaflet Sync Demo


またまたLeaflet.syncです.ブログへの埋め込み方がわかりました.

今回の北海道の地震液状化が生じた清田区里塚の地形を分析してました.

北海道札幌市清田区里塚の地震の影響の被害画像がヤバい!液状化で里塚ヤバすぎ,清田区おわった便りTwitterから届く,,, | What a story!in my Book

こう見ると谷底で人が住むところじゃないですね.

大きい版はこちら

github.com

Leaflet.Syncで真備町の写真を並べてみました

地理院真備町の被災直後の写真を公開していたいので,Leaflet.Syncで被災前,被災後,地形分類図を並べてみました.

いつもどおり,ブログに埋め込みたかったのですが,駄目だったので,githubに上げました.

Leaflet Sync Demo - with three maps listening

github.com

  • 被災箇所には旧川があったようです.
  • 破堤箇所は,特別水位が上がる箇所というよりも,越水で危険な状態にたまたまその箇所が切れた感じです.支川合流部で若干の弱部であったかもしれないですが,どこが切れてもおかしくない感じです.

小田川は,合流点付替えという大事業が今年度より行われる予定でした.確かに,高梁川の背水影響を受けづらくなりますが,抜本的な対策とは言えません.

今回の災害から,高い堤防は破堤時の浸水深が大きくなり,人的被害のリスクが高まることが理解できます.

あまり誰も言わないですが,築堤によるリスクは大きいな思いました.

私の改修案は,引き堤および旧川の復元です.これしかサスティナブルな開発はないと思います.


CDN(Content delivery network)を初めて使いました

私が知識不足なだけですが,CDNを初めて知りました. かなり便利ですね.

本番環境では禁止されているようですが,テスト環境だとサーバを立てずにnpmを使えるは非常に役立ちます.

2種類使ったので書いておきます.

UNPKG

npmのCDN化用です.jsやcssを使うためにコードに埋め込んでます.

RawGit

githubCDN化するのに使ってます.今回はhtmlですが,textやjsonCDN化に使えそうです.

Re:一次元河床変動計算

混合粒径モデル

ソースありがとうございます。私の方でもソースを読んでみて、計算回してみたいと思います。

個人的にはもう少し細かく粒径を切って計算を回して、粒度分布の形状をみたいかなと。

粒度分布の検討

中身を確認した上で、 上記のソースで表層の粒径と流砂量の分布を見てみたいと思います。

正規分布(粒径を対数変換)による粒度分布の推定(2個以上の混合正規分布の推定)のためのEMアルゴリズムは書けたので、 条件を決めていくつか試してみる予定です。言語はRなので、Pythonに変えるつもりです。

現在の河床変動解析のフレームワークで綺麗な対数正規分布になるようになるという結論なら、 それはそれで面白いかなと思ってます。

水位予測

記事について吟味させてください。 個人的には難しいところもありますが、甘えもある気がしています。

あと、計画や設計等々の絡みもあるので、一貫した態度を持たないといろいろ不整合が出てしまうのがなぁ、と。

一般の方々にも予測の難しさを知ってもらえるよう、あえていろいろディスクローズするのも 一つの手かなと思います。

国交省のとんでも水位予測について

あまりコラム的なことを書くのは趣旨が異なるので嫌なのですが思うことが多くあるのでまとめておきます.

水位予測の現状

去年の産経の記事がわかりやすいです.

https://www.sankei.com/premium/news/170711/prm1707110004-n1.html

要するに, - 精度が悪すぎて公開はできない. - 気象庁的には何とかしたいから指定区間は独自公開.

気象庁からみたら国交省の予測はありえない.基本がなってないので鼻で笑うレベルらしいです.

河川の予測は国交省しか出来ない的な印象ですが,測量データを公開しないので計算できないだけ.ほんとに誰かなんとかしてほしいです.

水位予測システムの改修

新システムに関する記事が以下です.

http://www.mlit.go.jp/common/001240498.pdf

https://www.sankei.com/region/news/180629/rgn1806290016-n1.html

水害リスクラインって...英語になってない...

モチベーションは正しいのですが,このモデルが相当まずい.

新たな水位予測システムの技術的なまず

このモデルは,京大 立川先生が提案したモデルです.

https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/194041/1/jscejhe.67.I_511.pdf

既存のモデルとの変更点は,「データ同化」と呼ばれる手法を用いて,観測値(水位)に基づき,計算条件を逐次修正(最適化)することにより予測精度を上げる点です. これ自体は問題ないのですが,問題は最適化している条件として,「粗度係数」と「流入量係数」としていることです.後者はまだ許される範囲(流入量係数じゃなくて流入量にしないと駄目)ですが,前者はありえない. そもそも,リアルタイム解析でパラメータを最適化するなんてありえない.こんなことしていたら,アポロは一生月に着きません. 気象予測をはじめとしたデータ同化を用いた数値予報システムはいくつもありますが,すべて状態量を同化しています.

なぜか国交省は,この論文一本を担保に今年度より2-3年,数十億かけて,この考え方によるシステムの入れ替えを行います.税金の無駄遣いもいいとこです.精度は期待できないでしょう.

急に立川さんがスタンスを変更

ところが,唯一の担保だった立川さんが今月号の土論で急にスタンスを変更してきました.

https://www.jstage.jst.go.jp/article/jscejhe/74/2/74_32/_pdf/-char/ja

最新なので有料です.(DBに)

結論は,パラメータの同化では予測精度がでないので,状態量(貯留関数なので貯留高)の同化しましょうということです.

国交省はもう走り出したのにどうするのでしょうか..

立川さん含めこの手の研究で気になること.

  • なぜか,水文関連では粒子フィルタが好まれているようですが,ほとんどの論文で誤差分布が示されていないので,非ガウシアン分布になっているのかが不明です.ざっくり見た感じでガウシアンっぽいので,メリットがないのでは,

  • 立川さん程度の計算ならバックプロパゲーションで解けるのでは?

今後の期待

  • 河川は気象と比較するとまだまだ観測データが少ないので同化が難しい.今後は観測データが増えるようなので精度向上が期待できる.
  • 降雨流出過程が概念モデルにならざる得ないので同化が難しい.ここはブレイクスルーが必要

愚痴っぽくなりましたが,一技術者として思うことです.