緯度経度間の距離 SciPy

緯度経度間の距離

地理関連の問題を解決する手法として、緯度経度間の距離を計算し、その結果をグラフ化することができます。

以下はその例です。

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
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

# 緯度経度のデータを定義
# 例: Tokyo, New York, London, Sydney
cities = {
'Tokyo': (35.6895, 139.6917),
'New York': (40.7128, -74.0060),
'London': (51.5074, -0.1278),
'Sydney': (-33.8688, 151.2093)
}

# 緯度経度をnumpy配列に変換
locations = list(cities.values())

# 距離を計算(ここではユークリッド距離を使用)
distances = cdist(locations, locations)

# 結果のグラフ化
plt.imshow(distances, cmap='viridis', interpolation='nearest')
plt.colorbar(label='Distance (degrees)')
plt.xticks(ticks=range(len(cities)), labels=cities.keys())
plt.yticks(ticks=range(len(cities)), labels=cities.keys())
plt.title('Distance between Cities')
plt.show()

このコードは、4つの都市(東京、ニューヨーク、ロンドン、シドニー)の緯度経度を使用し、それらの間の距離を計算し、結果をヒートマップとして可視化しています。

各都市間の距離が色で表され、より遠い場所ほど色が濃くなります。

こうした手法を用いて、地理的な位置関係を視覚的に表現することができます。

ソースコード解説

このコードは、PythonのSciPyライブラリを使用して、4つの都市の緯度経度からそれぞれの都市間の距離を計算し、その結果をヒートマップで可視化するものです。

詳細に見ていきましょう。

ライブラリのインポート

1
2
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

必要なライブラリをインポートしています。

  • cdist: SciPyの距離関数を使用して、配列内の要素間の距離を計算するために使用されます。
  • matplotlib.pyplot: グラフ描画のためのライブラリです。

都市の緯度経度の定義

1
2
3
4
5
6
cities = {
'Tokyo': (35.6895, 139.6917),
'New York': (40.7128, -74.0060),
'London': (51.5074, -0.1278),
'Sydney': (-33.8688, 151.2093)
}

4つの都市(東京、ニューヨーク、ロンドン、シドニー)の緯度経度データを辞書形式で定義しています。

緯度経度のnumpy配列への変換

1
locations = list(cities.values())

cities辞書から値(緯度経度)を取り出し、それをnumpy配列に変換しています。

都市間の距離の計算

1
distances = cdist(locations, locations)

cdist関数を使用して、各都市間の距離を計算しています。
ここではデフォルトでユークリッド距離が使用されています。

結果のグラフ化

1
2
3
4
5
6
plt.imshow(distances, cmap='viridis', interpolation='nearest')
plt.colorbar(label='Distance (degrees)')
plt.xticks(ticks=range(len(cities)), labels=cities.keys())
plt.yticks(ticks=range(len(cities)), labels=cities.keys())
plt.title('Distance between Cities')
plt.show()

imshow関数を使って都市間の距離をヒートマップとして表示しています。
colorbar関数でカラーバーを表示し、xticksおよびyticksでx軸とy軸の目盛りを設定し、都市名をラベルとして表示しています。
最後にtitle関数でグラフのタイトルを設定し、show関数でグラフを表示しています。

Rosenbrock関数 SciPy

Rosenbrock関数

非線形最適化問題の一例として、Rosenbrock関数の最小化問題を考えてみましょう。

Rosenbrock関数は、次のように定義されます:

$
f(x, y) = (a - x)^2 + b * (y - x^2)^2
$

ここで、$ a = 1, b = 100 $とします。

この関数の最小値は$ 0 $で、その位置は$ (a, a^2) = (1, 1) $です。

Pythonでは、Scipyの最適化関数を使用してこの最適化問題を解くことができます。

また、Matplotlibを使用して結果をグラフ化することも可能です。

以下に、PythonでRosenbrock関数の最小化問題を解くコードを示します。

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
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# Rosenbrock関数を定義
def rosenbrock(x):
return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# Rosenbrock関数の勾配を定義
def rosenbrock_der(x):
return np.array([-2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2), 200 * (x[1] - x[0]**2)])

# Rosenbrock関数のHessianを定義
def rosenbrock_hess(x):
return np.array([[2 + 1200 * x[0]**2 - 400 * x[1], -400 * x[0]],
[-400 * x[0], 200]])

# 初期値を設定
x0 = np.array([0.5, 0])

# 最適化問題を解く
res = minimize(rosenbrock, x0, method='trust-constr', jac=rosenbrock_der, hess=rosenbrock_hess)

# 結果を表示
print(res.x)

# 結果をグラフ化
x = np.linspace(-1, 2, 100)
y = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x, y)
Z = rosenbrock([X, Y])

plt.contour(X, Y, Z, levels=np.logspace(-2, 2, 20))
plt.xlabel('x')
plt.ylabel('y')
plt.title('Rosenbrock function')
plt.show()

このコードは、Rosenbrock関数の最小化問題を解くためのPythonのコードです。

まず、Rosenbrock関数とその勾配 Hessianを定義します。

次に、初期値を設定し、scipy.optimize.minimize関数を使用して最適化問題を解きます。

最後に、結果を表示し、結果をグラフ化します。

ソースコード解説

ソースコードを各部分を説明します。

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

1
2
3
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

必要なライブラリをインポートしています。

  • numpy: 数値計算用のライブラリ
  • scipy.optimize.minimize: 最適化問題を解くための関数
  • matplotlib.pyplot: グラフ描画のためのライブラリ

2. Rosenbrock関数の定義

1
2
def rosenbrock(x):
return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

Rosenbrock関数を定義しています。
この関数は最小化問題でよく使用される非線形関数で、$ (f(x, y) = (1 - x)^2 + 100(y - x^2)^2) $です。

3. Rosenbrock関数の勾配とHessianの定義

1
2
3
4
5
6
def rosenbrock_der(x):
return np.array([-2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2), 200 * (x[1] - x[0]**2)])

def rosenbrock_hess(x):
return np.array([[2 + 1200 * x[0]**2 - 400 * x[1], -400 * x[0]],
[-400 * x[0], 200]])

Rosenbrock関数の勾配とHessian(ヘシアン、2階偏微分行列)をそれぞれ定義しています。

4. 初期値の設定と最適化問題の解法

1
2
x0 = np.array([0.5, 0])
res = minimize(rosenbrock, x0, method='trust-constr', jac=rosenbrock_der, hess=rosenbrock_hess)

初期値を設定し、minimize関数を使ってRosenbrock関数の最小化問題を解いています。
method='trust-constr'は信頼領域法を用いて最適化を行います。
jac引数とhess引数はそれぞれ勾配とHessianの関数を指定しています。

5. 結果の表示

1
print(res.x)

最適化結果(最小値に達する$ x $の値)を出力しています。

6. 結果のグラフ化

1
2
3
4
5
6
7
8
9
10
x = np.linspace(-1, 2, 100)
y = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x, y)
Z = rosenbrock([X, Y])

plt.contour(X, Y, Z, levels=np.logspace(-2, 2, 20))
plt.xlabel('x')
plt.ylabel('y')
plt.title('Rosenbrock function')
plt.show()

Rosenbrock関数の等高線プロットを作成しています。
X軸とY軸の値を用意し、等高線の高さを関数値に対応させて関数を可視化しています。

結果解説

このグラフは、Rosenbrock関数の3次元表示を示しています。
Rosenbrock関数は、非線形最適化問題の一つで、最小値は$ 0 $で、その位置は$ (a, a^2) = (1, 1) $です。

グラフの$x軸$は、$x$の値を表し、$y軸$は$y$の値を表します。
Z軸は、Rosenbrock関数の値を表します。
この関数の最小値は$ 0 $で、その位置は$ (a, a^2) = (1, 1) $です。

グラフの中心部分は、関数の最小値を示しています。
この部分は、$ x=1 $と$ y=1 $の位置で、これが関数の最小値を示しています。

グラフの上部は、関数の値が大きい部分を示しています。
この部分は、$ x $や$ y $の値が大きいと、関数の値も大きくなります。

グラフの下部は、関数の値が小さい部分を示しています。
この部分は、$ x $や$ y $の値が小さいと、関数の値も小さくなります。

このグラフを見ると、Rosenbrock関数の最小値は、$ x=1 $と$ y=1 $の位置で、これが関数の最小値を示しています。

この位置は、$ x $や$ y $の値が最も小さいときに達します。
そのため、この位置は、最適化問題を解くための最良の解となります。

月と地球の運動 matplotlib

月と地球の運動

月と地球の運動をPythonを使って3Dで表現してみます。

このコードでは、地球が静止していると仮定し、月の動きのみを描画します。

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
68
69
70
71
72
73
74
75
76
77
78
79
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# 定数定義
G = 6.67430e-11 # 万有引力定数
m1 = 5.972e24 # 地球の質量(kg)
m2 = 7.35e22 # 月の質量(kg)
r = 384400000 # 地球から月までの距離(m)
v1 = 0 # 地球の初期速度(m/s)
v2 = 1022 # 月の初期速度(m/s)
timestep = 1000 # タイムステップ
num_steps = 10000 # ステップ数

# 初期条件
x1, y1, z1 = 0, 0, 0 # 地球の初期位置
x2, y2, z2 = r, 0, 0 # 月の初期位置

# 速度成分
vx1, vy1, vz1 = v1, 0, 0
vx2, vy2, vz2 = 0, v2, 0

# 位置、速度の更新
x1_values, y1_values, z1_values = [], [], []
x2_values, y2_values, z2_values = [], [], []

for _ in range(num_steps):
# 地球と月の距離
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
r = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)

# 万有引力の計算
F = G * m1 * m2 / r ** 2

# 重力の影響による加速度
ax1 = F * dx / r / m1
ay1 = F * dy / r / m1
az1 = F * dz / r / m1
ax2 = -F * dx / r / m2
ay2 = -F * dy / r / m2
az2 = -F * dz / r / m2

# 速度の更新
vx1 += ax1 * timestep
vy1 += ay1 * timestep
vz1 += az1 * timestep
vx2 += ax2 * timestep
vy2 += ay2 * timestep
vz2 += az2 * timestep

# 位置の更新
x1 += vx1 * timestep
y1 += vy1 * timestep
z1 += vz1 * timestep
x2 += vx2 * timestep
y2 += vy2 * timestep
z2 += vz2 * timestep

# 位置を保存
x1_values.append(x1)
y1_values.append(y1)
z1_values.append(z1)
x2_values.append(x2)
y2_values.append(y2)
z2_values.append(z2)

# 3D軌道のプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x2_values, y2_values, z2_values, label='Moon')
ax.plot([0], [0], [0], marker='o', markersize=8, color='blue', label='Earth')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Z coordinate (m)')
ax.set_title('Moon Orbit around Earth')
ax.legend()
plt.show()

このコードは、月の軌道を3Dで描画し、地球を静止して描画します。

月は地球の周りを楕円軌道で回る様子が視覚化されます。

ソースコード解説

コードの各部分を詳しく見ていきましょう:

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

  • import matplotlib.pyplot as plt: プロット用のMatplotlibライブラリをインポートします。
  • from mpl_toolkits.mplot3d import Axes3D: 3Dグラフを作成するためのAxes3Dをインポートします。
  • import numpy as np: 数値計算のためのNumPyライブラリをインポートします。

2. 定数の定義:

  • G = 6.67430e-11など、物理定数(万有引力定数など)や天体の質量、初期位置、速度などの値を定義します。

3. 初期条件:

  • x1, y1, z1など、地球と月の初期位置を設定します。

4. 速度成分:

  • vx1, vy1, vz1など、地球と月の初期速度を設定します。

5. 位置、速度の更新:

  • forループ内で、地球と月の位置および速度を時間ステップごとに計算して更新します。ニュートンの運動方程式と万有引力の法則に基づいて、加速度を計算し、それをもとに速度と位置を更新します。

6. 位置の保存:

  • 各ステップでの地球と月の位置を保存するためのリストを作成し、更新された位置をそれぞれのリストに追加します。

7. 3D軌道のプロット:

  • 地球と月の位置を元に、3Dグラフを作成します。
  • figとaxを使用して3Dグラフを作成し、plotメソッドで月の軌道を、plotメソッドとmarkerを使用して地球を示します。
  • 軸ラベル、タイトル、凡例を設定しています。

最後に、plt.show()で作成したグラフを表示します。

このコードは、地球と月の間の引力による運動をシミュレートし、月の楕円軌道を地球の周りにプロットしています。

結果解説

このグラフは、地球と月の間の万有引力を考慮した際の月の軌道を3Dで表現したものです。

青い点が地球を表し、月の軌道がその周りを回る様子が示されています。
地球は固定された位置にあり、月はその周りを楕円軌道で動いています。

月の軌道は、地球に引かれる重力によって決まります。
月が地球から遠ざかるときには引力が弱まり、近づくときには引力が強くなります。
その結果、月は楕円軌道を描きます。

この3Dグラフは、地球を中心に月が回転しながら地球の周りを周回する様子を視覚的に表現しています。
月の軌道が立体的に示されており、空間内での運動がわかりやすくなっています。

3D散布図/3Dサーフェスプロット Plotly

3D散布図/3Dサーフェスプロット

Plotlyを使って3Dグラフを描く方法はいくつかあります。

以下に、3D散布図と3Dサーフェスプロットの例を示します。

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

1
2
import plotly.graph_objects as go
import numpy as np

3D散布図

以下のコードは、3D散布図を描く例です。

ここでは、ヘリックスの方程式を使って3D空間に点を配置します。

1
2
3
4
5
6
# ヘリックスの方程式
t = np.linspace(0, 10, 50)
x, y, z = np.cos(t), np.sin(t), t

fig = go.Figure(data=[go.Scatter3d(x=x, y=y, z=z, mode='markers')])
fig.show()

3Dサーフェスプロット

以下のコードは、3Dサーフェスプロットを描く例です。

ここでは、3D空間にサーフェスを配置します。

1
2
3
4
5
6
x = np.outer(np.linspace(-2, 2, 30), np.ones(30))
y = x.copy().T
z = np.cos(x ** 2 + y ** 2)

fig = go.Figure(data=[go.Surface(x=x, y=y, z=z)])
fig.show()

これらのコードは、Plotlyを使って3D散布図と3Dサーフェスプロットを描きます。

ソースコード解説

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

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

1
2
import plotly.graph_objects as go
import numpy as np

3D散布図

次に、ヘリックスの方程式を使用して3D空間に点を配置します。

ヘリックスの方程式は、3次元空間でのヘリックスの形状を定義する数学的な式です。

ここでは、時間tに対してxとyの位置を計算し、zは時間t自体とします。

1
2
3
# ヘリックスの方程式
t = np.linspace(0, 10, 50)
x, y, z = np.cos(t), np.sin(t), t

これらの位置を使用して3D散布図を描画します。

go.Scatter3dは、3D空間での散布図を作成するためのPlotlyの関数です。

x=x, y=y, z=zは、散布図の座標を指定し、mode='markers'は、各点をマーカーとして表示することを指定します。

1
2
fig = go.Figure(data=[go.Scatter3d(x=x, y=y, z=z, mode='markers')])
fig.show()

3Dサーフェスプロット

次に、3Dサーフェスプロットを描画します。

3Dサーフェスプロットは、3次元空間でのサーフェスを表示するためのプロットです。

ここでは、3D空間でのサーフェスの形状を定義する数学的な式を使用します。

1
2
3
x = np.outer(np.linspace(-2, 2, 30), np.ones(30))
y = x.copy().T
z = np.cos(x ** 2 + y ** 2)

これらの位置を使用して3Dサーフェスプロットを描画します。

go.Surfaceは、3Dサーフェスプロットを作成するためのPlotlyの関数です。

x=x, y=y, z=zは、サーフェスの座標を指定します。

1
2
fig = go.Figure(data=[go.Surface(x=x, y=y, z=z)])
fig.show()

このコードは、Plotlyを使用して3D散布図と3Dサーフェスプロットを描画します。

カーネルメソッド matplotlib

カーネルメソッド

カーネルメソッドは、データを高次元空間にマッピングすることで、高次元空間上での計算を低次元空間上で行うことで計算効率を向上させる手法の一つです。

以下に、ガウスカーネルを用いたカーネルメソッドの具体的な例を示します。

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

1
2
import numpy as np
import matplotlib.pyplot as plt

次に、カーネル関数を定義します。
ここでは、ガウスカーネルを使用します。

1
2
def kernel(x, y):
return np.exp(-0.5*(x - y)**2 / 2)

次に、目的関数を定義します。

1
2
def f(x):
return 2 * x + 5 * np.cos(x)

次に、データを生成します。

1
2
3
4
5
6
7
np.random.seed(seed=4)
n = 13
m = 100
X = np.linspace(0, 10, m)
y_truth = f(X)
x = 0.5 + 8.5 * np.random.rand(n)
y = f(x)

次に、カーネルマトリックスを計算します。

1
K = np.array([[kernel(x[i], x[j]) for j in range(n)] for i in range(n)])

次に、カーネルメソッドによる予測を行います。

1
2
3
4
5
6
7
8
c = np.linalg.inv(K) @ y
y_pred = []
for i in range(m):
sum = 0
for j in range(n):
sum += c[j] * kernel(x[j], X[i])
y_pred.append(sum)
y_pred = np.array(y_pred)

最後に、予測結果をグラフに表示します。

1
2
3
4
5
6
fig, ax = plt.subplots()
ax.scatter(x, y, label="data")
ax.plot(X, y_pred, color="red", label="predict")
ax.plot(X, y_truth, label="ground truth")
ax.legend()
plt.show()

このコードは、ガウスカーネルを用いたカーネルメソッドの具体的な例を示しています。

結果解説

カーネルメソッドの例題で生成されるグラフは、データ点、予測値、真の値を表示したものです。

データ点:

青い点で表示されています。
これは、元のデータ点を表しています。

予測値:

赤い実線で表示されています。
これは、カーネルメソッドによる予測値を表しています。

真の値:

青い実線で表示されています。
これは、真の値(この場合は、カーネルメソッドが予測しようとしている関数の値)を表しています。

このグラフは、カーネルメソッドがどのようにデータを高次元空間にマッピングし、その高次元空間上での計算を低次元空間上で行っているかを視覚的に示しています。

また、カーネルメソッドがどの程度真の値を正確に予測できているかも示しています。

制約充足問題 python-constraint

制約充足問題

簡単な制約充足問題を解く例として、次の条件を満たす整数の組み合わせを見つける問題を考えましょう。

条件:

  1. 3つの整数$ ( x, y, z ) $を見つける。
  2. $ ( x, y, z ) $は 1 から 10 までの整数。
  3. $ ( x + y + z = 20 ) $の条件を満たす。

これをpython-constraintを使って解きます。

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
from constraint import Problem

def find_numbers():
problem = Problem()

# 変数の定義
problem.addVariable('x', range(1, 11))
problem.addVariable('y', range(1, 11))
problem.addVariable('z', range(1, 11))

# 制約の設定
def constraint_function(x, y, z):
if x + y + z == 20:
return True

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

solutions = problem.getSolutions()
return solutions

solutions = find_numbers()
if solutions:
print("整数の組み合わせ:")
for solution in solutions:
print(f"x: {solution['x']}, y: {solution['y']}, z: {solution['z']}")
else:
print("条件を満たす組み合わせは見つかりませんでした。")

このコードは、1から10までの整数$ (x, y, z) $の組み合わせを見つけて、条件$ (x + y + z = 20) $を満たす解を見つけます。

解が見つかれば、整数の組み合わせが出力されます。

[実行結果]

整数の組み合わせ:
x: 10, y: 9, z: 1
x: 10, y: 8, z: 2
x: 10, y: 7, z: 3
x: 10, y: 6, z: 4
x: 10, y: 5, z: 5
x: 10, y: 4, z: 6
x: 10, y: 3, z: 7
x: 10, y: 2, z: 8
x: 10, y: 1, z: 9
x: 9, y: 10, z: 1
x: 9, y: 9, z: 2
x: 9, y: 8, z: 3
x: 9, y: 7, z: 4
x: 9, y: 6, z: 5
x: 9, y: 5, z: 6
x: 9, y: 4, z: 7
x: 9, y: 3, z: 8
x: 9, y: 2, z: 9
x: 9, y: 1, z: 10
x: 8, y: 10, z: 2
x: 8, y: 9, z: 3
x: 8, y: 8, z: 4
x: 8, y: 7, z: 5
x: 8, y: 6, z: 6
x: 8, y: 5, z: 7
x: 8, y: 4, z: 8
x: 8, y: 3, z: 9
x: 8, y: 2, z: 10
x: 7, y: 10, z: 3
x: 7, y: 9, z: 4
x: 7, y: 8, z: 5
x: 7, y: 7, z: 6
x: 7, y: 6, z: 7
x: 7, y: 5, z: 8
x: 7, y: 4, z: 9
x: 7, y: 3, z: 10
x: 6, y: 10, z: 4
x: 6, y: 9, z: 5
x: 6, y: 8, z: 6
x: 6, y: 7, z: 7
x: 6, y: 6, z: 8
x: 6, y: 5, z: 9
x: 6, y: 4, z: 10
x: 5, y: 10, z: 5
x: 5, y: 9, z: 6
x: 5, y: 8, z: 7
x: 5, y: 7, z: 8
x: 5, y: 6, z: 9
x: 5, y: 5, z: 10
x: 4, y: 10, z: 6
x: 4, y: 9, z: 7
x: 4, y: 8, z: 8
x: 4, y: 7, z: 9
x: 4, y: 6, z: 10
x: 3, y: 10, z: 7
x: 3, y: 9, z: 8
x: 3, y: 8, z: 9
x: 3, y: 7, z: 10
x: 2, y: 10, z: 8
x: 2, y: 9, z: 9
x: 2, y: 8, z: 10
x: 1, y: 10, z: 9
x: 1, y: 9, z: 10

ソースコード解説

このコードは、Pythonのpython-constraintライブラリを使用して、特定の制約を満たす整数の組み合わせを見つける例です。

1. ライブラリのインポートと関数定義:

  • from constraint import Problem:constraintライブラリからProblemクラスをインポートします。
  • find_numbers()関数:制約充足問題を解くための関数を定義します。

2. 制約充足問題の定義:

  • problem = Problem():制約充足問題を定義するためのProblemオブジェクトを作成します。
  • problem.addVariable('x', range(1, 11)):変数 x を1から10の整数の範囲で定義します。
  • problem.addVariable('y', range(1, 11)):変数 y を1から10の整数の範囲で定義します。
  • problem.addVariable('z', range(1, 11)):変数 z を1から10の整数の範囲で定義します。

3. 制約の設定:

  • constraint_function()関数:x + y + z = 20 となるような条件を満たす制約関数を定義します。
  • problem.addConstraint(constraint_function, ['x', 'y', 'z']):constraint_functionを使用して、x + y + z = 20 の制約を追加します。

4. 問題の解決:

  • problem.getSolutions():Problemオブジェクトを解き、条件を満たす整数の組み合わせを探します。

5. 解の表示:

  • solutionsリストに解が含まれている場合は、各解の x, y, z の値を表示します。
  • 解が見つからなかった場合は、条件を満たす組み合わせが見つからなかったことを通知します。

線形最適化問題 CVXPY

線形最適化問題

CVXPYを使用して、簡単な線形最適化問題を解いてみましょう。

ここでは、2つの変数を持つ最小化の線形関数を扱います。例として、次のような問題を解いてみます。

$$
\text{minimize} \quad 3x + 4y
$$
$$
\text{subject to} \quad x + 2y \geq 10
$$
$$
\text{subject to} \quad 3x - y \leq 12
$$
$$
\text{subject to} \quad x, y \geq 0
$$

以下がそのコードです。

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
import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

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

# 目的関数と制約条件を定義
objective = cp.Minimize(3*x + 4*y)
constraints = [x + 2*y >= 10, 3*x - y <= 12, x >= 0, y >= 0]

# 最適化問題の定義
problem = cp.Problem(objective, constraints)

# 問題を解く
problem.solve()

# 最適な解の表示
print("Optimal value:", problem.value)
print("x =", x.value)
print("y =", y.value)

# グラフ描画
x_values = np.linspace(0, 10, 100)
y1_values = (10 - x_values) / 2
y2_values = 3 * x_values - 12

plt.figure(figsize=(6, 6))
plt.plot(x_values, y1_values, label='x + 2y >= 10')
plt.plot(x_values, y2_values, label='3x - y <= 12')
plt.fill_between(x_values, np.maximum(0, y1_values), np.minimum(10, (12 + y2_values) / 3), color='gray', alpha=0.3)
plt.plot(x.value, y.value, 'ro') # 最適解のプロット
plt.xlabel('x')
plt.ylabel('y')
plt.title('Optimization with CVXPY')
plt.legend()
plt.grid(True)
plt.show()

これにより、線形最適化問題が解かれ、最適解が計算されます。

また、制約条件と最適解が含まれたグラフも描画されます。

この場合、最適解は赤い点で示されます。

ソースコード解説

このコードは、CVXPYを使用して線形最適化問題を解決し、結果をグラフで視覚化するものです。

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

1
2
3
import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

必要なライブラリ(CVXPY、NumPy、Matplotlib)をインポートします。

2. 変数の定義

1
2
x = cp.Variable()
y = cp.Variable()

最適化問題で使用する変数 x と y を定義しています。

3. 目的関数と制約条件の設定

1
2
objective = cp.Minimize(3*x + 4*y)
constraints = [x + 2*y >= 10, 3*x - y <= 12, x >= 0, y >= 0]

最小化したい目的関数と制約条件を設定しています。

この場合、目的関数は$ (3x + 4y) $で、制約条件は不等式制約です。

4. 最適化問題の定義

1
problem = cp.Problem(objective, constraints)

CVXPYの Problem クラスを使って最適化問題を定義しています。

目的関数と制約条件を引数として渡しています。

5. 問題の解決

1
problem.solve()

定義した最適化問題を解いています。
CVXPYが問題を解くための適切な最適化アルゴリズムを選択し、最適な解を見つけようとします。

6. 最適な解の表示

1
2
3
print("Optimal value:", problem.value)
print("x =", x.value)
print("y =", y.value)

最適な目的関数の値とそれに対する最適な変数 x と y の値を表示しています。

7. グラフの描画

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
x_values = np.linspace(0, 10, 100)
y1_values = (10 - x_values) / 2
y2_values = 3 * x_values - 12

plt.figure(figsize=(6, 6))
plt.plot(x_values, y1_values, label='x + 2y >= 10')
plt.plot(x_values, y2_values, label='3x - y <= 12')
plt.fill_between(x_values, np.maximum(0, y1_values), np.minimum(10, (12 + y2_values) / 3), color='gray', alpha=0.3)
plt.plot(x.value, y.value, 'ro') # 最適解のプロット
plt.xlabel('x')
plt.ylabel('y')
plt.title('Optimization with CVXPY')
plt.legend()
plt.grid(True)
plt.show()

制約条件をグラフにプロットし、最適な解を赤い点として表示しています。

fill_between は制約条件の領域を色で塗りつぶし、plt.plot は最適解を赤い点でプロットしています。

結果解説

この結果は、与えられた線形最適化問題の最適な解を示しています。
最適な目的関数の値は約20です。

最適解は、$ (x \approx 0) $および$ (y \approx 5) $です。
つまり、制約条件を満たしながら最適な解が見つかりました。

グラフは、制約条件を表しています。
青い線はそれぞれ制約条件を示し、グレーの領域はすべての制約条件を満たす領域を示しています。
赤い点は最適解を表しており、制約条件の交差点に位置しています。
つまり、目的関数を最小化するための最適な$ (x) と (y) $の値を示しています。

最適解が$ (x \approx 0) $となる理由は、制約条件 $ (3x - y \leq 12) $により、$ (x) $が増加すると$ (y) $も増加する必要があるためです。
この条件を満たしつつ、目的関数を最小化するためには$ (x) $を極小化し$ (y) $を最大化するのが最適な戦略であり、その結果、$ (x \approx 0) $となります。

Bin Packing Problem(ビンパッキング問題) PuLP

Bin Packing Problem(ビンパッキング問題)

Bin Packing Problem(ビンパッキング問題)という問題をPythonで解く例を示します。

ビンパッキング問題とはビン(容器)に物を詰める最適化問題で、ここでは2つの容器に5つのアイテムを詰める例を示します。

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

1
pip install pulp

以下の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
from pulp import *

# ビンの数とアイテムの数を設定
bin_count = 2
item_count = 5

# 各ビンの容量、アイテムの重さを設定
bin_capacity = [10, 10]
item_weight = [2, 5, 4, 7, 1]

# 問題設定
prob = LpProblem("BinPacking", LpMinimize)

# アイテムをビンに入れるかどうかを示すノードを作成(バイナリ変数)
x = [[LpVariable("x(%s,%s)" % (i, j), cat="Binary") for j in range(item_count)] for i in range(bin_count)]

# 目的関数: 使用するビンの数を最小化する
prob += lpSum(x[i][j] for i in range(bin_count) for j in range(item_count))

# 制約条件: 各アイテムは1つのビンにしか入らない
for j in range(item_count):
prob += lpSum(x[i][j] for i in range(bin_count)) == 1

# 制約条件: 各ビンの容量を超えてアイテムを入れられない
for i in range(bin_count):
prob += lpSum(item_weight[j] * x[i][j] for j in range(item_count)) <= bin_capacity[i]

# 最適化計算
prob.solve()

# 結果表示
for i in range(bin_count):
print("*****Bin", i, "*****")
total_weight = 0
for j in range(item_count):
if value(x[i][j])==1:
print("Item", j)
total_weight+=item_weight[j]
print("Total weight:", total_weight)

このコードでは2つのビンに5つのアイテムを最適に配置しています。

出力結果は各ビンに収められたアイテムとその合計重量を表示します。

[実行結果]

*****Bin 0 *****
Item 1
Item 2
Total weight: 9
*****Bin 1 *****
Item 0
Item 3
Item 4
Total weight: 10

ソースコード解説

このコードは、Bin Packing Problem(ビンパッキング問題)を解くためにPuLP(Pythonの最適化モデリングライブラリ)を使用しています。

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

1
from pulp import *

PuLPからすべての関数をインポートしています。

2. 問題のパラメータ設定

1
2
3
4
bin_count = 2
item_count = 5
bin_capacity = [10, 10]
item_weight = [2, 5, 4, 7, 1]

ビンの数やアイテムの数、各ビンの容量、アイテムの重さを設定しています。

3. 問題の設定

1
prob = LpProblem("BinPacking", LpMinimize)

PuLPを使って「BinPacking」という名前の最小化問題を定義しています。

4. 変数の定義

1
x = [[LpVariable("x(%s,%s)" % (i, j), cat="Binary") for j in range(item_count)] for i in range(bin_count)]

ビンにアイテムを入れるかどうかを示す変数 x をバイナリ変数として定義しています。

5. 目的関数の設定

1
prob += lpSum(x[i][j] for i in range(bin_count) for j in range(item_count))

使用するビンの数を最小化する目的関数を定義しています。

6. 制約条件の設定

1
2
3
4
5
for j in range(item_count):
prob += lpSum(x[i][j] for i in range(bin_count)) == 1

for i in range(bin_count):
prob += lpSum(item_weight[j] * x[i][j] for j in range(item_count)) <= bin_capacity[i]

各アイテムが1つのビンにしか入らない制約条件と、各ビンの容量を超えない制約条件を設定しています。

7. 最適化計算

1
prob.solve()

PuLPを使って最適化問題を解いています。

8. 結果表示

1
2
3
4
5
6
7
8
for i in range(bin_count):
print("*****Bin", i, "*****")
total_weight = 0
for j in range(item_count):
if value(x[i][j])==1:
print("Item", j)
total_weight+=item_weight[j]
print("Total weight:", total_weight)

最適化計算の結果を表示しており、各ビンに入ったアイテムとその総重量を出力しています。

株価予測 Prophet

株価予測

株価の予測を行ってみましょう。

Prophetを使用して株価の未来予測を行うことが可能です。

以下は、仮想の株価データを用いてProphetを使った株価予測の例です。

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
import pandas as pd
import numpy as np
from prophet import Prophet
import matplotlib.pyplot as plt

# ランダムな株価データを生成
np.random.seed(42)
dates = pd.date_range(start='2023-01-01', periods=365, freq='D')
prices = np.cumsum(np.random.randn(365)) + 100 # ランダムな価格変動

data = {
'ds': dates,
'y': prices
}

df = pd.DataFrame(data)

# Prophetモデルの作成と学習
model = Prophet()
model.fit(df)

# 予測用データの作成
future = model.make_future_dataframe(periods=30) # 30日間の予測
forecast = model.predict(future)

# 予測結果をグラフ化
fig = model.plot(forecast)
plt.xlabel('Date')
plt.ylabel('Stock Price')
plt.title('Random Stock Price Forecast')
plt.show()

この例では、乱数を使って生成したランダムな価格変動を持つ株価データを作成し、Prophetモデルを使って未来30日間の株価を予測しています。

生成されたランダムなデータを用いているため、実際の株価とは関係ありませんが、Prophetがこのようなデータでも時系列の予測を行うことができることを示しています。

ソースコード解説

このソースコードは、PythonのPandas、NumPy、Prophet(Facebookが提供する時系列データの予測ライブラリ)、Matplotlibを使用して、ランダムな株価データを生成し、Prophetモデルを使って未来の株価を予測し、その結果をグラフ化しています。

1. ランダムな株価データの生成

1
2
3
4
5
6
7
8
9
10
np.random.seed(42)
dates = pd.date_range(start='2023-01-01', periods=365, freq='D')
prices = np.cumsum(np.random.randn(365)) + 100 # ランダムな価格変動

data = {
'ds': dates,
'y': prices
}

df = pd.DataFrame(data)

この部分では、np.random.seed(42)で乱数のシードを設定し、365日分の日付とランダムな株価のデータを生成しています。
これらのデータを'ds'(日付)と'y'(株価)のカラムを持つPandasのDataFrameに格納しています。

2. Prophetモデルの作成と学習

1
2
model = Prophet()
model.fit(df)

Prophetのモデルを作成し、fit()メソッドを使ってデータを学習させています。
Prophetは時系列データの予測に特化したモデルで、ここでは与えられた株価データを元にモデルを構築しています。

3. 予測用データの作成と予測

1
2
future = model.make_future_dataframe(periods=30)  # 30日間の予測
forecast = model.predict(future)

make_future_dataframe()メソッドを使って、元のデータから未来の30日間の日付を含むDataFrameを作成しています。
それを使ってpredict()メソッドで予測を行い、forecastに未来の予測結果を格納しています。

4. 予測結果をグラフ化

1
2
3
4
5
fig = model.plot(forecast)
plt.xlabel('Date')
plt.ylabel('Stock Price')
plt.title('Random Stock Price Forecast')
plt.show()

最後の部分では、plot()メソッドを使って予測結果をグラフ化しています。
横軸には日付、縦軸には株価が表示され、予測されたデータに対する実際のデータと予測の境界が示されています。
これにより、モデルが学習した過去のデータに対する予測と未来の予測の推移が視覚化されます。

結果解説

このグラフは、Prophetによる株価の予測を示しています。
横軸は日付を示し、縦軸は株価を表しています。

以下はグラフの要素についての詳細です。

1. 黒色の点:

グラフ上の黒色の点は、実際の株価データを表しています。
過去の実際のデータポイントがこの線上にプロットされています。

2. 青色の線:

青色の線は、過去の実際のデータポイントに対するProphetモデルの適合を表しています。
点線が黒色の点に密着しているほど、モデルが過去のデータをよく予測できていることを示します。

3. 青色の領域:

青色の線の周囲には、淡い青色の領域があります。
これは予測区間を表しており、将来の株価の変動がこの領域内に含まれる可能性を示しています。
通常、時間が経つにつれてこの領域が広がる傾向があります。

このグラフ全体を通じて、過去の実際のデータに対するモデルの適合度や未来の予測値が視覚的に表現されています。

モデルが過去のデータに適合し、予測区間が示されていることで、将来の株価の予測に対する確信度や不確実性を示すことができます。

プロジェクトスケジューリング問題 NetworkX

プロジェクトスケジューリング問題

プロジェクトスケジューリング問題は、タスクの完了時間やリソースの最適な割り当てを計画する問題です。

具体的な例として、4つのタスクとそれらの間の依存関係があるとしましょう。

  • Task A: 5日かかる
  • Task B: 3日かかる、Task A完了後に開始可能
  • Task C: 2日かかる、Task A完了後に開始可能
  • Task D: 4日かかる、Task BとTask C完了後に開始可能

これをPythonで解決してみましょう。

ここではnetworkxライブラリを使用してグラフを作成し、トポロジカルソートを行います。

まず、networkxをインストールします。

1
!pip install networkx

次に、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
import networkx as nx
import matplotlib.pyplot as plt

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

# タスクの追加
tasks = {
'A': {'duration': 5},
'B': {'duration': 3},
'C': {'duration': 2},
'D': {'duration': 4}
}

# タスク間の依存関係を追加
dependencies = [('A', 'B'), ('A', 'C'), ('B', 'D'), ('C', 'D')]

# グラフにタスクと依存関係を追加
for task, data in tasks.items():
G.add_node(task, duration=data['duration'])

for dependency in dependencies:
G.add_edge(dependency[0], dependency[1])

# トポロジカルソート
schedule = list(nx.topological_sort(G))

# タスクのスケジュール表示
print("スケジュール:")
for i, task in enumerate(schedule):
print(f"Day {i + 1}: Task {task}")

# グラフの描画
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels=True, node_size=800, node_color='skyblue', font_weight='bold', font_size=12, edge_color='gray', arrowsize=20)
plt.title('Project Scheduling')
plt.show()

このコードは、依存関係を持つタスクをグラフで表現し、トポロジカルソートを使ってスケジュールを作成します。

そして、タスクのスケジュールを表示し、グラフを描画します。

ソースコード解説

このコードは、Pythonのnetworkxとmatplotlibライブラリを使用してプロジェクトスケジュールを計画し、グラフで可視化するものです。

以下にコードの詳細を説明します:

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

1
2
import networkx as nx
import matplotlib.pyplot as plt
  • networkxはグラフの作成や操作、アルゴリズムの実行などをサポートするライブラリです。
  • matplotlib.pyplotはグラフを描画するためのライブラリです。

2. グラフの作成

1
G = nx.DiGraph()
  • DiGraph()関数を使って有向グラフを作成しています。

3. タスクの定義

1
2
3
4
5
6
tasks = {
'A': {'duration': 5},
'B': {'duration': 3},
'C': {'duration': 2},
'D': {'duration': 4}
}
  • タスクとその所要時間を辞書形式で定義しています。

4. タスク間の依存関係の定義

1
dependencies = [('A', 'B'), ('A', 'C'), ('B', 'D'), ('C', 'D')]
  • タスク間の依存関係をタプルのリストとして定義しています。

5. グラフにタスクと依存関係を追加

1
2
3
4
5
for task, data in tasks.items():
G.add_node(task, duration=data['duration'])

for dependency in dependencies:
G.add_edge(dependency[0], dependency[1])
  • タスクとその所要時間をグラフにノードとして追加し、タスク間の依存関係をエッジとして追加しています。

6. トポロジカルソート

1
schedule = list(nx.topological_sort(G))
  • topological_sort()関数を使って、トポロジカルソートを行い、タスクのスケジュールを計画します。

7. タスクのスケジュール表示

1
2
3
print("スケジュール:")
for i, task in enumerate(schedule):
print(f"Day {i + 1}: Task {task}")
  • タスクのスケジュールをコンソールに出力します。

8. グラフの描画

1
2
3
4
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels=True, node_size=800, node_color='skyblue', font_weight='bold', font_size=12, edge_color='gray', arrowsize=20)
plt.title('Project Scheduling')
plt.show()
  • spring_layout()関数でノードの位置を計算し、draw()関数でグラフを描画しています。
  • 各ノードはタスクを表し、エッジはタスク間の依存関係を示しています。

これにより、プロジェクトのタスクとその依存関係をグラフで視覚化し、トポロジカルソートによってスケジュールを立てています。

結果解説

この結果は、与えられたタスクとその依存関係に基づいて計画されたスケジュールを示しています。

それぞれのタスクが実行される日数と順序は以下のようになっています:

  • Day 1: Task A
  • Day 2: Task B
  • Day 3: Task C
  • Day 4: Task D

それぞれのタスクの日数に加え、依存関係がスケジュールに影響を与えています。
Task Aは最初に実行され、その後、Task BとTask Cがそれぞれ独立して実行されます。
最後に、Task BとTask Cの完了後にTask Dが実行されます。

グラフを見ると、各タスクがノードとして表され、依存関係が矢印で示されています。
このグラフは、タスク間の依存関係を視覚化し、トポロジカルソートによってタスクの実行順序が決定されました。
それにより、各タスクの開始と終了時期が特定されました。