compute-math/lab2/HookJeeves.py

import numpy as np
import AbsoluteMeasurementErrors as abserror
from itertools import product

result_log = ""


def log(*values, sep=" ", end="\n"):
    global result_log
    result_log += sep.join(list(map(str, values))) + end
    print(*values, sep=sep, end=end)


def get_errors_vector(x):
    return np.array(list(map(abserror.get_error, x)))


def sum_errors_vector(a, b):
    return np.array([abserror.sum_error(a[i], b[i]) for i in range(len(a))])


def frac_errors_vector(a, b):
    return np.array([abserror.frac_error(a[i], b) for i in range(len(a))])


def point_info(x, f, x_err):
    return f"x = {x}, x_err = {x_err}, f(x) = {f(x)}"


def hooke_jeeves(
    f,
    x0,
    delta0,
    epsilon,
    alpha,
    r=lambda x: True,
    fit=lambda a, b: a < b,
    x_err=None,
    delta_err=None,
):
    x = np.array(x0)
    delta = np.array(delta0)

    if x_err is None:
        x_err = get_errors_vector(x)

    if delta_err is None:
        delta_err = get_errors_vector(delta)

    iteration = 0

    while np.linalg.norm(delta) > epsilon:
        iteration += 1

        log()
        log("=" * 40)
        log("Итерация", iteration)

        log("Текущая базовая точка", point_info(x, f, x_err))

        sample_x, sample_x_err = exploratory_search(
            f, x, delta, r=r, fit=fit, x_err=x_err, delta_err=delta_err
        )

        if fit(f(sample_x), f(x)):
            log("Исследующий поиск УДАЧНЫЙ")

            x_p, x_p_err = sample_search(
                f, x, sample_x, r=r, fit=fit, x1_err=x_err, x2_err=sample_x_err
            )

            if fit(f(x_p), f(x)):
                log("Поиск по образцу УДАЧНЫЙ")
                x, x_err = x_p, x_p_err
            else:
                log("Поиск по образцу ПРОВАЛЕН")
                x, x_err = sample_x, sample_x_err

            log()
            log("Новая базовая точка", point_info(x, f, x_err))

        else:
            log("Исследующий поиск ПРОВАЛЕН")

            log("Уменьшаем шаг", delta, "->", delta / alpha)
            delta = delta / alpha
            delta_err = frac_errors_vector(delta, alpha)
            log("ε =", np.linalg.norm(delta))

    return x, f(x), x_err


# def exploratory_search(
#     f, x, delta, r=lambda x: True, fit=lambda a, b: a < b, x_err=None
# ):
#     log()
#     log("Выполняем исследующий поиск")

#     if x_err is None:
#         x_err = get_errors_vector(x)

#     x_new = np.array(x)
#     f_x_new = f(x_new)

#     for i in range(len(x)):
#         x_up = x_new.copy()
#         x_down = x_new.copy()

#         x_up[i] += delta[i]
#         x_down[i] -= delta[i]

#         f_x_up = f(x_up)
#         f_x_down = f(x_down)

#         if not r(x_up) or not r(x_down):
#             log("Ограничение ОДЗ")

#         if fit(f_x_up, f_x_new) and fit(f_x_up, f_x_down) and r(x_up):
#             x_new = x_up
#             f_x_new = f_x_up
#         elif fit(f_x_down, f_x_new) and r(x_down):
#             x_new = x_down
#             f_x_new = f_x_down

#     if any([x_new[i] != x[i] for i in range(len(x))]):
#         log("Найдена точка", x_new, f(x_new))

#     return x_new, x_err


def exploratory_search(
    f, x, delta, r=lambda x: True, fit=lambda a, b: a < b, x_err=None, delta_err=None
):
    log()
    log("Выполняем исследующий поиск")

    if x_err is None:
        x_err = get_errors_vector(x)

    if delta_err is None:
        delta_err = get_errors_vector(delta)

    dxs = product(*[[d, -d] for d in delta])

    x = np.array(x)
    t_x_err = None
    t_f = f(x)
    t_x = None

    for dx in dxs:
        if fit(f(x + dx), t_f) and r(x + dx):
            t_f = f(x + dx)
            t_x = x + dx
            t_x_err = x_err + delta_err

    if t_x is not None:
        log("Найдена точка", t_x, t_f)
        return t_x, t_x_err

    return x, x_err


def sample_search(
    f, x1, x2, r=lambda x: True, fit=lambda a, b: a < b, x1_err=None, x2_err=None
):
    log()
    log("Выполняем поиск по образцу")

    if x1_err is None:
        x1_err = get_errors_vector(x1)

    if x2_err is None:
        x2_err = get_errors_vector(x2)

    x2_last = x2.copy()
    x2_last_err = x2_err.copy()

    while fit(f(x2), f(x1)):
        x2_last = x2.copy()
        x2_last_err = x2_err.copy()

        if not r(x2 + (x2 - x1)):
            log("Ограничение ОДЗ")
            break

        x2, x1 = x2 + (x2 - x1), x2

        t_err = x2_err.copy()
        x2_err = x2_err + (x2_err + x1_err)
        x1_err = t_err.copy()

    log("Найдена точка", x2_last, f(x2_last))
    return x2_last, x2_last_err


def example():
    # начальная базовая точка
    x0 = [100, 100]

    # начальное значение приращения
    delta0 = [10, 10]

    # коэффициент приращения
    alpha = 2

    # условие окончания поиска
    epsilon = 1e-4

    # коэффициенты при неизвестных (порядок обратный)
    c = [[1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]

    def fit(a, b):
        """
        Функция сравнения значений оптимизируемой функции
        > - ищем максимум
        < - ищем минимум
        """
        return a > b

    def r(x):
        """
        Функция, описывающая область допустимых значений.
        Должна возвращать истину, если точка находится в ОДЗ.
        """
        return -10000 < x[0] < 10000 and -10000 < x[1] < 10000

    def f(x):
        """
        Функция, которую мы должны оптимизировать
        """
        return sum(
            [
                sum([c[i][j] * x[i] ** j for j in range(len(c[i]))])
                for i in range(len(c))
            ]
        )

    x, val, x_err = hooke_jeeves(f, x0, delta0, epsilon, alpha, r=r, fit=fit)

    log()
    log("=" * 40)
    log("Точка экстремумв:", x)
    log("Абсолютная погрешность:", x_err)
    log("Максимальное значение функции:", val)

    with open("result.log", "w", encoding="utf-8") as file:
        file.write(result_log)


def main():
    example()


if __name__ == "__main__":
    main()