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

Fortran Python Ruby bash

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

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

Fortran

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

杉下右京

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

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

Python

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

MSM25XX地形

Python 気象学

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

MSM25XXの地形 (対MSM2303差分)

*1:https://www.data.jma.go.jp/suishin/jyouhou/pdf/648.pdf

*2:https://www.data.jma.go.jp/suishin/jyouhou/pdf/606.pdf

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

Fortran

とあるプログラムのデバッグをしていたときのこと。浮動小数点例外が起きてエラーで止まってしまいました。当初はゼロ割を疑っていたのですが、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

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のクリーンインストールをする羽目になりました。

本の買取参考価格をスクレイピングで取得する

手持ちの本を売りにだそうとして、買取参考価格を調べられることに気づきました。ISBNコードを入れてWebページの結果を見てそれをエクセルに貼るということを手動でやっているうちに、だんだん作業が面倒くさくなってきて、これ自動でできないか、と考えました。ISBNコードから該当の本は一意に決まりますので、検索して得られるWebページをparseしてタイトルと買取参考価格をdivタグから取得できればよいわけです。実際、実装してみると案外すんなりできました。自分はRubyに割と慣れているのでRubyでやりましたが、最近の時流だとPythonでやるひとが多いのかもしれません。買取参考価格は変動するものなので、スクリプト実行するだけで更新してくれるという再現性持たせられたのも良い点と思っています。 ちなみに、参考価格を見ていると、高値になっている本は、古本市場に出回りにくい、所謂、良書なんだなと気づかされるわけですね。

値が不定のFortranプログラム

次のようなFortranプログラムを考えます。namelistから値を読み込むような仕組みです。ところが用意したnamelistに値が格納されておらず空白となっていました。さて実行するとどうなるでしょうか。

program main
  implicit none

  integer(4) :: isw_fss

  namelist/namprm/ isw_fss

  open(10,file="namelist.txt")
  read(10,namprm)
  close(10)

  write(6,namprm)
  print*, "isw_fss = ", isw_fss
end program main

$ cat namelist.txt
&NAMPRM
isw_fss=,
&END

$ gfortran test2.f90
$ ./a.out
&NAMPRM
 ISW_FSS=1206747284 ,
 /
 isw_fss =   1206747284

自分の環境 (WSL2 Debian12) では、実行のたびにでたらめな値が入ってしまいました。怖い挙動ですね。ネームリストのデフォルト値を内部で設定しておけば事故は防げます。

Nautilusからsmb接続ができなくなって困った話

職場で使っているLinux端末（Ubuntu22.04）で、あるときを境にNautilusからネットワーク共有（smb接続）ができなくなって困っていました。症状としてはsmb:// でアドレスを入れても赤字表示になり「接続」ボタンはグレーアウトしてしまい接続できないというものです。ちょうど別の不具合を解消するためにいろんなパッケージを入れたり削除したりしてたので、よくわからないうちに必要なパッケージも消えてしまったのかもしれないと想像します。 調べても原因が分からなかったのでchatGPTに訊いてみたところ、smbclientが入っているかどうかと回答されたので、それを足掛かりに改めてWeb検索してみると以下ページにたどり着きました。自分とまったく同じ症状だったので、これだと思いました。パッケージがインストールされているか確認してみたところ、smbclientもgvfs-backendsもインストールされていませんでした。前者をインストールしただけでは改善しなかったので後者もインストールするとsmb接続できるようになりました。めでたし、 forums.ubuntulinux.jp ちなみに sudo mount する方法で一時凌ぎした時期もありましたが、トンデモ事故を起こしてしまったので、二度とやりません。

e2pdf

プログラムのソースコードをPDF化するのに、e2psをちょっと改造してe2pdfとして自分は使っています。コードは以下のとおりです。実行権限をつけてPATHを通せばよいです。文字コードをEUC-JPとしないとe2psが動かないのでnkfの処理が間に入っています。

#!/bin/bash
set -ex

FILENAME=${1:?}

BASENAME=${FILENAME%.*}
EXTENTION=${FILENAME##*.}

OPT="-size 11"
nkf -e ${FILENAME:?} > ${BASENAME:?}.tmp
e2ps ${OPT:?} -p -head "${FILENAME:?}" ${BASENAME:?}.tmp > ${FILENAME:?}.ps
ps2pdf14 ${FILENAME:?}.ps

rm -f ${BASENAME:?}.tmp ${FILENAME:?}.ps
exit 0

概ねこれで満足できてはいるのですけれども、文末に文字化けと思われる記号が現われることがあって、原因が分からないでいます。PostScriptには何も出ていないのでps2pdf14で処理するときに何か起きているくらいの想像しかできていません。

Rubyで経過時間を任意の時刻形式にする

以下のように、ログに時刻の情報が記載されていたとしよう

2025/02/27 Thu 00:55:54       2025/02/27 Thu 01:25:37

この情報から所要時間を %H:%M:%S の形式で取り出すことを考える。なお、この時刻表記は date コマンドで

date +"%Y/%m/%d %a %H:%M:%S"

として得られる時刻情報である。Rubyで処理することを考える。うまく処理して、開始時刻と終了時刻を文字列に格納することができたとする。

st = "2025/02/27 Thu 00:55:54"
en = "2025/02/27 Thu 01:25:37"

Timeクラスのparseメソッドを使うと、Timeクラスのオブジェクトに変換してくれる。

require 'time'
pst = Time.parse(st)
#=> 2025-02-27 00:55:54 +0900

Timeクラスのオブジェクトどうしの差分をとると、秒数が返ってくる。

pst = Time.parse(st)
pen = Time.parse(en)
p pen-pst #=> 1783.0

この秒数から欲しい形式 %H:%M:%S に整形することをやってもよいが、うまい方法があり、以下のようにすれば実現できる。秒数を 1970/1/1 からの経過秒とみなし、Timeオブジェクトに対して Time#strftime で整形するというわけ。

Time.at(pen-pst).utc.strftime("%H:%M:%S")
#=> "00:29:43"