AI tutorial

Minmax – Schopenhauer wśród algorytmów

Rzadko kiedy mamy komfort braku przeciwników działających na naszą niekorzyść. A* spisywał się znakomicie w pustym labiryncie, ale co kiedy znajdziemy tam wrogie nam istoty? Albo jeśli świat gry jest niedeterministyczny? Zanim to, zastanówmy się nad prostszą grą – kółko i krzyżyk.

Zasady są wszystkim doskonale znane. Dwójka graczy na zmianę wykonuje swoje ruchy. Określony układ X i O na planszy nazwiemy stanem. W każdym stanie są dozwolone pewne akcje – różne dla różnych graczy (postawienie X oraz O na tym samym polu traktujemy jako dwie różne akcje). Niektóre stany są stanami końcowymi, tzn. takimi dla których nie ma dozwolonych akcji. Oczywiście są to sytuacje gdzie cała plansza jest już zapełniona oraz stany w których jeden z graczy odniósł zwycięstwo.

W jaki sposób należy podejść do problemu? Każdej naszej akcji będzie towarzyszyć reakcja ze strony przeciwnika, musimy ją przewidzieć i uwzględnić w procesie myślowym. Tradycyjnie można stworzyć drzewo. Korzeniem jest stan początkowy kiedy żaden ruch nie został wykonany. Załóżmy, że zaczyna gracz X. Na poziomie pierwszym znajdą się wszystkie możliwe stany po wykonaniu przez niego ruchu. Na poziomie drugim kontrolę przejmuje gracz O. I tak na zmianę aż dotrzemy do stanów końcowych.

Stany końcowe charakteryzują się posiadaniem pewnej liczby, którą będę nazywać oceną lub wartością. Mówi ona jak pożądany przez gracza jest taki stan. Łatwo zauważyć, że to co jest pożądane przez jednego gracza, będzie niepożądane przez drugiego. Niech wygrana gracza X będzie oznaczona jako +1, remis jako 0, a wygrana gracza O przez -1. Gracz X będzie dążył do maksymalizowania (max) a gracz O do minimalizowania (min).

Próbujemy podjąć decyzje, które doprowadzą do zmaksymalizowania wyniku (rozważamy przypadek gracza max). Spójrzmy na ilustrację przykładowego problemu (nie kółka i krzyżyka). Trójkąt skierowany w górę symbolizuje gracza max, w dół – min.

minmax1Analizę zacznijmy od lewego dolnego rogu. W tym punkcie gracz min ma do wyboru dwie decyzje. Jedna doprowadzi do końca gry z wynikiem +5, druga +7. Jako że chce zminimalizować ten wynik, wybierze +5. Po przejściu przez cały poziom wartości prezentują się następująco:

minmax2Poziom wyżej decyzje podejmuje max. Wie, że niezależnie co zrobi, przeciwnik zagra na jego niekorzyść. Spośród dostępnych opcji wybiera tą która będzie pomimo działań nieprzyjaciela doprowadzi do lepszego (z jego punktu widzenia) stanu. W ten sposób wypełniamy całe drzewo.

minmax3

Ścieżka pokazana na czerwono pokazuje ciąg decyzji. Max wybiera lewy węzeł. W odpowiedzi na to min wybierze prawy i tak do liścia z wartością. W ten sposób max uwzględnia co zrobi min aby mu przeszkodzić. Trudno jest tutaj mówić o maksymalizowaniu zysku, bliższe prawdy będzie minimalizowanie strat.

Całość da się zgrabnie zapisać w formie rekurencyjnej (nie jest tutaj uwzględnione zapamiętywanie decyzji do podjęcia – zależy to mocno od konkretnej implementacji.

minmax(stan):
1. jeśli stan jest stanem końcowym zwróć jego wartość.
2. Jeśli stan jest węzłem max
2.1. wywołaj na rzecz każdego potomka minmax, jako węzeł min
2.2. zwróć wartość maksymalną
3. w przeciwnym wypadku (węzeł min)
3.1. wywołaj na rzecz kazdego potomka minmax, jako węzeł max
3.2. zwróć wartość minimalną

Implementacja

Cały kod jest dostępny w repozytorium https://github.com/Xevaquor/aipac Przykład jest dla gry w kółko i krzyżyk, dla innych gier idea się nie zmienia.

Algorytm korzysta z obiektów state i game. Dostarczają one potrzebnych informacji na temat świata gry. Są to:
game.get_legal_moves – zwraca listę dozwolonych akcji w danym stanie gry
game.apply_move – zwraca stan po wykonaniu danej akcji
game.has_won/has_lose – informuje czy w danym stanie nastąpiło zwycięstwo/porażka gracza
game.is_terminal – czy dany stan jest stanem końcowym

Funkcja oceniająca liść:

def score(self, state, game):
	"""
	Computes state score.
	:type state: tictactoe.State
	:type game: tictactoe.TicTacToe
	:return: Score value, None if state is not terminal. (Is not leaf in tree)
	"""
	if game.has_won(state, self.symbol):
		return 1
	elif game.has_lose(state, self.symbol):
		return -1
	elif game.is_terminal(state):
		return 0
	else:
		return None

Operacje na węźle min (dla max jest analogicznie). Obliczamy wartości wszystkich potomków i wybieramy tego z najmniejszą.

def min(self, state, game):
	"""
	Processes min-node.
	:type state: tictactoe.State
	:type game: tictactoe.TicTacToe
	:return: tuple, node value, action taken
	"""
	# get legal moves for given state
	legal_moves = game.get_legal_moves(state, self.opponent_symbol)
	# compute score for current state
	state_score = self.score(state, game)
	# if state is unknown (None), continue. Return score otherwise.
	if state_score is not None:
		return state_score, None

	children = []
	# for every legal move
	for move in legal_moves:
		# get state after executing action
		next_state = game.apply_move(state, move)
		# get value of child node (for children of min node, it is a maximum).
		# action is irrelevant
		value, dumb = self.max(next_state, game)
		children.append((value, move))

	# select minimal successor
	return min(children)

I ostatecznie funkcja podejmująca decyzję. Przyjęliśmy że jesteśmy graczem max.

def make_move(self, state, game):
	"""
	Selects action to make in given state. Assumes max player.
	:type state: tictactoe.State
	:type game: tictactoe.TicTacToe
	:return: action to take
	"""
	score, move = self.max(state, game)
	return move

Własności

Minmax ma kilka ciekawych własności. Pierwszą z nich jest przeczesywanie całego drzewa. Nawet przy tak prostej grze decyzja jest podejmowana we względnie długim czasie. Proces można przyśpieszyć kończąc wywołania na pewnym poziomie zagnieżdżenia. Wtedy zamiast schodzić dalej, szacujemy wartość aktualnego stanu. Oczywiście za kolejnym wywołaniem zejdziemy o jeden poziom niżej poprawiając nasze szacunki, ale w ten sposób niestety tracimy gwarancję dokładnego rozwiązania. Aby ograniczyć wielkość przestrzeni można użyć alfa-beta cięć. Umożliwiają one ucięcie pewnych gałęzi drzewa, o których wiadomo, że na pewno nie wniosą nic nowego. Nie będę tego opisywać, gdyż nie przynosi to spektakularnej poprawy wydajności.

Drugą sprawą jest skrajny pesymizm (stąd tytuł posta). Weźmy takie drzewo:

minmax4

Wybierając lewą ścieżkę możliwe wyniki to +7 lub +6. Dla prawej są to +5 i +100. Minmax wybierze lewą, gdyż pesymistycznie założy optymalną grę przeciwnika, który nie pozwoli zdobyć +100. Jeżeli przeciwnik nie jest idealny mógłby się pomylić i pozwolić nam na pogrom. Gdybyśmy poszli w prawo, w najgorszym wypadku dostaniemy +5. Minmax pójdzie w lewo uzyskując +6 lub +7 (w zależności od decyzji przeciwnika). Cholera, może warto zaryzykować? Mamy w najgorszym razie 2 do stracenia (7 – 5), ale możemy zyskać 100 – 5 = 95. Tym problemem zajmiemy się w następnej kolejności – nie wiemy jak zachowają się wrogowie w labiryncie, przypuszczamy również, że nie zachowują się oni optymalnie.

Minmax można rozszerzyć dla większej ilości graczy. Np. jeden gracz max, dwóch graczy min (wtedy poziomy drzewa kształtują się: MAX-MIN-MIN-MAX-MIN-MIN…). Każdy z graczy może również maksymalizować/minimalizować inną wartość we wzajemnie powiązanej krotce danych. Przykład oraz ciekawy skutek uboczny – następnym razem.

Podsumowując: jeśli gra jest deterministyczna oraz wiemy że przeciwnik jest idealny, minmax pozwoli zachować się najlepiej jak to możliwe. Jeżeli nie – przygotuje się na najgorsze.

Reklamy

Chodź, pomaluj mój świat – kolorowanie mapy (grafu)

Ostatnim razem opracowaliśmy reprezentację niezbyt ogólną. Zanim zaczniemy ulepszać algorytm poprawmy implementację. Zamiast używać zmiennych i ograniczeń zamodelujemy całość w postaci grafu.

geo-map-australia-and-new-zealand-contour

Weźmy następujący problem. (Przykład pochodzi z książki Artificial Intelligence: A Modern Approach).
Mamy mapę Australii podzieloną na okręgi administracyjne. Należy ją pokolorować, oczywiście tak aby żadne dwa sąsiadujące nie używały tego samego koloru. Zgodnie z ostatnio używanymi oznaczeniami mamy 7 zmiennych:
\text{WA}, \text{NT}, \text{ST}, \text{Q}, \text{NSW}, \text{V}, \text{T} \in \{ \text{R}, \text{G}, \text{B} \}
gdzie R, G i B są kolorami. Ograniczenia prezentują się następująco:
\text{WA} \neq \text{ST} \newline  \text{NT} \neq \text{ST} \newline  \text{NSW} \neq \text{ST} \newline  \text{V} \neq \text{ST} \newline  \text{NT} \neq \text{Q} \newline  \text{NSA} \neq \text{V}\newline  \text{NSA} \neq \text{Q} \newline  \text{WE} \neq \text{NT}

Całość problemu można przedstawić jako graf. Zmienne będą wierzchołkami, ograniczenia – krawędziami między węzłami oznaczającymi związane zmienne.aus

Skutkiem ubocznym takiej reprezentacji jest natychmiastowe zauważenie, że Tasmania (T) nie ma narzuconych żadnych ograniczeń. Jeśli taki graf nie jest spójny (czyli da się go „rozdzielić” na kilka części) oznacza to, że poszczególne „kawałki” (składowe spójne) są osobnymi podproblemami, niezależnymi od siebie. Bardzo pożądana sytuacja.

Zaraz, wszystko super kiedy w ramach jednego ograniczenia występują dwie zmienne. A co jeśli jest ich więcej? Wtedy wprowadzamy inny typ wierzchołków do oznaczenia więzów. Dla przykładu z poprzedniego posta (kwadraty magiczne, wersja 3×3 aby było mniej bazgrołów):

magicsquaregraph3x3

Kolory nie wskazują na nic szczególnego, służą jedynie podniesieniu czytelności. Kwadraty są węzłami symbolizującymi ograniczenia miedzy zmiennymi.

Uaktualnijmy implementację.

Na początek ponumerujmy wierzchołki od zera do 6 (włącznie) w kolejności: WA, NT, ST, Q, NSW, V, T.

Wprowadzamy dziedzinę zmiennych jako obiekty klasy Color

class Color:
    def __init__(self):
        pass

    Red, Green, Blue = range(3)

Zmienne przyjmują teraz wartości z dziedziny {Red, Green, Blue} przez co zmienia się konstruktor klasy Variable.

class Variable(object):
    def __init__(self, fixed_value=None):
        self.domain =  [Color.Red, Color.Green, Color.Blue] if fixed_value is None else [fixed_value]
        self._value = fixed_value
	#...

W klasie Color będącej odpowiednikiem poprzedniej MagicSquare definiujemy graf. Robię to poprzez listę krawędzi, gdyż tak będzie mi później bardzo wygodnie.

class Color(object):
    def __init__(self):
        self.constraints = [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3), (2, 4), (2, 5), (3, 4), (4, 5)]
	#...

Modyfikuję początkowe przypisanie tak aby składało się z samych wartości None

def get_initial_assignment(self):
        return Assingment(self.get_variables_from_values([None] * 7))

Dodaję drobną metodę sprawdzającą czy dla pojedynczej krawędzi są spełnione warunki. Jeśli zmienne są sobie równe to nastąpiło naruszenie, jeżeli nie wszystko jest OK. Oczywiście sytuacja gdzie obie są równe None jest jak najbardziej dozwolona.

def edge_invalid(self, a, b):
	if not a.is_assigned or not b.is_assigned:
		return False

	return a.value == b.value

Na koniec uaktualniam is_violating_constraints

def is_violating_constraints(self, assignment):
	# alias co by mniej pisać 😉
	v = assignment.variables

	for a, b in self.constraints:
		if self.edge_invalid(v[a],v[b]):
			return True

        return False

Wynik wykonania:

geo-map-australia-and-new-zealand-contour

Mamy teraz bazę, gotową do dalszego szlifowania, o czym już wkrótce.

Całość kodu

Obrazek tytułowy pochodzi ze slajdów do kursu CS 188: Artificial Intelligence prowadzonego przez UC Berkeley.

A Kind of Magic – problemy spełnialności więzów (Constraint Satisfaction Problems) – część 1

Nie zawsze sposób osiągnięcia rozwiązania jest ważny, czasem liczy się tylko poprawna odpowiedź. Przykładem takiego problemu może być sudoku – trzeba wymyślić co wpisać w wolne pola, nie jest istotne w jakiej kolejności i w jaki sposób. Innym przykładem może być kolorowanie mapy – żadne dwa sąsiadujące kraje nie mogą używać tego samego koloru. Podobnie jak ustalenie planu zajęć: wykładowca i student może być tylko w jednym miejscu na raz, sala nie może mieścić dwóch wykładów jednocześnie, niektóre zajęcia wymagają określonego sprzętu etc. Tego typu problemy określa się tytułowym mianem Constraint Satisfaction Problems (CSP) lub w języku Mickiewicza: Problemy spełnialności więzów.

Jako przykładowy problem weźmy kwadraty magiczne. Chyba wszyscy spotkali się z nimi na wczesnym etapie edukacji, ale dla przypomnienia: jest to tabela (macierz) NxN z liczbami posiadająca następująca własność: suma każdego wiersza jest taka sama jak suma każdej kolumny jak również sumy obu przekątnych.

180px-Magicsquareexample.svgFormalizacja

Sformalizujmy nieco zagadnienie (przykład dla kwadratu 4×4 o sumie równej 34).
\text{Magic Square} =  \begin{pmatrix}  v_{0} & v_{1} & v_{2} & v_{3} \\  v_{4} & v_{5} & v_{6} & v_{7} \\  v_{8} & v_{9} & v_{10} & v_{11} \\  v_{12} & v_{13} & v_{14} & v_{15}  \end{pmatrix}
Mamy 16 zmiennych – po jednej na każde pole. Niektóre z nich mogą mieć już na początku ustalone wartości (mała dygresja: skoro są już ustalone to nie muszą być zmiennymi, ale takie podejście ułatwi w tym wypadku implementację).
V = (v_0, v_1, \dotsc v_{14}, v_{15})
Każda ze zmiennych przyjmuje wartości z określonej dziedziny. Ogólnie dziedzina może być dowolnym niepustym zbiorem, także ciągłym lub nieskończonym. W rozważanym przypadku jest to skończony zbiór liczb naturalnych powiedzmy od 1 do 20. (Dalsze rozważania są prowadzone dla problemów które mają skończoną dziedzinę). Każda zmienna może mieć różną dziedzinę – w szczególności gdy jej wartość jest ustalona może być to zbiór jednoelementowy.
\forall_{v \in V} \in \{1 \dotsc 20\}
Pomiędzy zmiennymi zachodzą pewne ograniczenia (więzy) mówiące o tym jak pewne zmienne mają się do siebie. Tutaj ograniczenia mówią o sumach poszczególnych wierszy/kolumn/przekątnych.
v_0 + v_1 + v_2 + v_3 = 34 \newline  v_4 + v_5 + v_6 + v_7 = 34 \newline  \vdots

Przypisaniem nazwiemy pewne przyporządkowanie (podstawienie) wartości do poszczególnych zmiennych (niekoniecznie wszystkich na raz). Przypisanie które spełnia ograniczenia będę nazywać spójnym.
Przypisanie spójne oraz kompletne (uwzględniające wszystkie zmienne) to rozwiązanie.

Jeśli ktoś pisał w Prologu to można nazwać to mniej więcej deklarowaniem ograniczeń, które wbudowany mechanizm wnioskowania będzie próbował spełnić 😉

Problem spełnialności więzów.

Tak sformułowane zagadnienie można określić jako Constraint Satisfaction Problem albo problem spełnialności więzów.

Zastanówmy się nad rozwiązaniem. Tradycyjnie zaczniemy od najprostszego a z czasem będziemy je ulepszać.

Brute force

Najprościej było by wygenerować wszystkie możliwe kombinacje. Można się jednak domyślić, że sprawdzenie tego zajmie wieki, ponadto nie ma to zbyt wiele wspólnego z jakąkolwiek inteligencją.

Backtracking search (przeszukiwanie z nawrotami)

Pomysł najlepiej będzie przedstawić na przykładzie. Załóżmy, że mamy 3 zmienne a,b,c \in \{1,2,3\} przy następujących ograniczeniach: a \ne b, b \ne c, a \ne c

csp1

Zaczynamy w korzeniu (z pustym przypisaniem). Schodzimy w dół (najpierw po lewej). Podstawienie a = 1 spełnia ograniczenia. Idziemy dalej otrzymując a = 1, b = 1. Nie spełnia ono warunków przez co cofamy się (wykonujemy nawrót). Bierzemy kolejnego potomka wierzchołka a = 1. Jest to b = 2. Podstawienie a = 1, b = 2 jest spójne. Przechodzimy w dół do węzła c = 1. Nie spełnia on warunków więc idziemy do c = 2. Tutaj podobnie: trafiamy w końcu na c = 3. Podstawienie a = 1, b = 2, c = 3 jest spójne i kompletne przez co jest rozwiązaniem.

Węzły oznaczone na czerwono to miejsca nawrotów, szare – wierzchołki nieodwiedzone.

Rekurencyjnie algorytm można zapisać w postaci:

Backtracking Search:
 1. Jeśli przypisanie jest rozwiązaniem, zakończ.
 2. v = pobierz nieprzypisaną zmienną
 3. Dla każdej wartości val z dziedziny zmiennej v:
 3.1. v = val
 3.2. Jeśli przypisanie nie spełnia ograniczeń:
 3.2.1 Oznacz v jako nieprzypisaną (lub równoważnie usuń z przypisania)
 3.2.2. Wróć do 4. (continue)
 3.3. wynik = wywołaj rekurencyjnie Backtracking Search
 3.4 Jeśli wynik jest różny od "rozwiązanie nie istnieje" zwróć wynik
 3.5 W przeciwnym wypadku ustaw v jako nieprzypisaną.
 3.6 Zwróć "rozwiązanie nie istnieje"

Implementacja

Zacznijmy od klasy reprezentującej zmienną. Przechowuje ona jej wartość lub informację, że jest nieustalona (_value = None) oraz domenę (domain). Jeśli jej wartość jest ustalona na stałe (fixed_value != None) jest to ustalane już na etapie konstrukcji obiektu. Metody __radd__, __repr__ oraz getter i setter dla value nie są istotne dla algorytmu – to szczegół implementacyjny.

class Variable(object):
    def __init__(self, fixed_value=None):
        # załpżono, że dziedzina dla wszystkich zmiennych jest taka sama [1..20]
        # ale równie dobrze może być przekazywana do konstruktora
        self.domain =  range(1,21) if fixed_value is None else [fixed_value]
        self._value = fixed_value

    # pozwala skorzystać z wbudowanej funkcji sum()
    def __radd__(self, other):
        if other is None: raise ValueError("Cannot sum None variable")
        return self.value + other

    # getter i setter dla value, rzuca wyjąktiem przy próbie przypisania
    # wartości spoza dziedziny
    @property
    def value(self):
        return self._value

    @value.setter
    def value(self, val):
        if val != None and not val in self.domain:
            raise DomainException(str(val))
        self._value = val

    def __repr__(self):
        return str(self.value)

    # informacja czy jest przypisana, tj czy value != None
    @property
    def is_assigned(self):
        return self.value != None

Klasa Assignment reprezentuje przypisanie, przechowuje listę wszystkich zmiennych. Pozwala sprawdzić kompletność oraz pobrać następną zmienną bez wartości.

class Assingment(object):
    # vars to lista zmiennych w problemie
    def __init__(self, vars):
        self.variables = vars

    # czy wszystkie zmienne mają wartość?
    def is_complete(self):
        not_assigned = [x for x in self.variables if x.value == None]
        return not not_assigned

    # weź pierwszą nieprzypisaną, jeśli takich nie ma rzuć wyjątek
    def get_unassigned_variable(self):
        for v in self.variables:
            if not v.is_assigned:
                return v
        raise AssertionError()

    # pokaż na ekranie w sensownym formacie
    def pretty_print(self):
        for i in range(0,4):
            print ' '.join(str(self.variables[i*4:i*4+4]))

Główna klasa problemu, pozwala odseparować ogólny algorytm rozwiązywania od szczegółów zagadnienia, co pozwoli łatwo korzystać z kodu dla różnych problemów. W podanej implementacji zakładamy znajomość sumy, jeżeli nie jest ona dostępna należy zmienić implementację funkcji is_violating_constraints tak aby nie przyrównywała do self.summary, ale sprawdzała czy wszystkie dają taki sam wynik.

    # sprawdza czy podstawienie narusza ograniczenia
    def is_violating_constraints(self, assigment):
        #alias coby mniej pisać 😉
        v = assigment.variables

        # spawdzenie wierszy
        for row_index in range(0,self.square_size):
            if not self.row_satisfies_constraint(row_index, v, self.summary):
                return True

        # i kolumn
        for col_index in range(0,self.square_size):
            if not self.col_satisfies_constraint(col_index, v, self.summary):
                return True

        # oraz przekątnych
        # tutaj oczywiście można zrobić to tego funkcje jak poprzednio ale już sobie odpuściłem

        if not assigment.is_complete():
            return False

        if v[0].value + v[5].value + v[10].value + v[15].value != self.summary:
            return True

        if v[3].value + v[6].value + v[9].value + v[12].value != self.summary:
            return True

        # sprawdziliśmy wszystkie warunki - dane przypisanie jest OK        
        return False

Na koniec właściwy algorytm. Warto zwrócić uwagę na fakt, że zmienna assignment jest pośrednio modyfikowana przy zmianach zmiennej var.

def backtracking_search(assignment, instance):
    # jeśli mamy rozwiązanie lub takowego nie ma - kończymy
    if assignment.is_complete() and not instance.is_violating_constraints(assignment):
        return assignment

    # pobierz zmienną bez wartości
    var = assignment.get_unassigned_variable()

    # dla wszystkich wartości z domeny zmiennej
    for val in var.domain:
        # przypisz tą wartość
        var.value = val
        # i sprawdź czy przypisanie spowodowało naruszenie więzów.
        if instance.is_violating_constraints(assignment):
            # jeśli tak cofnij przypisanie i spróbuj z następną wartością
            var.value = None
            continue
        # przypisanie było poprawne, powtórz dka kolejnej zmiennej (zejscie w dół drzewa)
        result = backtracking_search(assignment, instance)
        # mamy wynik, wróć w górę
        if result is not None:
            return result
        # niżej nie udało się spełnić ograniczeń, wyzeruj podstawienie
        else:
            var.value = None

    # nawrót
    return None

Przykładowe użycie:

magic = MagicSquare()
initial = magic.get_initial_assignment()
res = backtracking_search(initial, magic)
if res is None:
	print "Nope, nie da się."
else:
	res.pretty_print()

Sprawdźmy jak poradzi sobie z przykładowym problemem. S=34, D = \{1 \dotsc 20 \}

task

Rozwiązanie

solvedPodsumowanie

Zaprezentowany algorytm zawsze znajduje rozwiązanie (o ile oczywiście istnieje), niemniej można go znacząco ulepszyć. W następnej części sprawdzimy co jest nie tak (jak można się domyślić: wydajność) i zajmiemy się tym.

Cały kod jest dostępny na GitHubie.

Gwiazda wieczoru – algorytm A* (A star)

Algorytm UCS był przez nas stosowany ostatnio z dużym powodzeniem. Dlaczego więc mamy z niego rezygnować i szukać czegoś innego? Ponieważ nie jest on zbyt inteligentny. Spójrzmy na trywialną instancję problemu labiryntu.

00000Człowiek od razu zauważy, że wystarczy iść w prawo. A co zobaczy UCS?

00000

Czerwone kropki oznaczają rozwinięte wierzchołki, czyli te które były brane pod uwagę jako potencjalnie prowadzące do celu.

Czy nie można było po prostu iść prosto do celu zamiast marnować czas na szukanie dookoła a nawet wstecz? Poradziliśmy sobie lepiej niż komputer ponieważ intuicyjnie czujemy, że idąc w lewo oddalamy się od celu. Jak teraz przełożyć naszą intuicję na kod?

Heurystyka

Aby komputer wiedział w który kierunek poszukiwań jest dobry, a który zły potrzebujemy pewnej funkcji która będzie szacować odległość stanu od rozwiązania. Taką funkcję nazywamy heurystyką. Idealna heurystyka zwraca wartość będącą dokładnym kosztem przejścia z bieżącego stanu do rozwiązania. Pozwala to na bezbłędną pracę algorytmu, bez rozwijania zbędnych wierzchołków. Oczywiście uzyskanie dokładnej wartości może być skomplikowane i czasochłonne, przez co w praktyce niemożliwe do zastosowania.
Okazuje się że do znalezienia optymalnego rozwiązania wystarczy aby była ona dopuszczalna oraz spójna. Co oznaczają te pojęcia? Intuicyjnie, spójność oznacza, że funkcja nigdy nie zwraca wartości większej niż faktyczny koszt przejścia z danego wierzchołka do wyjścia. Formalnie
\forall_{v \in V} h(v) \le C(v)
gdzie V to zbiór wszystkich wierzchołków (stanów), h to heurystyka, a C to rzeczywisty koszt. Spójność natomiast można kojarzyć z nierównością trójkąta. Spójrzmy na schemat:

consistencyJeżeli zachodzi:
h(v) \le C(v,u) + h(u)
wtedy heurystyka jest spójna. Jak można to interpretować? Poprzedni warunek gwarantował, że „globalnie” nie przeszacujemy kosztu do celu. Spójność gwarantuje że to samo dzieje się lokalnie (dla każdego wierzchołka). Skoro h(v) podał nam pewną wartość oznacza to, że nie istnieje żadna droga o mniejszym koszcie – nie tylko bezpośrednio v->t ale również przez inne wierzchołki np. v->u->t. Dodatkowo spójność implikuje dopuszczalność

Wróćmy do problemu labiryntu (bez monet) i zastanówmy się nad heurystyką. Niezależnie od układu ścian nie da się przejść drogą krótszą niż tą wyznaczoną linią prostą (po prostej też może się nie dać ale jest to OK – warunek spójności jest spełniony). Taka heurystyka prezentuje się następująco (v to wierzchołek dla którego szacujemy koszt, g to wierzchołek docelowy).
h(v, g) = \sqrt{(v_x - g_x)^2 + (g_x - v_x)^2}
Jeśli poruszamy się jedynie w kierunkach NSWE możemy pójść o krok dalej i zastosować metrykę taksówkową.
h(v, g) = |v_x - g_x| + |v_y - g_y|

Algorytm A*

Mamy już sposób oceny potencjału wierzchołka. Jak to wykorzystać? Wystarczy do priorytetu dodać wartość heurystyki, to sprawi że bardziej obiecujące wierzchołki będą rozwijane wcześniej. Taka modyfikacja nosi nazwę algorytmu A*. Co ważne przy spójnej heurystyce nadal mamy gwarancję uzyskania optymalnego rozwiązania. Porównanie do UCS można przedstawić następująco:

astar vs ucs

A* przeszukuje przestrzeń stanów szybciej w kierunku rozwiązania niż w innych kierunkach. Sprawdźmy:

EsF688

so_good

Idealnie…

Implementacja

W związku z powyższym implementacja będzie wymagała jedynie niewielkich zmian.
Na początek funkcja heurystyczna

def heuristics(node, instance):
    node_pos = (node[0], node[1])
    target_pos = (instance.get_target_state()[0], instance.get_target_state()[1])
    return taxi(node_pos, target_pos)


def taxi(p1, p2):
    return abs(p1[0] - p2[0]) + abs(p1[1] - p2[1])

W obiekcie instance pojawiła się metoda get_target_state() która zwraca stan końcowy.

def astar(instance, h):
    closed = set()  # zbiór zamknięty
    # inicjujemy zbiór otwarty stanem początkowym. Decyzję ustawiamy na null
    # koszt na 0 - nie ma to znaczenia. Rodzicem jest również null (jest to
    # korzeń drzewa
    # fringe dla UCS jest kolejką priorytetową (niższy priorytet powoduje szybsze zdjęcie
    #elementu
    #enqueue - put
    #dequeue - get
    fringe = PriorityQueue()
    #format wierzchołka to:
    #(priorytet,[(stan, decyzja, koszt), rodzic])
    #jest to wymagane przez kolejkę.
    root_node = Node()
    root_node.parent = None
    root_node.cost = 0
    root_node.step = None
    root_node.state = instance.get_start_state()
    fringe.put((0, root_node))
    #znaleziony cel
    target = None

    while True:
        #jeśli zbiór otwarty jest pusty to nie istnieje droga do celu
        if fringe.empty():
            return []
        #pobierz kolejny węzeł z kolejki - ten o najniższym priorytecie
        #ignorujemy koszt pobrany z kolejki, zamiast niego używamy własności cost węzła
        node = fringe.get()[1]
        node_cost = node.cost

        #jeśli jesteśmy w stanie docelowym, ustaw cel i zakończ
        if instance.is_target_state(node.state):
            target = node
            break

        #jeśli węzeł nie był rozwijany
        if node.state not in closed:
            #dodaj go do zbioru zamkniętego (innymi słowy oznacz jako rozwinięty)
            closed.add(node.state)
            #rozwiń go
            children = instance.get_children(node.state)

            #print node.state, children
            #i dla każdego następnika
            for child in children:
                child_state, child_step, child_cost = child
                #dodaj informację o poprzedniku (node jest rodzicem child)
                #jako koszt ustaw sumę następnika i koszt dojścia do rodzica -
                #został on odczytany przy rozpakowywaniu krotki zwróconej przez
                #fringe.get()
                heuristic_cost = h(child_state, instance) #koszt do wyjścia oszacowany przez heurystykę
                vertex = Node()
                vertex.step = child_step
                vertex.cost = child_cost + node_cost
                vertex.parent = node
                vertex.state = child_state
                new_node = (vertex.cost + heuristic_cost, vertex) #priorytetem jest dotychczasowy
                #koszt + wartość heurystyki.
                #i wrzuć do kolejki (zbioru otwartego)
                fringe.put(new_node)

    #lista decyzji prowadzących do rozwiązania
    solution = []
    #zaczynamy od węzła z wynikiem
    i = target
    #dopóki ma rodzica (nie jesteśmy w korzeniu)
    while i.parent is not None:
        #dodaj decyzję która nas tutaj doprowadziła
        solution.append(i.step)
        #przejdź do rodzica
        i = i.parent
    #podążaliśmy od wyjścia do startu przez co trzeba odwrócić kolejność
    solution.reverse()

    return solution

Na koniec zobaczmy na porównanie UCS oraz A* dla kilku przykładowych problemów. Użyto heurystyki jak w powyższym kodzie.

Bez żadnych zmian możemy zastosować A* do problemu ze zbieraniem monet. (Heurystyka pozostała taka samo, dla lepszych efektów należało by w niej uwzględnić monety pozostałe na planszy). Tutaj niestety trudno jest pokazać dokładnie, które stany zostały odwiedzone.

Tak jak UCS, A* można wykorzystać do całkowicie innych problemów. „Wystarczy” opracować dobrą heurystykę.

Następnym razem zmusimy komputer do rozwiązania sudoku.

Wilk, koza i kapusta – o UCS raz jeszcze

Ostatnio udało nam się znaleźć optymalne wyjście z labiryntu, uwzględniając różne czasy przejścia przez różne miejsca. Skomplikujmy zadanie. Niech robot przed wyjściem zbierze wszystkie monety. Jako bazę wykorzystajmy kod z poprzedniej części (przykładowe wyniki będą prezentowane dla wersji z czterema kierunkami ruchu zamiast ośmioma ale to nic nie zmienia. Rodzaje pól również zostają bez zmian (i również nie jest to istotne).

Reprezentacja stanu

Do tej pory stan świata był przedstawiony za pomocą pary współrzędnych (x,y). Teraz musimy dodać informację o monetach (sposób w jaki ich pozycje są ustalane nie będzie omawiany – nie jest to istotne z punktu widzenia rozważanego problemu). Najprostszym pomysłem wydaje się tablica wartości logicznych z informacją czy na danym polu znajduje się moneta. Zobaczmy ile będzie takich stanów. Labirynt ma rozmiar NM . Dla każdego pola są możliwe dwie wartości. Ilość możliwych kombinacji ustawień monet to 2^{MN} . W efekcie dostajemy coś koło 2^{MN} \ cdot MN . Dużo, nawet bardzo dużo. Spróbujmy inaczej – do informacji o pozycji robota dodajmy listę pól z monetami. Monet jest K. Możliwych list jest 2^K . Dostajemy 2^K \cdot MN. Jeśli monet jest znacząco mniej niż kafli labiryntu (K \ll MN) będzie to bardzo duża oszczędność. Osiągnęliśmy nie najgorszy efekt – 2^n rośnie wykładniczo (każde zmniejszenie n jest wartościowe) więc taka zmiana może być warta zachodu (innym miłym dodatkiem jest zmniejszenie użycia pamięci na reprezentację stanu wraz z rozwojem sytuacji kiedy monet jest coraz mniej).

Podsumowując, stan jest reprezentowany przez krotkę (x, y, coins_left) gdzie coins_left jest krotką z monetami pozostałymi na planszy.

Stan końcowy musi spełnić teraz dwa warunki. Stoimy na kaflu wyjściowym oraz nie pozostała żadna moneta do zebrania. (nic nie stoi ma przeszkodzie aby wystarczało tylko zebrać monety).

# sprawdź czy podany stan (x,y, c) jest stanem końcowym
def is_target_state(self, state):
	x, y, c = state
        return self.get_tile(state) == Tile.Target and len(c) == 0

Pozostało jeszcze zaimplementować get_after_decision oraz get_children. Niewielkie modyfikacje służą jedynie sprawdzeniu czy wskutek działania została zebrana moneta. Jedyna magia w kodzie polega na zamianie krotki na listę i odwrotnie.

# pobierz stan po podjęciu decyzji d w stanie state
def get_after_decision(self, state, d):
	x, y, c = state
	cc = list(c)
	dx, dy = Directions[d]
	if (x, y) in cc:
		cc.remove((x, y))
	return x + dx, y + dy, tuple(cc)
	
def get_children(self, state):
	x, y, c = state
	children = []
	# dla wszystkich potencjalnie możliwych ruchów
	for d in Directions:
		cc = list(c)
		dx, dy = Directions[d]
		# w zależności od typu kafla ustaw koszt przejścia
		t = self.get_tile((x + dx, y + dy, c))
		cost = 2 if t == Tile.Swamp else 1
		# jeżeli kafel zawierał monetę ustaw ją jako zebraną
		if (x + dx, y + dy) in cc:
			cc.remove((x + dx, y + dy))

		# jeżeli kafel nie jest ścianą dodaj do następników
		if t != Tile.Wall:
			children.append(((x + dx, y + dy, tuple(cc)), d, cost))
	# zapamiętaj fakt rozwinięcia (póki co nieistotne)
	self.expanded[y][x] = True
	if (x, y) in self.coins:
		self.coins.remove((x, y))

	return children

Zobaczmy jak to działa w praktyce.

dzxMWf

Zebranie monet

CRiKgn

Zebranie monet oraz wyjście

Voilà!

Wilk, koza i kapusta

Ariadna straciła pracę na rzecz komputera – Tezeuszowi do znalezienia Minotaura wystarczy UCS, ale we wstępie było o decyzjach w ogólności. Okazuje się, że jest to metoda na tyle uniwersalna aby poradzić sobie także z innymi problemami. Weźmy wszystkim znaną zagadkę z przewiezieniem przez rzekę wilka, kozy i kapusty. Brzmi ona następująco:

Pewien kupiec ma wilka, kozę i kapustę. Musi się on przeprawić przez rzekę i ma do dyspozycji jedną łódkę, do której może wejść w jednym momencie tylko on sam i jeden z przewożonych towarów. Wiadomo, że koza zje kapustę, a wilk kozę, gdy tylko któraś z tych par zostanie bez opieki. W jaki sposób kupiec ma pokonać rzekę nie tracąc przy tym żadnego z przewożonych towarów?

Kiedy człowiek zacznie się nad nią zastanawiać może pomyśleć następująco: „Ok, za pierwszym razem nie mogę zabrać ani wilka ani kapusty – musi to być koza. Kiedy wrócę mogę zabrać dowolnego z nich ale nie mogę zostawić ich razem…”
W poprzednim problemie nie było stanów które powodowały by porażkę które teraz się pojawiły. Jak możemy je obsłużyć? Na przynajmniej 3 sposoby.

1. Wykrywać stan porażki i informować o tym algorytm żeby go unikał.
2. Nie dopuszczać do pojawienia się takiego stanu w liście następników po rozwinięciu wierzchołka.
3. Przy próbie rozwinięcia stanu porażki zwrócić pustą listę następników – algorytm będzie musiał poszukać innej drogi gdyż ta jest ślepą uliczką.

Zdecydowałem się na 3 sposób gdyż jest on wg mnie najłatwiejszy w implementacji.

W jaki sposób należy pamiętać stan? Korzystam z z krotki 4 wartości logicznych mówiących czy dany osobnik jest po właściwej stronie rzeki. Dla wygody opakowałem to w klasę State.

class State:
    def __init__(self):
        self.human = False
        self.wolf = False
        self.goat = False
        self.cabbage = False

    def to_tuple(self):
        return self.human, self.wolf, self.goat, self.cabbage

    def from_tuple(self, t):
        self.human = t[0]
        self.wolf = t[1]
        self.goat = t[2]
        self.cabbage = t[3]

    def __str__(self):
        return 'Human: ' + str(self.human) + ' Wolf: ' + str(self.wolf) + ' Goat: ' + str(
            self.goat) + ' Cabbage: ' + str(self.cabbage)

Stan startowy oznacza że wszyscy są po złej stronie (4xFalse).

self.startState = State()

Stan końcowy wskazuje, że wszyscy są po właściwej stronie.

# sprawdź czy podany stan jest stanem końcowym
def is_target_state(self, state):
	s = State()
	s.from_tuple(state)
	return s.human and s.wolf and s.goat and s.cabbage

Sprawdzenie czy stan oznacza porażkę opera się prostym warunku.

def is_lose_state(self, s):
	return s.wolf == s.goat and s.wolf != s.human or s.goat == s.cabbage and s.goat != s.human

Rozwijając wierzchołek sprawdzamy czy stan nie oznacza porażki, jeżeli nie generujemy listę dozwolonych przejść.

def get_children(self, state):
	s = State()
	s.from_tuple(state)
	children = []

	if self.is_lose_state(s):
		return []

	#zawsze można przepłynąć pustą łódką
	no_transport = State()
	no_transport.human = not s.human
	no_transport.wolf = s.wolf
	no_transport.goat = s.goat
	no_transport.cabbage = s.cabbage

	children.append((no_transport.to_tuple(), "MOVE HUMAN", 1))

	if s.human == s.wolf:
		xd = State()
		xd.human = not s.human
		xd.wolf = not s.wolf
		xd.goat = s.goat
		xd.cabbage = s.cabbage
		children.append((xd.to_tuple(), "MOVE WOLF", 1))
	if s.human == s.goat:
		xd = State()
		xd.human = not s.human
		xd.goat = not s.goat
		xd.wolf = s.wolf
		xd.cabbage = s.cabbage
		children.append((xd.to_tuple(), "MOVE GOAT", 1))
	if s.human == s.cabbage:
		xd = State()
		xd.human = not s.human
		xd.cabbage = not s.cabbage
		xd.goat = s.goat
		xd.wolf = s.wolf
		children.append((xd.to_tuple(), "MOVE CABBAGE", 1))

	return children
	

I to wszystko. UCS nawet nie dotykaliśmy, wystarczyło zapewnić implementację kilku metod dla obiektu instance w sposób charakterystyczny dla problemu.

Panaceum?

Jeśli metoda jednolitego kosztu gwarantuje nam optymalne rozwiązanie to czy potrzebujemy czegoś więcej? Skoro to pytanie zostało postawione można się domyślić, że tak. O tym co możemy zrobić lepiej już niedługo.

Wiódł ślepy kulawego – Uniform-Cost Search [Część 2] – AI Tutorial 6

Część pierwsza

Rozwiązanie zaproponowane przez BFS nie jest optymalne – nie uwzględnia wyższego kosztu przechodzenia przez bagna. Potrzebujemy czegoś lepszego.

Uniform-Cost Search (metoda jednolitego kosztu)

Bazą dla dalszej pracy będzie algorytm BFS. Jak się okazuje jest to jego przypadek szczególny. Będziemy potrzebowali kolejki priorytetowej jako zbioru otwartego. W skrócie jest to struktura do której wrzucamy elementy skojarzone z pewnym priorytetem. Wyciągamy zawsze element o najniższym priorytecie.
Ustalmy czym jest priorytet w naszym problemie. Szukamy ścieżki najkrótszej, to jest o najniższym koszcie przejścia. Wynika z tego że najpierw chcemy zbadać ścieżki tańsze. Koszt przejścia między stanami nada się doskonale jako wartość priorytetu. Wróćmy do najprostszego przykładu:

4x4

Oryginalny graf wyglądał następująco:

4x4

Zmodyfikujmy sytuację tak aby wejście na pole (2,1) było bardziej kosztowne. Uaktualniony graf:

4x4(1)

Liczby w nawiasach oznaczają koszt przejścia.

Teraz drzewo. Zwróćmy uwagę, że każdy węzeł drzewa musi pamiętać całkowity koszt dojścia do niego – chcemy przecież zminimalizować koszt całej ścieżki. Gdybyśmy brali pod uwagę tylko koszt podjęcia danej decyzji dostalibyśmy algorytm który wybierał by zawsze najpierw tańszą decyzję – co globalnie wcale nie musi przekładać się na jakość rozwiązania.

4x4 ucs(1)Przeanalizujmy sytuację. W pierwszym kroku do kolejki trafiają dwa węzły na drugim poziomie drzewa (fioletowe). Zgodnie z regułą priorytetów do rozwinięcia wybierany jest lewy (ma niższy priorytet). Prawy pozostaje w kolejce (zbiorze otwartym). Po rozwinięciu dodawane są zielony i niebieski.
Teraz w kolejce są trzy węzły – prawy fioletowy, niebieski i zielony. Wszystkie mają priorytet równy dwa. Warto jednak zauważyć że fioletowy ma 0+2 (koszt rodzica albo koszt wsteczny + koszt własny) kiedy zielony i niebieski mają po 1+1.
Załóżmy że wybrany do rozwinięcia został fioletowy. Dodajemy czerwony i szary do kolejki. Stan kolejki:
zielony – 1+1=2
niebieski – 1+1=2
szary – 2+1=3
czerwony – 2+1=3
Do wyboru jest zielony i niebieski. Złożyło się pechowo i trafiło na niebieski. Po jego rozwinięciu powstaną węzły o koszcie co najmniej 3. Stan kolejki:
zielony – 1+1=2
szary – 2+1=3
czerwony – 2+1=3
kilka innych co najmniej 3 (koszt wsteczny = 2 + koszt własny co najmniej 1)
Wyciągamy z kolejki zielony (ma najniższy priorytet), jest to rozwiązanie. Znaleźliśmy drogę SE omijającą drogie pole (2,1).

Wcześniej było powiedziane, że BFS jest szczególnym przypadkiem UCS. Dlatego implementacja będzie bardzo podobna.

# coding=utf-8
from Queue import PriorityQueue


def ucs(instance):
    closed = set()  # zbiór zamknięty
    # inicjujemy zbiór otwarty stanem początkowym. Decyzję ustawiamy na null
    #koszt na 0 - nie ma to znaczenia. Rodzicem jest również null (jest to
    #korzeń drzewa
    #fringe dla UCS jest kolejką priorytetową (niższy priorytet powoduje szybsze zdjęcie
    #elementu
    #enqueue - put
    #dequeue - get
    fringe = PriorityQueue()
    #format wierzchołka jest nieco inny teraz to:
    #(priorytet,[((x,y), decyzja, koszt), rodzic])
    #jeżeli przez 'node' oznaczyć strukturę wierzchołka używaną poprzednio (w bfs i dfs)
    #to teraz jest to: (priorytet, node)
    #jest to wymagane przez kolejkę.
    fringe.put((0, [(instance.get_start_state(), None, 0), None]))
    #znaleziony cel
    target = None

    while True:
        #jeśli zbiór otwarty jest pusty to nie istnieje droga do celu
        if fringe.empty():
            return []
        #pobierz kolejny węzeł z kolejki - ten o najniższym priorytecie
        #zmiena node jest takie samego typu jak poprzednio (pobrany wierzchołek
        #rozpakowaliśmy)
        node_cost, node = fringe.get()

        #jeśli jesteśmy w stanie docelowym, ustaw cel i zakończ
        if instance.is_target_state(node[0][0]):
            target = node
            break

        #jeśli węzeł nie był rozwijany
        if node[0][0] not in closed:
            #dodaj go do zbioru zamkniętego (innymi słowy oznacz jako rozwinięty)
            closed.add(node[0][0])
            #rozwiń go
            children = instance.get_children(node[0][0])
            #i dla każdego następnika
            for child in children:
                #sprawdź koszt następnika
                child_cost = child[2]
                #dodaj informację o poprzedniku (node jest rodzicem child)
                #jako koszt ustaw sumę następnika i koszt dojścia do rodzica -
                #został on odczytany przy rozpakowywaniu krotki zwróconej przez
                #fringe.get()
                vertex = (node_cost + child_cost, [child, node])
                #i wrzuć do kolejki (zbioru otwartego)
                fringe.put(vertex)

    #lista decyzji prowadzących do rozwiązania
    solution = []
    #zaczynamy od węzła z wynikiem
    i = target
    #dopóki ma rodzica (nie jesteśmy w korzeniu)
    while i[1] is not None:
        #dodaj decyzję która nas tutaj doprowadziła
        solution.append(i[0][1])
        #przejdź do rodzica
        i = i[1]
    #podążaliśmy od wyjścia do startu przez co trzeba odwrócić kolejność
    solution.reverse()

    return solution

Spójrzmy na poprzednio używaną instancję i sprawdźmy jak poradzi sobie UCS.

3MdZLMUdało się! Zobaczmy jeszcze jeden przykład gdzie będzie lepiej widać, że komputer omija kosztowne pola jeśli jest inna opcja. Jeśli jednak prowadzi przez nie krótsza droga, wykorzystuje je.

06ORB7

Obiekt instance

Przeanalizujmy omijany do tej pory obiekt instance. Z poprzedniej części wiemy, że dostarcza metod:

get_start_state() – zwraca stan początkowy (x,y)
is_target_state((x,y)) – sprawdza czy podany stan jest stanem docelowym
get_children((x,y)) – pobiera następniki podanego stanu w formacie ((x,y), Decyzja(kierunek), Koszt).

Zobaczmy implementację na potrzeby labiryntu. (nieistotne fragmenty zostały wycięte – cały plik będzie wkrótce dostępny w repozytorium).

# Możliwe decyzje. Wartość to ruch w osi X i Y.
Directions = {
    'East': (1, 0),
    'North': (0, -1),
    'South': (0, 1),
    'West': (-1, 0)}


# Typy pól (kafli). Kolejno: puste, ściana, bagno, pole startowe, wyjście
class Tile:
    def __init__(self):
        pass

    Blank, Wall, Swamp, Start, Target = range(5)


# Klasa reprezentująca problem znajdowania drogi w labiryncie.
class Maze:
    def __init__(self):
        # układ pól na planszy
        self.layout = []
        # miejsce startowe (ustawiany przy wczytywaniu mapy)
        self.startState = None

    # pobierz kafel pod daną pozycją. state jest krotką (x,y)
    def get_tile(self, state):
        x, y = state
        return self.layout[y][x]

    # sprawdź czy podany stan (x,y) jest stanem końcowym
    def is_target_state(self, state):
        return self.get_tile(state) == Tile.Target

    # pobierz stan początkowy
    def get_start_state(self):
        return self.startState

    # pobierz stan po podjęciu decyzji d w stanie state
    def get_after_decision(self, state, d):
        x, y = state
        dx, dy = Directions[d]
        return x + dx, y + dy

    # rozwiń wierzchołek
    def get_children(self, state):
        x, y = state
        children = []
        # dla wszystkich potencjalnie możliwych ruchów
        for d in Directions:
            dx, dy = Directions[d]
            # w zależności od typu kafla ustaw koszt przejścia
            t = self.get_tile((x + dx, y + dy))
            cost = 2 if t == Tile.Swamp else 1
            # jeżeli kafel nie jest ścianą dodaj do następników
            if t != Tile.Wall:
                children.append(((x + dx, y + dy), d, cost))
        # zapamiętaj fakt rozwinięcia (póki co nieistotne)
        self.expanded[y][x] = True

        return children

Ruch po skosie

Wprowadźmy jeszcze jedną drobną zmianę – umożliwienie ruchu po skosie.

Dodajemy 4 nowe kierunki:

Directions = {
    'East': (1, 0),
    'North': (0, -1),
    'South': (0, 1),
    'West': (-1, 0),
    'NorthEast': (1, -1),
    'NorthWest': (-1, -1),
    'SouthEast': (1, 1),
    'SouthWest': (-1, 1)
}

Jeśli poruszamy się po skosie mnożymy koszt razy \sqrt{2} \approx 1,41 – tyle razy dłuższe jest przejście po skosie względem przejścia po prostej pionowej/poziomej.

# [...]
cost = 2 if t == Tile.Swamp else 1
if dx != 0 and dy != 0:
	cost *= 1.41
# jeżeli kafel nie jest ścianą dodaj do następników
# [...]

Glnt5w

Podsumowanie

Potrafimy już znaleźć drogę do wyjścia. Następnym razem poszerzymy zagadnienie o zebranie wszystkich monet przed wyjściem, a także zastanowimy się czy pokazaną metodę można wykorzystać w innych sytuacjach, niezwiązanych z szukaniem drogi.

qweqwe

Wiódł ślepy kulawego czyli o wyszukiwaniu drogi [część 1] – AI TUTORIAL 5

Do tej pory skupiliśmy się na zbieraniu wiedzy i poszukiwaniu pewnych wzorców, które pozwolą nam zamodelować rzeczywistość. Dzisiaj zastanówmy się nad podejmowaniem decyzji. Od naszej AI oczekujemy, że będzie podejmować racjonalne decyzje. Np. robot stojący przed jezdnią ma za zadanie przejść bezpiecznie na drugą stronę. W tym celu raczej nie wejdzie prosto pod koła samochodu, ale poczeka na wolną drogę. AI musi zdecydować czy już iść czy jeszcze nie. Podany przykład jest oczywiście bardzo prosty, ale idąc tym tropem rozszerzmy problem.

Weźmy klasyczne zagadnienie znajdowania najkrótszej drogi. Jak się potem okaże, nasze rozumowanie będzie bardzo przydatne w innych problemach. Rozważmy postać znajdującą się w labiryncie. Wiemy gdzie jesteśmy, dokąd zmierzamy oraz znamy układ ścian. Chcemy się stąd wydostać, najlepiej jak najszybciej. Aby nieco urozmaicić problem załóżmy, że koszt poruszania się zależy od typu podłoża.

0

Przykładowa instancja problemu. Niebieski – ściany, żółty – miejsce startowe, zielony – wyjście, brązowy – bagno w którym poruszamy się wolniej.

 

Za każdym razem możemy podjąć jedną z kilku decyzji – ruch w kierunku N/S/E/W. (Płn, Płd, Wsch, Zach. – dalej będę używał tych oznaczeń).
Oczekujemy znalezienia ciągu decyzji które zaprowadzą nas do wyjścia (np. NEESSENNNNNWWWNNE).

Potrzebujemy pewnego modelu świata. Zauważmy, że rozważane uniwersum może być w jednym z wielu stanów. Taki stan jest opisany jednoznacznie przez pozycję bohatera. Ile jest tych stanów? Dla labiryntu MxN jest ich również MN. Oczywiście pewne z nich są niepoprawne – robot nie może być wewnątrz ściany ale interesuje nas rząd wielkości a nie konkretna liczba. Gdybyśmy mieli dwójkę bohaterów było by to (MN)^2.

Pomiędzy tymi stanami możemy poruszać się w określony sposób. Można to doskonale przedstawić w postaci grafu. Wierzchołki reprezentują poszczególne stany (dozwolone/poprawne). A krawędzie/łuki między nimi określone decyzje. W takim grafie jak najbardziej mogą pojawić się cykle, choć intuicyjnie czujemy że błądzenie NSNSNSNSNS nie jest racjonalne. Spójrzmy na prosty labirynt i graf dla niego.

4x4 4x4

W tym momencie problem sprowadził się do znalezienia pewnej ścieżki, od stanu (wierzchołka) początkowego (żółtego) do stanu końcowego (zielonego). Jeżeli uwzględnimy różny koszt przejścia pomiędzy stanami, możemy opisać nim poszczególne krawędzie. Dodatkowo jeśli można przejść od stanu A do B ale nie od B do A to zamiast krawędzi umieszczamy łuk. (Inny przykład: podejmując decyzję „wciśnij czerwony guzik” przechodzimy ze stanu „bomba uzbrojona” do stanu „bomba zdetonowana” bez możliwości powrotu.)

bomb(1)

Liczby w nawiasach oznaczają koszt przejścia między stanami.

 

Drzewo decyzji.

Kolejno podejmowane decyzje możemy przedstawić w postaci drzewa. Korzeniem jest sytuacja początkowa kiedy jeszcze nic nie zdecydowaliśmy. Z każdego wierzchołka możemy przejść do tylu innych ile jest dozwolonych decyzji. Liśćmi tego drzewa są miejsca gdzie nie ma żadnej dozwolonej decyzji, lub miejsca gdzie cel jest spełniony. Jak możemy przypuszczać, omawiane grafy i drzewa mogą być BARDZO duże – z tego powodu będziemy generować je leniwie tj. kolejne części pojawią się dopiero gdy są niezbędne. Jak się za chwile okaże bardzo często będą one nieskończenie duże. Łuki między wierzchołkami oczywiście mogą mieć wagi podobnie jak w grafie. Zobaczmy jak będą wyglądało drzewo dla przedstawionego wcześniej labiryntu.

4x4tree

Nawet dla tak prostego problemu jest ono nieskończenie duże! Wynika to bezpośrednio z faktu istnienia cyklu w grafie. Nie będzie to jednak trudne do ominięcia.

Przeszukiwanie grafu

Pod pojęciem rozwinięcia wierzchołka rozumiem pobranie jego następników.

Wciąż szukamy wyjścia (najkrótszego) z labiryntu. Dla podanego wyżej przykładu oczywistym jest SE lub ES. Teraz pora na komputer.
Przypomnijmy, że w przedstawionym modelu świata szukamy ścieżki która pozwoli przejść od stanu startowego do stanu końcowego. Po drodze mamy cykle których chcielibyśmy uniknąć. Jak to zrobić? Wystarczy, że będziemy zapamiętywać już odwiedzone wierzchołki w tzw. zbiorze zamkniętym.
Idea jest następująca. Przechowujemy zbiór kandydatów do rozwinięcia(zbiór otwarty) oraz listę już rozwiniętych wierzchołków (zbiór zamknięty). Zgodnie z pewną strategią decydujemy który element zbioru otwartego wybrać do badania. Sprawdzamy czy jest on stanem końcowym, jeżeli nie – rozwijamy. Sprawdzamy czy był już rozwinięty. Jeżeli nie – rozwijamy a wszystkie następniki dodajemy zbioru otwartego zgodnie ze strategią. Procedurę kontynuujemy tak długo aż znajdziemy rozwiązanie lub skończą się elementy w zbiorze otwartym. Bardziej systematycznie:

Graph Search
S – stan/wierzchołek początkowy
T – stan/wierzchołek docelowy
1. zbiór otwarty = {S}
2. zbiór zamknięty = {}
3. Jeżeli zbiór otwarty jest pusty zwróć błąd.
4. akt = wybierz element ze zbioru otwartego zgodnie ze strategią.
5. jeśli akt == S zwróć sukces
6. jeśli akt nie jest w zbiorze zamkniętym
6.1. dodaj v do zbioru zamkniętego
6.2. dla każdego v w rozwiń(akt):
6.1.1 dodaj v do zbioru otwartego
7. Wróć do 3

Gdzie zbiór zamknięty jest zbiorem w sensie zbioru, a zbiór otwarty może być listą, stosem, etc lub w ogólności kolejką priorytetową.

No dobrze, w ten sposób dojdziemy do wierzchołka z rozwiązaniem – dowiedzieliśmy się że taka droga istnieje. Nie jest to zbyt satysfakcjonujące gdyż potrzebujemy konkretnych, kolejno podejmowanych decyzji. Możemy to zrobić na przynajmniej dwa sposoby.
1. Dla każdego węzła pamiętać decyzję która nas do niego doprowadziła oraz wykorzystanego poprzednika przechowywane w dedykowanym słowniku/mapie.
2. Zmienić nieco format wierzchołka tak aby pamiętał ścieżkę która do niego doprowadziła. Np:

4x4 full

W pokazanej implementacji używam drugiego sposobu.

Do pozyskiwania informacji o problemie (miejsce startu, cel, następniki, etc) używam obiektu instance który zapewnia metody:

get_start_state – zwraca stan początkowy (x,y)
is_target_state((x,y)) – sprawdza czy podany stan jest stanem docelowym
get_children((x,y)) – pobiera następniki podanego stanu w formacie ((x,y), Decyzja(kierunek), Koszt).

Dodatkowo fringe to zbiór otwarty. Wierzchołki są pamiętane w formacie:
[((x,y), Decyzja, Koszt), Poprzednik]

Zmieniając funkcję strategię modyfikujemy sposób przeszukiwania drzewa. Zastanówmy się teraz nad możliwościami.

DFS & BFS

Algorytmów Deep First Search i Breadth First Search nie będę szeroko omawiał – zostało to już zrobione milion razy w innych miejscach w sieci. Pokażę tylko jak wyglądają strategie dla nich.
DFS
W tym wypadku zbiorem otwartym jest stos, w punkcie 4.1 zdejmujemy element z wierzchu stosu, w 6.1.2 odkładamy na szczyt.
Implementacja:

# coding=utf-8
def dfs(instance):
    closed = set()  # zbiór zamknięty
    #inicjujemy zbiór otwarty stanem początkowym. Decyzję ustawiamy na null
    #koszt na 0 - nie ma to znaczenia. Rodzicem jest również null (jest to
    #korzeń drzewa
    #formalnie jest to lista, ale jest używana jako stos
    #pop - pop
    #push - append
    fringe = [[(instance.get_start_state(), None, 0), None]]
    #znaleziony cel
    target = None

    while True:
        #jeśli zbiór otwarty jest pusty to nie istnieje droga do celu
        if len(fringe) == 0:
            return []
        #pobierz kolejny węzeł z wierchu stosu
        node = fringe.pop()

        #jeśli jesteśmy w stanie docelowym, ustaw cel i zakończ
        if instance.is_target_state(node[0][0]):
            target = node
            break

        #jeśli węzeł nie był rozwijany
        if node[0][0] not in closed:
            #dodaj go do zbioru zamkniętego (innymi słowy oznacz jako rozwinięty)
            closed.add(node[0][0])
            #rozwiń go
            children = instance.get_children(node[0][0])
            #i dla każdego następnika
            for child in children:
                #dodaj informację o poprzedniku (node jest rodzicem child)
                vertex = [child, node]
                #i wrzuć na szczyt stosu (zbioru otwartego)
                fringe.append(vertex)

    #lista decyzji prowadzących do rozwiązania
    solution = []
    #zaczynamy od węzła z wynikiem
    i = target
    #dopóki ma rodzica (nie jesteśmy w korzeniu)
    while i[1] is not None:
        #dodaj decyzję która nas tutaj doprowadziła
        solution.append(i[0][1])
        #przejdź do rodzica
        i = i[1]
    #podążaliśmy od wyjścia do startu przez co trzeba odwrócić kolejność
    solution.reverse()

    return solution

Spójrzmy jak to działa. Mamy przykładowy labirynt:

0

Białe koło to nasz robot. Jak poradzi sobie DFS?

nf0VJw

Wyjście znalazł ale trudno nazwać je optymalnym. Nie ma nawet pozorów optymalności.

BFS

Tutaj używamy kolejki – dodanie elementu to zakolejkowanie na końcu, element jest pobierany z początku.
Kod jest bardzo podobny do poprzedniego – zmieniła się tylko strategia dodawania i wybierania ze zbioru otwartego.

# coding=utf-8
from Queue import Queue

def bfs(instance):
    closed = set()  # zbiór zamknięty
    #inicjujemy zbiór otwarty stanem początkowym. Decyzję ustawiamy na null
    #koszt na 0 - nie ma to znaczenia. Rodzicem jest również null (jest to
    #korzeń drzewa
    #fringe dla BFS jest kolejką
    #enqueue - put
    #dequeue - get
    fringe = Queue()
    fringe.put([(instance.get_start_state(), None, 0), None])
    #znaleziony cel
    target = None

    while True:
        #jeśli zbiór otwarty jest pusty to nie istnieje droga do celu
        if fringe.empty():
            return []
        #pobierz kolejny węzeł z wierzchu stosu
        node = fringe.get()

        #jeśli jesteśmy w stanie docelowym, ustaw cel i zakończ
        if instance.is_target_state(node[0][0]):
            target = node
            break

        #jeśli węzeł nie był rozwijany
        if node[0][0] not in closed:
            #dodaj go do zbioru zamkniętego (innymi słowy oznacz jako rozwinięty)
            closed.add(node[0][0])
            #rozwiń go
            children = instance.get_children(node[0][0])
            #i dla każdego następnika
            for child in children:
                #dodaj informację o poprzedniku (node jest rodzicem child)
                vertex = [child, node]
                #i wrzuć na szczyt stosu (zbioru otwartego)
                fringe.put(vertex)

    #lista decyzji prowadzących do rozwiązania
    solution = []
    #zaczynamy od węzła z wynikiem
    i = target
    #dopóki ma rodzica (nie jesteśmy w korzeniu)
    while i[1] is not None:
        #dodaj decyzję która nas tutaj doprowadziła
        solution.append(i[0][1])
        #przejdź do rodzica
        i = i[1]
    #podążaliśmy od wyjścia do startu przez co trzeba odwrócić kolejność
    solution.reverse()

    return solution

aUZOyU

Wynik jest znacznie lepszy niż poprzednio. Ściślej mówiąc: jeżeli zignorować dodatkowy koszt przechodzenia przez bagno jest optymalny.

Ciąg dalszy nastąpi

W następnym odcinku zajmiemy się uwzględnieniem różnego kosztu przechodzenia różnych pól.

Broyden–Fletcher–Goldfarb–Shanno – AI Tutorial 4

Poprzedni post z cyku.

Miało być o racjonalnych decyzjach ale podjąłem takową dzisiaj i będzie o czymś innym.

Wróćmy jeszcze do uczenia maszynowego i wprowadźmy kilka usprawnień.
Na początek zobaczmy dane na jakich będziemy pracować. (X – 0, gwiazdka – 1).

raw

Pierwszą zmianą będzie użycie innej metody minimalizowania. Gradient prosty jest prosty i tu w zasadzie kończą się jego zalety. Jednakże zamiast samemu implementować inny mechanizm skorzystajmy z tego co dostarcza biblioteka do obliczeń symbolicznych SymPy.
Wybór padł na algorytm Broyden–Fletcher–Goldfarb–Shanno. Nie będę go opisywać gdyż nie jest to blog o analizie numerycznej (jeszcze nie ( ͡° ͜ʖ ͡°) ).
Przegląd innych metod minimalizacji jest dostępny w dokumentacji.
Użycie jest bardzo proste, wystarczy zaimportować odpowiednią funkcję.

from scipy.optimize import fmin_bfgs

A następnie wywołać ją, przekazując jej funkcję do zoptymalizowania oraz miejsce od którego zacząć:

theta = fmin_bfgs(J, initialTheta)

Znajomość pochodnych nie jest konieczna, wystarczy nam sama funkcja J

#ostatecznie uznałem, że tak będzie czytelniej
def single_penalty(Theta, x, y):
    return -np.log(h(Theta, x)) if y == 1 else -np.log(1-h(Theta, x))                
        
def J(Theta):
    sum = 0.
    for i in range(m):
        sum += single_penalty(Theta, TrainingData[i,:], Y[i])
    return sum/m

Funkcja fmin_bfgs wypisze nam jak przebiegła optymalizacja oraz zwróci nam wynik (poszukiwane przez nas \theta )Np:

Optimization terminated successfully.
Current function value: 0.000002
Iterations: 30
Function evaluations: 487
Gradient evaluations: 37

Takie podejście do sprawy daje nam komfort swobodnych eksperymentów z funkcjami g_\theta(x) i J(\theta). Zacznijmy od najprostszej wersji hipotezy:
g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2

Możemy się spodziewać, że dla testowych danych trudno będzie poprowadzić sensownie prostą linię. I tak jest w istocie.

1

Wygląda źle, ale nie oczekiwaliśmy niczego innego. Dodajmy jeden wyraz do hipotezy.

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2

2

Eeee nie. Coś się zaczęło co prawda „wyginać” ale daleko nam do celu. Kolejne dwa wyrazy:

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2

3Znacznie lepiej. 3 kolejne okręgi (od środka) oznaczają kolejno P(y=1) = 1%, P(y=1) = 50%, P(y=1) = 99%. Im będą bliżej siebie, tym bardziej wyraźna granica decyzji i mniej obszarów niepewności. Czy dodanie kolejnego wyrazu może coś poprawić? Sprawdźmy

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2 + \theta_6 x_1^2x_2^2

4

Kształt zmienił się zdecydowanie, obszar niepewności również. Czy jest to lepsze dopasowanie niż poprzednio to już kwestia dyskusyjna. Faktem jest jednak że algorytm wytrenował się poprawnie do podanych danych. Dodajmy jeszcze coś.

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2 + \theta_6 x_1^2x_2^2 + \theta_7 x_1^2 x_2 + \theta_8 x_1 x_2^2

5

Sytuacja jest podobna. Ostatni krzyk rozpaczy, dodajemy trzecie potęgi.

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2 + \theta_6 x_1^2x_2^2 + \theta_7 x_1^2 x_2 + \theta_8 x_1 x_2^2 + \theta_9 x_1^3 + \theta_{10} x_2^3

6
Wystarczy. Dodawanie kolejnych potęg nic nie wnosi, a wręcz może zaszkodzić. W przyszłości poznamy mechanizm zniechęcania do wykorzystywania zbyt skomplikowanych krzywych.

Podsumowanie

Teraz już na pewno odstawiamy na bok (tymczasowo) uczenie maszynowe i przechodzimy do czegoś innego. O tym jak znaleźć drogę w labiryncie już niedługo.

Całość kodu

Regresja logistyczna (logistic regression) – AI tutorial 3

Pierwszy post z cyklu
Drugi post z cyklu

Zostawmy na jakiś czas egzaminy, przejdźmy do innych problemów.
Idąc do lekarza mówimy co nas boli, robimy badania (jeśli NFZ przypadkiem jeszcze nie zbankrutował) a wyniku pracy inteligentnej osoby dostajemy odpowiedź czy jesteśmy chorzy. Załóżmy, że chcemy zastąpić lekarzy maszynami. Doktor na podstawie zdobytej wiedzy oraz doświadczenia (czyli de facto również wiedzy) podejmuje decyzję, czy jesteśmy zdrowi czy chorzy. Jak to się ma do poprzedniego problemu? Trochę podobnie, trochę inaczej. Na wejściu jest to samo – n zmiennych wejściowych (po ang. zwane features). Wyjście jest inne – zamiast dostać wartość z ciągłego zbioru (punkty) dostajemy binarną odpowiedź chory/zdrowy. Mamy tutaj do czynienia z problemem klasyfikacji – efekt pracy może należeć tylko do pewnego skończonego zbioru. Wymaga to nieco innego podejścia. Inne przykłady takich problemów:

X oznacza wektor zmiennych wejściowych

  • Email klasyfikowany jako spam X → {0, 1}
  • Email klasyfikowany do pewnej kategorii (praca, prywatne, newslettery) X → {0, 1, 2}
  • Kategoria aut (sedan, pickup etc) na podstawie zdjęcia X → {0, 1, 2, 3, 4…}
  • Klasyfikacja płci X → {0, 1}

Testowany obiekt musi być w dokładnie jednej kategorii.

Na razie ograniczmy się do odpowiedzi tak/nie: y \in \{0,1\}. Postawmy taki problem: mamy chorego na określoną chorobę i chcemy określić czy może zostać wyleczony czy też nie. Do dyspozycji mamy takie dane jak jego wiek, waga, wyniki badań, płeć etc. Na początku załóżmy że rozważamy tylko wiek. Przykładowe dane mogą wyglądać następująco:

pic1

Kółko oznacza 0, gwiazdka 1 – będę dalej używać tej notacji.

Widzimy, że tendencja jest bardzo prosta – im pacjent starszy tym gorsze rokowania. Gdybyśmy chcieli wyznaczyć przedziały dla y=0 (wyzdrowienie) oraz y=1 (śmierć) mogło by wyglądać to tak:

pic2

Widać tutaj pewien problem – w okolicach miejsca podziału pewnie obiekty zostaną źle sklasyfikowane, nie da się tego całkowicie uniknąć. Co możemy z tym zrobić? Zamiast mówić „Hej, umrzesz” lepiej będzie powiedzieć „Umrzesz z prawdopodobieństwem 87%”. Oczywistą zaletą jest, że w niepewnej okolicy będziemy dostawać wyniki oscylujące wokół 50% co wskaże fakt niepewności klasyfikacji. Rzućmy okiem na przykład z dwiema zmiennymi wejściowymi:

db1

Jedna zmienna jest na osi X, druga na Y. Zaznaczona linia wskazuje miejsce w którym tniemy płaszczyznę (czy to zawsze i wszędzie muszą być linie proste? Nie, pokażemy wkrótce inne warianty).

Hipoteza h(x)

Ostatnio szukaliśmy takiej funkcji h która dla pewnego zestawu danych zwróci nam przewidywany wynik egzaminu. Teraz szukamy takiej która zwróci nam prawdopodobieństwo tego, że y = 1. W związku z tym funkcja musi być ograniczona 0 \le h_\theta(x) \le 1 . Konstrukcja takiej funkcji może być bardzo problematyczna.
Prościej było by uznać że jeśli h(x) > 0 – wynik pozytywny (y = 1). W takim rozumieniu tracimy jednak informację o prawdopodobieństwie chociaż intuicyjnie czujemy że im większa wartość tym większa pewność wyniku pozytywnego, a im mniejsza negatywnego. Potrzebujemy więc pewnej funkcji (nazwijmy ją S) która zamieni wygodną dla nas reprezentację na tą bardziej użyteczną.

S: (-\infty; +\infty) \mapsto (0;1)

Taka funkcja istnieje, nazywa się ona sigmoid function (tłumaczenia się nie podejmę) i jest zdefiniowana następująco:

S(t) = \frac{1}{1+e^{-t}}

Spójrzmy na wykres

 

Nada się znakomicie – dla t = 0 zwraca dokładnie 0,5. im dalej w kierunku nieskończoności tym bliżej 0 lub 1.

Uporządkujmy nieco oznaczenia. Funkcja hipotezy h_\theta(x) będzie teraz korzystać z funkcji S aby „przetłumaczyć” wynik z formy wygodnej dla nas na tą użyteczną. Do kreślenia prostych używaliśmy ostatnio h_\theta(\boldsymbol{X}) = \boldsymbol{\theta}^T \boldsymbol{X} . Teraz też będzie odpowiednia. Oznaczmy ją tym razem jako g(X):
g(\boldsymbol {X}) = \boldsymbol{\theta}^T \boldsymbol{X}
pamiętając że mamy do dyspozycji funkcję S:
S(t) = \frac{1}{1+e^{-t}}

Otrzymujemy ostatecznie:

h_\theta(\boldsymbol{X}) = S(g(\boldsymbol{X}))
lub równoważnie:
h_\theta(\boldsymbol{X}) = \frac{1}{1+e^{-\boldsymbol{\theta}^T \boldsymbol{X}}}

Kiedy będziemy chcieli zmienić zachowanie funkcji zmodyfikujemy tylko g(X)

Właśnie, jak mamy interpretować wartość zwracane przez g? Przy regresji liniowej była to wartość oczekiwana dla danego wejścia (co zresztą było bardzo dobrze widać na wykresie). Do tej pory ustaliliśmy: funkcja h zwraca prawdopodobieństwo wyniku pozytywnego (y = 1). Jeżeli jest ono większe niż 50% klasyfikujemy wynik jako pozytywny, w przeciwnym razie jako negatywny. Progiem „przełączania” jest wartość 0.5 funkcji h czyli wartość 0 funkcji g. Wynika z tego, że równanie g(x) = 0 pokazuje linię (powierzchnię, hiperpłaszczyznę etc) po której jednej stronie wartości są klasyfikowane jako pozytywne, a po drugiej – negatywne. Spójrzmy na przykłady

db2

 

Dotychczas używana przez nas metoda (konkretnie funkcja g) pozwalała tylko na oddzielanie elementów linią prostą. Nic jednak nie stoi na przeszkodzie aby była nieco bardziej skomplikowana, np:

g(x) = \theta_0 x_0 +\theta_1 x_1+\theta_2 x_2+ \theta_3 x_1^2 +\theta_4 x_2^2 +\theta_5 x_1 x_2 + \dots

Spowoduje to dopasowanie pewną krzywą, im dłuższe wyrażenie tym potencjalnie bardziej skomplikowaną. Rozwiniemy temat przy okazji sieci neuronowych (przykład na coś innego niż prosta pojawi się jeszcze w tym poście).

Funkcja kosztu J

Funkcję hipotezy mamy za sobą, co z funkcją kosztu J? Poprzednio miała formę:

J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Wyizolujmy z niej część odpowiedzialną za obliczenie kosztu pojedynczego zestawu danych. Będzie to:
h_\theta(x^{(i)}) - y^{(i)}

Zastanówmy się jak powinna wyglądać obecnie. Z całą pewnością chcemy, aby poprawny wynik był obarczony zerowym kosztem. Z drugiej strony jeżeli nasza hipoteza stwierdzi z dużą pewnością nieprawidłowy wynik, penalizacja powinna być wysoka. Pamiętajmy że poruszamy się w przedziale (0;1).

Na początek załóżmy, że y = 0. Jak ma kształtować się kara dla różnych wyników h(x)? Na przykład tak:

j0

Dla prawidłowego rezultatu kara jest równa zero. Im bardziej algorytm był pewny złej odpowiedzi tym większą karę ponosi (y dąży do nieskończoności wraz z x dążącym do 1).

Dla y = 1 sytuacja jest odwrotna:

j1

Wyrażenie opisujące pojedynczy zestaw danych będzie więc wyglądać tak:

\text{single penalty}(\boldsymbol{X}) =  \begin{cases}  -\log(h_\theta(\boldsymbol{X})) & \text{if } y = 1 \\  -\log(1-h_\theta(\boldsymbol{X})) & \text{if } y = 0  \end{cases}

Możemy pozbyć się niepożądanego warunku zapisując całość jako: (ze względu na y \in \{0,1\})

\text{single penalty}(\boldsymbol{X}) =-y \log(h_\theta(\boldsymbol{X})) - (1-y) \log(1-h_\theta(\boldsymbol{X}))

Podstawiając to do poprzedniego wyrażenia na wartość funkcji J(\boldsymbol{\theta}) otrzymujemy ostatecznie:

J(\boldsymbol{\theta}) = -\frac{1}{m} \left( \sum\limits_{i=1}^m \left (  y \log(h_\theta(\boldsymbol{X})) - (1-y) \log(1-h_\theta(\boldsymbol{X})    \right ) \right )

Optymalizacja

Tradycyjnie, kiedy mamy już opisną funkcję J problem jest sprowadzony do jej minimalizacji. Wyznaczmy pochodne kierunkowe względem wektora \boldsymbol{\theta}.

\frac{\partial}{\partial \theta_j} J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_j^{(i)}

???

Nie, to nie jest błąd. Czy dostaliśmy to samo (z dokładnością do mnożenia przez 2)? Prawie. Pamiętajmy że obecne h_\theta(x) wygląda inaczej niż w regresji liniowej.

Całość algorytmu gradientu prostego jest identyczna, więc nie będę go opisywał ponownie.

Implementacja

Zaimplementujmy wersję w której funkcja g(x) będzie prezentować się następująco:

g(x) = \theta_0 x_0 + \theta_1 x_1^2 + \theta_2 x_2^2

Jest to nic innego jak elipsa o środku w O=(0,0).

Funkcje S oraz h w kodzie Pythona:

def sigmoid(x):
    return 1/(1 + np.e**(-x))

def h(Theta, X):
    z = Theta[0] * X[0] + Theta[1] * X[1] **2 + Theta[2] * X[2]**2
    return sigmoid(z)

Obliczenia pochodnych cząstkowych są identyczne jak poprzednio podobnie jak i sam gradient. Spójrzmy na przykładowy wynik działania. Mając dostępne następujące dane:

raw

Algorytm wytyczył krzywą dla której prawdopodobieństwo obu wyników jest takie samo (P(y=1)=0,5 w sposób następujący:

line

Dodatkowo krzywe dla 75% prawdopodobieństwa y=1 oraz y=0.

imgNastępnym razem zostawimy uczenie się na pewien czas i zastanowimy się nad podejmowaniem rozsądnych decyzji.

Całość kodu

Regresja liniowa wielu zmiennych (AI tutorial 2)

Poprzedni post z cyklu

Założenie, że wynik egzaminu zależy tylko od obecności na wykładach jest bardzo naiwne. Bez problemu możemy wskazać inne czynniki – ocena z ćwiczeń, IQ, nastrój prowadzącego (im lepszy tym łagodniej ocenia), ilość godzin spędzonych na nauce, etc. Uogólnijmy więc metodę regresji liniowej na wiele zmiennych. Przykłady będę prowadził na dwóch, ale podane algorytmy działają dla dowolnej liczby (w granicach mocy obliczeniowej).

Powróćmy do przykładu z egzaminem. Tym razem do dyspozycji mamy również ocenę z laboratoriów. Dostępne dane treningowe prezentują się następująco:

training

Hipoteza h(x)

Trend wydaje się następujący: im więcej nieobecności tym gorszy wynik, im lepsza ocena tym więcej punktów.
Metoda postępowania jest bardzo podobna do poprzedniej, właściwie jedyną rzeczą którą będziemy robić będzie uogólnianie. Na początek przyjrzyjmy się funkcji h_\theta(x) . Poprzednio miała postać:

h_\theta(x) = \theta_0 + \theta_1 x

Oznaczmy opuszczone wykłady przez x_1 a ocenę z laborek przez x_2. W ogólnym przypadku x_k oznacza k-tą zmienną wejściową. Funkcja będzie wyglądać następująco:

h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2

W przypadku ogólnym (dla n zmiennych) będzie to:

h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

Spróbójmy uprościć nieco zapis (nie mam tu na myśli wstawienia znaku \sum ). Zmiennych wejściowych jest wiele, zapiszmy je jako wektor.

\boldsymbol{X} = \left [ x_1, x_2 \dots x_n \right ] n-elementowy wektor wierszowy \boldsymbol{X} \in \mathbb{R}^n

Parametrów θ jest również wiele – zróbmy z nich wektor \boldsymbol{\theta} \in \mathbb{R}^{n+1}. n + 1 ponieważ mamy \theta_0 nie powiązane z żadnym x.

\boldsymbol{\theta} = \left [ \theta_0, \theta_1 \dots \theta_n \right ]

Na razie nie pomogło nam to zbyt mocno, ilosć elementów jest różna. Wystarczy jeden prosty trick aby ją unormować. Dodajmy „fałszywą” zmienną x_0 = 1 i przypiszmy ją do pary razem z \theta_0

h_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

Mnożenie przez jeden oczywiście nic nie zmieni, zostaje to co było ale w wygodniejszej formie. X wygląda teraz tak:

\boldsymbol{X} = \left [ x_0 ,x_1, x_2 \dots x_n \right ], \quad \boldsymbol{X} \in \mathbb{R}^{n+1}, x_0 = 1

Dochodzimy do konkluzji. Przypominając sobie coś algebry całą funkcje zapisujemy w postaci:

h_\theta(\boldsymbol{X}) = \boldsymbol{\theta}^T \boldsymbol{X}

Zapis matematyczny jest zwięzły, ale czy przekłada się to na kod? Tak (zakładając że korzystamy z sensownego środowiska). Używam Pythona w wersji 3.4 oraz NumPy, funkcja h prezentuje się następująco:

def h(Theta, X):
    return Theta.transpose().dot(X)

Gdzie Theta oraz X są wektorami wierszowymi. transpose() zwraca transpozycję macierzy (pamiętajmy, że wektor to szczególny przypadek macierzy gdzie jeden z wymiarów jest równy 1). dot() dokonuje mnożenia macierzowego.

Funkcja kosztu J

Nie wiem czy funkcja kosztu to poprawne tłumaczenie ale uważam je za całkiem zręczne – koszt podobnie jak funkcję J chcemy zminimalizować.

Co z funkcją J? Nic. Mamy jedynie kosmetyczną zmianę, zamiast dwóch parametrów \theta_0 oraz \theta_1 mamy teraz jeden – wektor \boldsymbol{\theta}.

J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Również pochodne cząstkowe są praktycznie identyczne. Kontynuując nasz przykład dla dwóch zmiennych mamy trzy pochodne kierunkowe (pamiętajmy że \boldsymbol{\theta} to wektor (n+1)-elementowy gdzie n to ilość zmiennych wejściowych).

\frac{\partial}{\partial \theta_0} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right)
\frac{\partial}{\partial \theta_1} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_1^{(i)}
\frac{\partial}{\partial \theta_2} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_2^{(i)}

Pojawiła się nowa notacja x_j^{(i)} Oznacza ona wartość zmiennej x_j w i-tym zestawie danych. Spójrzmy na przykład. Weźmy i = 8, j = 2. Na podstawie tabelki z początku postu odczytujemy że jest to 4,5 (wiersz 8, druga zmienna, dla przypomnienia wartości w kolumnie „punkty na egzaminie” to nasze y).

Pochodne po \theta_1 i \theta_2 możemy swobodnie uogólnić.
Czy to samo możemy zrobić z pochodną po \theta_0? Tak, x_0 jest przecież zawsze równe jeden. Ostatecznie wszystkie pochodne wyznaczamy ze wzoru:
\frac{\partial}{\partial \theta_j} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_j^{(i)}

Kod funkcji obliczajacej wartość pochodnej funkcji J po j-tej zmiennej dla danego \boldsymbol{\theta}

def Jpartial(j,Theta):
    sum = 0
    for i in range(m):
        sum += (h(Theta, TrainingData[i]) - Y[i]) * TrainingData[i][j]
    return 2 * sum / m

Pętlą for przechodzimy po wszystkich zestawach danych. TrainingData jest macierzą z danymi treningowymi. Dla rozważanych danych (tabela) wygląda ona tak:

TrainingData = np.array([     
	[1.,1.	,4.5],[1.,1.	,3. ],[1.,1.	,4. ],
	[1.,2.	,3.5],[1.,2.	,5. ],[1.,2.	,3. ],
	[1.,3.	,3. ],[1.,3.	,4.5],[1.,3.	,4. ],
	[1.,4.	,3. ],[1.,4.	,3.5],[1.,4.	,5. ],
	[1.,5.	,4. ],[1.,5.	,3. ],[1.,5.	,3.5],
	[1.,6.	,4.5],[1.,6.	,4.5],[1.,6.	,4.5],
	[1.,7.	,4. ],[1.,7.	,3. ],[1.,7.	,5. ],
	[1.,8.	,3. ],[1.,8.	,4. ],[1.,8.	,3.5],
	[1.,9.	,3. ],[1.,9.	,4. ],[1.,10.	,4.5],
	[1.,10.	,5. ],[1.,10.	,3. ]])
m,n = TrainingData.shape

Jest to macierz (w terminologii NumPy tablica) o m wierszach i n + 1 kolumnach. Kolumna 0 to oczywiście sztuczne x_0.
Zapis

TrainingData[i][j]

Oznacza wybranie wartości z i-tego wiersza i j-tej kolumny.

TrainingData[i]

to wybór całego i-tego wiersza.
Y to wektor z wynikami (trzecia kolumna w tabelce).

Metoda gradientu prostego dla wielu zmiennych

Już poprzednio szukanie minimum wymagało modyfikowania dwóch zmiennych. Dwie pętle for i jesteśmy w domu. Zanim to zrobimy zastanówmy się nad warunkiem stopu. Poprzednio podałem bardzo ogólne „skończ kiedy uznasz za stosowne”. Teraz zastanówmy się kiedy stosownie będzie zakończyć.

 

Analiza matematyczna mówi nam, że w optimum lokalnym wartości pochodnych są równe zero. Szansa na to, że dostaniemy *dokładnie* zero jest niewielka (a właściwie to zerowa – im bliżej jesteśmy tym wolniej idziemy, do zera zbliżamy się asymptotycznie). Pamiętajmy że obliczenia numeryczne są mniej lub bardziej przybliżone. Wybierzmy zatem pewną liczbę ε poniżej której uznamy, że funkcja się zbiegła. Jeżeli wszystkie pochodne są mniejsze od niej stwierdzamy osiągnięcie minimum i kończymy.

converge

Zbieżność do zera. X – iteracje, Y = J(θ)

notcoverage

Rozbieżność. Wartość α jest zbyt duża. X – iteracje, Y = J(θ)

 

def gradient_descent(startTheta):
    theta = startTheta
    converge = False
    i = 0
    while not converge:
        differentials = np.zeros(n)
        for j in range(n):
            differentials[j] = Jpartial(j, theta)
        for j in range(n):
            theta[j] -= alpha * differentials[j]
			
        if i % 500 == 0:
            print (np.amax(differentials))
        converge = np.amax(differentials) < epsylon
    return theta

startTheta to punkt początkowy. Dopóki nie stwierdzimy zbieżności obliczamy wartości wszystkich pochodnych. np.zeros(n) tworzy tablicę w której będziemy zapisywać ich wartości. Po ich obliczeniu przesuwamy wszystkie θ. Jeśli wszystkie są mniejsze od ε kończymy.
Dodatkowo co 500 iteracji wypisujemy informację która pozwoli nam ocenić jak długo musimy jeszcze czekać.

Wyniki

Zobaczmy jakich wartości θ wyuczył się algorytm.
\theta_0 = -0,173
\theta_1 = -0,396
\theta_2 = 1,935

Opuszczone wykłady są skorelowane z modyfikatorem -0,4. Im więcej opuszczonych tym mniej punktów. Ocena z laborek ma mnożnik ~2. Widać że komputer doszedł do „rozsądnych” wniosków (na podstawie zadanych wartości treningowych). Zauważyć można również fakt, że ocena jest ~4x bardziej istotna przy przewidywaniu wyniku. Zapytajmy teraz o kilku studentów i sprawdźmy co na temat ich egzaminu powie algorytm.

results

Pamiętajmy, że jakość otrzymywanych wyników zależy w kluczowym stopniu od podanego zestawu treningowego!

Całość kodu