Skip to content

Commit 35dc965

Browse files
bartzbeielsteinbartzbeielstein
committed
tests for Kriging
1 parent 9e67062 commit 35dc965

5 files changed

Lines changed: 267 additions & 35 deletions

File tree

‎notebooks/testKriging.ipynb‎

Lines changed: 190 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Test Functions for Kriging"
8+
]
9+
},
10+
{
11+
"cell_type": "code",
12+
"execution_count": null,
13+
"metadata": {},
14+
"outputs": [],
15+
"source": [
16+
"from spotPython.build.kriging import Kriging\n",
17+
"import numpy as np\n",
18+
"import matplotlib.pyplot as plt\n",
19+
"from numpy import linspace, arange\n",
20+
"rng = np.random.RandomState(1)\n",
21+
"X = linspace(start=0, stop=10, num=1_000).reshape(-1, 1)\n",
22+
"y = np.squeeze(X * np.sin(X))\n",
23+
"training_indices = rng.choice(arange(y.size), size=6, replace=False)\n",
24+
"X_train, y_train = X[training_indices], y[training_indices]\n",
25+
"S = Kriging(name='kriging', seed=124)\n",
26+
"S.fit(X_train, y_train)\n",
27+
"mean_prediction, std_prediction, s_ei = S.predict(X, return_val=\"all\")\n",
28+
"plt.plot(X, y, label=r\"$f(x)$\", linestyle=\"dotted\")\n",
29+
"plt.scatter(X_train, y_train, label=\"Observations\")\n",
30+
"plt.plot(X, mean_prediction, label=\"Mean prediction\")\n",
31+
"plt.fill_between(\n",
32+
" X.ravel(),\n",
33+
" mean_prediction - 1.96 * std_prediction,\n",
34+
" mean_prediction + 1.96 * std_prediction,\n",
35+
" alpha=0.5,\n",
36+
" label=r\"95% confidence interval\",\n",
37+
" )\n",
38+
"plt.legend()\n",
39+
"plt.xlabel(\"$x$\")\n",
40+
"plt.ylabel(\"$f(x)$\")\n",
41+
"_ = plt.title(\"Gaussian process regression on noise-free dataset\")\n",
42+
"plt.show()"
43+
]
44+
},
45+
{
46+
"cell_type": "code",
47+
"execution_count": null,
48+
"metadata": {},
49+
"outputs": [],
50+
"source": [
51+
"from spotPython.build.kriging import Kriging\n",
52+
"import numpy as np\n",
53+
"import matplotlib.pyplot as plt\n",
54+
"from numpy import linspace, arange\n",
55+
"rng = np.random.RandomState(1)\n",
56+
"X = linspace(start=0, stop=10, num=1_0).reshape(-1, 1)\n",
57+
"y = np.squeeze(X * np.sin(X))\n",
58+
"training_indices = rng.choice(arange(y.size), size=6, replace=False)\n",
59+
"X_train, y_train = X[training_indices], y[training_indices]\n",
60+
"S = Kriging(name='kriging', seed=124)\n",
61+
"S.fit(X_train, y_train)\n",
62+
"mean_prediction, std_prediction, s_ei = S.predict(X, return_val=\"all\")\n",
63+
"# Kriging is a interpolator, so the mean prediction should be equal to the training points:\n",
64+
"# check if the difference between the mean prediction and the true value in the training points is smaller than 1e-6\n",
65+
"assert np.allclose(mean_prediction[training_indices], y[training_indices], atol=1e-6)\n"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"## Expected Improvement"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {},
79+
"outputs": [],
80+
"source": [
81+
"from spotPython.build.kriging import Kriging\n",
82+
"from math import erf\n",
83+
"import numpy as np\n",
84+
"S = Kriging(name='kriging', seed=124)\n",
85+
"S.mean_cod_y = [0.0, 0.0, 0.0, 0.0, 0.0]\n",
86+
"# asset that the S.exp_imp(1.0, 0.0) is equal to 0.0\n",
87+
"assert 0.0 == S.exp_imp(1.0, 0.0)\n",
88+
"# assert that the S.exp_imp(0.0, 1.0) is equal to 1/sqrt(2 pi)\n",
89+
"# assert S.exp_imp(0.0, 1.0) == 1/np.sqrt(2*np.pi)\n",
90+
"# play safe and use np.allclose\n",
91+
"assert np.allclose(S.exp_imp(0.0, 1.0), 1/np.sqrt(2*np.pi), atol=1e-6)\n",
92+
"assert np.allclose(S.exp_imp(1.0, 1.0), -(0.5 + 0.5*erf(-1/np.sqrt(2))) + 1/np.sqrt(2*np.pi)*np.exp(-1/2), atol=1e-6)"
93+
]
94+
},
95+
{
96+
"cell_type": "markdown",
97+
"metadata": {},
98+
"source": [
99+
"# set_de_bounds"
100+
]
101+
},
102+
{
103+
"cell_type": "code",
104+
"execution_count": null,
105+
"metadata": {},
106+
"outputs": [],
107+
"source": [
108+
"from spotPython.build.kriging import Kriging\n",
109+
"S = Kriging(name='kriging', seed=124)\n",
110+
"S.set_de_bounds()\n",
111+
"assert S.de_bounds == [[-3, 2]]\n",
112+
"from spotPython.build.kriging import Kriging\n",
113+
"n = 10\n",
114+
"S = Kriging(name='kriging', seed=124, n_theta=n)\n",
115+
"S.set_de_bounds()\n",
116+
"assert len(S.de_bounds) == n\n",
117+
"n=2\n",
118+
"p=4\n",
119+
"S = Kriging(name='kriging', seed=124, n_theta=n, n_p=p, optim_p=True)\n",
120+
"S.set_de_bounds()\n",
121+
"assert len(S.de_bounds) == n+p\n",
122+
"S = Kriging(name='kriging', seed=124, n_theta=n, n_p=p, optim_p=False)\n",
123+
"S.set_de_bounds()\n",
124+
"assert len(S.de_bounds) == n"
125+
]
126+
},
127+
{
128+
"cell_type": "markdown",
129+
"metadata": {},
130+
"source": [
131+
"## extract_from_bounds"
132+
]
133+
},
134+
{
135+
"cell_type": "code",
136+
"execution_count": 2,
137+
"metadata": {},
138+
"outputs": [
139+
{
140+
"name": "stdout",
141+
"output_type": "stream",
142+
"text": [
143+
"[1 2]\n",
144+
"[3]\n"
145+
]
146+
}
147+
],
148+
"source": [
149+
"import numpy as np\n",
150+
"from spotPython.build.kriging import Kriging\n",
151+
"n=2\n",
152+
"p=4\n",
153+
"S = Kriging(name='kriging', seed=124, n_theta=n, n_p=p, optim_p=True, noise=False)\n",
154+
"S.extract_from_bounds(np.array([1, 2, 3]))\n",
155+
"print(S.theta)\n",
156+
"print(S.p)\n",
157+
"\n",
158+
"\n"
159+
]
160+
},
161+
{
162+
"cell_type": "code",
163+
"execution_count": null,
164+
"metadata": {},
165+
"outputs": [],
166+
"source": []
167+
}
168+
],
169+
"metadata": {
170+
"kernelspec": {
171+
"display_name": "spotCondaEnv",
172+
"language": "python",
173+
"name": "python3"
174+
},
175+
"language_info": {
176+
"codemirror_mode": {
177+
"name": "ipython",
178+
"version": 3
179+
},
180+
"file_extension": ".py",
181+
"mimetype": "text/x-python",
182+
"name": "python",
183+
"nbconvert_exporter": "python",
184+
"pygments_lexer": "ipython3",
185+
"version": "3.11.6"
186+
}
187+
},
188+
"nbformat": 4,
189+
"nbformat_minor": 2
190+
}

‎pyproject.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
77

88
[project]
99
name = "spotPython"
10-
version = "0.7.6"
10+
version = "0.7.7"
1111
authors = [
1212
{ name="T. Bartz-Beielstein", email="tbb@bartzundbartz.de" }
1313
]

‎src/spotPython/build/kriging.py‎

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -127,19 +127,18 @@ def __init__(
127127
counter : Counter.
128128
129129
Examples:
130-
Surrogate of the x*sin(x) function, see [1].
131-
132130
>>> from spotPython.build.kriging import Kriging
133131
import numpy as np
134132
import matplotlib.pyplot as plt
133+
from numpy import linspace, arange
135134
rng = np.random.RandomState(1)
136135
X = linspace(start=0, stop=10, num=1_000).reshape(-1, 1)
137136
y = np.squeeze(X * np.sin(X))
138137
training_indices = rng.choice(arange(y.size), size=6, replace=False)
139138
X_train, y_train = X[training_indices], y[training_indices]
140139
S = Kriging(name='kriging', seed=124)
141140
S.fit(X_train, y_train)
142-
mean_prediction, std_prediction = S.predict(X)
141+
mean_prediction, std_prediction, s_ei = S.predict(X, return_val="all")
143142
plt.plot(X, y, label=r"$f(x)$", linestyle="dotted")
144143
plt.scatter(X_train, y_train, label="Observations")
145144
plt.plot(X, mean_prediction, label="Mean prediction")
@@ -225,14 +224,18 @@ def exp_imp(self, y0: float, s0: float) -> float:
225224
float: The expected improvement value.
226225
227226
Examples:
228-
229227
>>> from spotPython.build.kriging import Kriging
230-
>>> S = Kriging(name='kriging', seed=124)
231-
>>> S.cod_y = [0.0, 0.0, 0.0, 0.0, 0.0]
232-
>>> S.mean_cod_y = [0.0, 0.0, 0.0, 0.0, 0.0]
233-
>>> S.exp_imp(1.0, 2.0)
234-
0.0
235-
228+
S = Kriging(name='kriging', seed=124)
229+
S.mean_cod_y = [0.0, 0.0, 0.0, 0.0, 0.0]
230+
S.exp_imp(1.0, 0.0)
231+
0.0
232+
>>> from spotPython.build.kriging import Kriging
233+
S = Kriging(name='kriging', seed=124)
234+
S.mean_cod_y = [0.0, 0.0, 0.0, 0.0, 0.0]
235+
# assert S.exp_imp(0.0, 1.0) == 1/np.sqrt(2*np.pi)
236+
# which is approx. 0.3989422804014327
237+
S.exp_imp(0.0, 1.0)
238+
0.3989422804014327
236239
"""
237240
# y_min = min(self.cod_y)
238241
y_min = min(self.mean_cod_y)
@@ -262,13 +265,10 @@ def set_de_bounds(self) -> None:
262265
self (object): The Kriging object.
263266
264267
Examples:
265-
266268
>>> from spotPython.build.kriging import Kriging
267-
>>> MyClass = Kriging(name='kriging', seed=124)
268-
>>> obj = MyClass()
269-
>>> obj.set_de_bounds()
270-
>>> print(obj.de_bounds)
271-
[[min_theta, max_theta], [min_theta, max_theta], ..., [min_p, max_p], [min_Lambda, max_Lambda]]
269+
S = Kriging(name='kriging', seed=124)
270+
S.set_de_bounds()
271+
print(S.de_bounds)
272272
273273
Returns:
274274
None

‎src/spotPython/fun/objectivefunctions.py‎

Lines changed: 0 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -683,21 +683,3 @@ def fun_random_error(self, X: np.ndarray, fun_control: Optional[Dict] = None) ->
683683
else:
684684
print(y)
685685
return y
686-
687-
def fun_hcf(X: np.ndarray, fun_control: Optional[Dict] = None) -> np.ndarray:
688-
def oneMass(params, Oerr) -> np.ndarray:
689-
Oeig, D, F = params
690-
nue = Oerr / Oeig
691-
V = F / (np.sqrt((1 - nue**2) ** 2 + (4 * D**2 * nue**2)))
692-
return V
693-
694-
def objective_function(params, Oerr, amplitudes) -> np.ndarray:
695-
return np.sum((amplitudes - oneMass(params, Oerr)) ** 2)
696-
697-
frequencies_filtered = fun_control["frequency_filtered"]
698-
amplitudes_filtered = fun_control["amplitude_filtered"]
699-
X = np.atleast_2d(X)
700-
y = np.array([], dtype=float)
701-
for i in range(X.shape[0]):
702-
y = np.append(y, objective_function(X[i], frequencies_filtered, amplitudes_filtered))
703-
return y

‎test/test_kriging.py‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
from spotPython.build.kriging import Kriging
2+
import numpy as np
3+
from math import erf
4+
5+
def test_interpolation_property():
6+
"""
7+
Kriging is a interpolator, so the mean prediction
8+
should be equal to the training points:
9+
check if the difference between the mean prediction
10+
and the true value in the training points is smaller than 1e-6.
11+
"""
12+
from spotPython.build.kriging import Kriging
13+
import numpy as np
14+
import matplotlib.pyplot as plt
15+
from numpy import linspace, arange
16+
rng = np.random.RandomState(1)
17+
X = linspace(start=0, stop=10, num=1_0).reshape(-1, 1)
18+
y = np.squeeze(X * np.sin(X))
19+
training_indices = rng.choice(arange(y.size), size=6, replace=False)
20+
X_train, y_train = X[training_indices], y[training_indices]
21+
S = Kriging(name='kriging', seed=124)
22+
S.fit(X_train, y_train)
23+
mean_prediction, std_prediction, s_ei = S.predict(X, return_val="all")
24+
assert np.allclose(mean_prediction[training_indices], y[training_indices], atol=1e-6)
25+
26+
def test_ei():
27+
"""
28+
Test computation of expected improvement based on (3.8) in Forrester et al. (2008).
29+
"""
30+
S = Kriging(name='kriging', seed=124)
31+
S.mean_cod_y = [0.0, 0.0, 0.0, 0.0, 0.0]
32+
# assert that the S.exp_imp(1.0, 0.0) is equal to 0.0,
33+
# because EI is zero when the std is zero
34+
assert 0.0 == S.exp_imp(1.0, 0.0)
35+
# assert that the S.exp_imp(0.0, 1.0) is equal to 1/sqrt(2 pi)
36+
# assert S.exp_imp(0.0, 1.0) == 1/np.sqrt(2*np.pi)
37+
# play safe and use np.allclose
38+
assert np.allclose(S.exp_imp(0.0, 1.0), 1/np.sqrt(2*np.pi), atol=1e-6)
39+
# assert S.exp_imp(1.0, 1.0) is correct based on (3.8) in Forrester et al. (2008)
40+
assert np.allclose(S.exp_imp(1.0, 1.0), -(0.5 + 0.5*erf(-1/np.sqrt(2))) + 1/np.sqrt(2*np.pi)*np.exp(-1/2), atol=1e-6)
41+
42+
def test_de_bounds():
43+
"""
44+
Test if the bounds for the DE algorithm are set correctly.
45+
"""
46+
S = Kriging(name='kriging', seed=124)
47+
S.set_de_bounds()
48+
assert S.de_bounds == [[-3, 2]]
49+
n = 10
50+
S = Kriging(name='kriging', seed=124, n_theta=n)
51+
S.set_de_bounds()
52+
assert len(S.de_bounds) == n
53+
n=2
54+
p=4
55+
S = Kriging(name='kriging', seed=124, n_theta=n, n_p=p, optim_p=True)
56+
S.set_de_bounds()
57+
assert len(S.de_bounds) == n+p
58+
S = Kriging(name='kriging', seed=124, n_theta=n, n_p=p, optim_p=False)
59+
S.set_de_bounds()
60+
assert len(S.de_bounds) == n

0 commit comments

Comments
 (0)