blackbox.py

import sys
import multiprocessing as mp
import numpy as np
import scipy.optimize as op
from scipy.interpolate import Rbf as rbf
import datetime


def get_default_executor():
    """
    Provide a default executor (a context manager
    returning an object with a map method).

    This is the multiprocessing Pool object () for python3.

    The multiprocessing Pool in python2 does not have an __enter__
    and __exit__ method, this function provides a backport of the python3 Pool
    context manager.

    Returns
    -------
    Pool : executor-like object
        An object with context manager (__enter__, __exit__) and map method.
    """
    if (sys.version_info > (3, 0)):
        Pool = mp.Pool
        return Pool
    else:
        from contextlib import contextmanager
        from functools import wraps

        @wraps(mp.Pool)
        @contextmanager
        def Pool(*args, **kwargs):
            pool = mp.Pool(*args, **kwargs)
            yield pool
            pool.terminate()
        return Pool


def search_min(f, domain, budget, batch, resfile,
               rho0=0.5, p=1.0,
               executor=get_default_executor()):
    """
    Minimize given expensive black-box function and save results into text file.

    Parameters
    ----------
    f : callable
        The objective function to be minimized.
    domain : list of lists
        List of ranges for each parameter.
    budget : int
        Total number of function calls available.
    batch : int
        Number of function calls evaluated simultaneously (in parallel).
    resfile : str
        Text file to save results.
    rho0 : float, optional
        Initial "balls density".
    p : float, optional
        Rate of "balls density" decay (p=1 - linear, p>1 - faster, 0<p<1 - slower).
    executor : callable, optional
        Should have a map method and behave as a context manager.
        Allows the user to use various parallelisation tools
        as dask.distributed or pathos.
    """
    # space size
    d = len(domain)

    # adjusting the budget to the batch size
    if budget % batch != 0:
        budget = budget - budget % batch + batch
        print('[blackbox] FYI: budget was adjusted to be ' + str(budget))

    # default global-vs-local assumption (50-50)
    n = budget//2
    if n % batch != 0:
        n = n - n % batch + batch
    m = budget-n

    # n has to be greater than d
    if n <= d:
        print('[blackbox] ERROR: budget is not sufficient')
        return

    # go from normalized values (unit cube) to absolute values (box)
    def cubetobox(x):
        return [domain[i][0]+(domain[i][1]-domain[i][0])*x[i] for i in range(d)]

    # generating R-sequence
    points = np.zeros((n, d+1))
    points[:, 0:-1] = rseq(n, d)

    # initial sampling
    for i in range(n//batch):
        print('[blackbox] evaluating batch %s/%s (samples %s..%s/%s) @ ' % (i+1, (n+m)//batch, i*batch+1, (i+1)*batch, n+m) + \
        str(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")) + ' ...')

        with executor() as e:
            points[batch*i:batch*(i+1), -1] = list(e.map(f, list(map(cubetobox, points[batch*i:batch*(i+1), 0:-1]))))

    # normalizing function values
    fmax = max(abs(points[:, -1]))
    points[:, -1] = points[:, -1]/fmax

    # volume of d-dimensional ball (r = 1)
    if d % 2 == 0:
        v1 = np.pi**(d/2)/np.math.factorial(d/2)
    else:
        v1 = 2*(4*np.pi)**((d-1)/2)*np.math.factorial((d-1)/2)/np.math.factorial(d)

    # subsequent iterations (current subsequent iteration = i*batch+j)

    for i in range(m//batch):
        print('[blackbox] evaluating batch %s/%s (samples %s..%s/%s) @ ' % (n//batch+i+1, (n+m)//batch, n+i*batch+1, n+(i+1)*batch, n+m) + \
        str(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")) + ' ...')

        # sampling next batch of points
        fit = rbf(*np.transpose(points), function='cubic')
        points = np.append(points, np.zeros((batch, d+1)), axis=0)

        for j in range(batch):
            r = ((rho0*((m-1.-(i*batch+j))/(m-1.))**p)/(v1*(n+i*batch+j)))**(1./d)
            cons = [{'type': 'ineq', 'fun': lambda x, localk=k: np.linalg.norm(np.subtract(x, points[localk, 0:-1])) - r}
                    for k in range(n+i*batch+j)]
            while True:
                minfit = op.minimize(fit, np.random.rand(d), method='SLSQP', bounds=[[0., 1.]]*d, constraints=cons)
                if np.isnan(minfit.x)[0] == False:
                    break
            points[n+i*batch+j, 0:-1] = np.copy(minfit.x)

        with executor() as e:
            points[n+batch*i:n+batch*(i+1), -1] = list(e.map(f, list(map(cubetobox, points[n+batch*i:n+batch*(i+1), 0:-1]))))/fmax

    # saving results into text file
    points[:, 0:-1] = list(map(cubetobox, points[:, 0:-1]))
    points[:, -1] = points[:, -1]*fmax
    points = points[points[:, -1].argsort()]

    labels = [' par_'+str(i+1)+(7-len(str(i+1)))*' '+',' for i in range(d)]+[' f_value    ']
    np.savetxt(resfile, points, delimiter=',', fmt=' %+1.4e', header=''.join(labels), comments='')

    print('[blackbox] DONE: see results in ' + resfile + ' @ ' + str(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")))


def rseq(n, d):
    """
    Build R-sequence (http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/).

    Parameters
    ----------
    n : int
        Number of points.
    d : int
        Size of space.

    Returns
    -------
    points : ndarray
        Array of points uniformly placed in d-dimensional unit cube.
    """
    phi = 2
    for i in range(10):
        phi = pow(1+phi, 1./(d+1))

    alpha = np.array([pow(1./phi, i+1) for i in range(d)])

    points = np.array([(0.5 + alpha*(i+1)) % 1 for i in range(n)])

    return points