算額（その445）

算額（その445）

十二 群馬県高崎市石原町 清水寺 享和3年(1803)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

享和3年癸亥4月(1803) 上州中里邑 黒崎九兵祀則
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

キーワード：円4個，外円，菱形
#Julia, #SymPy, #算額, #和算

外円内に甲円 1 個，乙円 2 個，菱形 1 個が入っている。外円と甲円の直径がそれぞれ 7 寸，3 寸のとき，乙円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r0 - r1, 0)
乙円の半径と中心座標を r2, (x2, y2)
とおき，連立方程式の解を求める。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive,
     r2::positive, x2::positive, y2::positive;

(r0, r1) = (7, 3) .// 2
eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = distance(0, r0 - 2r1, sqrt(r0^2 - r1^2), -r1, x2, y2) - r2^2 
res = solve([eq1, eq2, eq3], (r2, x2, y2));

解は 1 組得られるが，係数が分数の長い式になるので，evalf() した値を示しておく。

乙円の直径は 2 * 1.2 = 2寸4分である。

(res[1][1].evalf(), res[1][2].evalf(), res[1][3].evalf())

   (1.20000000000000, 2.24499443206436, 0.500000000000000)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (7, 3) .// 2
   (r2, x2, y2) = (res[1][1].evalf(), res[1][2].evalf(), res[1][3].evalf())
   @printf("r2 = %g;  x2 = %g;  y2 = %g\n", r2, x2, y2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :red)
   circle(0, r0 - r1, r1, :blue)
   circle2(x2, y2, r2, :brown)
   b = sqrt(r0^2 - r1^2)
   plot!([0, b, 0, -b, 0], [-r0, -r1, r0 - 2r1, -r1, -r0], color=:orange, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r0 - r1, " 甲円:r1, (0,r0-r1)", :black, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :brown, :center, delta=-delta)
       point(0, r0 - 2r1, " r0-2r1")
       point(0, -r1, " -r1")
       point(0, r0, " r0", :red, :left, :bottom)
       point(sqrt(r0^2 - r1^2), -r1, "(√(r0^2-r1^2),-r1) ", :orange, :right, :vcenter)
   end
end;

算額（その444）

算額（その444）

十二 群馬県高崎市石原町 清水寺 享和3年(1803)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

享和3年癸亥4月(1803) 上州棟高邑 坂本丹治亮春
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

キーワード：円5個，長方形

長方形内に大円 1 個，中円，小円がそれぞれ 2 個入っており，互いに内接・外接しあっている。長方形の長辺が 40寸8分のとき，それぞれの円の直径を求めよ。

原点を図のようにおき，長方形の長辺を a，短辺を 2b とする。
大円の半径と中心座標を r1, (a - r1, 0)
中円の半径と中心座標を r2, (r2, b - r2)
小円の半径と中心座標を r3, (r3, r3)
とおき，連立方程式の解を求める。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive;

b = r1
eq1 = (a - r1 - r2)^2 + (b - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + r3^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + (b - r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (a*(9 - 4*sqrt(2))/8, a/8, a*(3/2 - sqrt(2)))
    (a*(4*sqrt(2) + 9)/8, a/8, a*(sqrt(2) + 3/2))

最初のものが適解である。

a = 40寸8分 のとき各円の直径は 34寸1分あまり, 10寸2分, 7寸あまりである。

2(a*(9 - 4*sqrt(2))/8)(a => 40.8).evalf() |> println
2(a/8)(a => 40.8).evalf() |> println
2(a*(3/2 - sqrt(2)))(a => 40.8).evalf() |> println

   34.1000866551777
   10.2000000000000
   7.00017331035544

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 408//10
   (r1, r2, r3) = (a*(9 - 4*sqrt(2))/8, a/8, a*(3/2 - sqrt(2)))
   b = r1
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   plot([0, a, a, 0, 0], [-b, -b, b, b, -b], color=:black, lw=0.5)
   circle(a - r1, 0, r1, :red)
   circle(r2, b - r2, r2, :blue)
   circle(r2, r2 - b, r2, :blue)
   circle(r3, r3, r3, :brown)
   circle(r3, -r3, r3, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a - r1, 0, " 大円:r1,(a-r1,0)", :red, delta=-delta)
       point(r2, b - r2, " 中円:r2,(r2,b-r2)", :black)
       point(r3, r3, " 小円:r3,(r3,r3)", :black)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その443）

算額（その443）

天明五年乙巳秋九月(1785) 藤田貞資門人 筑後州久留米藩 清水與市道香

藤田貞資(1789)：神壁算法巻下
http://www.wasan.jp/jinpeki/jinpekisanpo2.pdf

求めるものと条件として与えられる数値が異なるだけで，算額（その31）と基本的に同じである。

図のように，直線の上に大円，中円，小円，甲円，乙円がそれぞれ外接しあっている。甲円，乙円の径を98寸1厘，121寸 としたとき，中円の径を求めよ。

それぞれの円の中心座標と半径を (0, r1, r1), (x2, r2, r2), (x3, r3, r3), (x4, r4, r4), (x5, r5, r5) とする。

include("julia-source.txt")

using SymPy

@syms r1, r2, r3::positive, r4::positive, r5::positive
@syms x2::positive, x3::positive, x4::positive, x5::positive;

(r4, r5) = (981, 1210) .// 20;
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2;
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2;
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2;
eq4 = x5^2 + (r1 - r5)^2 - (r1 + r5)^2;
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2;
eq6 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2;
eq7 = (x5 - x3)^2 + (r3 - r5)^2 - (r3 + r5)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x2, x3, x4, x5))

   2-element Vector{NTuple{7, Sym}}:
    (3917133*sqrt(1090)/3682898 + 304705467/7365796, 121/2, 3917133*sqrt(1090)/52441 + 260073891/104882, 3993*sqrt(1090)/2714 + 118701/1357, 188259786/310753 + 11751399*sqrt(1090)/621506, 118701/1357 + 32373*sqrt(1090)/6785, 3993*sqrt(1090)/2714 + 118701/1357)
    (35254197*sqrt(1090)/3682898 + 2742349203/7365796, 70508394*sqrt(1090)/14891881 + 6218626689/29783762, 2340665019/104882 - 35254197*sqrt(1090)/52441, 2340665019/5236663 + 176270985*sqrt(1090)/10473326, -401684184/310753 + 35254197*sqrt(1090)/621506, 356103/1357 + 97119*sqrt(1090)/6785, 11979*sqrt(1090)/2714 + 356103/1357)

2番めのものが適解である。

name = ["r1", "r2", "r3", "x2", "x3", "x4", "x5"]
j = 2
for i in 1:7
   println("$(name[i]) = $(res[j][i].evalf())")
end

   r1 = 688.343020749222
   r2 = 365.108908025960
   r3 = 122.232157448001
   x2 = 1002.63686078867
   x3 = 580.129821644955
   x4 = 734.990886123080
   x5 = 408.140920542540

中円の直径は 730.21781605192 である。算額の答えでは 729 となっており，かなり違いがある。

2res[2][2].evalf()

   730.21781605192

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, r5) = (981, 1210) ./ 20
   r1, r2, r3, x2, x3, x4, x5 = (35254197*sqrt(1090)/3682898 + 2742349203/7365796, 70508394*sqrt(1090)/14891881 + 6218626689/29783762, 2340665019/104882 - 35254197*sqrt(1090)/52441, 2340665019/5236663 + 176270985*sqrt(1090)/10473326, -401684184/310753 + 35254197*sqrt(1090)/621506, 356103/1357 + 97119*sqrt(1090)/6785, 11979*sqrt(1090)/2714 + 356103/1357)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  r5 = %g\n", r1, r2, r3, r4, r5)
   @printf("中円の直径 = %g\n", 2r2)
   plot()
   circle(0, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :brown)
   circle(x4, r4, r4, :green)
   circle(x5, r5, r5, :red)
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, r1, "大円 r1")
       point(x2, r2,"中円 r2")
       point(x3, r3,"小")
       point(x4, r4,"乙")
       point(x5, r5,"甲")
   end
end;

Plotlyを使用した対話型データ可視化の作成

Plotlyを使用した対話型データ可視化の作成

データの可視化は，データ分析と情報の共有において極めて重要である。Plotlyは，対話的なデータ可視化を作成するための強力なPythonライブラリで，グラフや図表を対話的に操作し，データのパターンを探索するのに役立つ。この記事では，Plotlyを使って対話的なデータ可視化を作成する基本的な手法について説明する。

Plotlyのインストール

まず最初に，Plotlyをインストールする。

pip install plotly

折れ線グラフの作成

Plotlyを使用して折れ線グラフを作成する基本的なステップを見てみよう。

import plotly.express as px

# データ
x = [1, 2, 3, 4, 5]
y = [10, 16, 5, 12, 7]

# 折れ線グラフを作成
fig = px.line(x=x, y=y, title='折れ線グラフの例')
fig.show()

このコードでは，Plotly Expressを使用して簡単に折れ線グラフを作成し，対話的に操作することができる。

散布図の作成

散布図はデータの分布や関係性を視覚化するのに役立つ。

# 散布図を作成
fig = px.scatter(x=x, y=y, title='散布図の例')
fig.show()

ヒートマップの作成

データのパターンや相関をヒートマップで可視化する。

# ヒートマップを作成
data = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
fig = px.imshow(data, title='ヒートマップの例')
fig.show()

対話的な機能の活用

Plotlyを使えば，グラフをズームイン，ズームアウト，ドラッグ，ツールチップ表示などの対話的な機能を簡単に追加できる。また，Dashというライブラリを使用すれば，Plotlyの可視化をウェブアプリケーションに統合することもできる。

まとめ

Plotlyはデータ科学者やデータエンジニアにとって優れたツールであり，対話的なデータ可視化を作成し，データの洞察を得るのに役立つ。この記事では基本的な操作を紹介しましたが，Plotlyにはさらに多くの機能がある。詳細なカスタマイズや高度なプロットの作成を学ぶことで，データの探索と解釈を強化できる。

NumPyとPandasを組み合わせて高速なデータ操作を実現する方法

NumPyとPandasを組み合わせて高速なデータ操作を実現する方法

データ操作はデータサイエンスと機械学習の中で不可欠なスキルである。NumPyとPandasはPythonのデータ操作において非常に強力なツールであるが，それらを組み合わせて利用することで，さらに高速かつ効率的なデータ処理が可能である。この記事では，NumPyとPandasを連携させてデータ操作を最適化する方法について説明する。

1. NumPyの配列をPandas DataFrameに変換

NumPyの高速な配列操作を利用しながら，Pandas DataFrameの柔軟性を享受するために，NumPyの配列をPandas DataFrameに変換する。

import numpy as np
import pandas as pd

# NumPy配列を作成
data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

# NumPy配列をPandas DataFrameに変換
df = pd.DataFrame(data, columns=['A', 'B', 'C'])

# DataFrameを表示
print(df)

2. データの選択とフィルタリング

Pandas DataFrame内でデータを選択およびフィルタリングする際に，NumPyのブールインデックスを活用する。

# 条件に合致する行を選択
selected_data = df[df['A'] > 3]

# 特定の列を選択
selected_column = df['B']

# 条件に合致するデータを置き換え
df.loc[df['A'] > 3, 'B'] = 0

3. データの計算と集約

NumPyのベクトル化演算を使用して，Pandas DataFrame内で効率的な計算と集約を行う。

# 列ごとの合計を計算
column_sum = df.sum()

# 条件に合致する行の合計を計算
conditional_sum = df[df['A'] > 3].sum()

# 列ごとの統計情報を取得
column_stats = df.describe()

4. グループ化と集約

Pandasのgroupbyメソッドを使用してデータをグループ化し，NumPyの関数を適用して集約する。

# グループ化して平均値を計算
grouped_data = df.groupby('A').mean()

# グループごとのデータ数を計算
group_sizes = df.groupby('A').size()

これらのテクニックを組み合わせることで，NumPyとPandasを利用して高速かつ効率的なデータ操作を実現できる。データサイエンスプロジェクトにおいて，データのクリーニング，前処理，分析，モデリングなどのステップにおいて，この連携は非常に役立つ。