ads2_2022/code/python/src/algorithms/random_walk/algorithms.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# IMPORTS
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

from src.thirdparty.maths import *;
from src.thirdparty.plots import *;
from src.thirdparty.types import *;

from models.generated.config import *;
from models.generated.commands import *;
from src.core.log import *;
from src.core.utils import *;
from src.models.random_walk import *;
from src.algorithms.random_walk.display import *;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# EXPORTS
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

__all__ = [
    'adaptive_walk_algorithm',
    'gradient_walk_algorithm',
    'metropolis_walk_algorithm',
];

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# CONSTANTS
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MAX_ITERATIONS = 1000; # um endlose Schleifen zu verhindern

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# METHOD adaptive walk
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

def adaptive_walk_algorithm(
    landscape: Landscape,
    r: float,
    coords_init: tuple,
    optimise: EnumOptimiseMode,
    verbose: bool,
):
    '''
    Führt den Adapative-Walk Algorithmus aus, um ein lokales Minimum zu bestimmen.
    '''

    # lege Fitness- und Umgebungsfunktionen fest:
    match optimise:
        case EnumOptimiseMode.max:
            f = lambda x: -landscape.fitness(*x);
        case _:
            f = lambda x: landscape.fitness(*x);
    nbhd = lambda x: landscape.neighbourhood(*x, r=r, strict=True);
    label = lambda x: landscape.label(*x);

    # initialisiere
    steps = [];
    x = coords_init;
    fx = f(x);
    fy = fx;
    N = nbhd(x);

    # führe walk aus:
    k = 0;
    while k < MAX_ITERATIONS:
        # Wähle zufälligen Punkt und berechne fitness-Wert:
        y = uniform_random_choice(N);
        fy = f(y);

        # Nur dann aktualisieren, wenn sich f-Wert verbessert:
        if fy < fx:
            # Punkt + Umgebung + f-Wert aktualisieren
            x = y;
            fx = fy;
            N = nbhd(x);
            step = Step(coords=x, label=label(x), improved=True, changed=True);
        else:
            # Nichts (außer logging) machen!
            step = Step(coords=x, label=label(x));

        # Nur dann (erfolgreich) abbrechen, wenn f-Wert lokal Min:
        if fx <= min([f(y) for y in N], default=fx):
            step.stopped = True;
            steps.append(step);
            break;

        steps.append(step);
        k += 1;

    return x;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# METHOD gradient walk
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

def gradient_walk_algorithm(
    landscape: Landscape,
    r: float,
    coords_init: tuple,
    optimise: EnumOptimiseMode,
    verbose: bool,
):
    '''
    Führt den Gradient-Descent (bzw. Ascent) Algorithmus aus, um ein lokales Minimum zu bestimmen.
    '''

    # lege Fitness- und Umgebungsfunktionen fest:
    match optimise:
        case EnumOptimiseMode.max:
            f = lambda x: -landscape.fitness(*x);
        case _:
            f = lambda x: landscape.fitness(*x);
    nbhd = lambda x: landscape.neighbourhood(*x, r=r, strict=True);
    label = lambda x: landscape.label(*x);

    # initialisiere
    steps = [];
    x = coords_init;
    fx = landscape.fitness(*x);
    fy = fx;
    N = nbhd(x);
    f_values = [f(y) for y in N];
    fmin = min(f_values);
    Z = [y for y, fy in zip(N, f_values) if fy == fmin];

    # führe walk aus:
    k = 0;
    while k < MAX_ITERATIONS:
        # Wähle zufälligen Punkt mit steilstem Abstieg und berechne fitness-Wert:
        y = uniform_random_choice(Z);
        fy = fmin;

        # Nur dann aktualisieren, wenn sich f-Wert verbessert:
        if fy < fx:
            # Punkt + Umgebung + f-Wert aktualisieren
            x = y;
            fx = fy;
            N = nbhd(y);
            f_values = [f(y) for y in N];
            fmin = min(f_values);
            Z = [y for y, fy in zip(N, f_values) if fy == fmin];
            step = Step(coords=x, label=label(x), improved=True, changed=True);
        else:
            # Nichts (außer logging) machen!
            step = Step(coords=x, label=label(x));

        # Nur dann (erfolgreich) abbrechen, wenn f-Wert lokal Min:
        if fx <= min([f(y) for y in N], default=fx):
            step.stopped = True;
            steps.append(step);
            break;

        steps.append(step);
        k += 1;

    return x;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# METHOD metropolis walk
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

def metropolis_walk_algorithm(
    landscape: Landscape,
    r: float,
    coords_init: tuple,
    T: float,
    annealing: bool,
    optimise: EnumOptimiseMode,
    verbose: bool,
):
    '''
    Führt den Metropolis-Walk Algorithmus aus, um ein lokales Minimum zu bestimmen.
    '''

    # lege Fitness- und Umgebungsfunktionen fest:
    match optimise:
        case EnumOptimiseMode.max:
            f = lambda x: -landscape.fitness(*x);
        case _:
            f = lambda x: landscape.fitness(*x);
    nbhd = lambda x: landscape.neighbourhood(*x, r=r, strict=True);
    label = lambda x: landscape.label(*x);

    # definiere anzahl der hinreichenden Schritt für Stabilität:
    n_stable = 2*(3**(landscape.dim) - 1);

    # initialisiere
    x = coords_init;
    fx = f(x);
    fy = fx;
    nbhd_x = nbhd(x);
    steps = [];
    step = Step(coords=x, label=label(x));

    # führe walk aus:
    k = 0;
    n_unchanged = 0;
    while k < MAX_ITERATIONS:
        # Wähle zufälligen Punkt und berechne fitness-Wert:
        y = uniform_random_choice(nbhd_x);
        fy = f(y);
        p = math.exp(-abs(fy-fx)/T);
        u = random_binary(p);

        # Aktualisieren, wenn sich f-Wert verbessert
        # oder mit einer Wahrscheinlichkeit von p:
        if fy < fx or u:
            # Punkt + Umgebung + f-Wert aktualisieren
            x = y;
            fx = fy;
            nbhd_x = nbhd(x);
            n_unchanged = 0;
            step = Step(coords=x, label=label(x), improved=(fy < fx), chance=u, probability=p, changed=True);
        else:
            # Nichts (außer logging) machen!
            n_unchanged += 1;
            step = Step(coords=x, label=label(x));

        # »Temperatur« ggf. abkühlen:
        if annealing:
            T = cool_temperature(T, k);

        # Nur dann (erfolgreich) abbrechen, wenn f-Wert lokal Min:
        if n_unchanged >= n_stable:
            step.stopped = True;
            steps.append(step);
            break;

        steps.append(step);
        k += 1;

    if verbose:
        for step in steps:
            print(step);

    return x;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# AUXILIARY METHODS
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

def cool_temperature(T: float, k: int, const: float = 2.) -> float:
    harm = const*(k + 1);
    return T/(1 + T/harm);