専門書の自炊

学生時代の蒐集癖で気象の専門書をたくさん持っているのですが、ごくたまにしか読まないため本棚に長らく飾られていました。いまとなっては絶版の本もあるので紙媒体に手を入れることに抵抗があったのですが、自らの血肉になることで初めて意味を持つので、意を決して自炊することにしました。iPadのようなタブレットで見ることができて快適です。早くやればよかったですね。転勤が2,3年くらいであるので、引っ越しの都度、重い専門書を一緒に運ばなくてもよくなるメリットもありますね。デメリットとしては、自炊したことを忘れてしまうことでしょうか

私の環境は、ScanSnap iX1500, HFS(R) ペーパーカッター (A4 ホワイト) 大型裁断機、です。裁断機はすごくパワフルで自炊に役立っていますが、如何せんでかいし重いので保管するのが難点です

Amazon.co.jp: HFS(R) ペーパーカッター (A4 ホワイト) 大型裁断機 400枚まで 事務 : 文房具・オフィス用品

ガベージコレクション

べつに全世界のユーザに見える形で投稿することもないのですが、日々アウトプットしていないと頭が強烈になまることが経験的にわかってきたので、とびとびになりそうな気がしますが、習慣にできるまではなんとか続けていけたらと思っています。あんまりプライベートなことは書くつもりはなく、基本的にはプログラミングにおけるハック的なことの紹介ができればそれでいいかなと思います。思い通りのコードがかけると楽しいと思っている程度の人間です。気が向いたら自分が作ったツールとかをまとめて紹介できればいいかなと思っています (大したものを作っているわけではありませんが)。

プログラミング言語には、garbage collection という動的に確保したメモリのうち不要なものを解放して再利用する仕組みがあります。いつまでも確保していては資源がいくらあっても足りないわけですね。自分も、いらなくなったものは定期的に処分しなければ、と思いつつ、処分の決断ができずに溜め込んでしまうのでした。

キャッシュ

中間ファイルをかますことでデータ読み込みの時間を短縮できることがあります。NuSDaS形式だと、ある水平面の要素のデータを取り出すnusdas_read関数で読み込むことが基本で、メンバーごとや高度面ごとのデータを読もうとすると、その分都度ループをまわしてデータを構築する必要があり、高解像度のデータになるとデータ読み込み処理がネックとなってきます。解析スクリプトの初期作成段階では割と試行錯誤がつきもので、何もいきなり本番相当のデータでやらなくてもいいのですが、実際のデータで読み込んでどう見えるかはやはり気になるものです。なので、このデータ読み込み処理がなんとか短縮できないかなと思うわけです。PythonにはNumpy配列をバイナリに保存できる機能があるのでそれを使うと、

import os
import numpy as np

isw_use_cache = True
datfile = "data.npz"

if isw_use_cache and os.path.exists(datfile):
    with np.load(datfile) as data:
        dat = data["rain"]
else:
    dat = ....
    np.savez(datfile, rain=dat)

のようになり、バイナリデータを1回でごそっと読むので読み込み処理を短縮できます。実際、手許の、データ読み込み+matplotlibによるお絵かきする解析スクリプトでは、何もしないと4分弱かかりますが、一度保存しておいた中間キャッシュデータを読むことで40秒程度で済みます。何回も試行錯誤するなら待ち時間の短縮はありがたいですよね。ただ、中間データを作っている分は当然容量食いますのでストレージに余裕があることが必要です。実行時間を取るかストレージを取るか、まさにトレードオフの関係ですね。np.savez_compressed() という圧縮保存方法もあるようなので、大容量データの場合はこちらの方がよいかもしれません。

敵がいなけりゃ

 自分が過去に書いた記事 (このブログに限らない) を見返していると、昔の自分は何て良いことを書いているんだ!(それに引き換え今の自分と来たら過去に書いたことなど全く覚えておらず、過去の教えも忘れて我流に走り、成長も維持もできずむしろ退化する始末、という一瞬の自省をするものの、またすぐに怠惰なサイクルに戻ってしまうのであった)。観測データが取り込まれず、数値予報でいうところの予報-予報サイクルに陥っているかのようである。連鎖を断ち切ってリセットするか、現実を取り込んで少しずつ取り戻すかということになるだろうか。たまには外の空気を吸うのがよいと言われる所以はこういうことなのだろう
 内向型は自分の意見を発信するのがよいとのアドバイスが巷にあるらしいので、できるかどうかはさておき、発信頻度を少し増やせたらいいなとぼんやり思っています (デッドストックはまああるんだけど)。何考えてるか分かんないって面と向かって言われることないけど、たぶんそう思ってるよねという

少しのことにも先達はあらまほしきことなり

パワーポイントで挿入するテキストボックス、フォントをMS Pゴシックから変更できないだろうかとずっと思っていながらも特に調べてこなかった自分がいたが、ようやく調べる気になった。Web検索するとトップに出てくる。変更したテキストボックスを右クリックして出てくるメニューから「既定のテキストボックスとスタイルとして設定」を選択すればOK。なんだ簡単にできるじゃないか。

検索すれば全て解決するわけではないが、検索しなければ解決策が見つかる類のものでもない。

パワーポイントついでで言うと数式を挿入するショートカットキーってないのかなと常々思っていた。これも検索してみよう。手許で動作確認した結果、どうやら Alt+「;」で行けるようだ。ほんとうはどこでショートカットキーを割り当てているかまで突き止めておくと安心できるのだが、そこまでは分からず。

matplotlibでの1次元プロットで欠損域にあたる背景を塗りつぶす

Copilotに教えてもらいながら作業。欠損域を塗りつぶす。最初に描画範囲を灰色で塗りつぶしておく。mask配列がTrueのところが欠損扱いになっているので、(False,False)区間だけ白で塗り直すという処理になる。サンプルスクリプトと描画例を示す。実行のたびに結果が変わるのを防ぐために、乱数のseedは固定しておく。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

def main():
    # random seed
    np.random.seed(5489)

    xx = np.linspace(0,1,101)
    yy1 = np.random.randn(xx.shape[0])
    yy2 = np.ma.masked_where(yy1 <= 0, yy1)

    fig = plt.figure(figsize=(8,4))
    ax1 = fig.add_subplot(1,1,1)
    ax1.plot(xx,yy1,linestyle=":",linewidth=1)
    ax1.plot(xx,yy2)

    # mask
    mask = yy2.mask
    ax1.axvspan(xx[0], xx[-1], color="lightgray", alpha=0.5)
    for i in range(len(xx) - 1):
        if (not mask[i]) and (not mask[i+1]):
            ax1.axvspan(xx[i], xx[i+1], color="white", alpha=1.0)

    plt.savefig("test11.png")
    plt.close()

if __name__ == "__main__":
    main()

共通モジュール作ればいいんだろう

毎回毎回似たようなコードをなんか書いているような気がして、雛形らしきものを作ってそれをコピペできるようにしているけれども、汎用性を高める意味ではメソッドとかモジュールにしたほうがいいんだろうなと思いつつ、やり方を調べる気になれなくて進んでいない。やってみると案外すんなりできてしまう予感はする。自分でハードルを上げてしまっている典型かも。

ログはちゃんとみましょう

とあるプログラムでインライン展開のオプションを付けると実行速度が遅くなってしまい、てっきり、今回の改修で用意したサブルーチン内の処理で最適化がうまくいっていないのかと思いましたが、chatGPTに訊いていろいろオプションを付けたり試す中でログをちゃんと見た方がいいとふと思い、インライン展開のオプション有無で差分がどこに出ているのかを見てみると、インライン展開が別のサブルーチンで起きていると分かりました。そのサブルーチン内のループにOMP指示行を付け足したところ高速化し、インライン展開のオプション有無にかかわらず速度が同等になりました。ログ確認大事。

杉下右京

コードの改良を加えたとき、既存の実行結果を変えるものではないことがもし期待されるのであれば、出力結果が一致することを確かめる必要があります。結果が変わってしまう場合、バグが混入している可能性があるからです。「あるはずなのにないもの」「ないはずなのにあるもの」があった場合、それは何か不可解なことが起こっているということです。違和感に気づく感性を身に付けたいものです。

お絵描きにはまずサンプルデータで

いきなり実データを使ってお絵描きを試みるのもよいですが、お絵描きする際に渡すデータのフォーマットを決めておいて、それに合わせた適当なサンプルデータでお絵描きしてみるとよいと思います。サンプルデータは np.random で作成すればよいでしょう。サンプルデータでお絵描きの細部にわたる調整を実施し、完成したところで実データのお絵描きをすれば効率的かと思います。インタフェースを整えれば、メソッドとして使いまわすこともできるでしょう。

MSM25XX地形

MSM25XXの改良では標高データの作成元がMERIT DEMに変わります*1。ではどれくらい標高が変わっているのでしょうか。可視化して調べてみましょう。地形データは気象業務支援センターのWebページに置いてありますので、そこから取得します。TOPO.MSM_5K というバイナリファイルがありますので、これを読み込みます。比較するデータは2023年3月の改良時に更新した地形データです。側面境界値に使うGSMにおいて地形データの作成元をMERIT DEMにしたことから、ほんのごくわずかですが、緩和領域での標高に変化がありました*2

MSM25XXの地形 (対MSM2303差分)

ループを抜けるのは exit?cycle?

とあるプログラムのデバッグをしていたときのこと。浮動小数点例外が起きてエラーで止まってしまいました。当初はゼロ割を疑っていたのですが、printデバッグしてみたところ、どうも違うらしく、じゃあ原因はなんだろうといろいろ試していたところ、コードに仕込んだprint文が標準出力されないことに気づき、バグの原因が判明しました。結論から言うと、cycle x_loop と書くべきところが exit x_loop になっており、まるっとループを抜ける仕様になっていました。違いを以下のサンプルコードで見てみましょう。

program main
  implicit none
  integer(4), parameter :: ixmax = 5
  integer(4), parameter :: ix_exit = 2
  real(4), parameter :: undef_s = huge(0.0e0)
  integer(4) :: ix1, ix2
  real(4) :: rx1(ixmax), rx2(ixmax)

  print*, "start x_loop1"
  x_loop1: do ix1 = 1, ixmax
     print*, 'ix1 = ', ix1
     rx1(ix1) = real(ix1)
     if ( ix1 == ix_exit ) then
        rx1(ix1) = undef_s
        print*, "exit x_loop1"
        exit x_loop1
     end if
  end do x_loop1

  print*, "start x_loop2"
  x_loop2: do ix2 = 1, ixmax
     print*, 'ix2 = ', ix2
     rx2(ix2) = real(ix2)
     if ( ix2 == ix_exit ) then
        rx2(ix2) = undef_s
        print*, "cycle x_loop2"
        cycle x_loop2
     end if
  end do x_loop2

  print*, "rx1 = ", rx1
  print*, "rx2 = ", rx2

end program main

これを実行すると私の環境では以下のようになりました。意図しているのは後者のほうですが、前者では途中でループを抜けてしまっているので、実数配列の後半が不定になってしまっており、特に4番目の要素には意図しない値が混入してしまっています。

$ ./a.out
 start x_loop1
 ix1 =            1
 ix1 =            2
 exit x_loop1
 start x_loop2
 ix2 =            1
 ix2 =            2
 cycle x_loop2
 ix2 =            3
 ix2 =            4
 ix2 =            5
 rx1 =    1.00000000       3.40282347E+38   0.00000000      -2.35710200E+21   0.00000000
 rx2 =    1.00000000       3.40282347E+38   3.00000000       4.00000000       5.00000000

exit なのか cycle なのかは、どういうプログラムにするかに依るので、よくよく考えて使い分けましょう。

Ruby+gnuplot

学生時代にはRubyDCL (GPhys/GGraph) を使った解析をしていたことから Ruby の方が慣れているのですが、最近の可視化は Python (matplotlib) による方法が主流で、Python での解析はやや不慣れです。できれば Ruby だけで完結させたいところなのですが、可視化のライブラリがないのがネックです (ないことはないはずだが)。 Ruby + gnuplot という方法が簡単なようですが、自分は gnuplot は不慣れです。gnuplot は学生時代に最初に教わった可視化の方法なんですが、どうにも身に付かずに終わってしまいました。gnuplot が使えるようになるとクイックルックに非常に便利なんだろうなと思っています。

Win11 Proにアップグレードしてリモートデスクトップ接続をする

昨年の夏にPCを新調したのですが、古いWin11はHomeからProにアップグレードして使っています。Win11 Proにするとリモートデスクトップ接続ができるようになり、ホストPCでVPN接続してやれば、リモートPCからあたかもVPN接続しているかのように扱えます。 なお、Win11 Homeのままリモートデスクトップ接続ができるというRDPWrapを試していた時期もありましたが、Windows Updateが走るとなんか動作しなくなるし、そもそもレジストリを破壊してしまうのでオススメしません。実際、もとに戻そうとしたら戻らなくて、やむなくWindowsクリーンインストールをする羽目になりました。