import math

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate

global l, r, x, y


def f(x, y):
    return x * y * beta / (x + y)


def RoundAndZeros(X, Y):
    X = round(X, rou)
    Y = round(Y, rou)
    X = 0 if X < 0 else X
    Y = 0 if Y < 0 else Y
    return X, Y


def system1(cond1, itter):
    X, Y = cond1[0], cond1[1]
    X, Y = RoundAndZeros(X, Y)
    dXdT = round(-f(X, Y), rou)
    dYdT = round(f(X, Y), rou) - gamma * Y
    return [dXdT, dYdT]


def system2(cond1, itter):
    X, Y = cond1[0], cond1[1]
    X, Y = RoundAndZeros(X, Y)
    dXdT = round(-f(X, Y), rou) - u if X > 0 else round(-f(X, Y), rou)
    dYdT = round(f(X, Y), rou) - gamma * Y
    return [dXdT, dYdT]


# вычисление значения y
def f_fun(l, r, e):
    global x, y
    X = np.linspace(l, r, n)  # for forward integration
    cond1 = [x, y]
    if e == 1:
        sol1 = integrate.odeint(system1, cond1, X)
        Psi, Psi1 = sol1[:, 0], sol1[:, 1]
        x, y = Psi[n - 1], Psi1[n - 1]
        x, y = RoundAndZeros(x, y)
        return Psi, Psi1
    if e == 2:
        sol2 = integrate.odeint(system2, cond1, X)
        Fi, Fi1 = sol2[:, 0], sol2[:, 1]
        return Fi, Fi1


def J(l, a, r):
    global x, y
    # решение задачи коши вперед
    x1, y1 = f_fun(l, a, 1) if l > a else 0, 0
    x2, y2 = f_fun(a, r, 2)
    x, y = x0, y0
    leftIntegral = integrate.quad(integrateLeft, l, a, args=y1[n - 1]) if l > a else [0]
    rightIntegral = integrate.quad(integrateRight, a, r, args=y2[n - 1])
    return leftIntegral[0] + rightIntegral[0]


def integrateLeft(t, y1):
    return q * y1


def integrateRight(t, y2):
    return q * y2 + c * u


def goldenSectionSearch(eps):
    global l, r
    l1, r1 = l, r
    phi = (1 + math.sqrt(5)) / 2
    resphi = 2 - phi
    a1 = l + resphi * (r - l)
    a2 = r - resphi * (r - l)
    f1 = J(l, a1, r)
    f2 = J(l, a2, r)
    while abs(f1 - f2) > eps:
        if f1 < f2:
            r = a2
            a2 = a1
            f2 = f1
            a1 = l + resphi * (r - l)
            f1 = J(l, a1, r)
        else:
            l = a1
            a1 = a2
            f1 = f2
            a2 = r - resphi * (r - l)
            f2 = J(l, a2, r)
    return (a1 + a2) / 2


# Левый конец
l = 0
# Правый конец
r = 60
t = r
eps = 0.01
n = t + 1  # odd integer number
x0 = 10000
y0 = 4000
x, y = x0, y0
u = 350
c = 1
q = 1
beta = 0.2
gamma = 0.1
X = np.linspace(0, n - 1, n)

# округление чисел после запятой
rou = 6

zolotoe = goldenSectionSearch(eps)
print("Дата рекомендованной вакцинации", zolotoe)
print("Затраты на вакцинацию при стоимости вакцины 20 по 35 привитых", J(0, zolotoe, t))

x1_1, y1_1 = f_fun(0, t, 1)
# plt.plot(X[0:n - 2], x1_1[0:n - 2], color="red", zorder=3)
# plt.plot(X[0:n - 2], y1_1[0:n - 2], color="green", zorder=3)

x, y = x0, y0
x1_2, y1_2 = f_fun(0, t, 2)
# plt.plot(X[0:n - 2], x1_2[0:n - 2], zorder=2)
# plt.plot(X[0:n - 2], y1_2[0:n - 2], color="purple", zorder=2)

n1 = n
# zolotoe = 61
zolotoe = round(zolotoe)
n = int(zolotoe)

# y2_1 = f_fun(0, zolotoe, True)
n1 = t
Jgraf = np.zeros(n1)
nj = 0
while nj < n1:
    Jgraf[nj] = J(0, nj, t)
    nj = nj + 1

plt.plot(X[0:n1 - 2], Jgraf[0:n1 - 2], color="gold", zorder=4)
#
# y2_1 = np.zeros(n)
#
# i = 0
# while i < n:
#     y2_1[i] = y1_1[i]
#     i = i + 1
#
# if i != 0:
#     x, y = x1_1[i - 1], y1_1[i - 1]
#
# n = n1 - n
# x2_2, y2_2 = f_fun(zolotoe, t, 2)
#
# n = n1
# y211 = np.zeros(n)
#
# i = 0
# while i < int(zolotoe) - 1:
#     y211[i] = y2_1[i]
#     i = i + 1
#
# j = 0
# while j < (n - zolotoe) - 1:
#     y211[i] = y2_2[j]
#     j = j + 1
#     i = i + 1
#
# ii = 0
# while ii < n:
#     if (y211[ii] <= 0):
#         y211[ii] == None
#     ii = ii + 1
#
# yy = [*y211]
# plt.plot(X[0:n - 2], yy[0:n - 2], color="green", zorder=1)
plt.show()