DIFFERENTIAL_EVOLUTION

=DIFFERENTIAL_EVOLUTION(func_expr, bounds, [strategy], [maxiter], [popsize], [tol], [mutation], [recombination], [seed], [polish])

=DIFFERENTIAL_EVOLUTION("x[0]**2 + x[1]**2", {-5,5;-5,5}, "best1bin", 500, 20, 0.0001, 0.6, 0.8, 123, TRUE)

=DIFFERENTIAL_EVOLUTION("(x[0]-2)**2 + (x[1]+1)**2", {-6,6;-6,6}, "best2bin", , , , {0.5,1.2}, 0.9, 321)

=DIFFERENTIAL_EVOLUTION("abs(x[0]) + abs(x[1])", {-3,3;-3,3}, , , , 0.000001, , , 999, FALSE)

=DIFFERENTIAL_EVOLUTION("(x[0]-1)**2 + (x[1]-0.5)**2 + (x[2]+2)**2", {-5,5;-5,5;-5,5}, , 200, 10, , , , 2024)

import math
import re
from typing import List, Optional, Sequence, Union
 
import numpy as np
from scipy.optimize import differential_evolution as scipy_differential_evolution
 
 
def differential_evolution(
    func_expr: str,
    bounds: List[List[float]],
    strategy: str = 'best1bin',
    maxiter: int = 1000,
    popsize: int = 15,
    tol: float = 0.01,
    mutation: Union[float, Sequence[float]] = 0.5,
    recombination: float = 0.7,
    seed: Optional[int] = None,
    polish: bool = True,
):
    """Minimize a multivariate function using differential evolution.
 
    This function wraps scipy.optimize.differential_evolution. For more information, see:
    https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html
 
    Args:
        func_expr: Objective expression using x[i] for variable access.
        bounds: 2D list of [min, max] pairs for each variable.
        strategy: Differential evolution strategy identifier. Defaults to 'best1bin'.
        maxiter: Maximum number of generations. Defaults to 1000.
        popsize: Population size multiplier. Defaults to 15.
        tol: Relative convergence tolerance. Defaults to 0.01.
        mutation: Mutation constant or [min, max] pair. Defaults to 0.5.
        recombination: Crossover probability. Defaults to 0.7.
        seed: Optional integer seed for reproducibility.
        polish: Whether to perform local search polishing. Defaults to True.
 
    Returns:
        [[x1, x2, ..., objective]] on success, or an error message string.
 
    This example function is provided as-is without any representation of accuracy.
    """
    if not isinstance(func_expr, str):
        return "Invalid input: func_expr must be a string."
    if func_expr.strip() == "":
        return "Invalid input: func_expr must be a non-empty string."
 
    func_expr = re.sub(r'\^', '**', func_expr)
 
    if "x" not in func_expr:
        return "Invalid input: func_expr must reference variable x (e.g., x[0])."
 
    if not isinstance(bounds, list) or len(bounds) == 0:
        return "Invalid input: bounds must be a 2D list of [min, max] pairs."
 
    processed_bounds = []
    for idx, bound in enumerate(bounds):
        if not isinstance(bound, list) or len(bound) != 2:
            return "Invalid input: each bound must be a [min, max] pair."
        try:
            lower = float(bound[0])
            upper = float(bound[1])
        except (TypeError, ValueError):
            return "Invalid input: bounds must contain numeric values."
        if lower > upper:
            return f"Invalid input: lower bound must not exceed upper bound for variable index {idx}."
        processed_bounds.append((lower, upper))
 
    dimension = len(processed_bounds)
 
    valid_strategies = {
        'best1bin', 'best1exp', 'rand1exp', 'randtobest1exp', 'currenttobest1exp',
        'best2exp', 'rand2exp', 'randtobest1bin', 'best2bin', 'rand2bin',
        'rand1bin', 'rand1binwithmem', 'rand1expwithmem'
    }
    if strategy not in valid_strategies:
        return f"Invalid input: strategy must be one of: {', '.join(sorted(valid_strategies))}"
 
    try:
        maxiter = int(maxiter)
        popsize = int(popsize)
        tol = float(tol)
        recombination = float(recombination)
    except (TypeError, ValueError):
        return "Invalid input: maxiter, popsize must be integers; tol, recombination must be floats."
 
    if maxiter <= 0:
        return "Invalid input: maxiter must be positive."
    if popsize <= 0:
        return "Invalid input: popsize must be positive."
    if tol <= 0:
        return "Invalid input: tol must be positive."
    if not (0.0 < recombination <= 1.0):
        return "Invalid input: recombination must be in (0, 1]."
 
    mutation_param: Union[float, Sequence[float]]
    if isinstance(mutation, (int, float)):
        try:
            mutation_value = float(mutation)
        except (TypeError, ValueError):
            return "Invalid input: mutation must be numeric or a pair of floats."
        if not (0 < mutation_value <= 2):
            return "Invalid input: mutation scalar must be in (0, 2]."
        mutation_param = mutation_value
    elif isinstance(mutation, (list, tuple)):
        if len(mutation) != 2:
            return "Invalid input: mutation range must contain exactly two values."
        try:
            mutation_low = float(mutation[0])
            mutation_high = float(mutation[1])
        except (TypeError, ValueError):
            return "Invalid input: mutation range must be numeric."
        if not (0 <= mutation_low <= mutation_high <= 2):
            return "Invalid input: mutation range must satisfy 0 <= low <= high <= 2."
        mutation_param = (mutation_low, mutation_high)
    else:
        return "Invalid input: mutation must be a float or length-2 list/tuple."
 
    rng_seed = None
    if seed is not None:
        try:
            rng_seed = int(seed)
        except (TypeError, ValueError):
            return "Invalid input: seed must be an integer."
 
    if not isinstance(polish, bool):
        return "Invalid input: polish must be a boolean."
 
    safe_globals = {
        "math": math,
        "np": np,
        "numpy": np,
        "__builtins__": {},
    }
    safe_globals.update({
        name: getattr(math, name)
        for name in dir(math)
        if not name.startswith("_")
    })
    safe_globals.update({
        "sin": np.sin,
        "cos": np.cos,
        "tan": np.tan,
        "asin": np.arcsin,
        "arcsin": np.arcsin,
        "acos": np.arccos,
        "arccos": np.arccos,
        "atan": np.arctan,
        "arctan": np.arctan,
        "sinh": np.sinh,
        "cosh": np.cosh,
        "tanh": np.tanh,
        "exp": np.exp,
        "log": np.log,
        "ln": np.log,
        "log10": np.log10,
        "sqrt": np.sqrt,
        "abs": np.abs,
        "pow": np.power,
        "pi": math.pi,
        "e": math.e,
    })
 
    x0_vector = [np.mean([b[0], b[1]]) for b in processed_bounds]
 
    try:
        initial_eval = eval(func_expr, safe_globals, {"x": x0_vector})
        float(initial_eval)
    except Exception as exc:
        return f"Error: Invalid expression at initial guess: {exc}"
 
    def objective(x_vector):
        try:
            local_x = [float(val) for val in np.atleast_1d(x_vector)]
            result = eval(func_expr, safe_globals, {"x": local_x})
            numeric_value = float(result)
        except Exception as exc:
            raise ValueError(f"Error evaluating func_expr: {exc}")
        if not math.isfinite(numeric_value):
            raise ValueError("Objective evaluation produced NaN or infinity.")
        return numeric_value
 
    try:
        result = scipy_differential_evolution(
            objective,
            bounds=processed_bounds,
            strategy=strategy,
            maxiter=maxiter,
            popsize=popsize,
            tol=tol,
            mutation=mutation_param,
            recombination=recombination,
            seed=rng_seed,
            polish=polish,
        )
    except ValueError as exc:
        return f"Error during differential_evolution: {exc}"
    except Exception as exc:
        return f"Error during differential_evolution: {exc}"
 
    if not result.success:
        return f"Error during differential_evolution: {result.message}"
 
    try:
        solution_vector = [float(val) for val in np.atleast_1d(result.x)]
    except (TypeError, ValueError):
        return "Error during differential_evolution: could not convert solution to floats."
 
    objective_value = float(result.fun)
    return [solution_vector + [objective_value]]