138 lines
5.1 KiB
Python
138 lines
5.1 KiB
Python
# core/doe.py
|
|
"""
|
|
Модуль с функциями для генерации планов эксперимента.
|
|
"""
|
|
import random
|
|
import numpy as np
|
|
|
|
def ispoz(N):
|
|
"""Проверяет, является ли N положительным целым числом (для использования в ze)."""
|
|
return isinstance(N, (int, float)) and not isinstance(N, bool) and N > 0
|
|
|
|
def ze(k, n):
|
|
"""Генерирует матрицу размера n x k, заполненную нулями (центральные точки)."""
|
|
if n is None or k is None:
|
|
return []
|
|
elif ispoz(k) and ispoz(n):
|
|
return [[0] * k for _ in range(n)]
|
|
return []
|
|
|
|
def ffe(k):
|
|
"""Полный факторный план 2^k. Возвращает список списков размером 2^k x k со значениями ±1."""
|
|
if k < 2:
|
|
raise ValueError("Factor must be at least 2")
|
|
N = 2 ** k
|
|
out = []
|
|
for i in range(N):
|
|
row = []
|
|
for j in range(k):
|
|
bit = ((i >> j) & 1) * 2 - 1
|
|
row.append(bit)
|
|
out.append(row)
|
|
return out
|
|
|
|
def ccp(k, alpha=1.0, center_points=1, fractional=None):
|
|
"""
|
|
Центральный композиционный план (ЦКП) для k факторов.
|
|
alpha – звёздное плечо.
|
|
center_points – количество центральных точек.
|
|
fractional – степень дробления ядра (если None, используется полный факторный план).
|
|
"""
|
|
if k < 2:
|
|
raise ValueError("Для ЦКП необходимо минимум 2 фактора.")
|
|
if fractional is not None:
|
|
if fractional <= 0 or fractional >= k:
|
|
raise ValueError("Степень дробления p должна быть от 1 до k-1")
|
|
base_k = k - fractional
|
|
if base_k < 2:
|
|
raise ValueError(f"Ядро дробного факторного плана должно иметь минимум 2 фактора, а получается {base_k}.")
|
|
if fractional is not None:
|
|
from core.doe import fractional_factorial
|
|
core_plan = fractional_factorial(k, fractional)
|
|
else:
|
|
core_plan = ffe(k)
|
|
star_points = []
|
|
for i in range(k):
|
|
row = [0.0] * k
|
|
row[i] = alpha
|
|
star_points.append(row)
|
|
row_neg = [0.0] * k
|
|
row_neg[i] = -alpha
|
|
star_points.append(row_neg)
|
|
center_rows = [[0.0] * k for _ in range(center_points)]
|
|
plan = core_plan + star_points + center_rows
|
|
return plan
|
|
|
|
def fractional_factorial(k, p=1):
|
|
"""Дробный факторный план 2^(k-p). p – степень дробления."""
|
|
if k < 2 or p <= 0 or p >= k:
|
|
raise ValueError("Некорректные параметры: k>=2, 0<p<k")
|
|
base_k = k - p
|
|
if base_k < 2:
|
|
raise ValueError(f"Ядро должно содержать минимум 2 фактора, а получается {base_k}.")
|
|
N = 2 ** base_k
|
|
base_plan = ffe(base_k)
|
|
plan = []
|
|
for row in base_plan:
|
|
new_row = row.copy()
|
|
for i in range(p):
|
|
prod = 1
|
|
for j in range(base_k):
|
|
prod *= row[j]
|
|
new_row.append(prod)
|
|
plan.append(new_row)
|
|
return plan
|
|
|
|
def plackett_burman(k):
|
|
"""
|
|
План Плакетта-Бермана для k факторов.
|
|
Число опытов N – ближайшее кратное 4, >= k+1.
|
|
Для N = 8, 16, 32 используется рекурсивное построение матрицы Адамара.
|
|
"""
|
|
if k < 2:
|
|
raise ValueError("Требуется минимум 2 фактора.")
|
|
N = 4 * ((k + 3) // 4)
|
|
|
|
def hadamard_rec(n):
|
|
if n == 1:
|
|
return np.array([[1]])
|
|
H = hadamard_rec(n//2)
|
|
top = np.hstack([H, H])
|
|
bottom = np.hstack([H, -H])
|
|
return np.vstack([top, bottom])
|
|
|
|
if (N & (N - 1)) == 0: # N — степень двойки
|
|
H = hadamard_rec(N)
|
|
else:
|
|
try:
|
|
from scipy.linalg import hadamard
|
|
H = hadamard(N)
|
|
except ImportError:
|
|
print("Предупреждение: scipy не установлена. Генерируем случайный план.")
|
|
plan = []
|
|
for _ in range(N):
|
|
row = [random.choice([-1, 1]) for _ in range(k)]
|
|
plan.append(row)
|
|
return plan
|
|
plan = H[:, 1:k+1].tolist()
|
|
return plan
|
|
|
|
def box_behnken(k):
|
|
"""
|
|
План Бокса-Бенкена для k факторов (k >= 3).
|
|
Содержит комбинации на серединах рёбер и центральные точки.
|
|
"""
|
|
if k < 3:
|
|
raise ValueError("Требуется минимум 3 фактора.")
|
|
plan = []
|
|
for i in range(k):
|
|
for j in range(i+1, k):
|
|
for a, b in [(1,1), (1,-1), (-1,1), (-1,-1)]:
|
|
row = [0.0] * k
|
|
row[i] = a
|
|
row[j] = b
|
|
plan.append(row)
|
|
for _ in range(3):
|
|
plan.append([0.0] * k)
|
|
return plan
|