ackley_mc.py

#!/usr/bin/env python3

import re,sys
import numpy as np
import random
import matplotlib.pyplot as plt


def ackley(x, a=20, b=0.2, c=2 * np.pi):
    x = np.array(x)
    n = len(x)
    x_squared_term = np.sum(x*x)
    x_cos_term = np.sum(np.cos(c * x))
    term1 = -a * np.exp(-b * np.sqrt(x_squared_term/ n))
    term2 = -np.exp(x_cos_term/n)
    y = term1 + term2 + a + np.exp(1)
    return y

def get_trial_x(x, x_delta):
    n = len(x)
    idim = random.randint(0,n-1)
    x_trial = x.copy()
    step = x_delta * (2 * random.uniform(0,1)-1)
    x_trial[idim]+= step
    return x_trial

def monte_carlo(n_step, x, x_delta, E, temp, fout, c_mult):
    energies, x_history = [],[]
    for step in range(n_step):
        x_trial = get_trial_x(x,x_delta)
        E_trial = ackley(x_trial)
        flag = False
        if (E_trial <= E):
            flag = True
        elif (E_trial > E):
            delta_E = E - E_trial
            if (random.uniform(0,1) < np.exp(delta_E/temp)):
                flag = True
        if (flag):
            x = x_trial
            E = E_trial
            fout.write(f"{step},{E},{','.join(map(str, x))}\n")
            energies.append(E)
            x_history.append(x.copy())
        temp = temp * c_mult
    return energies,x_history


def ackley_mc(filename):

    # try to read parameter file
    try:
        stream = open(filename)
    except IOError as err:
        sys.stderr.write('{}: {}\n'.format(sys.argv[0],err))
        exit(1)
    
    # assign variables
    for line in stream: 
        line = line.strip()
        if line.startswith('ini_temp'):
            ini_temp = float(re.findall(r'ini_temp\s+([+-]?\d+(?:\.\d+)?)', line)[0])
        elif line.startswith('final_temp'):
            final_temp = float(re.findall(r'final_temp\s+([+-]?\d+(?:\.\d+)?)', line)[0])
        elif line.startswith('n_step'):
            n_step = int(re.findall(r'n_step\s+(\d+)', line)[0])
        elif line.startswith('x_delta'):
            x_delta = float(re.findall(r'x_delta\s+([+-]?\d+(?:\.\d+)?)', line)[0])
        elif line.startswith('seed'):
            seed = int(re.findall(r'seed\s+(\d+)', line)[0])
        elif line.startswith('foutname'):
            foutname = re.findall(r'foutname\s+(.+)', line)[0] #gibt str anstatt list zurück
        elif line.startswith('x_ini'):
            x_ini = [float(x) for x in re.findall(r'-?\d+(?:\.\d+)?', line)]

    # set random seed number
    random.seed(seed)

    # set starting values
    x = x_ini.copy()
    E = ackley(x)
    c = -(1/n_step)*np.log(final_temp/ini_temp)
    c_mult = np.exp(-c)

    # open output file
    with open(foutname, 'w') as fout:
        fout.write("Step,Energy,dim\n")
        energies, x_history = monte_carlo(n_step, x, x_delta, E, ini_temp, fout, c_mult)

    #plotting
    steps = list(range(len(energies)))
    x_history = np.array(x_history)

    #plot energy over steps
    plt.figure(figsize=(10, 4))
    plt.subplot(1, 2, 1)
    plt.plot(steps, energies, label="Ackley Energy")
    plt.xlabel("Step")
    plt.ylabel("Energy")
    plt.title("Energy developement")
    plt.grid(True)
    #plot x_val over steps
    plt.subplot(1, 2, 2)
    for i in range(x_history.shape[1]):
        plt.plot(steps, x_history[:, i], label=f"x[{i}]")
    plt.xlabel("Step")
    plt.ylabel("x value")
    plt.title("X-coordinates developement")
    plt.legend()
    plt.grid(True)

    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    filename = sys.argv[1]
    ackley_mc(filename)