compute-math/lab2/HookJeeves.py

# 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, 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()

    error_growth_coeff = 0.05  # КОЭФФИЦИЕНТ РОСТА ОШИБКИ

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

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

        x2, x1 = x2 + step, x2

        # Новая погрешность
        step_error = error_growth_coeff * np.abs(step)
        x2_err = x2_err + step_error
        x1_err = x2_err.copy()

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


class TARGET:
    MIN = 0
    MAX = 1


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

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

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

    # условие окончания поиска
    epsilon = 0.0001

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

    x1r = [-10000, 10000]
    x2r = [-10000, 10000]

    target = TARGET.MIN

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

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

    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()