制約充足問題

$Python$の制約充足問題を解くためにpython-constraintライブラリを使用したサンプルコードを示します。

この例では、$3$つの変数$（x、y、z）$に対する制約条件を設定し、制約条件を満たす解を見つけます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

from constraint import Problem

# 問題を作成
problem = Problem()

# 変数の定義
problem.addVariable('x', [1, 2, 3])
problem.addVariable('y', [4, 5, 6])
problem.addVariable('z', [7, 8, 9])

# 制約条件の定義
def custom_constraint(x, y, z):
    if x + y > z:
        return True

problem.addConstraint(custom_constraint, ['x', 'y', 'z'])

# 問題を解く
solutions = problem.getSolutions()

# 解の表示
for solution in solutions:
    print(solution)

このコードでは、Problemオブジェクトを作成し、addVariableメソッドを使用して変数を定義します。

次に、addConstraintメソッドを使用して制約条件を定義します。

最後に、getSolutionsメソッドを使用して制約条件を満たす解を見つけます。

制約条件は、$3$つの変数$（x、y、z）$の値が与えられ、それらの値が特定の条件を満たすかどうかを評価します。

この例では、$x$と$y$の合計が$z$より大きい場合に制約条件を満たします。

[実行結果]

{'x': 3, 'y': 6, 'z': 8}
{'x': 3, 'y': 6, 'z': 7}
{'x': 3, 'y': 5, 'z': 7}
{'x': 2, 'y': 6, 'z': 7}

最小全域木（Minimum Spanning Tree, MST）問題

最小全域木（Minimum Spanning Tree, MST）問題は、グラフ理論における基本的な最適化問題の一つです。

与えられた重み付きの無向グラフにおいて、すべてのノードをつなぐために必要なエッジの集合であって、その総重みが最小となるような部分グラフを見つける問題です。

具体的には、無向グラフ $ ( G = (V, E) ) $が与えられたとき、その部分グラフ$ ( T ) $で以下の条件を満たすものを最小全域木と呼びます：

1. $( T ) $は$ ( G ) $の全域木である（すべてのノードが連結されている）。
2. $( T ) $のエッジ数は$ ( V - 1 ) $である。
3. $( T ) $のエッジの重みの合計が最小である。

最小全域木問題は、通信ネットワーク、電気回路、道路ネットワークなどの多くの応用分野で重要な役割を果たしています。

具体的な例としては、都市間の最短経路を確立する通信ネットワークや、最小コストですべての地域を結ぶ電力網などがあります。

最小全域木問題は、PrimのアルゴリズムやKruskalのアルゴリズムなどのアルゴリズムを用いて解くことができます。
これらのアルゴリズムは、グラフのエッジの重みを考慮して、最小の重みを持つエッジを選択しながら部分グラフを構築していきます。

要するに、最小全域木問題は与えられたグラフ内のすべてのノードを最小のコストでつなぐ木を見つける問題であり、多くの応用分野で効果的に利用されています。

プログラム例

NetworkXを使用して最小全域木問題を解決し、その結果をグラフ化する方法を以下に示します。

まず、グラフを作成し、最小全域木を計算し、その結果を描画します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import networkx as nx
import matplotlib.pyplot as plt

# グラフの作成
G = nx.Graph()
G.add_weighted_edges_from([
    (0, 1, 1),
    (0, 2, 3),
    (1, 2, 1),
    (1, 3, 4),
    (2, 3, 2),
    (3, 4, 5),
    (4, 5, 1),
    (5, 0, 2)
])

# 最小全域木の取得
mst = nx.minimum_spanning_tree(G)

# グラフの描画
pos = nx.spring_layout(G)  # ノードの配置を計算
plt.figure(figsize=(12, 8))

# 元のグラフを描画（灰色のエッジ）
nx.draw(G, pos, with_labels=True, node_color='lightblue', edge_color='gray', node_size=500, font_size=15)

# 最小全域木を描画（赤色のエッジ）
nx.draw_networkx_edges(mst, pos, edge_color='red', width=2)

# エッジの重みを表示
edge_labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)

plt.title("Minimum Spanning Tree in the Graph")
plt.show()

【コードの説明】

1. networkxとmatplotlib.pyplotをインポートします。
2. グラフ G を作成し、ノード間の重み付きエッジを追加します。
3. nx.minimum_spanning_tree関数を使用して最小全域木 mst を取得します。
4. ノードの配置を計算するために nx.spring_layout を使用します。
5. plt.figure を使用して描画領域のサイズを設定します。
6. 元のグラフ G を灰色のエッジで描画します。
7. 最小全域木 mst を赤色のエッジで描画します。
8. エッジの重みを表示します。
9. グラフにタイトルを設定し、表示します。

このコードを実行すると、与えられたグラフの最小全域木が赤色のエッジで強調表示されたグラフが表示されます。

最小全域木は、グラフ内のすべてのノードを最小コストで連結する部分グラフです。

[実行結果]

ソースコード解説

このソースコードは、NetworkXを使用してグラフを作成し、最小全域木を計算し、結果を可視化するためのものです。

以下では、コードの各部分を詳しく説明します。

1. ライブラリのインポート

1
2

import networkx as nx
import matplotlib.pyplot as plt

• networkxはグラフ理論やネットワーク解析のためのPythonライブラリです。
• matplotlib.pyplotはグラフを描画するためのPythonライブラリです。

2. グラフの作成

1
2
3
4
5
6
7
8
9
10
11

G = nx.Graph()
G.add_weighted_edges_from([
    (0, 1, 1),
    (0, 2, 3),
    (1, 2, 1),
    (1, 3, 4),
    (2, 3, 2),
    (3, 4, 5),
    (4, 5, 1),
    (5, 0, 2)
])

• 無向グラフ G を作成します。
• add_weighted_edges_fromメソッドを使用して、重み付きのエッジを追加します。
  エッジはタプルの形で指定され、最初の2つの要素はノードのペアで、3番目の要素はエッジの重みです。

3. 最小全域木の取得

1

mst = nx.minimum_spanning_tree(G)

• nx.minimum_spanning_tree関数を使用して、与えられたグラフの最小全域木を計算します。

4. グラフの描画

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

pos = nx.spring_layout(G)  # ノードの配置を計算
plt.figure(figsize=(12, 8))

# 元のグラフを描画（灰色のエッジ）
nx.draw(G, pos, with_labels=True, node_color='lightblue', edge_color='gray', node_size=500, font_size=15)

# 最小全域木を描画（赤色のエッジ）
nx.draw_networkx_edges(mst, pos, edge_color='red', width=2)

# エッジの重みを表示
edge_labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)

plt.title("Minimum Spanning Tree in the Graph")
plt.show()

• nx.spring_layout関数を使用して、グラフ内のノードの配置を計算します。
• plt.figure(figsize=(12, 8))で描画領域のサイズを設定します。
• 元のグラフ G を灰色のエッジで描画します。
• 最小全域木 mst を赤色のエッジで描画します。
• nx.get_edge_attributes関数を使用して、グラフのエッジの属性（ここでは重み ）を取得します。
• nx.draw_networkx_edge_labels関数を使用して、エッジの重みを表示します。
• plt.title関数を使用して、グラフにタイトルを設定します。

以上が、与えられたソースコードの詳細な説明です。

これにより、グラフがどのように構成され、最小全域木がどのように計算されて可視化されるかが理解できるでしょう。

グラフ解説

[実行結果]

グラフに表示される内容を詳しく説明します。

1. ノード（頂点）:

   • グラフ内の各ノードは、円で表現されています。
     円の中には、ノードの名前が表示されています。
     この例では、ノードは$0$から$5$までの整数で表されています。

2. エッジ（辺）:

   • グラフ内の各エッジは、ノード間の接続関係を表しています。
     エッジは、ノード間の線で表現されます。
   • エッジの上に表示される数字は、エッジの重み（コストや距離など）を示しています。
     この例では、エッジの重みは整数です。

3. 最小全域木:

   • 最小全域木は、グラフ内のすべてのノードを含み、かつ全てのノードが直接または間接的に接続されている部分グラフです。
   • 最小全域木は赤色のエッジで強調表示されています。
   • 最小全域木は、元のグラフ内のすべてのノードを最小コストで連結するように構築されます。

4. エッジの重み:

   • エッジの上に表示される数字は、エッジの重みを示しています。
     最小全域木のエッジには重みが表示されますが、元のグラフのすべてのエッジには表示されません。

5. グラフのタイトル:

   • グラフのタイトルは、「Minimum Spanning Tree in the Graph」です。
     これは、描画されているグラフが元のグラフ内の最小全域木を示していることを示しています。

これらの要素が組み合わさって、グラフ全体が構成されています。

これにより、最小全域木がどのように元のグラフから選ばれるかや、各エッジの重みが最小全域木の構築にどのように影響するかがわかります。

Helicoid（ヘリコイド）

Helicoid（ヘリコイド）を描画してみましょう。

Helicoidは、微分幾何学や応用数学でよく知られている曲面の一つで、螺旋状の形状を持っています。

Helicoidは次のようなパラメトリック方程式で表されます。

$$
x(u, v) = u \cos(v)
$$
$$
y(u, v) = u \sin(v)
$$
$$
z(u, v) = c v
$$

ここで、$( c ) $は螺旋の巻きの間隔を調整する定数です。

プログラム例

以下に、PythonコードでHelicoidを描画する方法を示します。

必要なライブラリのインストール

1

pip install matplotlib numpy

PythonコードでHelicoidを描画

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメトリック方程式のパラメータ範囲を設定
u = np.linspace(-10, 10, 100)
v = np.linspace(0, 4 * np.pi, 100)
u, v = np.meshgrid(u, v)

# Helicoidの定数
c = 1

# Helicoidの方程式
x = u * np.cos(v)
y = u * np.sin(v)
z = c * v

# 3Dプロットを作成
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# サーフェスプロット
surf = ax.plot_surface(x, y, z, cmap='viridis')

# カラーバー
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

# タイトルと軸のラベルを設定
ax.set_title('Helicoid')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')

# プロットを表示
plt.show()

説明

1. numpyを使って$u$と$v$の範囲を設定し、グリッドを作成します。
2. Helicoidのパラメトリック方程式を使用して$x$, $y$, $z$の値を計算します。
3. matplotlibとmpl_toolkits.mplot3dを使って3Dグラフを描画します。
4. グラフには軸ラベル、タイトル、カラーバーが含まれています。

このコードを実行すると、螺旋状の形状を持つHelicoidを描画した3Dグラフが表示されます。

[実行結果]

ソースコード解説

ソースコードを詳しく説明します。

1. 必要なライブラリのインポート

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

説明

• numpy: 数値計算ライブラリであり、配列操作や数学関数の計算に使用します。
• matplotlib.pyplot: グラフ描画ライブラリであり、データの可視化に使用します。
• mpl_toolkits.mplot3d.Axes3D: matplotlibの3Dプロット機能を使用するためのモジュールです。

2. パラメトリック方程式のパラメータ範囲を設定

1
2
3

u = np.linspace(-10, 10, 100)
v = np.linspace(0, 4 * np.pi, 100)
u, v = np.meshgrid(u, v)

説明

• uとvは、Helicoidをパラメトリックに表現するためのパラメータです。
• np.linspace(start, stop, num)は、startからstopまでの範囲をnum個の等間隔で分割した値を生成します。
  • uは$-10$から$10$までの範囲で$100$個の値を持ちます。
  • vは$0$から$4π$までの範囲で$100$個の値を持ちます。
• np.meshgrid(u, v)は、uとvのグリッドを作成し、各組み合わせに対して対応するx, y, z座標を計算できるようにします。

3. Helicoidの定数

1

c = 1

説明

• cはHelicoidの螺旋の巻きの間隔を調整する定数です。
  この例では、cは$1$に設定されています。

4. Helicoidの方程式

1
2
3

x = u * np.cos(v)
y = u * np.sin(v)
z = c * v

説明

• Helicoidのパラメトリック方程式を用いて、各パラメータに対するx, y, zの値を計算します。
  • x = u * np.cos(v):
    uとvの値を使ってx座標を計算します。
  • y = u * np.sin(v):
    uとvの値を使ってy座標を計算します。
  • z = c * v:
    vの値に定数cを掛けてz座標を計算します。

5. 3Dプロットを作成

1
2

fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

説明

• fig = plt.figure(figsize=(10, 7)):
  図のサイズを指定して新しい図を作成します。
  figsizeは幅と高さをインチ単位で指定します。
• ax = fig.add_subplot(111, projection='3d'):
  3Dプロット用のサブプロットを追加します。
  • 111は$1$行$1$列のサブプロットグリッドを作成し、その$1$番目のサブプロットを指します。
  • projection='3d'は3Dプロットを作成することを指定します。

6. サーフェスプロット

1

surf = ax.plot_surface(x, y, z, cmap='viridis')

説明

• ax.plot_surface(x, y, z, cmap='viridis'):
  x, y, zのデータを使って3Dサーフェスプロットを作成します。
  • cmap='viridis'はカラーマップを指定し、サーフェスの色をz軸の値に応じて変化させます。

7. カラーバー

1

fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

説明

• fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5):
  サーフェスプロットのカラーバーを追加します。
  • shrink=0.5はカラーバーの高さを縮小して、全体の高さの$50%$にします。
  • aspect=5はカラーバーのアスペクト比を設定します。

8. タイトルと軸のラベルを設定

1
2
3
4

ax.set_title('Helicoid')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')

説明

• ax.set_title('Helicoid'): プロットにタイトル「Helicoid」を設定します。
• ax.set_xlabel('X-axis'): X軸にラベル「X-axis」を設定します。
• ax.set_ylabel('Y-axis'): Y軸にラベル「Y-axis」を設定します。
• ax.set_zlabel('Z-axis'): Z軸にラベル「Z-axis」を設定します。

9. プロットを表示

1

plt.show()

説明

• plt.show(): 作成した図を表示します。
  このコマンドが実行されると、Helicoidの3Dサーフェスプロットがウィンドウに表示されます。

以上が、コードの各部分の詳しい説明です。

このコードは、Helicoidという螺旋状の形状を持つ曲面を3Dで描画するためのものです。

グラフ解説

[実行結果]

グラフに表示される内容は以下の通りです：

1. グラフの形状

• Helicoid（ヘリコイド）の形状が描画されます。
  これは螺旋状の形状を持つ曲面であり、中心から外側に向かって広がる形をしています。

2. カラーマップ

• グラフはカラーマップ（viridis）を使用して色付けされています。
  このカラーマップは高さ（$z$軸の値）に応じて色が変わります。
  これにより、$z$軸方向の変化を視覚的に容易に捉えることができます。

3. カラーバー

• カラーバーがグラフの右側に表示され、色と高さの対応を示します。
  カラーバーは、サーフェスの高さがどの色に対応するかを示し、$z$軸の値の範囲を視覚化します。

4. タイトルと軸ラベル

• タイトル: グラフの上部に「Helicoid」というタイトルが表示されます。
  これにより、描画されている曲面がHelicoidであることがわかります。
• 軸ラベル: 各軸にラベルが付けられています。
  • X軸には「X-axis」
  • Y軸には「Y-axis」
  • Z軸には「Z-axis」

これにより、各軸が何を表しているかが明示されます。

グラフの表示

• 3Dサーフェスプロットとして、$x$, $y$, $z$の値から計算されたHelicoidの曲面が表示されます。
  曲面は、$u$と$v$という2つのパラメータに基づいて生成される点の集合から形成されています。

これらの要素が組み合わさって、Helicoidの3D曲面が視覚的にわかりやすく表示されています。

アニメーション３Dグラフ

Google Colab環境を使って、アニメーションする３Dグラフを描画します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

# Google Colabでの設定
from IPython.display import HTML
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation

# データの生成
x = np.linspace(-2 * np.pi, 2 * np.pi, 200)
y = np.linspace(-2 * np.pi, 2 * np.pi, 200)
x, y = np.meshgrid(x, y)
z = np.sin(np.sqrt(x**2 + y**2))

# プロットの設定
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 初期状態のプロット
line = ax.plot_surface(x, y, z, cmap='viridis')

# アニメーションの更新関数
def update(frame):
    global z, line
    ax.clear()
    z = np.sin(np.sqrt(x**2 + y**2) + frame / 10.0)
    line = ax.plot_surface(x, y, z, cmap='viridis')
    ax.set_zlim(-1, 1)
    return line,

# アニメーションの設定
ani = FuncAnimation(fig, update, frames=200, interval=50, blit=False)

# アニメーションをHTMLとして表示
HTML(ani.to_jshtml())

このコードでは、np.sqrt(x**2 + y**2)のようなより複雑な関数を使い、カラーマップviridisを適用しています。

これにより、アニメーションがカラフルで視覚的に魅力的になります。

【実行手順】

1. Google Colabを開きます。
2. 新しいノートブックを作成します。
3. 上記のコードをセルにコピーして貼り付けます。
4. セルを実行します。

アニメーション3Dグラフがカラフルで複雑な形状で表示されるはずです。

[実行結果]

ソースコード解説

ソースコードを詳しく説明します。

1. 必要なライブラリのインポート

1
2
3
4
5

from IPython.display import HTML
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation

• IPython.displayのHTML：Google ColabやJupyter NotebookでHTML形式の出力を表示するためのモジュール。
• numpy (np): 数値計算ライブラリ。数値データの生成と操作に使用。
• matplotlib.pyplot (plt): データの可視化ライブラリ。プロットの生成に使用。
• mpl_toolkits.mplot3dのAxes3D: 3Dプロットを生成するためのツールキット。
• matplotlib.animationのFuncAnimation: アニメーションを生成するためのモジュール。

2. データの生成

1
2
3
4

x = np.linspace(-2 * np.pi, 2 * np.pi, 200)
y = np.linspace(-2 * np.pi, 2 * np.pi, 200)
x, y = np.meshgrid(x, y)
z = np.sin(np.sqrt(x**2 + y**2))

• np.linspace(-2 * np.pi, 2 * np.pi, 200): -2πから2πまでの範囲を200等分した配列を生成。
• np.meshgrid(x, y): 2次元グリッドを生成。xとyの全ての組み合わせを含む。
• np.sqrt(x**2 + y**2): 各グリッドポイントのxとyの二乗和の平方根を計算。
• np.sin(...): 平方根の値をサイン関数に適用し、zの値を計算。

3. プロットの設定

1
2

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

• plt.figure(): 新しい図（figure）を作成。
• fig.add_subplot(111, projection='3d'): 3Dプロットを作成するための軸（axes）を追加。

4. 初期状態のプロット

1

line = ax.plot_surface(x, y, z, cmap='viridis')

• ax.plot_surface(x, y, z, cmap='viridis'): x, y, zのデータを用いて3Dサーフェスプロットを作成。
  cmap='viridis'はカラーマップを指定。

5. アニメーションの更新関数

1
2
3
4
5
6
7

def update(frame):
    global z, line
    ax.clear()
    z = np.sin(np.sqrt(x**2 + y**2) + frame / 10.0)
    line = ax.plot_surface(x, y, z, cmap='viridis')
    ax.set_zlim(-1, 1)
    return line,

• update(frame): アニメーションの各フレームで呼び出される関数。
  • ax.clear(): 現在のプロットをクリア。
  • z = np.sin(np.sqrt(x**2 + y**2) + frame / 10.0): zの値をフレームごとに更新。
  • line = ax.plot_surface(x, y, z, cmap='viridis'): 更新されたzの値で新しいサーフェスプロットを作成。
  • ax.set_zlim(-1, 1): z軸の範囲を設定。
  • return line,: 更新されたプロットを返す。

6. アニメーションの設定

1

ani = FuncAnimation(fig, update, frames=200, interval=50, blit=False)

• FuncAnimation(fig, update, frames=200, interval=50, blit=False): アニメーションを作成するための関数。
  • fig: アニメーションを表示する図。
  • update: フレームごとに呼び出される更新関数。
  • frames=200: フレームの数。
  • interval=50: フレーム間の時間間隔（ミリ秒）。
  • blit=False: フル再描画を行うかどうか。

7. アニメーションをHTMLとして表示

1

HTML(ani.to_jshtml())

• HTML(ani.to_jshtml()): アニメーションをHTML形式に変換し、Google ColabやJupyter Notebookで表示。

1. はじめに

当アプリは、ユーザーのプライバシーを重視しています。本プライバシーポリシーは、当アプリがユーザーの個人情報やその他のデータを一切収集しないことを明示するものです。

2. 収集する情報

当アプリは、ユーザーの個人情報やその他のデータを一切収集しません。

3. 情報の使用

当アプリは、ユーザーの個人情報やその他のデータを一切収集しないため、これらの情報を使用することはありません。

4. 情報の共有

当アプリは、ユーザーの個人情報やその他のデータを一切収集しないため、これらの情報を第三者と共有することはありません。

5. 情報の保護

当アプリは、ユーザーの個人情報やその他のデータを一切収集しません。

6. クッキー

当アプリは、クッキーを使用しません。

7. 子供のプライバシー

当アプリは、ユーザーの個人情報やその他のデータを一切収集しません。

8. プライバシーポリシーの変更

当アプリは、本プライバシーポリシーを随時更新することがあります。変更があった場合、当ページにてお知らせします。

9. お問い合わせ

プライバシーポリシーに関する質問や懸念がある場合は、以下の連絡先までお問い合わせください。

• メールアドレス: angularity@msn.com

混合整数計画法（Mixed-Integer Programming, MIP)

混合整数計画法（Mixed-Integer Programming, MIP）は、連続変数と整数変数の両方を含む最適化問題を解くための手法です。

以下に具体例を示し、それをPythonで解く方法を紹介します。

例題

ある工場では、2種類の製$品 (A) $と$ (B) $を生産しています。

工場の目標は、利益を最大化することです。以下の情報が与えられています。

• 製品$ (A) $の単位当たりの利益は$ 40 $ドル。
• 製品$ (B) $の単位当たりの利益は$ 30 $ドル。
• 製品$ (A) $の製造には$ 2 $時間の機械時間と$ 1 $単位の材料が必要。
• 製品$ (B) $の製造には$ 1 $時間の機械時間と$ 1 $単位の材料が必要。
• 使用可能な機械時間は$ 8 $時間。
• 使用可能な材料は$ 5 $単位。
• 製品$ (A) $と$ (B) $の生産数は整数でなければならない。

この問題を解くために、次のように定式化します。

数学的定式化

【最大化する目的関数】
　$\text{Maximize } 40x_A + 30x_B$

【制約条件】
　$2x_A + x_B \leq 8$ (機械時間の制約)
　$x_A + x_B \leq 5$ (材料の制約)
　$x_A, x_B \geq 0$
　$x_A, x_B \in \mathbb{Z} (整数条件)$

Pythonによる解法

ここでは、Pythonのcvxpyライブラリを使ってこの問題を解きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

import cvxpy as cp

# 変数の定義
x_A = cp.Variable(integer=True)
x_B = cp.Variable(integer=True)

# 目的関数の定義
objective = cp.Maximize(40 * x_A + 30 * x_B)

# 制約の定義
constraints = [
    2 * x_A + x_B <= 8,  # 機械時間の制約
    x_A + x_B <= 5,      # 材料の制約
    x_A >= 0,            # 非負の制約
    x_B >= 0             # 非負の制約
]

# 問題の定義
prob = cp.Problem(objective, constraints)

# 問題を解く
result = prob.solve()

# 結果の表示
print(f"Optimal value: {result}")
print(f"Optimal solution: x_A = {x_A.value}, x_B = {x_B.value}")

結果の解説

このコードを実行すると、以下のような結果が得られます。

1
2

Optimal value: 180.0
Optimal solution: x_A = 3.0, x_B = 2.0

この結果は、製品$ (A) $を$ 3 $単位、製品$ (B) $を$ 2 $単位生産することが利益を最大化する最適な解であることを示しています。

このときの最大利益は$ 180 $ドルです。

解説

1. 変数の定義:

1
2

x_A = cp.Variable(integer=True)
x_B = cp.Variable(integer=True)

ここでは、製品$ (A) $と$ (B) $の生産数を整数変数として定義します。

2. 目的関数の定義:

1

objective = cp.Maximize(40 * x_A + 30 * x_B)

目的関数は、製品$ (A) $と$ (B) $の生産による利益を最大化するものです。

3. 制約の定義:

1
2
3
4
5
6

constraints = [
   2 * x_A + x_B <= 8,  # 機械時間の制約
   x_A + x_B <= 5,      # 材料の制約
   x_A >= 0,            # 非負の制約
   x_B >= 0             # 非負の制約
]

ここでは、機械時間と材料の使用量に関する制約を定義しています。

4. 問題の定義と解法:

1
2

prob = cp.Problem(objective, constraints)
result = prob.solve()

問題を定義し、solveメソッドで解きます。

5. 結果の表示:

1
2

print(f"Optimal value: {result}")
print(f"Optimal solution: x_A = {x_A.value}, x_B = {x_B.value}")

最適な値と解を表示します。

このようにして、混合整数計画法を用いて最適な生産計画を立てることができます。

二次計画問題 (Quadratic Programming, QP)

凸最適化問題の具体例として、二次計画問題 (Quadratic Programming, QP) を挙げます。
この問題は、凸二次関数を最小化するもので、以下の形式を持ちます：

$$
\text{minimize} \quad \frac{1}{2}x^T Q x + c^T x
$$
$$
\text{subject to} \quad Ax \leq b
$$

ここで、$( Q ) $は対称かつ正定値行列、$( c ) $はベクトル、$( A ) $は行列、$( b ) $はベクトルです。

具体例：
$$ Q = \begin{bmatrix} 2 & 0 \ 0 & 2 \end{bmatrix}, \quad c = \begin{bmatrix} -4 \ -6 \end{bmatrix}, \quad A = \begin{bmatrix} 1 & 1 \ -1 & 2 \ 2 & 1 \end{bmatrix}, \quad b = \begin{bmatrix} 3 \ 2 \ 4 \end{bmatrix} $$

この例では、二次関数$ (\frac{1}{2}x^T Q x + c^T x) $を最小化し、線形制約$ (Ax \leq b) $を満たす$ (x) $を見つけます。

以下は、Pythonでこの問題を解くコードです。

cvxpyライブラリを使用して問題を解きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import cvxpy as cp
import numpy as np

# 定義されたパラメータ
Q = np.array([[2, 0], [0, 2]])
c = np.array([-4, -6])
A = np.array([[1, 1], [-1, 2], [2, 1]])
b = np.array([3, 2, 4])

# 変数の定義
x = cp.Variable(2)

# 目的関数の定義
objective = cp.Minimize(0.5 * cp.quad_form(x, Q) + c.T @ x)

# 制約の定義
constraints = [A @ x <= b]

# 問題の定義
prob = cp.Problem(objective, constraints)

# 問題を解く
result = prob.solve()

# 結果の表示
print(f"Optimal value: {result}")
print(f"Optimal solution: x = {x.value}")

このコードでは、次の手順を踏みます：

1. cvxpyライブラリをインポートします。
2. パラメータ$ (Q)$、$(c)$、$(A)$、および$ (b) $を定義します。
3. 最適化変数$ (x) $を定義します。
4. 目的関数を定義します。
5. 制約を定義します。
6. 問題を定義し、解きます。
7. 最適解と最適値を表示します。

実行すると、最適な解$ (x) $とそのときの最小値を得ることができます。

[実行結果]

Optimal value: -10.399999999999999
Optimal solution: x = [1.2 1.6]

ソースコード解説

コードの各部分を詳しく説明します。

1. 必要なライブラリのインポート

1
2

import cvxpy as cp
import numpy as np

最初に、最適化ライブラリであるCVXPYと数値計算ライブラリであるNumPyをインポートしています。
CVXPYは凸最適化問題を簡潔に定義し解くためのライブラリです。

2. パラメータの定義

1
2
3
4
5

# 定義されたパラメータ
Q = np.array([[2, 0], [0, 2]])
c = np.array([-4, -6])
A = np.array([[1, 1], [-1, 2], [2, 1]])
b = np.array([3, 2, 4])

ここでは、最適化問題に必要なパラメータを定義しています。

• Q: 二次形式の係数行列（正定値行列）
• c: 線形項の係数ベクトル
• A: 制約条件の係数行列
• b: 制約条件の右辺ベクトル

3. 変数の定義

1
2

# 変数の定義
x = cp.Variable(2)

次に、最適化変数 x を定義しています。
この変数は2次元のベクトルです。

4. 目的関数の定義

1
2

# 目的関数の定義
objective = cp.Minimize(0.5 * cp.quad_form(x, Q) + c.T @ x)

目的関数を定義しています。
この問題の目的は、以下の形式の関数を最小化することです：

$$
\frac{1}{2} x^T Q x + c^T x
$$

cp.quad_form(x, Q) は$ (x^T Q x) $を計算し、c.T @ x は$ (c^T x) $を計算します。
これらを足し合わせて目的関数を作成し、cp.Minimize を使って最小化します。

5. 制約の定義

1
2

# 制約の定義
constraints = [A @ x <= b]

制約条件を定義しています。
ここでは、以下の形式の線形不等式制約を設定しています：

$$
A x \leq b
$$

6. 問題の定義

1
2

# 問題の定義
prob = cp.Problem(objective, constraints)

最適化問題全体を定義しています。
cp.Problem 関数は、目的関数と制約条件を受け取り、最適化問題を作成します。

7. 問題を解く

1
2

# 問題を解く
result = prob.solve()

定義した最適化問題を解きます。
prob.solve() メソッドは、問題を解いて最適値を返します。

8. 結果の表示

1
2
3

# 結果の表示
print(f"Optimal value: {result}")
print(f"Optimal solution: x = {x.value}")

最適化の結果を表示します。

• result は最適値（目的関数の最小値）です。
• x.value は最適解（変数 x の値）です。

全体の流れ

1. ライブラリをインポートし、パラメータを定義。
2. 最適化変数を定義。
3. 目的関数を定義。
4. 制約条件を定義。
5. 目的関数と制約条件を組み合わせて最適化問題を定義。
6. 問題を解く。
7. 結果を表示する。

このコードは、二次計画問題を解くための一般的な方法を示しており、CVXPYライブラリを使うことで問題の定義と解決を簡潔に行えることを示しています。

結果解説

Optimal value: -10.399999999999999

最適値（目的関数の最小値）は$ (-10.4) $です。
これは、与えられた二次計画問題において、制約条件 $(Ax \leq b) $を満たす$ (x) $に対して目的関数が最小となる値です。

Optimal solution: x = [1.2, 1.6]

最適解$ (x = [1.2, 1.6]) $です。
これは、目的関数$ (\frac{1}{2}x^T Q x + c^T x) $が最小値をとる$ (x) $の値です。
この$ (x) $に対して、制約条件$ (Ax \leq b) $がすべて満たされます。

計算の詳細

目的関数の計算

与えられた$ (x = [1.2, 1.6]) $に対して目的関数を計算します。

$$
\frac{1}{2}x^T Q x + c^T x
$$

まず、二次項を計算します：
$$
x^T Q x = \begin{bmatrix} 1.2 & 1.6 \end{bmatrix} \begin{bmatrix} 2 & 0 \ 0 & 2 \end{bmatrix} \begin{bmatrix} 1.2 \ 1.6 \end{bmatrix} = 1.2 \times 2 \times 1.2 + 1.6 \times 2 \times 1.6 = 2.88 + 5.12 = 8.0
$$

次に、線形項を計算します：
$$
c^T x = \begin{bmatrix} -4 & -6 \end{bmatrix} \begin{bmatrix} 1.2 \ 1.6 \end{bmatrix} = -4 \times 1.2 + -6 \times 1.6 = -4.8 - 9.6 = -14.4
$$

したがって、目的関数の値は：
$$
\frac{1}{2} \times 8.0 + (-14.4) = 4.0 - 14.4 = -10.4
$$

これは、表示された最適値と一致します。

制約条件の確認

次に、与えられた最適解 $(x = [1.2, 1.6]) $が制約を満たすことを確認します：

$$
Ax = \begin{bmatrix} 1 & 1 \ -1 & 2 \ 2 & 1 \end{bmatrix} \begin{bmatrix} 1.2 \ 1.6 \end{bmatrix} = \begin{bmatrix} 1 \times 1.2 + 1 \times 1.6 \ -1 \times 1.2 + 2 \times 1.6 \ 2 \times 1.2 + 1 \times 1.6 \end{bmatrix} = \begin{bmatrix} 2.8 \ 2.0 \ 4.0 \end{bmatrix}
$$

これに対して、制約条件 $(b = [3, 2, 4]) $と比較します：
$$
\begin{bmatrix} 2.8 \ 2.0 \ 4.0 \end{bmatrix} \leq \begin{bmatrix} 3 \ 2 \ 4 \end{bmatrix}
$$

すべての不等式が満たされていることが確認できます。

結論

この問題に対する最適解は$ (x = [1.2, 1.6]) $であり、最適値は$ (-10.4) $です。
与えられた制約条件を満たしつつ、目的関数の値を最小化する解が求められました。

トラフィックエンジニアリング

問題:

ISP（インターネットサービスプロバイダー）は、特定の時間帯におけるネットワークトラフィックを効率的に管理し、ボトルネックを防ぎたいと考えています。
各リンクのトラフィック容量と予測されるトラフィック量が与えられたとき、トラフィックをどのように分散すれば全体の遅延を最小化できるかを求めてください。

入力:

• ノードとエッジのリスト。
• 各エッジのトラフィック容量。
• 各ノード間の予測トラフィック量。

出力:

• トラフィックの分配方法（各経路にどれだけのトラフィックを流すか）。

解法:

この問題は「多品種フロー問題」として知られ、線形計画法やネットワークシンプレックス法を用いて解決できます。

これらの例題は、通信網最適化の典型的な問題であり、現実のネットワーク設計や管理において非常に重要です。
各問題に対する解法は、適切なアルゴリズムとデータ構造を利用して実装できます。

例題

具体例として、簡単なトラフィックエンジニアリング問題を設定し、Pythonで解決し、グラフ化します。

具体例では、リンクのトラフィック容量と予測されるトラフィック量を考慮して、トラフィックを最適に分配します。

具体例

以下のようなネットワークがあるとします：

• ノード: A, B, C, D
• エッジとその容量:
  • A-B: $10$
  • A-C: $15$
  • B-D: $10$
  • C-D: $10$
  • B-C: $5$

予測されるトラフィック量は以下の通りです：

• AからDへのトラフィック: $10$
• BからCへのトラフィック: $5$

解法

この問題を線形計画法を使って解決します。

Pythonのscipy.optimize.linprogを使用します。

Pythonコード

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

1

pip install scipy networkx matplotlib

次に、トラフィックエンジニアリングの具体例を解くPythonコードを示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from scipy.optimize import linprog

# 定義されたエッジと容量
edges = [('A', 'B', 10), ('A', 'C', 15), ('B', 'D', 10), ('C', 'D', 10), ('B', 'C', 5)]
nodes = set(sum(([u, v] for u, v, _ in edges), []))
node_list = sorted(nodes)

# ノード間の予測トラフィック量
traffic = {('A', 'D'): 10, ('B', 'C'): 5}

# エッジの数とノードの数
num_edges = len(edges)
num_nodes = len(node_list)

# 線形計画法のセットアップ
c = np.ones(num_edges)  # 最小化するコスト関数の係数（単純に全て1）
A_eq = np.zeros((num_nodes, num_edges))
b_eq = np.zeros(num_nodes)

# 各エッジのインデックスを保持
edge_indices = {edge[:2]: idx for idx, edge in enumerate(edges)}

# ノードのインデックスを保持
node_indices = {node: idx for idx, node in enumerate(node_list)}

# 各エッジの制約をセットアップ
for (u, v, capacity) in edges:
    idx = edge_indices[(u, v)]
    A_eq[node_indices[u], idx] = 1
    A_eq[node_indices[v], idx] = -1

# トラフィックの制約をセットアップ
for (src, dst), demand in traffic.items():
    b_eq[node_indices[src]] += demand
    b_eq[node_indices[dst]] -= demand

# 各エッジの容量制約をセットアップ
bounds = [(0, capacity) for u, v, capacity in edges]

# 線形計画法を解く
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')

# 結果の表示
if res.success:
    print("Optimal traffic distribution:")
    for (u, v, capacity), flow in zip(edges, res.x):
        print(f"Flow on edge {u}-{v}: {flow:.2f}")
else:
    print("Optimization failed.")

# グラフの描画
G = nx.DiGraph()
for (u, v, capacity), flow in zip(edges, res.x):
    G.add_edge(u, v, capacity=capacity, flow=flow)

pos = nx.spring_layout(G)
edge_labels = {(u, v): f"{flow:.2f}/{capacity}" for (u, v, capacity), flow in zip(edges, res.x)}

plt.figure(figsize=(10, 8))
nx.draw(G, pos, with_labels=True, node_size=3000, node_color="lightblue", font_size=16, font_weight="bold", arrowsize=20)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=14)
nx.draw_networkx_edges(G, pos, width=2, edge_color='black', arrows=True)
plt.title("Optimized Traffic Distribution in the Network")
plt.show()

ソースコード解説

1. ライブラリのインポート: numpy, networkx, matplotlib, scipy.optimizeをインポートします。
2. ネットワークの定義: ノード、エッジ、容量、予測トラフィック量を定義します。
3. 線形計画法のセットアップ: コスト関数、制約行列、制約ベクトル、境界を設定します。
4. 線形計画法の解決: linprogを使用して問題を解きます。
5. 結果の表示: 各エッジに対する最適なトラフィック分配を表示します。
6. グラフの描画: networkxを使用してネットワークグラフを描画し、各エッジのフローを表示します。

このコードを実行すると、最適なトラフィック分配が計算され、その結果がグラフとして視覚化されます。

結果解説

表示される結果は、最適なトラフィック分配を示しています。

[実行結果]

以下に各エッジのフローについて説明します。

フロー結果

1. Flow on edge A-B: 5.00

   • AからBへのフローが$5.00$です。
     これは、BからCへのトラフィック量を処理するためにAからBに送られるトラフィックです。

2. Flow on edge A-C: 5.00

   • AからCへのフローが$5.00$です。
     これは、AからDへのトラフィック量の一部がAからCを経由してDに送られることを示しています。

3. Flow on edge B-D: 10.00

   • BからDへのフローが$10.00$です。
     これは、AからDへのトラフィック量がA-B-D経由で送られていることを示しています。
     つまり、AからBに送られたトラフィックのうち、$5.00$はBからDに直接送られています。
     また、B-Cへのトラフィック量も含まれています。

4. Flow on edge C-D: 0.00

   • CからDへのフローは$0.00$です。
     これは、AからDへのトラフィックがすべてA-B-D経由で送られているため、C-D経由でのトラフィックは発生していないことを示しています。

5. Flow on edge B-C: 0.00

   • BからCへのフローは$0.00$です。
     BからCへのトラフィック量$5.00$は、AからCを経由して送られているため、B-C間にはトラフィックが発生していません。

全体の流れ

• AからDへのトラフィック（10）:

  • AからBを経由して$5.00$が送られ、残りの$5.00$はAからCを経由して送られています。
  • BからDへのフローは$10.00$で、これはAからDへのトラフィック全体を運ぶために使われています。

• BからCへのトラフィック（5）:

  • AからBに送られた$5.00$がBからDに送られていますが、最終的にはAからCを経由してCに到達する形で処理されています。

この結果は、全体のネットワークトラフィックが最適に分配され、エッジの容量制限内でトラフィックが処理されていることを示しています。

二次計画法（Quadratic Programming, QP）

CVXPYを使用して二次計画法（Quadratic Programming, QP）の問題を解く方法を説明します。

二次計画法の一般的な形式は以下のようになります。

目的関数：
$$
\text{minimize} \quad \frac{1}{2} x^T P x + q^T x
$$

制約：
$$
Gx \leq h
$$
$$
Ax = b
$$

ここで、$(x)$は決定変数のベクトル、$(P)$は正定値対称行列、$(q)$はベクトル、$(G)$、$(h)$、$(A)$、および$(b)$はそれぞれ行列またはベクトルです。

プログラム例

以下はPythonのコードで、CVXPYを使用してこのような問題を解く方法の一例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

import cvxpy as cp
import numpy as np

# 定義する二次計画法の問題
P = np.array([[2, 0], [0, 2]])    # 目的関数の2次項の係数
q = np.array([-2, -5])            # 目的関数の1次項の係数
G = np.array([[1, 2], [1, -4], [-1, 0], [0, -1]]) # 不等式制約の係数
h = np.array([3, 3, 0, 0])        # 不等式制約の右辺
A = np.array([[1, 1]])            # 等式制約の係数
b = np.array([1])                 # 等式制約の右辺

# 変数を定義
x = cp.Variable(2)

# 目的関数を定義
objective = cp.Minimize((1/2) * cp.quad_form(x, P) + q @ x)

# 制約を定義
constraints = [G @ x <= h, A @ x == b]

# 問題を定義
prob = cp.Problem(objective, constraints)

# 問題を解く
result = prob.solve()

# 結果を出力
print(f"最適値: {result}")
print(f"最適解: {x.value}")

説明

1. 行列およびベクトルの定義:

   • P、q、G、h、A、b それぞれ問題の行列およびベクトルを定義します。

2. 変数の定義:

   • x = cp.Variable(2) で、最適化変数 $( x ) $を定義します。
     ここでは2次元の変数を考えています。

3. 目的関数の設定:

   • objective = cp.Minimize((1/2) * cp.quad_form(x, P) + q @ x) で、目的関数を定義しています。

4. 制約の設定:

   • constraints = [G @ x <= h, A @ x == b] で不等式制約と等式制約を定義しています。

5. 問題の解決:

   • prob = cp.Problem(objective, constraints) で問題を作成し、result = prob.solve() で解を求めます。

6. 結果の出力:

   • 最適値と最適解を出力します。

このようにして、CVXPYを使用して二次計画法の問題をシンプルに解くことができます。

CVXPYは非常に強力であり、もっと複雑な制約や目的関数にも対応することができます。

[実行結果]

最適値: -4.0
最適解: [1.78280298e-24 1.00000000e+00]

結果解説

CVXPYを使用した二次計画法の解の結果について説明します。

解の結果

• 最適値 (Objective Value): $-4.0$
• 最適解 (Optimal Solution): $[1.78280298e-24, 1.00000000e+00]$

1. 最適値 (Objective Value)

最適値: -4.0 という結果は、設定した二次計画法の目的関数を最小化したときの最小値が$ -4.0 $になることを示しています。
つまり、与えられた制約の範囲内で目的関数の評価が最も低くなる点（最適解）での値が$ -4.0 $です。

2. 最適解 (Optimal Solution)

最適解: [1.78280298e-24, 1.00000000e+00] という結果は、$( x ) $の最適値が$ [1.78280298e-24, 1.00000000e+00] $であることを示しています。
これが、与えられた目的関数を最小にする$ ( x ) $の値です。

結果の詳細

• 最適解$ ( x )$:
  • $( x_1 ) = 1.78280298e-24$
  • $( x_2 ) = 1.00000000e+00$

ここで、$( x_1 ) $が非常に小さい値（ほぼゼロ）になっていますが、これは最適解の計算過程での丸め誤差や数値計算の性質からくるものです。

実質的には$ ( x_1 \approx 0 ) $と解釈しても問題ありません。

制約の確認

この結果が制約を満たしていることを確認してみましょう。

【不等式制約】
$[ Gx \leq h ]$

ここで$ ( G ) $と$ ( h ) $は次のように定義されます。
$[ G = \begin{pmatrix} 1 & 2 \ 1 & -4 \ -1 & 0 \ 0 & -1 \end{pmatrix}, \quad h = \begin{pmatrix} 3 \ 3 \ 0 \ 0 \end{pmatrix} ]$

最適解$ ( x ) $が $[0, 1] $付近であることを考慮します。

• 第1の制約: $( 1 \cdot 0 + 2 \cdot 1 \leq 3 ) $ は$ ( 2 \leq 3 ) $で成立。
• 第2の制約: $( 1 \cdot 0 - 4 \cdot 1 \leq 3 ) $は$ ( -4 \leq 3 ) $で成立。
• 第3の制約: $( -1 \cdot 0 \leq 0 ) $は$ ( 0 \leq 0 ) $で成立。
• 第4の制約: $( 0 - 1 \cdot 1 \leq 0 ) $は$ ( -1 \leq 0 ) $で成立。

これらの不等式制約はすべて満たされています。

【等式制約】
$[ Ax = b ]$

ここで $( A ) $と$ ( b ) $は次のように定義されます。
$[ A = \begin{pmatrix} 1 & 1 \end{pmatrix}, \quad b = 1 ]$

• 制約: $( 1 \cdot 0 + 1 \cdot 1 = 1 )$

こちらの等式制約も満たされています。

まとめ

今回求めた最適解と最適値は、与えられた二次計画法の目的関数および制約条件に対して正しく計算されています。

最適解$ [1.78280298e-24, 1.00000000e+00] $が正しい最適解であり、これにより目的関数の最小値 $ -4.0 $が得られました。

• 最適値 (目的関数の最小値) は$ -4.0$。
• 最適解は$ ( x_1 \approx 0 )$、$( x_2 = 1 )$。

バーンズリーのシダ（Barnsley's fern）

バーンズリーのシダ（Barnsley’s fern）は、確率的な方法で生成されるフラクタル構造の一種であり、数学的な美しさと興味深さで知られています。

このシダは、数学者であるMichael Barnsleyによって導入されました。

バーンズリーのシダは、次のような確率的な方法で生成されます。

まず、初期条件として $(0, 0) $の点を選びます。

その後、一連の確率的な変換を適用して新しい点を生成します。

これらの変換は、特定の確率で選択され、シダの形を生成します。

バーンズリーのシダの生成には、以下のような変換が使用されます：

1. $1%$の確率で、点$ (x, y) $は$ (0, 0.16y) $に変換されます。
   これにより、葉の部分が生成されます。

2. $85%$の確率で、点$ (x, y) $は$ (0.85x + 0.04y, -0.04x + 0.85y + 1.6) $に変換されます。
   これにより、茎と枝が生成されます。

3. $7%$の確率で、点$ (x, y) $は$ (0.20x - 0.26y, 0.23x + 0.22y + 1.6) $に変換されます。
   これにより、さらに複雑な枝が生成されます。

4. $7%$の確率で、点$ (x, y) $は$ (-0.15x + 0.28y, 0.26x + 0.24y + 0.44) $に変換されます。
   これにより、より細かい枝が生成されます。

これらの変換を多数回繰り返すことで、バーンズリーのシダの構造が形成されます。

このフラクタル構造は、自然界に見られるシダの葉や枝の形状に類似しており、その美しさと複雑さから数学的な芸術の一形態として広く愛されています。

プログラム例

以下は、NetworkXを使用してバーンズリーのシダを生成し、グラフで表現するサンプルコードです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

import random
import networkx as nx
import matplotlib.pyplot as plt

# バーンズリーのシダ生成関数
def barnsley_fern(n):
    # 初期条件
    x, y = 0, 0
    points = [(x, y)]
    
    # バーンズリーのシダ生成
    for _ in range(n):
        r = random.random()
        if r <= 0.01:
            x, y = 0, 0.16 * y
        elif r <= 0.86:
            x, y = 0.85 * x + 0.04 * y, -0.04 * x + 0.85 * y + 1.6
        elif r <= 0.93:
            x, y = 0.20 * x - 0.26 * y, 0.23 * x + 0.22 * y + 1.6
        else:
            x, y = -0.15 * x + 0.28 * y, 0.26 * x + 0.24 * y + 0.44
        points.append((x, y))
    
    return points

# バーンズリーのシダを生成
fern_points = barnsley_fern(10000)

# グラフの作成
G = nx.Graph()

# ノードの追加
for i, (x, y) in enumerate(fern_points):
    G.add_node(i, pos=(x, y))

# エッジの追加
for i in range(len(fern_points) - 1):
    G.add_edge(i, i+1)

# グラフの描画
pos = nx.get_node_attributes(G, 'pos')
nx.draw(G, pos, node_size=0, edge_color='green', alpha=0.5, linewidths=0.5)
plt.title("Barnsley Fern using NetworkX")
plt.show()

このコードでは、バーンズリーのシダを生成するための関数 barnsley_fern を定義し、その結果をNetworkXのグラフに変換して描画しています。

生成されたグラフは、バーンズリーのシダの形を模倣した美しい構造を持っています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用してバーンズリーのシダ（Barnsley’s fern）を生成し、NetworkXを使用してそのグラフを描画するものです。

以下では、コードの構造と各部分の役割について詳しく説明します。

1. バーンズリーのシダ生成関数の定義

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

def barnsley_fern(n):
    x, y = 0, 0  # 初期条件
    points = [(x, y)]  # 座標を格納するリスト

    # バーンズリーのシダ生成
    for _ in range(n):
        r = random.random()  # 0から1の乱数を生成
        if r <= 0.01:
            x, y = 0, 0.16 * y
        elif r <= 0.86:
            x, y = 0.85 * x + 0.04 * y, -0.04 * x + 0.85 * y + 1.6
        elif r <= 0.93:
            x, y = 0.20 * x - 0.26 * y, 0.23 * x + 0.22 * y + 1.6
        else:
            x, y = -0.15 * x + 0.28 * y, 0.26 * x + 0.24 * y + 0.44
        points.append((x, y))  # 新しい座標をリストに追加

    return points

この部分では、バーンズリーのシダを生成する関数 barnsley_fern が定義されています。

この関数は、生成する点の数 n を引数として受け取ります。

バーンズリーのシダは、一連の確率的な変換を繰り返すことで生成されます。

各ステップで、ランダムに選択された変換に基づいて新しい座標が計算され、リスト points に追加されます。

2. バーンズリーのシダを生成

1

fern_points = barnsley_fern(10000)

この行では、barnsley_fern 関数を呼び出して、$10000$個の点からなるバーンズリーのシダを生成し、fern_points に格納しています。

3. グラフの作成

1

G = nx.Graph()

この行では、NetworkXのグラフオブジェクト G を作成しています。

4. ノードの追加

1
2

for i, (x, y) in enumerate(fern_points):
    G.add_node(i, pos=(x, y))

この部分では、バーンズリーのシダで生成された点をグラフのノードとして追加しています。

各ノードには、座標を属性として持つことができます。

ここでは、各ノードに座標情報を pos という名前の属性として追加しています。

5. エッジの追加

1
2

for i in range(len(fern_points) - 1):
    G.add_edge(i, i+1)

この部分では、連続する点の間にエッジを追加しています。

これにより、シダの枝の構造が表現されます。

6. グラフの描画

1
2
3
4

pos = nx.get_node_attributes(G, 'pos')
nx.draw(G, pos, node_size=0, edge_color='green', alpha=0.5, linewidths=0.5)
plt.title("Barnsley Fern using NetworkX")
plt.show()

最後に、NetworkXの draw 関数を使用してグラフを描画しています。

ノードの座標は pos から取得し、それに基づいてノードとエッジが描画されます。

結果として得られるグラフは、バーンズリーのシダの形状を美しく表現しています。