合計特殊出生率（３）

合計特殊出生率（２）からの続き

物理学者の白石さんが，合計特殊出生率の件で強い主張をしていたので確認してみる。ようは，合計特殊出生率を決めているのは有配偶率であって，人口学的に確かめられている。それ以外の恣意的な原因を持ち出しているのはみんな出鱈目だ。ということだった。

仮にそうであっても，有配偶率自身を規定する原因は考えられるのだから，そこまで強く主張する必要はまったくないのに，と思う。まあいつものことだ。

都道府県別の合計特殊出生率と直近の国勢調査で得られた都道府県別の15-49歳女性の有配偶率を散布図にしてみる。1953年から2022年までの両者の関係も出せないかと調べるのだけれど，必要な情報はデータの海に埋もれてしまっていて，簡単には見つかりそうになかった。どの生成AIも簡単には探しだしてはくれない。

とりあえず，Juliaでコーディングしてみる。

using Gadfly

using Compose

using DataFrames

X = [46.6,47.4,49.4,47.7,49.9,51.2,50.0,50.4,50.9,49.9,49.3,49.0,44.2,49.0,50.4,52.0,50.4,52.3,49.5,51.7,51.4,51.6,51.9,51.4,52.2,44.4,43.5,48.8,47.0,48.5,49.2,51.9,48.9,50.7,49.9,49.3,51.1,49.6,46.1,46.2,48.5,48.6,48.5,49.5,49.4,49.2,45.8]

Y = [1.06,1.23,1.16,1.07,1.10,1.22,1.21,1.22,1.19,1.25,1.14,1.14,0.99,1.13,1.23,1.35,1.34,1.46,1.32,1.34,1.31,1.25,1.29,1.29,1.38,1.11,1.19,1.29,1.21,1.33,1.44,1.46,1.32,1.33,1.40,1.36,1.40,1.31,1.30,1.26,1.46,1.49,1.47,1.39,1.49,1.48,1.60]

Labels = ["北海道","青森","岩手","宮城","秋田","山形","福島","茨城","栃木","群馬","埼玉","千葉","東京","神奈川","新潟","富山","石川","福井","山梨","長野","岐阜","静岡","愛知","三重","滋賀","京都","大阪","兵庫","奈良","和歌山","鳥取","島根","岡山","広島","山口","徳島","香川","愛媛","高知","福岡","佐賀","長崎","熊本","大分","宮崎","鹿児島","沖縄"]

plot(x=X, y=Y, label=Labels, Geom.point, Geom.label, Theme(major_label_font="CMU Serif",minor_label_font="CMU Serif",major_label_font_size=12pt,minor_label_font_size=12pt), Guide.xlabel("女性有配偶率"),Guide.ylabel("合計特殊出生率"))

沖縄は特異点になっているが，合計特殊出生率とこれに対応する年代の女性の有配偶率には正の相関がみられる。

図：女性の有配偶率と合計特殊出生率の関係

P. S.  x軸とy軸の日本語ラベルのつけ方は，ChatGPT-4oに教えてもらった。

合計特殊出生率（２）

合計特殊出生率（１）からの続き

去年の今ごろ，奥村晴彦さんの合計特殊出生率のRプログラムを写経していた。すっかり忘れていたけれど，2023年の新しい合計特殊出生率が過去最低の1.20人を記録したいうニュースで思い出すことになった。

都道府県別でみれば，東京都の0.99や北海道の1.06をはじめ，大都市部や東北地方で低く，沖縄の1.60を筆頭に，九州，山陰，北陸などで高い。

今年は，コマンドラインで実行すべく，Rプログラムを修正した。あいかわらず日本語表示でつまづいたので，ChatGPT-4oに聞きながら試行錯誤した結果，何とか解決できた。なお，スタートを自分が生まれた1953年にしてみた。

修正したRプログラムのtfr.rの内容（注：remove \ before copy and paste）は

library(showtext)

font_add("Hiragino", "/System/Library/Fonts/ヒラギノ角ゴシック W3.ttc")

showtext_auto()

X <- read.table('/Users/koshi/Desktop/birth.txt', header = TRUE)

par(family = "Hiragino")  # 利用可能なフォントを指定

plot(X\$年, X\$合計特殊出生率, type="o", pch=16, xlab="", ylab="")

t = c(range(X\$合計特殊出生率), X\$合計特殊出生率[length(X[,1])])

axis(4, t, t)

title("合計特殊出生率", line=0.5)

#title("Total Fertility Rate", line=0.5)

# x軸とy軸のメモリの位置を手動で指定

x_ticks <- pretty(X$年, 15)  # 15個のメモリを生成

y_ticks <- pretty(X$合計特殊出生率, 16)  # 16個のメモリを生成

# x軸に対してグリッド線を描画

for (i in x_ticks) {

  abline(v = i, col = "lightgray", lty = "dotted")

}

# y軸に対してグリッド線を描画

for (i in y_ticks) {

  abline(h = i, col = "lightgray", lty = "dotted")

}

これで，Rscript tfr.r とすると，pdfの出力，Rplots.pdf が得られる。なお，showtextというライブラリは，コマンドラインで R として立ち上げたプロンプトに，install.packages("extrafont") としてインストールしておく。また，データは birth.txt にあらかじめ入れてある。

年    出生数 合計特殊出生率

1953 1868040           2.70

1954 1769580           2.48

1955 1730692           2.37

1956 1665278           2.22

1957 1566713           2.04

1958 1653469           2.11

1959 1626088           2.04

1960 1606041           2.00

1961 1589372           1.96

1962 1618616           1.98

1963 1659521           2.01

1964 1716761           2.05

1965 1823697           2.14

1966 1360974           1.58

1967 1935647           2.23

1968 1871839           2.13

1969 1889815           2.13

1970 1934239           2.13

1971 2000973           2.16

1972 2038682           2.14

1973 2091983           2.14

1974 2029989           2.05

1975 1901440           1.91

1976 1832617           1.85

1977 1755100           1.80

1978 1708643           1.79

1979 1642580           1.77

1980 1576889           1.75

1981 1529455           1.74

1982 1515392           1.77

1983 1508687           1.80

1984 1489780           1.81

1985 1431577           1.76

1986 1382946           1.72

1987 1346658           1.69

1988 1314006           1.66

1989 1246802           1.57

1990 1221585           1.54

1991 1223245           1.53

1992 1208989           1.50

1993 1188282           1.46

1994 1238328           1.50

1995 1187064           1.42

1996 1206555           1.43

1997 1191665           1.39

1998 1203147           1.38

1999 1177669           1.34

2000 1190547           1.36

2001 1170662           1.33

2002 1153855           1.32

2003 1123610           1.29

2004 1110721           1.29

2005 1062530           1.26

2006 1092674           1.32

2007 1089818           1.34

2008 1091156           1.37

2009 1070036           1.37

2010 1071305           1.39

2011 1050807           1.39

2012 1037232           1.41

2013 1029817           1.43

2014 1003609           1.42

2015 1005721           1.45

2016  977242           1.44

2017  946146           1.43

2018  918400           1.42

2019  865239           1.36

2020  840835           1.33

2021  811622           1.30

2022  770759           1.26

2023  727277           1.20

図：1953-2023  の合計特殊出生率の推移

［１］令和５年(2023)人口動態統計月報年計（概数）の概況

散歩地図

チョイソコてんり（２）からの続き

朝の散歩，今日の目標は北東にあるチョイソコてんりの小田中町公民館だったので，片道35分ほどだった。グーグルマップを見誤ったせいでちょっと遠回りになり40分かかった。

そんなわけで，自宅を中心として地図上に一定距離の半径の円を描いたものがあるとなんとなくうれしいような気がした。こういうときは，地理情報システムGISに対応したアプリとかデータとかプログラムを駆使することになる。Mathematicaの該当部分の関数を眺めたところ，それだけでおなかが一杯になって挫折した。

さあ，ChatGPT-4oの出番である。早速お願いしたところ一発で目的が達成できた。ゴイゴイスー（津田さんは中田先生の女婿だって知ってました？）。

Ｑ（koshix）：

　天理市の地図のうえで，自宅を中心として半径が500ｍの倍数となる円を1倍から6倍まで描画したいのです。なるべく簡単で具体的な手順や方法あるいはプログラムを一つ教えてください。言語は，JuliaでもMathematicaでもPythonでもRでも構いません

Ａ（GPT-4o）：

Pythonを使用して天理市の地図上に自宅を中心として半径が500mの倍数となる円を描画する方法を示します。この手順では、以下のライブラリを使用します。

・folium: インタラクティブな地図の作成

・geopy: 緯度経度の計算

以下の手順に従ってください。

１．必要なライブラリをインストールします。

(sh)

pip install folium geopy

２．次に、Pythonスクリプトを作成します。以下のコードをコピーして実行してください。

(python)

import folium

from geopy.distance import geodesic

# 自宅の緯度と経度を設定（天理市の例として）

home_location = (34.596047, 135.819249)

# フォリウムマップの作成

m = folium.Map(location=home_location, zoom_start=15)

# 半径500mの倍数となる円を1倍から6倍まで描画

for i in range(1, 7):

    radius = 500 * i

    folium.Circle(

        location=home_location,

        radius=radius,

        color='blue',

        fill=True,

        fill_opacity=0.1

    ).add_to(m)

# 自宅のマーカーを追加

folium.Marker(location=home_location, popup="Home").add_to(m)

# 地図を保存

m.save('tenri_map.html')

このコードは、以下の手順で動作します。

１．home_locationに自宅の緯度と経度を設定します。この例では天理市の中心付近の座標を使用していますが、実際の自宅の座標に置き換えてください。

２．folium.Mapを使用して地図を作成し、初期の表示位置とズームレベルを設定します。

３．半径500mの倍数となる円を1倍から6倍まで描画します。folium.Circleを使用して円を描画し、radius引数で半径をメートル単位で指定します。

４．自宅の位置にマーカーを追加します。

５．作成した地図をHTMLファイルとして保存します。

保存されたtenri_map.htmlファイルをブラウザで開くと、自宅を中心とした円が描かれた地図が表示されます。

ということで，説明も十分読まずに，pythonプログラムをコピペしてjupyter環境で実行すると，アラ不思議，あっというまにインタラクティブな地図上に同心円が書けてしまった。見栄えのために修正したのは透明度だけだ。

図：散歩地図（GPT-4oがつくったpythonプログラムの出力）

purl

Zennのあたりに五万度が一つ落ちていたのでインストールしてみた。

普段は，brew もしくは pip もしくはパッケージでしかインストールしないので，makeするのは久しぶりだ。プログラミング言語は，go でインストールしていたはずだったが見つからないのでbrewで再インストール。Github のcatasuy/purlからソースコードを引っ張ってきて，makeするだけ。.gitがないといわれたので，空の.gitを作ってやり直しているうちに出来ていた。

purl というコマンドだ。sed や grep では正規表現がperl互換ではなくて使いにくいことや，perlのone linerでもいろいろ不都合があるということで開発されたものだ。特徴は，perl 正規表現互換で，OS非依存で，軽量で，オプションがわかりやすいということ。もちろん日本語OKだ。

コマンドの実行例は

(1) purl -replace "@search@replace@" sample.txt 置換確認

(2) purl -overwrite -replace "@search@replace@" sample.txt 置換実行

(3) cat sample.txt | purl -replace "@search@replace@" パイプライン

(4) purl -filter "error" -filter "warning" sample.txt 複数ファイルフィルタ

(5) purl -exclude "debug" -exclude "info" sample.txt 複数ファイル除去

普段の作業は，正規表現を用いて次のように実現できた。

purl -replace "@^.*staff@@" spr.txt > spr0.txt

astropy

「Pythonで太陽系の天体位置を描画してみる」というのを見つけた。python上の天文学統合データ処理環境であるastropyというライブラリを使うものだ。

jupyter lab を起動して，さっそく試してみると，いきなりmatplotlibが見つからないというところでつまづいた。以前に試したpythonのコードでは，matplotlibを読み込んでちゃんと作図がでてている。その古いコードでもエラーがでてきた。

macbook air の環境は，基本的にはhomebrew で更新しているので安心できるのだが，問題は，python環境である。python は3.9〜3.12までの複数のバージョンが同居することになっている。これは基幹ソフトなので，下手に古いものを消してしまうとややこしいのだ。そこで，一般には，python環境を指定するpyenvとかvenvで仮想環境を設定した上で，様々なライブラリを加えていくことになる。

が，ですね，それはそれで面倒なので，アドホックにバージョンアップを繰り返しながら，pip install ホゲホゲをつづけているために，わけわかめ状態になっている今日この頃なのだ。jupyterlabもバージョンが上がるたびに，見えなくなってしまうので，毎回 brew link jupyterlabを繰り返す始末だ。70歳になると，自分の体の中のDNA情報システムだけでなく，外の電子情報システムまで老化が進んでくるということだ。

まあ，macOSの場合は，機種更新の際に新規インストールすればそれなりに，なんとかなるかもしれないけれど，複雑な記録や認証情報を引き継ぐのは厄介だ。特に，iOSでは，指紋認証を含め，現金が動く認証アプリがウヨウヨいるので，次期機種更新のことを考えただけで気が遠くなりそうだ。

話を戻して，astropyについて。jupyter環境でpython kernelを立ち上げているところから，!エスケープで!pip listをやってみると，matplotlibは存在している。!pip install matplotlib でも致命的エラーはでない。ChatGPTに相談してみると，カーネルがどうなっているか確認せよとのこと。

 jupyter kernelspec list

  python3              /opt/homebrew/Cellar/jupyterlab/4.0.9_2/libexec/lib/python3.12/site-packages/ipykernel/resources

  julia-1.9            /Users/koshi/Library/Jupyter/kernels/julia-1.9

  maxima               /Users/koshi/Library/Jupyter/kernels/maxima

  wolframlanguage12    /Users/koshi/Library/Jupyter/kernels/wolframlanguage12

はいはい，表に見えていたのは python3.11.6だったけれど，実は，python.12.1が使われており，ここにはmatplotlibが入っていませんでした。さっそく，このポイントで，matplotlibとかscipyとかastropyとかpytest-astropyとかをpipでインストールした。別の参考資料に指示があった，astropyが正しく導入できているかどうかのテストも，試行錯誤の後にOKとなった。

　　　cd /opt/homebrew/Cellar/jupyterlab/4.0.9_2/libexec/bin

　　　./pip install astropy

　　　./pip install ephem

import astropy as ap

>>> ap.test()

platform darwin -- Python 3.12.1, pytest-7.4.3, pluggy-1.3.0

Running tests with Astropy version 6.0.0.

Running tests in lib/python3.12/site-packages/astropy.

Date: 2023-12-19T22:26:53

Platform: macOS-14.2-arm64-arm-64bit

Executable: /opt/homebrew/Cellar/jupyterlab/4.0.9_2/libexec/bin/python

Full Python Version: 

3.12.1 (main, Dec  7 2023, 20:45:44) [Clang 15.0.0 (clang-1500.0.40.1)]

encodings: sys: utf-8, locale: UTF-8, filesystem: utf-8

byteorder: little

float info: dig: 15, mant_dig: 15

Package versions: 

Numpy: 1.26.2

Scipy: 1.11.4

Matplotlib: 3.8.2

h5py: not available

Pandas: not available

PyERFA: 2.0.1.1

Cython: not available

Scikit-image: not available

asdf-astropy: not available

Using Astropy options: remote_data: none.

rootdir: /opt/homebrew/Cellar/jupyterlab/4.0.9_2/libexec

plugins: astropy-0.11.0, remotedata-0.4.1, hypothesis-6.92.1, cov-4.1.0, filter-subpackage-0.1.2, astropy-header-0.2.2, anyio-4.1.0, mock-3.12.0, doctestplus-1.1.0, arraydiff-0.6.1

collected 28085 items / 5 skipped  

-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html

========== 27425 passed, 425 skipped, 240 xfailed, 80 warnings in 184.56s (0:03:04) ===========

<ExitCode.OK: 0>

ところが，やっぱりもとの太陽系の天体位置描画が動かない。最後に，Animation関数中のblit=False をblit=Trueにしてようやく何か表示されてほっとしたが，そこまでで内部にエラーがあるとのメッセージが出て異常終了してしまった。チーン。

図：ここまでは出たけれどアニメーションにはなっていなかった残念な図

astropyはイメージ処理やデータ処理を本業とした，プロ用の本格的なシステムだった。

久々のJulia

プログラミングはAIに任せる時代になってしまったので，自分でコードを書くのは棋士の卵が詰め将棋の練習をしているようなことかもしれない。

Homebrewでmacにインストールされているソフトウェアを定期的に更新しているが，jupyterlabは，毎回，手動により brew link jupyterlab をすることを強いるのでうっとうしい。そんなこともあって，しばらくJuliaから足が遠のいていた。

鈴木貫太郎のYouTubeの数学の問題では，しばしば，整数の冪を比較するような問題が登場する。そこで，2^m と3^n が非常に近くなるような整数を探してみようと思った。m,nを1から順に条件を判定しながら増やしていけばよいかと思ったけれど，実際にコードを書いてみると，面倒なロジックは組まずに虱潰しに調べたほうが簡単そうであった。

久々にプログラミングすると，浮動小数点数から整数への変換がどうなっているかなど，基本的なところを忘れてしまっている。困ったものだ。これはconvert(T, x)だと書いてあったのだが，エラーになる。ChatGPTに相談してみたら，floor(T,x)とceil(T,x)を使いなさいとのことだった。その結果が次のコードである。

function dt(N,eps)

# N=1000000

# eps = 1.0e-6

ratio = log(2)/log(3)

    for i in 1:N

        jmin = floor(Int64,i*ratio)

        jmax = ceil(Int64,i*ratio)

        for j in jmin:jmax

            sol = i*log(2) - j*log(3)

            if abs(sol) < eps

                println("  2**",i,"/","3**",j," = ",exp(sol))

            end

        end

    end

end

for k = 1:8

    N=10^k

    eps=1/N

    println("N=",N," eps=",eps)

    @time(dt(N,eps))

end

N=10 eps=0.1

  2**8/3**5 = 1.0534979423868305

  0.000176 seconds (98 allocations: 2.750 KiB)

N=100 eps=0.01

  2**84/3**53 = 0.9979140462573083

  0.000160 seconds (97 allocations: 2.562 KiB)

N=1000 eps=0.001

  0.000003 seconds

N=10000 eps=0.0001

  2**1054/3**665 = 0.9999563468421858

  2**2108/3**1330 = 0.9999126955899699

  0.000351 seconds (207 allocations: 5.609 KiB)

N=100000 eps=1.0e-5

  2**50508/3**31867 = 0.9999927350845753

  0.000408 seconds (103 allocations: 2.656 KiB)

N=1000000 eps=1.0e-6

  2**301994/3**190537 = 1.0000000644940903

  2**603988/3**381074 = 1.0000001289881848

  2**905982/3**571611 = 1.0000001934822835

  0.003297 seconds (311 allocations: 8.562 KiB)

N=10000000 eps=1.0e-7

  2**301994/3**190537 = 1.0000000644940903

  0.027930 seconds (105 allocations: 4.938 KiB)

N=100000000 eps=1.0e-8

  2**85137581/3**53715833 = 0.9999999925494194

  0.264204 seconds (229 allocations: 17.953 KiB)

 いちおうできたことにしておこう。

Bluesky（２）

Bluesky（１）からの続き

いよいよTwitterがややこしくなってきた。7月1日，イーロン・マスク（1971-）はTwitterアカウントの種類に応じて一時的に閲覧回数の制限を設けた。認証済みアカウントは6000件（1万件），未認証アカウントで600件（1000件），作成されたばかりの未認証アカウントで300件（500件）というわけだ（括弧内は7月2日の泥縄的修正値）。制限の理由については，インターネット上のデータ収集や処理における不正な操作の急増に対応するためとしている。

実際には，TwitterがAPI課金を急騰させたことを原因とするウェブスクレーピングの増加だとか，サーバーを運用しているクラウドとの契約の関係だとか，Twitterコードのバグが自分自身に対するDOS攻撃をしたのが原因だとか，諸説がある。いずれにせよTwitterのコントロールがいよいよ効いていないことが明らかになった。

このため，再びTwitterから他のSNSへの流失が進んでいるらしい。その候補の一つはジャック・ドーシー（1976-）がはじめたBlueskyだけれど，招待制の登録者はまだ20万人だ（Twitterは4.5億人）。他に T2 とか MetaのThreads が話題になっている。どうなることか。現在のところ，Twitterに流れる情報の方がまだまだ有用度が高い。

そのBlueskyの歩き方というドキュメントがmd形式で流れていた。pdfに変換しようと考えた。たぶんpandocが使えるはずなのだが，調べてみると一手間必要だった。これをまとめると次のようになる。

$ brew install pandoc ✓

$ brew cask install mactex ✓

~/.zshrc に次の行を追加

export PATH=/usr/local/texlive/2019/bin/x86_64-darwin/:$PATH

$ source ~/.zshrc

$ pandoc bluesky.md -o bluesky.pdf --pdf-engine=xelatex -V documentclass=bxjsarticle -V classoption=pandoc

これをシェルスクリプトにすると

$ cat md2pdf.sh

#!/bin/zsh

# usage mdspdf.sh test

# input test.md output test.pdf

#

pandoc $1.md -o $1.pdf --pdf-engine=xelatex -V documentclass=bxjsarticle -V classoption=pandoc

これで無事にmdファイルをpdfファイルに変換できた。

P. S. Blueskyの招待コードが若干余っているので必要な方はメール，TwitterのDM，Facebook Messengerなどでご連絡ください。

合計特殊出生率（１）

奥村晴彦先生が，Rで合計特殊出生率のグラフを描いていたのでまね（写経）をしてみた。

厚生労働省の発表では，2022年の合計特殊出生率が1.26となって（前年は1.30），2005年に並んで過去最低を記録した，というのがニュースになったからかな。

写経は成功したが，日本語タイトルが文字化けしていた。ネットで検索すると解決できた。次に，グリッドラインを入れるべく検索したが要領を得ないし，なんだかずれてしまう。GPT-4に相談したところ2回目で正解にたどりつけた。最後に，奥村先生はクリップボードからデータを入力するという高等技術を使っていたので，低レベルのテキストファイル入力に切り替えた。自分で考える力がどんどん喪失されていく・・・

その結果が次のRプログラムと図である。

# X = read.table(pipe("pbpaste"), header=TRUE) （奥村先生のテクニック for macOS）

X <- read.table('/Users/koshi/Desktop/birth.txt', header =TRUE)

par(family = "HiraKakuProN-W3")

plot(X$年, X\$合計特殊出生率, type="o", pch=16, xlab="", ylab="")

t = c(range(X\$合計特殊出生率), X$合計特殊出生率[length(X[,1])])

axis(4, t, t)

title("合計特殊出生率", line=0.5)

# x軸とy軸のメモリの位置を手動で指定

x_ticks <- pretty(X\$年, 11)  # 11個のメモリを生成

y_ticks <- pretty(X\$合計特殊出生率, 9)  # 9個のメモリを生成

# x軸に対してグリッド線を描画

for (i in x_ticks) {

  abline(v = i, col = "lightgray", lty = "dotted")

}

# y軸に対してグリッド線を描画

for (i in y_ticks) {

  abline(h = i, col = "lightgray", lty = "dotted")

}

図：合計特殊出生率の推移（1970-2022）

数独

日本経済新聞の土曜版のクイズ欄では，数独（ナンプレ）が毎週出題されていた。以前は赤で配置可能な数字をメモしながら解くものだと信じていたのだが，ある日電車の隣に座ったおじさんが，メモなしでスイスイ解いているのを見てしまった。それ以降は，簡単なものなら補助メモなしで解けるようになった。なにごとも信念が重要である。

その土曜版のクイズ欄は，4月から漢字クイズに変わってしまった。えぇーっ。せっかく惚け防止のために頭の体操（10-20分）をしていたのに・・・。まあ，数独を含む各種クイズがセットになった日経脳活動クイズのiPhoneアプリ版があるといえばあったけれど消してしまった。紙に書かないと，落ち着かない旧世代なのだ。

それならば，自分で数独の問題を作るプログラムを書けばよい。さっそくGPT-4に聞いてみたら，pythonプログラムを出してきた（ありふれているので学習済みだろう）。さらに，第2案として，sudokuモジュールをインポートする2-3行ですむ簡単な方法が提案された。pip でsudokuをインストールせよとのことだが，それではだめで，py-sudokuならうまくいった。

GPT-4と相談しながら，修正を繰り返し，結果をcsvファイルに格納するプログラムができた。これをexcelで読み込んで出力すればOKだ。

#!/Users/koshi/bin/python

import csv

from sudoku import Sudoku

# Initializes a Sudoku puzzle with 3 x 3 sub-grid and

# generates a puzzle with 75% of the cells empty

puzzle = Sudoku(3).difficulty(0.65)

puzzle.show()

# Convert the puzzle to a list of lists

puzzle_data = []

for row in puzzle.board:

    puzzle_data.append(list(row))

# Write the puzzle data to a CSV file

with open('puzzle.csv', 'w', newline='') as csvfile:

    csv_writer = csv.writer(csvfile)

    for row in puzzle_data:

        # Replace 0 or None with an empty string to represent empty cells

        formatted_row = [str(cell) if cell not in (0, None) else '' for cell in row]

        csv_writer.writerow(formatted_row)

 その結果わかったこと。日経のクイズは超難問の場合でも解けるようにかなり配慮されている。難易度を易しめに設定しても，単純にランダムに消すことで数独問題をつくるとかなり解くのが難しい極値にはまり込むことがあるような気がする。

図：GPT-4と協力して作成したコードの出力例

［１］py-sudoku 1.0.3（pypi.org）