diff --git a/ackley_mc.py b/ackley_mc.py new file mode 100644 index 0000000000000000000000000000000000000000..e0a123459e2f51de5341517eb9217857b7aba535 --- /dev/null +++ b/ackley_mc.py @@ -0,0 +1,123 @@ +#!/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) + + + + + \ No newline at end of file