Adaptive Scenario Subset Selection CMA-ES for minmax problem

as3_cma_es.py

import math
import numpy as np
import random
from scipy.stats import bernoulli
from scipy import stats
from scipy.stats import chi2
import ddcma
"""To use AS3-CMA-ES, dd-cma.py should be imported from 
https://gist.github.com/youheiakimoto/1180b67b5a0b1265c204cba991fa8518"""


def proj_p(p, EPS):
    """Projetion of P
    Parameters
    ----------
    p : ndarray (1D)
        probability vector
    EPS : float, >= 0
        smallest possible probability
        
    Returns
    -------
    probability vector in [EPS/(M-1), 1 - EPS] with sum(p) = 1
    """
    m = len(p)
    for i in range (m):
        p[i] = min(1, max(EPS, p[i]))
    pp = p
    return pp


def as3cma(f, init_x, init_p, init_sigma, xb, cp, epe, lam, mode, gamma, cmaterm, 
             MAXEVAL=10000, MAXITR=10000, MINSIGMA=1e-10, MAXBI=1, log='test.txt'):
    
    """min_x max_y f(x, y)

    Parameters
    ----------
    f : callable
        f(x, y) 
        x: design vector (1D array of size n)
        y: scenario index (integer between 0 and m-1)
    init_x : 1D array of size n
        initial x vector
    init_p : 1D array of size m
        initial p vector, whose elements are positive and sum to one
    init_sigma : 1D array, positive
        coordinate-wise std for x update
    cp : float, 0 < cp <= 1
        learning rate for p update
    eps : float, 0 <= eps < 1
        minimum p-value
    lam : float, 4 + int(3 * log(n))
        population size for the CMA-ES
    mode : 
        if mode = 0 : sampling all scenarios
        if 0 < mode < 1 : Adaptive support scenario subset and mode works as eta
        if 1 <= mode : Adaptive support scenario subset with fixed lambda_s = mode 
    gamma :
        probability where the cumulative density function of a chi-square distribution with n degrees of freedom for the CMA-ES
    cmaterm : termination criteria for the CMA-ES, 1D array of size 4
        cmaterm[0] : minimum sigma
        cmaterm[1] : maximum ratio between maximum and minimum eigen value of covariance matrix
        cmaterm[2] : maximum ratio between maximum and minimum standard deviasion
        cmaterm[3] : minimum value for maximum standard deviasion
    MAXEVAL : integer
        maximum number of f-calls
    MAXITR : integer
        maximum number of iterations
    MAXBI : number of restarts
        MAXBI = 0 : non restart
        MAXBI = 1 : simple restart with changing mean vector
        if number of restart is less than MAXBI, population size is double, otherwize population size is fixed. 
        
    Returns
    -------
    xbest : best solution candidate
    fbest : max_y f(xbest, y) 
    """
    # Initialization
    # Use DD-CMA 
    cma = ddcma.DdCma(init_x, init_sigma, lam, flg_variance_update=False)
    p = np.array(init_p, copy=True)
    n = cma.N
    m = len(p)
    
    neval = 0
    conv = 0
    
    nominalx=[]
    res = np.zeros(1+1+1+n+n+n+m)
    binum = 0
    # Main Loop
    for t in range(MAXITR):         
        S = np.zeros(len(p), dtype=int)
        tS = np.zeros(len(p), dtype=int)

        # all select
        if 0==mode:
            S = np.arange(0,m)
            lambda_s = m
            
#       scenario samplin
        if 0<mode and mode<1:
            lambda_s = 0
            for i in range (m):
                ss = np.random.binomial(1, p=p[i])
                if (ss==1):
                    S[lambda_s] = i
                    lambda_s += 1
                
            if(lambda_s==0):
                S[lambda_s] = random.randint(0,m)
                lambda_s=1
        
#       fixed number of scenarios
        if mode>=1:
            pp = p/sum(p)
            S = np.random.choice(m, size=int(np.floor(mode)), replace=False, p=pp)
            lambda_s=int(np.floor(mode))
            
        # CMA Sampling
        arx, ary, arz = cma.sample()
        
        # Evaluation
        fx_arr = np.array([[f(x, S[i]) for i in range(lambda_s)] for x in arx])
        fx_max_idx = np.argmax(fx_arr, axis=1)
        fx_max = np.max(fx_arr, axis=1)
        neval += lambda_s * cma.lam

        # X Update
        idx = np.argsort(fx_max)
        cma.update(idx, arx, ary, arz)
        
        # compute dp
        # positive update
        dp = np.zeros(m)
        ct = 1
        for i in range(cma.lam):
            z=np.sum(arz[i,:]**2)
            if stats.chi2.cdf(z, n)<=gamma:
                dp[S[fx_max_idx[i]]] += cp
                if ct==1:
                    cpidx=np.ones(1, dtype='int') * fx_max_idx[i]
                    ct +=1
                cpidx=np.append(cpidx, fx_max_idx[i])

        ## compute negative learning rate cn    
        idxp =np.unique(fx_max_idx)
        cn = cp
        if 0<mode and mode<1:
            cn =  cp*(mode*cma.lam)/max(m - mode*cma.lam - 1,  mode*cma.lam)
        if 1<=mode:
            cn = cp*cma.lam/(lambda_s-len(idxp))

        tS[:] = S[:]
        for i in range(len(idxp)):
            tS = tS[tS!=S[idxp[i]]]
            
        # negative update 
        for i in range(lambda_s-len(idxp)):
            dp[tS[i]] -= cn
            
        # p update
        p+=dp
        
        ## cat p
        p = proj_p(p, eps)


        ## output to log files
        idx = 0
        res[idx] = neval; idx += 1
        res[idx] = min(fx_max); idx += 1
        res[idx] = cma.sigma; idx += 1
        res[idx:idx+n] = cma.xmean; idx += n
        res[idx:idx+n] = cma.coordinate_std; idx += n
        res[idx:idx+n] = cma.S; idx += n
        res[idx:idx+m] = p ; idx += m
        with open(log, 'ba') as flog:
            np.savetxt(flog, res, newline=' ')
            flog.write(b"\n")
        

#     restart
        if cma.sigma < cmaterm[0]:
            conv = 2
        if np.max(cma.S) / np.min(cma.S) > cmaterm[1]:
            conv = 3
        if np.max(cma.coordinate_std) / np.min(cma.coordinate_std) > cmaterm[2]:
            conv = 4
        if np.max(cma.coordinate_std) < cmaterm[3]:
            conv = 5
        
        if conv > 1:
            bid = np.argmin(fx_max)
            nominalx.append(arx[bid])
            print ('restart  '+str(neval)+'  '+str(conv))
            init_x = np.random.rand(n)*np.random.rand(n)*(xb[1] - xb[0]) + xb[0]
            init_sigma = np.ones(n)* np.ones(n)*(xb[1] - xb[0])/4
            p = np.ones(m) * init_p
            binum += 1
            if binum < MAXBI:
                lam = cma.lam * 2
            else:
                lam = cma.lam
            # restart DD-CMA
            cma = ddcma.DdCma(init_x, init_sigma, lam, flg_variance_update=False)
            conv = 0
            
# Termination
        if neval >= MXEVAL:
            break 
            
# check best solution candidate
    bid = np.argmin(fx_max)
    nominalx.append(arx[bid])   
    fsize=len(nominalx)
    flist=np.zeros(fsize)
    for i in range(fsize):
        f_arr =  np.array([f(nominalx[i], S[j]) for j in range(m)])
        flist[i] = np.max(f_arr)
    fbest = np.min(flist)
    bid = np.argmin(flist)
    xbest = nominalx[bid]
    return fbest, xbest


if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt

    m =30
    n = 10
    nsupport = 20 ## number of support scenarios
    xb = [-4, 4]
    def ftest(x, y):
        yy = range(1, m+1, 1)
        N = len (x)
        omega = np.pi/nsupport
        romega = np.zeros(N)
        romega[0] = (np.cos(omega*yy[y]))
        romega[1] = (np.sin(omega*yy[y]))
        alpha = np.tan(omega)**-2
        fxy = (np.dot(x, x) - (1 + alpha) * np.dot(x, romega)**2)
        if (nsupport<yy[y]):
            omega = 2*np.pi/(m-nsupport)
            romega[0] = (np.cos(omega*(yy[y]-nsupport)))
            romega[1] = (np.sin(omega*(yy[y]-nsupport)))            
            fxy = 2*np.dot(x-romega, x-romega) -8 
        return fxy
    
    ## Experiment
    f = ftest

    ## initial parameters for the CMA-ES
    init_x = np.random.rand(n)*(xb[1] - xb[0]) + xb[0]
    cmalam = (4 + int(3 * math.log(len(init_x))))
    init_p = np.ones(m) * 0.1
    init_sigma = np.ones(n)*(xb[1] - xb[0])/4
    cmaterm = [1e-8, 1e+8, 1e+8, 1e-8]

    ## initial parameters for AS3
    gamma = 0.99
    cp = 0.3
    eps= 1/m
    eta = 0.3
    mode = eta

    MXEVAL = int(1e+05) 
    MXITR = int(MXEVAL/cmalam)
    MAXBI = 1
    log = 'AS3-CMA-test.txt'
    with open(log, 'w'): 
        pass
    
    fbest, xbest = as3cma( f, init_x, init_p, init_sigma, xb, cp, eps, cmalam, mode, gamma, cmaterm,
        MAXEVAL=MXEVAL, MAXITR=MXITR, MAXBI=MAXBI, log=log)

    print ('best function value = ', fbest)
    print ('best x = ', xbest)

    dat = np.loadtxt(log)
    ## Plot
    plt.figure(figsize=(15,10))
    ax1 = plt.subplot(231)
    ax1.plot(dat[:,0], dat[:,1], label='f values')
    ax1.plot(dat[:,0], dat[:,2], label='$\sigma$')
    ax1.set_yscale("log")
    plt.xlabel("#f-calls")
    plt.legend()
    plt.grid()

    ax2 = plt.subplot(232)   
    plt.plot(dat[:,0], np.sum(dat[:,-m:], axis=1))
    plt.xlabel("#f-calls")
    plt.ylabel("$\sum p_t$")
    plt.grid()

    ax3 = plt.subplot(233)
    ax3.plot(dat[:,0], dat[:,3:3+n])
    plt.xlabel("#f-calls")
    plt.ylabel("$m^t$")
    plt.grid()

    ax4 = plt.subplot(234)
    ax4.plot(dat[:,0], dat[:,3+n:3+2*n])
    ax4.set_yscale("log")
    plt.xlabel("#f-calls")
    plt.ylabel("std")
    plt.grid()

    ax5 = plt.subplot(235)
    ax5.plot(dat[:,0], dat[:,3+2*n:3+3*n])
    plt.xlabel("#f-calls")
    plt.ylabel("Eigen value")
    plt.grid()

    plt.tight_layout()
    plt.savefig('AS3-CMA.pdf')