"""
Random Lines Algorithm finds the optimal minimum of a function.

Sahin, I. (2013). Minimization over randomly selected lines.  An
International  Journal Of Optimization And Control: Theories &
Applications (IJOCTA), 3(2), 111-119.
http://dx.doi.org/10.11121/ijocta.01.2013.00167
"""

# Author : Ismet Sahin

__all__ = ["random_lines", "particle_swarm"]

from itertools import count

from numpy import zeros, ones, asarray, sqrt, arange, isfinite
from numpy.random import rand, random_integers


def print_every_five(step, x, fx, k):
    if step % 5 == 0:
        print(step, ":", fx[k], x[k])


def random_lines(cfo, NP, CR=0.9, epsilon=1e-10, abort_test=None, maxiter=1000):
    """
    Random lines is a population based optimizer which using quadratic
    fits along randomly oriented directions.

    *cfo* is the cost function object.  This is a dictionary which contains
    the following keys:

        *cost* is the function to be optimized.  If *parallel_cost* exists,
        it should accept a list of points, not just a single point on each
        evaluation.

        *n* is the problem dimension

        *x0* is the initial point

        *x1* and *x2* are lower and upper bounds for each parameter

        *monitor* is a callable which is called each iteration using
        *callback(step, x, fx, k)*, where *step* is the iteration number,
        *x* is the population, *fx* is value of the cost function for each
        member of the population and *k* is the index of the best point in
        the population.

        *f_opt* is the target value of the optimization

    *NP* is the number of fit parameters

    *CR* is the cross-over ratio, which is the proportion of dimensions
    that participate in any random orientation vector.

    *epsilon* is the convergence criterion.

    *abort_test* is a callable which indicates whether an external processes
    requests the fit to stop.

    *maxiter* is the maximum number of generations

    Returns success, num_evals, f(x_best), x_best.
    """
    if "parallel_cost" in cfo:
        mapper = lambda v: asarray(cfo["parallel_cost"](v.T), "d")
    else:
        mapper = lambda v: asarray(list(map(cfo["cost"], v.T)), "d")
    monitor = cfo.get("monitor", print_every_five)

    n = cfo["n"]

    X = rand(n, NP)  # will hold original vectors

    # CREATE FIRST GENERATION WITH LEGAL PARAMETER VALUES AND EVALUATE COSTS
    # m th member of the population
    for m in range(0, NP):
        X[:, m] = cfo["x1"] + (cfo["x2"] - cfo["x1"]) * X[:, m]
    if "x0" in cfo:
        X[:, 0] = cfo["x0"]
    f = mapper(X)

    n_feval = NP
    f_best, i_best = min(zip(f, count()))

    # CHECK INITIAL STOPPING CRITERIA
    if abs(cfo["f_opt"] - f_best) < epsilon:
        satisfied_sc = 1
        x_best = X[:, i_best]
        return satisfied_sc, n_feval, f_best, x_best

    for L in range(1, maxiter + 1):
        # finding destination vector
        i_Xj = random_integers(0, NP - 2, NP)
        i_ge = i_Xj >= arange(0, NP)
        i_Xj[i_ge] += 1

        # choosing muk
        muk = 0.01 + 0.49 * rand(NP)
        inx = rand(NP) < 0.5
        muk[inx] = -muk[inx]

        # find xk and fk s
        Xi = X
        Xj = X[:, i_Xj]
        P = Xj - Xi
        Xk = Xi + (ones((n, 1)) * muk) * P
        fk = mapper(Xk)
        n_feval = n_feval + NP

        # find quadratic models
        if any(muk == 0) or any(muk == 1):
            satisfied_sc = 0
            x_best = X[:, i_best]
            print("muk cannot be zero or one !!!")
            return satisfied_sc, n_feval, f_best, x_best

        fi = f
        fj = f[i_Xj]
        b = (muk / (muk - 1)) * fj - ((muk + 1) / muk) * fi - (1 / (muk * (muk - 1))) * fk
        a = fj - fi - b

        crossovers = []
        for k in range(0, NP):
            if abs(a[k]) < 1e-30 or (a[k] < 0 and fk[k] > fi[k] and fk[k] > fj[k]) or not isfinite(a[k]):
                # xi survives
                continue
            else:
                # xi may not survive
                mustar = -b[k] / (2 * a[k])
                xstar = Xi[:, k] + mustar * P[:, k]

                # choosing random numbers for crossover
                rn = rand(n)
                indi = rn < 0.5 * (1 - CR)
                indj = rn > 0.5 * (1 + CR)
                xstar[indi] = Xi[indi, k]
                xstar[indj] = Xj[indj, k]

                # map into feasible set
                inx = xstar < cfo["x1"]
                xstar[inx] = cfo["x1"][inx]
                inx = xstar > cfo["x2"]
                xstar[inx] = cfo["x2"][inx]

                crossovers.append((k, xstar))

        if len(crossovers) > 0:
            idx, xstar = [asarray(v) for v in zip(*crossovers)]
            fstar = mapper(xstar.T)
            n_feval += len(crossovers)

            # xi does not survive, xstar replaces it
            update = fstar < fi[idx]
            f[idx[update]] = fstar[update]
            X[:, idx[update]] = xstar[update, :].T

        # CHECKING STOPPING CRITERIA
        f_best, i_best = min(zip(f, count()))
        if abs(cfo["f_opt"] - f_best) < epsilon:
            satisfied_sc = 1
            x_best = X[:, i_best]
            return satisfied_sc, n_feval, f_best, x_best
        if abort_test():
            break

        monitor(L, X, f, i_best)

    return 1, n_feval, f_best, X[:, i_best]


def particle_swarm(cfo, NP, epsilon=1e-10, maxiter=1000):
    """
    Particle swarm is a population based optimizer which uses force and
    momentum to select candidate points.

    *cfo* is the cost function object.  This is a dictionary which contains
    the following keys:

        *cost* is the function to be optimized.  If *parallel_cost* exists,
        it should accept a list of points, not just a single point on each
        evaluation.

        *n* is the problem dimension

        *x0* is the initial point

        *x1* and *x2* are lower and upper bounds for each parameter

        *monitor* is a callable which is called each iteration using
        *callback(step, x, fx, k)*, where *step* is the iteration number,
        *x* is the population, *fx* is value of the cost function for each
        member of the population and *k* is the index of the best point in
        the population.

        *f_opt* is the target value of the optimization

    *NP* is the number of fit parameters

    *epsilon* is the convergence criterion.

    *abort_test* is a callable which indicates whether an external processes
    requests the fit to stop.

    *maxiter* is the maximum number of generations

    Returns success, num_evals, f(x_best), x_best.
    """

    if "parallel_cost" in cfo:
        mapper = lambda v: asarray(cfo["parallel_cost"](v.T), "d")
    else:
        mapper = lambda v: asarray(list(map(cfo["cost"], v.T)), "d")
    monitor = cfo.get("monitor", print_every_five)

    n = cfo["n"]
    c1 = 2.8
    c2 = 1.3
    phi = c1 + c2
    K = 2 / abs(2 - phi - sqrt(phi * phi - 4 * phi))

    X = rand(n, NP)  # will hold original vectors
    V = zeros((n, NP))

    # CREATE FIRST GENERATION WITH LEGAL PARAMETER VALUES AND EVALUATE COSTS
    rn1 = rand(n, NP)
    # m th member of the population
    for m in range(0, NP):
        extend = cfo["x2"] - cfo["x1"]
        X[:, m] = cfo["x1"] + extend * X[:, m]
        V[:, m] = 2 * rn1[:, m] * extend - extend

    if "x0" in cfo:
        X[:, 0] = cfo["x0"]
    f = mapper(X)

    n_feval = NP
    P = X[:]

    f_best, i_best = min(zip(f, count()))
    for L in range(2, maxiter + 1):
        rn2 = rand(n, NP)
        for i in range(0, NP):
            # r = rand(2)
            r = rn2[:, i]
            V[:, i] = V[:, i] + r[0] * c1 * (P[:, i] - X[:, i]) + r[1] * c2 * (P[:, i_best] - X[:, i])
            V[:, i] = K * V[:, i]

            X[:, i] = X[:, i] + V[:, i]

        f_temp = mapper(X)
        idx = f_temp < f
        f[idx] = f_temp[idx]
        P[:, idx] = X[:, idx]

        n_feval = n_feval + NP

        # CHECKING STOPPING CRITERIA
        f_best, i_best = min(zip(f, count()))
        if abs(cfo["f_opt"] - f_best) < epsilon:
            satisfied_sc = 1
            return satisfied_sc, n_feval, f_best, X[:, i_best]

        monitor(L, X, f, i_best)

    satisfied_sc = 0
    return satisfied_sc, n_feval, f_best, X[:, i_best]


def example_call(optimizer=random_lines):
    from numpy.random import seed

    seed(1)
    cost = lambda x: x[0] ** 2 + x[1] ** 2
    n = 2
    x1 = -5 * ones(n)
    x2 = 5 * ones(n)
    f_opt = 0
    cfo = {"cost": cost, "n": n, "x1": x1, "x2": x2, "f_opt": f_opt}

    NP = 10 * n
    satisfied_sc, n_feval, f_best, x_best = optimizer(cfo, NP)
    print(satisfied_sc, "n:%d" % n_feval, f_best, x_best)


def main():
    print("=== Random Lines")
    example_call(random_lines)
    print("=== Particle Swarm")
    example_call(particle_swarm)


if __name__ == "__main__":
    main()