PythonWorkshop/main.py at master · Parowicz/PythonWorkshop · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
from typing import List, Tuple

import matplotlib.pyplot as plt

from solver import Solver, Chromosome
from timeseries import holt, holt_winters, mape
from data import load_wb_xml


FORECAST_HORIZON = 10


def f(x: float) -> float:
    return (14 * x ** 2) + (6 * x) - 8


def find_root():
    """
    Find root of 14x^2 + 6x - 8
    """

    def fit(chromosome: Chromosome) -> float:
        return abs(f(chromosome[0]))

    solver = Solver(generations=10)

    model = solver.minimize(fit, [(-1000, 1000)])

    print("Example 1".center(80, "*"))
    print("f(%f)=%f" % (model[0], f(model[0])))
    print("*" * 80)


def find_roots():
    """
    Find all roots of 14x^2 + 6x - 8
    We know that this function does have exactly 2 roots in a real domains
    """

    def fit(chromosome: Chromosome) -> float:
        root1 = abs(f(chromosome[0]))
        root2 = abs(f(chromosome[1]))

        return root1 + root2 - abs(chromosome[0] - chromosome[1])

    solver = Solver(generations=150)

    model = solver.minimize(fit, [(-1000, 1000), (-1000, 1000)])

    print("Example 2".center(80, "*"))
    print("x1=%f, x2=%f" % (model[0], model[1]))
    print("f(x1)=%f, f(x2)=%f" % (f(model[0]), f(model[1])))
    print("*" * 80)


def optimize_holt(data: List[Tuple[int, float]], horizon: int):
    """
    Optimize Holt parameters

    :param data: Time series to forecast
    :param horizon: How many step forward we want to forecast
    """

    years, measurements = zip(*data)
    estimated_years = years + tuple(range(years[-1], years[-1] + horizon))

    def fit(chromosome: Chromosome) -> float:
        smoothing, _ = holt(
            measurements,
            chromosome[0],
            chromosome[1],
            chromosome[2],
            chromosome[3],
            horizon=horizon,
        )

        return mape(measurements, smoothing)

    solver = Solver()
    model = solver.minimize(fit, [(0, 1), (0, 1), (0, 100), (0, 100)])
    print("Example 3".center(80, "*"))
    print("MAPE = %f" % fit(model))
    smoothing, forecast = holt(
        measurements, model[0], model[1], model[2], model[3], horizon=horizon
    )
    plt.plot(years, measurements, label="Population, total")
    plt.plot(estimated_years, smoothing + forecast, label="Smoothing")
    plt.title(f"Population, total - Benin. Holt smoothing")
    plt.legend()

    plt.savefig("total_population.png")
    plt.clf()
    print("*" * 80)


def optimize_holt_winters(data: List[Tuple[int, float]], horizon: int):
    """
    Optimize Holt-Winters parameters

    :param data: Time series to forecast
    :param horizon: How many step forward we want to forecast
    """

    years, measurements = zip(*data)
    estimated_years = years + tuple(range(years[-1], years[-1] + horizon))

    def fit(chromosome: Chromosome) -> float:
        smoothing, _ = holt_winters(
            measurements, chromosome[0], chromosome[1], chromosome[2], 12, horizon
        )

        return mape(measurements, smoothing)

    solver = Solver()
    model = solver.minimize(fit, [(0, 1), (0, 1), (0, 1)])

    print("Example 4".center(80, "*"))
    print("MAPE: %f" % fit(model))
    smoothing, forecast = holt_winters(
        measurements, model[0], model[1], model[2], 12, horizon
    )
    plt.plot(years, measurements, label="Population, total")
    plt.plot(estimated_years, smoothing + forecast, label="Smoothing")
    plt.title("Population, total - Benin. Holt-Winters smoothing")
    plt.legend()

    plt.savefig("total_population2.png")
    plt.clf()
    print("*" * 80)


if __name__ == "__main__":
    find_root()
    find_roots()

    persistence_data = load_wb_xml(r"Persistence_to_last_grade_of_primary.xml", "WLD")

    optimize_holt(persistence_data, FORECAST_HORIZON)
    optimize_holt_winters(persistence_data, FORECAST_HORIZON)