snp.py

# encoding: utf-8
import math
from typing import Any
import numpy as np
import matplotlib.pyplot as plt
import argparse

from sho import make, algo, iters, plot, num, bit, pb

# Dimension of the search space.
d = 2

########################################################################
# Interface
########################################################################


def get_args_parser():
    can = argparse.ArgumentParser()

    can.add_argument(
        "-n",
        "--nb-sensors",
        metavar="NB",
        default=3,
        type=int,
        help="Number of sensors",
    )

    can.add_argument(
        "-r",
        "--sensor-range",
        metavar="RATIO",
        default=0.3,
        type=float,
        help="Sensors' range (as a fraction of domain width, max is √2)",
    )

    can.add_argument(
        "-w",
        "--domain-width",
        metavar="NB",
        default=30,
        type=int,
        help="Domain width (a number of cells). If you change this you will probably need to update `--target` accordingly",
    )

    can.add_argument(
        "-i",
        "--iters",
        metavar="NB",
        default=100,
        type=int,
        help="Maximum number of iterations",
    )

    can.add_argument(
        "-s",
        "--seed",
        metavar="VAL",
        default=None,
        type=int,
        help="Random pseudo-generator seed (none for current epoch)",
    )

    solvers = [
        "num_greedy",
        "bit_greedy",
        "num_sim_anneal",
        "num_random",
        "num_evolutionary",
        "bit_sim_anneal"
    ]
    can.add_argument(
        "-m",
        "--solver",
        metavar="NAME",
        choices=solvers,
        default="num_evolutionary",
        help="Solver to use, among: " + ", ".join(solvers),
    )

    can.add_argument(
        "-t",
        "--target",
        metavar="VAL",
        default=30 * 30,
        type=float,
        help="Objective function value target",
    )

    can.add_argument(
        "-y",
        "--steady-delta",
        metavar="NB",
        default=50,
        type=float,
        help="Stop if no improvement after NB iterations",
    )

    can.add_argument(
        "-e",
        "--steady-epsilon",
        metavar="DVAL",
        default=0,
        type=float,
        help="Stop if the improvement of the objective function value is lesser than DVAL",
    )

    can.add_argument(
        "-a",
        "--variation-scale",
        metavar="RATIO",
        default=0.3,
        type=float,
        help="Scale of the variation operators (as a ration of the domain width)",
    )

    can.add_argument(
        "--nb-population",
        metavar="NB",
        default=100,
        type=int,
        help="Size of the initial population for evolutionary algorithm.",
    )

    can.add_argument(
        "--nb-selected",
        metavar="NB",
        type=int,
        default=10,
        help="Number of selected individuals before each generation.",
    )

    can.add_argument(
        "--nb-offspring",
        metavar="NB",
        default=100,
        type=int,
        help="Number of offspring for each generation.",
    )

    can.add_argument(
        "--quality-threshold",
        metavar="DVAL",
        default=660,
        type=float,
        help="Quality threshold. Used to plot the probability of a run being under the quality threshold.",
    )

    return can


def run_algorithm(the, iters_func):
    val, sol, sensors = None, None, None
    if the.solver == "num_greedy":
        val, sol = algo.greedy(
            make.func(
                num.cover_sum,
                domain_width=the.domain_width,
                sensor_range=the.sensor_range,
                dim=d * the.nb_sensors,
            ),
            make.init(
                num.rand,
                dim=d * the.nb_sensors,
                scale=the.domain_width,
            ),
            make.neig(
                num.neighb_square,
                scale=the.variation_scale,
                domain_width=the.domain_width,
            ),
            iters_func,
        )
        sensors = num.to_sensors(sol)

    elif the.solver == "bit_greedy":
        val, sol = algo.greedy(
            make.func(
                bit.cover_sum,
                domain_width=the.domain_width,
                sensor_range=the.sensor_range,
                dim=d * the.nb_sensors,
            ),
            make.init(
                bit.rand,
                domain_width=the.domain_width,
                nb_sensors=the.nb_sensors,
            ),
            make.neig(
                bit.neighb_square,
                scale=the.variation_scale,
                domain_width=the.domain_width,
            ),
            iters_func,
        )
        sensors = bit.to_sensors(sol)
    elif the.solver == "num_sim_anneal":
        val, sol = algo.simulated_annealing(
            make.func(
                num.cover_sum,
                domain_width=the.domain_width,
                sensor_range=the.sensor_range,
                dim=d * the.nb_sensors,
            ),
            make.init(
                num.init_with_t,
                t=2000,
                dim=d * the.nb_sensors,
                scale=the.domain_width,
            ),
            make.neig(
                num.neighb_square,
                scale=the.variation_scale,
                domain_width=the.domain_width,
            ),
            make.temp(0.99),
            make.proba(),
            make.rand(),
            iters_func,
        )
        sensors = num.to_sensors(sol)

    elif the.solver == "bit_sim_anneal":
        val, sol = algo.simulated_annealing(
            make.func(
                bit.cover_sum,
                domain_width=the.domain_width,
                sensor_range=the.sensor_range,
                dim=d * the.nb_sensors,
            ),
            make.init(
                bit.init_with_t,
                t=2000,
                domain_width=the.domain_width,
                nb_sensors=the.nb_sensors,
            ),
            make.neig(
                bit.neighb_square,
                scale=the.variation_scale,
                domain_width=the.domain_width,
            ),
            make.temp(0.99),
            make.proba(),
            make.rand(),
            iters_func,
        )
        sensors = bit.to_sensors(sol)
    elif the.solver == "num_random":
        val, sol = algo.random(
            make.func(
                num.cover_sum,
                domain_width=the.domain_width,
                sensor_range=the.sensor_range,
                dim=d * the.nb_sensors,
            ),
            make.init(
                num.rand,
                dim=d * the.nb_sensors,
                scale=the.domain_width,
            ),
            iters_func,
        )
        sensors = num.to_sensors(sol)
    elif the.solver == "num_evolutionary":
        val, sol = algo.evolutionary(
            make.func(
                num.cover_sum,
                domain_width=the.domain_width,
                sensor_range=the.sensor_range,
                dim=d * the.nb_sensors,
            ),
            make.init(
                num.init_evolutionary,
                dim=d * the.nb_sensors,
                scale=the.domain_width,
                nb_population=the.nb_population,
            ),
            num.best,
            make.selection(num.selection, nb_selected=the.nb_selected),
            num.evaluate,
            make.variation(
                num.variation,
                nb_offspring=the.nb_offspring,
                domain_width=the.domain_width,
                variation_scale=the.variation_scale,
                nb_sensors=the.nb_sensors,
            ),
            make.replacement(
                num.replacement,
                nb_next_generation=the.nb_selected,
                strat=None,
            ),
            iters_func,
        )
        sensors = num.to_sensors(sol)
    return val, sol, sensors


def plot_sol(the, history, sensors, shape):
    fig = plt.figure()

    if the.nb_sensors == 1 and the.domain_width <= 50:
        ax1 = fig.add_subplot(121, projection="3d")
        ax2 = fig.add_subplot(122)

        f = make.func(
            num.cover_sum,
            domain_width=the.domain_width,
            sensor_range=the.sensor_range,
            dim=d * the.nb_sensors,
        )
        plot.surface(ax1, shape, f)
        plot.path(ax1, shape, history)
    else:
        ax2 = fig.add_subplot(121)

    domain = np.zeros(shape)
    domain = pb.coverage(domain, sensors, the.sensor_range * the.domain_width)
    domain = plot.highlight_sensors(domain, sensors)
    ax2.imshow(domain)

    plt.show()


def main(args=None):
    can = get_args_parser()

    the = can.parse_args(args)

    # Minimum checks.
    assert 0 < the.nb_sensors
    assert 0 < the.sensor_range <= math.sqrt(2)
    assert 0 < the.domain_width
    assert 0 < the.iters

    # Do not forget the seed option,
    # in case you would start "runs" in parallel.
    np.random.seed(the.seed)

    # Weird numpy way to ensure single line print of array.
    np.set_printoptions(linewidth=np.inf)  # type: ignore

    # Common termination and checkpointing.
    history: list[Any] = []
    iters_func = make.iter(
        iters.several,
        agains=[
            make.iter(iters.max, nb_it=the.iters),
            make.iter(
                iters.save, filename=the.solver + ".csv", fmt="{it} ; {val} ; {sol}\n"
            ),
            make.iter(iters.log, fmt="\r{it} {val}"),
            make.iter(iters.history, history=history),
            make.iter(iters.target, target=the.target),
            iters.steady(the.steady_delta, the.steady_epsilon),
        ],
    )

    # Erase the previous file.
    with open(the.solver + ".csv", "w") as fd:
        fd.write("# {} {}\n".format(the.solver, the.domain_width))

    val, sol, sensors = run_algorithm(the, iters_func)

    # Fancy output.
    print("\n{} : {}".format(val, sensors))

    shape = (the.domain_width, the.domain_width)

    should_plot = False
    if should_plot:
        plot_sol(the, history, sensors, shape)

    # coefficient when objective function is called multiples times in one epoch
    coef = get_coef_objective_func(the)

    # constant when objective function is called multiples times in one epoch
    initial = get_constant_objectiv_func(the)
    return val, len(history) * coef + initial, sol, sensors


def get_constant_objectiv_func(the):
    return the.nb_population if the.solver == "num_evolutionary" else 0


def get_coef_objective_func(the):
    return the.nb_offspring if the.solver == "num_evolutionary" else 1


if __name__ == "__main__":

    main()