Skip to content
Snippets Groups Projects
Commit 856e954b authored by Jiří Kalvoda's avatar Jiří Kalvoda
Browse files

PaSt zápočťák

parent bc1bc9b8
No related branches found
No related tags found
No related merge requests found
Showing
with 611 additions and 0 deletions
#!/usr/bin/env python3
import main
main.web.run()
print(main.timer)
../data_lib.py
\ No newline at end of file
import plotly
import plotly.subplots
import plotly.graph_objects as go
import numpy as np
import data_lib
algo_to_name = {
"greedy": "Hladové ř.",
"rg": "Rekurzivní ř.",
"rsg": "Hvězdičkové r. ř.",
"semidef_prog(10)": "Semidefinitní ř.",
}
algo_to_color = {
"greedy": '#8a00d4',
"rg": '#d527b7',
"rsg": '#f782c2',
"semidef_prog(10)": '#f9c46b',
}
def draw_algo_graph(fig, data, algo, n, val_getter=lambda x:x.score, name=None, color=None):
d = data_lib.group_by_n(data.pipelines[algo])[n]
y = [0 for i in range(2*n)]
for i in d:
y[val_getter(i)] += 1
fig.add_trace(go.Histogram(x=[val_getter(i) for i in d], name=name or algo_to_name[algo], xbins=dict(size=1), marker=dict(color=color or algo_to_color[algo])))
fig.add_trace(go.Box(x=[val_getter(i) for i in d], name=name or algo_to_name[algo], showlegend=False, marker=dict(color=color or algo_to_color[algo]), boxmean=True), row=2, col=1)
fig.update_xaxes(range=[0, 2*n])
# xaxis=np.arange(0, 2*n)
fig.update_layout(yaxis_title="Počet řešení")
fig.update_xaxes(title_text="Skóre", row=2, col=1)
def intro_graph(algo):
fig = plotly.subplots.make_subplots(rows=2, cols=1, row_heights=[0.8, 0.2], shared_xaxes=True, vertical_spacing = 0.05)
draw_algo_graph(fig, data_lib.Data("log-intr", validate_versions=False), algo, 200)
return fig
from typing import List, Optional, Union, Tuple
import types
tag_names = ["a", "abbr", "acronym", "address", "applet", "area", "article", "aside", "audio", "b", "base", "basefont", "bdi", "bdo", "big", "blockquote", "body", "br", "button", "canvas", "caption", "center", "cite", "code", "col", "colgroup", "data", "datalist", "dd", "del", "details", "dfn", "dialog", "dir", "div", "dl", "dt", "em", "embed", "fieldset", "figcaption", "figure", "font", "footer", "form", "frame", "frameset", "head", "header", "hgroup", "h1", "h2", "h3", "h4", "h5", "h6", "hr", "html", "i", "iframe", "img", "input", "ins", "kbd", "keygen", "label", "legend", "li", "link", "main", "map", "mark", "menu", "menuitem", "meta", "meter", "nav", "noframes", "noscript", "object", "ol", "optgroup", "option", "output", "p", "param", "picture", "pre", "progress", "q", "rp", "rt", "ruby", "s", "samp", "script", "section", "select", "small", "source", "span", "strike", "strong", "style", "sub", "summary", "sup", "svg", "table", "tbody", "td", "template", "textarea", "tfoot", "th", "thead", "time", "title", "tr", "track", "tt", "u", "ul", "var", "video", "wbr"]
class Markup(str):
pass
class Html():
def __init__(html, Markup=Markup):
html.Markup = Markup
html.doctype_header = "<!DOCTYPE html>\n"
def public(f):
setattr(html, f.__name__, f)
return f
@public
class EscapeError(RuntimeError):
pass
@public
def escape(a):
return a.replace("&", "&amp;").replace("<", "&lt;").replace(">", "&gt;").replace("'", "&#39;").replace("\\", '#34;')
@public
def escape_attribute(x: str) -> str:
return str(x).replace("&", "&amp;").replace('"', "&quot;")
@public
def escape_attribute_name(x: str) -> str:
for c in "<>='\" ":
if c in x:
raise EscapeError
return x
def escape_tag_name(x: str) -> str:
return escape_attribute_name(x)
Element = Union[str, Markup, 'Bucket']
html.Element = Element
@public
def indent_str(indent, line_mode, disable_indent, s):
if disable_indent: indent = 0
r = ("\n" + "\t"*indent).join(s.split("\n"))
if line_mode:
return r
else:
return "\t"*indent + r + "\n"
@public
class Bucket:
content: List[Element]
is_paired: bool
builder: 'Builder'
def __init__(self, builder):
self.builder = builder
self.content = []
self.before_tag = None
def add(self, x: Element):
if isinstance(x, types.FunctionType):
before_current_tag = self.builder.current_tag
self.builder.current_tag = self
try:
x(self)
finally:
self.builder.current_tag = before_current_tag
else:
if isinstance(self, Tag):
self.is_paired = True
self.content.append(x)
def __call__(self, *arg, **kvarg):
for i in arg:
self.add(i)
for k, v in kvarg.items():
self.add_attribute(remove_leading_underscore(k), v)
return self
def add_tag(self, name: str, attributes: List[Tuple[str, str]]):
t = Tag(self.builder, name, attributes)
self.add(t)
return t
def wrap(self):
pass
def print(self):
out = []
self.serialize_append_to_list(out, 0, False, False)
return Markup("".join(out))
def print_file(self):
out = [html.doctype_header+"\n"]
self.serialize_append_to_list(out, 0, False, False)
return Markup("".join(out))
def __enter__(self):
if self.before_tag is not None:
raise RuntimeError("Duplicit __enter__")
if isinstance(self, Tag):
self.is_paired = True
self.before_tag = self.builder.current_tag
self.builder.current_tag = self
return self
def __exit__(self, exc_type, exc_value, exc_traceback):
if self.before_tag is None:
raise RuntimeError("__exit__ before __enter__")
self.builder.current_tag = self.before_tag
self.before_tag = None
def serialize_append_to_list(self, out, indent, line_mode, disable_indent):
for i in self.content:
if isinstance(i, Bucket):
i.serialize_append_to_list(out, indent, line_mode, disable_indent)
elif isinstance(i, Markup):
out.append(indent_str(indent, line_mode, disable_indent, i))
else:
out.append(indent_str(indent, line_mode, disable_indent, escape(str(i))))
@public
class Tag(Bucket):
name: str
attributes: List[Tuple[str, str]]
def __init__(self, builder, name: str, attributes: List[Tuple[str, str]]):
super().__init__(builder)
self.name = name
self.attributes = attributes
self.is_paired = False
def add_attribute(self, k, v):
self.attributes.append((k, v))
def format_attributes(self):
return " ".join(f'{escape_attribute_name(i[0])}="{escape_attribute(i[1])}"' for i in self.attributes)
def serialize_append_to_list(self, out, indent, line_mode, disable_indent):
if self.is_paired:
out.append(indent_str(indent, line_mode, disable_indent, f"<{escape_tag_name(self.name)} {self.format_attributes()}>"))
super().serialize_append_to_list(out, indent + (0 if line_mode else 1), line_mode, disable_indent or self.name == "pre")
if self.name == "pre":
print(self)
out.append(indent_str(indent, line_mode, disable_indent, f"</{escape_tag_name(self.name)}>"))
else:
out.append(indent_str(indent, line_mode, disable_indent, f"<{escape_tag_name(self.name)} {self.format_attributes()} />"))
def set_paired(self, val=True):
self.is_paired = val
return self
@public
class Line(Bucket):
def serialize_append_to_list(self, out, indent, line_mode, disable_indent):
if line_mode:
super().serialize_append_to_list(out, indent, True, disable_indent)
else:
out.append(indent_str(indent, line_mode, disable_indent, "")[:-1])
super().serialize_append_to_list(out, indent, True, disable_indent)
out.append("\n")
@public
class Builder:
current_tag: Bucket
root_tag: Bucket
def __init__(self, tag: Bucket = None):
if tag is None:
tag = Bucket(self)
self.root_tag = tag
self.current_tag = tag
def add_tag(self, name: str, attributes: List[Tuple[str, str]] = []):
return self.current_tag.add_tag(name, attributes)
def __call__(self, *arg, **kvarg):
self.current_tag(*arg, **kvarg)
return self
def print(self):
return self.root_tag.print()
def print_file(self):
return self.root_tag.print_file()
@public
class HtmlBuilder(Builder):
def __init__(self):
super().__init__(Tag(self, "html", []))
def remove_leading_underscore(s):
if s == "":
return s
if s[0] == "_":
return s[1:]
return s
@public
class WrapAfterBuilder(Builder):
def __init__(self, f):
super().__init__(Bucket(self))
self.wrap_done = False
self.wrap_function = f
def wrap(self, *arg, **kvarg):
if self.wrap_done:
return
self.wrap_done = True
content = self.root_tag.content
self.root_tag = None
self.current_tag = None
self.root_tag = self.wrap_function(self, content, *arg, **kvarg) or self.root_tag
return self
def print(self, *arg, **kvarg):
self.wrap()
return super().print(*arg, **kvarg)
def print_file(self, *arg, **kvarg):
self.wrap()
return super().print_file(*arg, **kvarg)
@public
def WrapAfterBuilder_decorator(f):
def l():
return WrapAfterBuilder(f)
return l
@public
def add_taglike(name, f):
assert hasattr(Builder, name) is False
def l1(self, *args, **kvarg):
x = f(self.builder, *args, **kvarg)
self.add(x)
return x
setattr(Bucket, name, l1)
setattr(Builder, f"_{name}", f)
def l3(self, *args, **kvarg):
return getattr(self.current_tag, name)(*args, **kvarg)
setattr(Builder, name, l3)
def l4(self, *args, **kvarg):
return getattr(self.builder, f"_{name}")(*args, **kvarg)
setattr(Bucket, f"_{name}", l4)
@public
def add_taglike_decorator(f):
add_taglike(f.__name__, f)
for name in tag_names:
def run(name):
def l1(builder, *args, **kvarg):
return Tag(builder, name, [(remove_leading_underscore(k), v) for k, v in kvarg.items()])(*args)
add_taglike(name, l1)
run(name)
@add_taglike_decorator
def line(builder, *args):
return Line(builder)(*args)
@add_taglike_decorator
def markup(builder, *args):
return Markup(*args)
@add_taglike_decorator
def bucket(builder, *args):
return Bucket(builder)(*args)
@add_taglike_decorator
def parse(b, html_str):
# return Markup(html_str)
import AdvancedHTMLParser
import html as python_html
parser = AdvancedHTMLParser.AdvancedHTMLParser()
parser.parseStr(html_str)
def parse(nd):
if type(nd)==str:
b(Markup(nd))
elif isinstance(nd, AdvancedHTMLParser.Tags.AdvancedTag):
t = Tag(b, nd.tagName, [(k, v) for k, v in nd.attributesList])
b(t)
if not nd.isSelfClosing:
with t:
for x in nd.childBlocks:
parse(x)
else:
raise RuntimeError(type(nd))
with b._line() as root:
for nd in parser.getRootNodes():
parse(nd)
return root
@public
def customtag_decorator(f):
def l(b, *args, **kvarg):
with b._bucket() as b:
return f(b, *args, **kvarg) or b
add_taglike(f.__name__, l)
return f
---
title: "PaSt: Zápočtová úloha"
lang: "cs"
highlight-style: "dracula"
---
# Problém
V této práci si představíme binary paintshop problem.
Jedná se o úlohu, kterou neumíme efektivně řešit a ani aproximovat [@computing],
takže mimo jiné probíhá aktivní výzkum snažící se najít algoritmus,
který je dobrý v průměrném případě (pro náhodný vstup).
V této zápočtové práci bych se rád podíval na některé z těchto algoritmů
a pomocí statistických metod pro zvolenou velikost vstupu určil, který z nich
je lepší.
# Zadání
Nejprve si pojďme představit zadání binary paint shop problému:
V řadě je $2n$ aut $n$ různých typů -- od každého typu 2.
Chtěli bychom od každého typu nabarvit jedno auto červeně a druhé modře.
Auto však na barvící linky vjíždí v pořadí, v jakém jsou v řadě.
Barvící linka je optimalizovaná na barvení velkého počtu aut jednou barvou.
Tedy měnit barvu, kterou se barví, je složitá a drahá záležitost.
Chceme tedy nají obarvení aut tak, aby od každého typu bylo jedno červené a jedno modré, přitom počet změn barev v řadě byl co nejmenší.
Počet změn barev značíme jako skóre algoritmu.
# Programy
## Hladové řešení
Jednoduché řešení je následující:
Půjdeme postupně po řadě.
První auto v řadě obarvíme řekněme červeně.
Dále budeme vždy první auto daného typu barvit stejně jako předcházející
auto.
Druhé auto daného typu obarvíme vždy zbývající barvou.
Snažit se nějak statisticky analyzovat a porovnávat tento algoritmus není moc
zajímavé, protože je známo, že střední hodnota počtu změn hladového algoritmu je $\sum_{k=0}^{n-1} {2k^2-1 \over 4k^2-1}$ [@glim_gr], což je pro velká $n$ zhruba ${1\over 2} n$ (formálně řečeno $\lim_{n\to\infty} {\mathbb{E}_n({\rm g}) \over n} = {1\over 2}$).
Následuje graf skóre $100$ běhů hladového algoritmu pro $n=200$.
```python {c=plotly}
import g, data_lib
fig = g.intro_graph("greedy")
```
## Rekurzivní hladové řešení [@glim_gr]
Rekurzivní hladové řešení postupuje následovně:
Z řady aut odstraní první auto a auto stejného typu.
Zbytek aut rekurzivně obarví a pak do řady přidá odebranou dvojici tak,
aby byl počet změn co nejmenší možný.
Střední hodnotu počtu změn barev tohoto algoritmu lze odhadnout pomocí:
${2\over 5}\ n - {8\over 15} \le \mathbb{E}_n({\rm rg}) \le {2\over 5}\ n + {7\over 10}$.
```python {c=plotly}
import g, data_lib
fig = g.intro_graph("rg")
```
## Hvězdičkové rekurzivní řešení [@docwp]
Při běhu rekurzivního hladového řešení se občas stane,
že některé dvojice aut lze obarvit oběma způsoby a počet
změn bude stále stejný. Takovéto auta označíme novou barvou `*`.
Přitom budeme dodržovat invariant, že `*` nikdy není na okrajích a nejsou dvě vedle sebe.
Pokud budeme dále přidávat auto poblíž `*` a jedno nastavení `*` na barvu bude lepší, tak `*` přenastavíme na danou barvu.
Existuje domněnka, která tvrdí, že $\mathbb{E}_n({\rm srg}) \le 0.361 n$.
```python {c=plotly}
import g, data_lib
fig = g.intro_graph("rsg")
```
## Řešení pomocí semidefinitního programování (Jiří Kalvoda, doposud nepublikováno)
Řešení se částečně podobá aproximačnímu algoritmu na maximální řez pomocí semidefinitního programování [@semidef].
V něm si můžeme představit vrcholy jako body $n-1$ rozměrné koule v $n$ dimenzionálním prostoru
a semidefinitní programování nám je rozmístí tak,
aby součet skalárních součinů bodů spojených hranou
byl co možná největší.
Pak náhodně vybereme nadrovinu procházející počátkem a vrcholy rozdělíme do partit podle toho na které straně nadroviny je daný bod.
V našem problému také necháme rozmisťovat body po $n-1$ rozměrné sféře.
Body nyní budou reprezentovat typy aut. Přesněji řečeno každý
bod bude reprezentovat první auto daného typu a druhé
auto daného typu bude reprezentovat virtuální bod k němu
středově symetrický dle počátku.
Když pak tedy rozřízneme body nadrovinou (která žádným bodem neprochází)
a autům přidělíme barvy dle toho, na které straně skončí, bude od každého typu právě jedno auto každé barvy.
Semidefinitní programování se bude snažit maximalizovat
součet skalárních součinů sousedních aut,
tedy sousední auta budou blízko, což znamená, že je malá pravděpodobnost,
že budou mít různou barvu.
Pro účely této práce vždy vyzkoušíme 10 náhodných řezů koule nadrovinu a z nich vybereme nejlepší možný.
Je známo, že tento algoritmus na každé instanci ve střední hodnotě přes nadroviny vrátí řešení, co je nejhůře o $0.22n$ horší než optimum na daném vstupu.
Navíc pomocí hodnoty optimalizační funkce semidefinitního programu
pro každou instanci můžeme získat i dolní odhad na skóre optima.
```python {c=plotly}
import g, data_lib, plotly.subplots
data = data_lib.Data("log-intr", validate_versions=False)
fig = plotly.subplots.make_subplots(rows=2, cols=1, row_heights=[0.8, 0.3], shared_xaxes=True, vertical_spacing = 0.05)
g.draw_algo_graph(fig, data, "semidef_prog(10)", 200, name="Řešení")
g.draw_algo_graph(fig, data, "semidef_prog(10)", 200, val_getter=lambda x: int(x.data["lower_bound"]+0.99), name="Dolní odhad", color="green")
```
# Statistická práce
Na níže uvedeném grafu jsou výsledky všech výše uvedených algoritmů pohromadě (pozor na to, že v grafu je zobrazena jen zajímavá část vertikální osy mezi $0.25 n$ a $0.6n$.
```python {c=plotly}
import g, data_lib, plotly.subplots
data = data_lib.Data("log-intr", validate_versions=False)
fig = plotly.subplots.make_subplots(rows=2, cols=1, row_heights=[0.8, 0.5], shared_xaxes=True, vertical_spacing = 0.05)
g.draw_algo_graph(fig, data, "greedy", 200)
g.draw_algo_graph(fig, data, "rg", 200)
g.draw_algo_graph(fig, data, "rsg", 200)
g.draw_algo_graph(fig, data, "semidef_prog(10)", 200)
fig.update_xaxes(range=[0.25*200, 0.6*200])
```
Z grafu můžeme odhadnout, že jako nejlepší algoritmus se jeví semidefinitní programování a po něm hvězdičkové rekurzivní řešení.
V následující části bych rád ukázal, že tomu tak skutečně je.
To vypadá jako jednoduchá úloha pro [Welchův t-test](https://en.wikipedia.org/wiki/Welch%27s_t-test).
Vygeneruji dalších 30 běhů každého algoritmu a na nich spustím:
```python {c=code_and_output highlight=True}
import data_lib
import scipy.stats
data = data_lib.Data("log-t-test", validate_versions=False)
rsg = [ x.score for x in data.pipelines["rsg"] ]
semidef_prog = [ x.score for x in data.pipelines["semidef_prog(10)"] ]
print(scipy.stats.ttest_ind(rsg, semidef_prog, equal_var=False))
```
Z tohoto bychom mohli usuzovat, že řešení pomocí semidefinitního programování je na $99\%$ lepší,
protože `pvalue` je dokonce o několik řádů menší než $0.01$.
V čem je tedy problém?
Problém t-testu je že nesamplovaná data musí pocházet z normální distribuce
a já vůbec netuším, jak bych dokazoval, že naše data se jí alespoň blíží (pokud to vůbec platí).
Pojďme tedy využít trošku více dřevorubecké řešení, kterým to dokážeme.
Budeme se snažit dokázat, že algoritmus pomocí semidefinitního programování je lepší než hvězdičkové rekurzivní řešení.
Tedy:\
$H_0:$ Hvězdičkové řešení je stejně dobré nebo dokonce lepší.\
$H_1:$ Řešení pomocí semidefinitního programování je lepší.
Víme, že rozptyl náhodné veličiny skóre řešení je nejvýše $n$, protože se jedná o náhodnou veličinu z intervalu $0$ až $2n$.
Když víme, že $\sigma^2(X) \le n$, tak $\sigma^2(X_1 + \cdots + X_k) \le kn$
(pro nezávislé $X_1,\dots,X_n$, což v našem případě jsou, protože vybíráme vstupy nezávisle).
Tedy $\sigma^2(\overline{X_k}) = \sigma^2({X_1 + \cdots + X_k \over k}) \le {kn \over k^2} = {n \over k}$.
Dle Cebyševovy nerovnosti tedy (pro $t>0$):
$$P\left(\mathbb{E}(X) \ge \overline{X_k} + t {n \over k}\right) \le {1 \over 2 t^2}$$
Což můžeme využít následovně:
Nechť $X$ a $Y$ jsou náhodné veličiny skóre dvou algoritmů (semidefinitního programování a hvězdičkového rekurzivního), kde naměříme $\overline{X_k} < \overline{Y_k}$.
$$
P(H_0) = P\left( \mathbb{E}(X) \ge \mathbb{E}(Y) \right)
= 1 - P\left( \mathbb{E}(X) < \mathbb{E}(Y) \right)
\le 1 - P\left( \mathbb{E}(X) < { \overline{X_k}+\overline{Y_k} \over 2} < \mathbb{E}(Y) \right)
\le
$$
$$
\le 1 - P\left( \mathbb{E}(X) < { \overline{X_k}+\overline{Y_k} \over 2} \right) + 1 - P\left( { \overline{X_k}+\overline{Y_k} \over 2} < \mathbb{E}(Y) \right)
\le P\left( \mathbb{E}(X) \ge \overline{X_k} + { \overline{Y_k}-\overline{X_k} \over 2} \right) + P\left( \overline{Y_k} - { \overline{Y_k}-\overline{X_k} \over 2} \ge \mathbb{E}(Y) \right)
$$
Nyní použijeme výše uvedenou nerovnost pro $t{n\over k} = { \overline{Y_k}-\overline{X_k} \over 2}$, tedy $t = k {\overline{Y_k}-\overline{X_k} \over 2n}$.
Symetrická nerovnost funguje i pro $Y$ a otočený směr odhadu.
$$
P\left( H_0 \right) \le
2 {1 \over 2 \left({k(\overline{Y_k}-\overline{X_k}) \over 2n}\right)^2} =
{4n^2 \over k^2 \left({\overline{Y_k}-\overline{X_k} }\right)^2} =
$$
Pokud chceme, abychom věděli, že algoritmus $X$ je lepší než algoritmus $Y$ (pro danou velikost vstupu) s jistotou $99\%$, musí platit, že:
$$
{4n^2 \over k^2\left({\overline{Y_k}-\overline{X_k} }\right)^2} < 0.01
$$
Z výše uvedených grafů můžeme doufat, že $\overline{Y_k}-\overline{X_k}$ bude alespoň zhruba $10$.
Pro $k = 500$ by tedy mohla být výsledná hodnota mohla dost malá.
Navrhneme tedy následující pokus:
Spustíme oba dva algoritmy na $500$ náhodně vybraných vstupů a spočítáme
jejich průměry.
Pokud výsledek výše uvedené formule bude menší než $0.01$, můžeme zamítnout $H_0$ s pravděpodobností $99\%$.
```python {c=code_and_output highlight=True}
import data_lib
import scipy.stats
n = 200
data = data_lib.Data("log-correct-test", validate_versions=False)
rsg = [ x.score for x in data.pipelines["rsg"] ]
semidef_prog = [ x.score for x in data.pipelines["semidef_prog(10)"] ]
k = len(rsg)
assert k == len(semidef_prog)
rsg_avg = sum(rsg) / k
semidef_prog_avg = sum(semidef_prog) / k
print(f"n: {n} k: {k}")
print("rsg: ", rsg_avg)
print("semidef_prog:", semidef_prog_avg)
assert semidef_prog_avg < rsg_avg
print(4*n**2 / k**2 / (rsg_avg - semidef_prog_avg)**2)
```
Nulovou hypotézu se nám tedy podařilo zamítnout a tudíž řešení pomocí semidefinitního programování je lepší.
Všimněte si, že pokud by platilo, že se jedná o normální rozdělení, tak nám stačí mnohem méně dat a máme mnohem větší jistotu.
References
==========
::: {#refs}
:::
#!/bin/bash
set -ueo pipefail
./build.py "$@"
rsync --recursive --info=progress2 --info=name1 --stats --delete out/ jirikalvoda@kam.mff.cuni.cz:WWW/past_zapoctak
File added
File added
File added
File added
File added
File added
File added
File added
File added
File added
File added
File added
File added
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment