OccamBayesian/occam_bayesian.py at master · mortele/OccamBayesian · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
# /usr/local/bin/python3
import os
import shlex
import subprocess
import matplotlib.pyplot as plt
from bayes_opt import BayesianOptimization
from file_read_backwards import FileReadBackwards
from visualize_progress import plot_gp


def _replace_in_file(file_name, key, value):
    """Create a new file with one change from the given file_name
    """
    tmp_name = file_name + '_tmp'
    old_name = file_name + '_old'
    with open(file_name, 'r') as in_file, open(tmp_name, 'w') as out_file:
        next_ = False
        for line in in_file:
            if next_:
                out_file.write(str(value) + '\n')
                next_ = False
            else:
                if key in line:
                    next_ = True
                out_file.write(line)
    os.rename(file_name, old_name)
    os.rename(tmp_name, file_name)


def _read_total_pressure():
    with FileReadBackwards('fort.7') as in_file:
        for _ in range(50):
            line = in_file.readline()
            if line:
                line = line.split()
                if len(line) > 2:
                    if line[1] == 'total' and line[2] == 'press':
                        pressure = float(line[0])
            else:
                break
    return pressure


def occam_parameters(steps=None, kappa=None):
    """Search replace numbers in OCCAM fort.1 and fort.3 files to specify
    simulation details.
    """
    _replace_in_file('fort.1', 'number_of_steps', steps)
    _replace_in_file('fort.3', 'compressibility', kappa)


def occam_function(path, executable_name='occamcg'):
    """The 'black-box' function we optimize, w.r.t. the parameters adjusted in the
    occam_parameters function.
    """
    subprocess.call(shlex.quote(os.path.join(path, executable_name)),
                    stdout=subprocess.DEVNULL,
                    shell=False)
    return _read_total_pressure()


def occam_optimize(path, steps, kappa, init_points=10):
    """Wrapper function for performing the entire optimization procedure.
    """
    target_pressure = 26.1514

    def opt_target(x):
        """Input for the BayesianOptimization procedure from the bayes_opt
        package.

        The Bayesian optimization is a maximization procedure, so we have to
        define a cost function which takes a maximum where the difference
        between the target pressure and the measured pressure is smallest.
        """
        occam_parameters(steps=steps, kappa=x)
        pressure = occam_function(path)
        # cost = 1.0 / np.sqrt(np.abs(pressure - target_pressure) + 0.5)
        # cost = -abs(pressure - target_pressure)
        # cost = 1.0 / ((pressure - target_pressure)**2 + 0.5)
        cost = - (pressure - target_pressure)**2
        return cost

    p_bounds = {'x': (kappa[0], kappa[1])}
    opt = BayesianOptimization(f=opt_target,
                               pbounds=p_bounds,
                               verbose=2,
                               random_state=92898)
    opt.maximize(init_points=init_points, n_iter=0, kappa=5)

    for _ in range(20):
        opt.maximize(init_points=0, n_iter=5)
        # plot_gp(opt, x, y, set_xlim=(kappa[0], kappa[1]))
        plot_gp(opt,
                set_xlim=(kappa[0], kappa[1]))
        plt.show()
    print(opt.max)


if __name__ == '__main__':
    OCCAM_PATH = os.path.join('..', 'OCCAM', 'bin')
    occam_optimize(OCCAM_PATH, 1000, (0.01, 1.0), init_points=10)