from utils import joukowski
    
# Instantiate the airfoil generator class, joukowski.WingGenerator
wing_generator = joukowski.WingGenerator()
    
# Display integer index table of airfoil generator parameters
wing_generator.show_index_table()

# Perform airfoil generation (mapping transformation). xi0, eta0, and alpha are specified by integer indices. See the output of wing_generator.show_index_table() for the index table.
wing = wing_generator.generate_wing(10, 15, 20)
    
# Show the constructed airfoil profile
wing.draw()
    
# Show the mapping parameters for transformation
wing_generator.show_parameters("My wing")

from utils import lbm
    
# Simulation grid size (nx: horizontal, ny: vertical)
nx, ny = 400, 160
    
# Instantiate the lbm.Solver class.
solver = lbm.Solver(nx, ny)
    
# Add the airfoil, wing, constructed above as a boundary condition (profile model)
solver.add_model(wing)
    
# Run 3000-time-step simulation
solver.integrate(steps=3000)
    
# Calculate the negative value of the lift-drag ratio, the objective function, from the simulation results (for this FMQA, minimize this value)
cost = -solver.fy_average / solver.fx_average
    
# Print the evaluated value of the objective function
print(f"{cost=:.2f}")

import os
import torch
import numpy as np
    
    
def seed_everything(seed=0):
    os.environ["PYTHONHASHSEED"] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.deterministic = True

from amplify import FixstarsClient
from datetime import timedelta
    
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=1000)  # timeout is 1000 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # If you use Amplify in a local environment, enter the Amplify API token.

import torch.nn as nn
    
    
class TorchFM(nn.Module):
    def __init__(self, d: int, k: int):
        super().__init__()
        self.V = nn.Parameter(torch.randn(d, k), requires_grad=True)
        self.lin = nn.Linear(
            d, 1
        )  # The first and second terms on the right-hand side are fully connected network
    
    def forward(self, x):
        out_1 = torch.matmul(x, self.V).pow(2).sum(1, keepdim=True)
        out_2 = torch.matmul(x.pow(2), self.V.pow(2)).sum(1, keepdim=True)
        out_inter = 0.5 * (out_1 - out_2)
        out_lin = self.lin(x)
        out = out_inter + out_lin
        return out

from torch.utils.data import TensorDataset, DataLoader
from sklearn.model_selection import train_test_split
from typing import Type
    
import copy
    
    
def train(
    X: np.ndarray,
    y: np.ndarray,
    model_class: Type[nn.Module],
    model_params: dict[str, int | float],
    batch_size=1024,
    epochs=3000,
    criterion=nn.MSELoss(),
    optimizer_class=torch.optim.AdamW,
    opt_params={"lr": 1},
    lr_sche_class=None,
    lr_sche_params=None,
):
    X_tensor, y_tensor = (
        torch.from_numpy(X).float(),
        torch.from_numpy(y).float(),
    )
    indices = np.array(range(X.shape[0]))
    indices_train, indices_valid = train_test_split(
        indices, test_size=0.2, random_state=42
    )
    
    train_set = TensorDataset(X_tensor[indices_train], y_tensor[indices_train])
    valid_set = TensorDataset(X_tensor[indices_valid], y_tensor[indices_valid])
    loaders = {
        "train": DataLoader(train_set, batch_size=batch_size, shuffle=True),
        "valid": DataLoader(valid_set, batch_size=batch_size, shuffle=False),
    }
    
    model = model_class(**model_params)
    best_model_wts = copy.deepcopy(model.state_dict())
    optimizer = optimizer_class(model.parameters(), **opt_params)
    scheduler = None
    if lr_sche_class is not None:
        scheduler = lr_sche_class(optimizer, **lr_sche_params)
    best_score = 1e18
    for _ in range(epochs):
        losses = {"train": 0.0, "valid": 0.0}
    
        for phase in ["train", "valid"]:
            if phase == "train":
                model.train()
            else:
                model.eval()
    
            for batch_x, batch_y in loaders[phase]:
                optimizer.zero_grad()
                out = model(batch_x).T[0]
                loss = criterion(out, batch_y)
                losses[phase] += loss.item() * batch_x.size(0)
    
                with torch.set_grad_enabled(phase == "train"):
                    if phase == "train":
                        loss.backward()
                        optimizer.step()
    
            losses[phase] /= len(loaders[phase].dataset)
    
        with torch.no_grad():
            model.eval()
            if best_score > losses["valid"]:
                best_model_wts = copy.deepcopy(model.state_dict())
                best_score = losses["valid"]
        if scheduler is not None:
            scheduler.step()
    
    with torch.no_grad():
        model.load_state_dict(best_model_wts)
        model.eval()
    return model

def gen_training_data_flow_onehot(D: int, N0: int, true_func, input_generator):
    X = np.zeros((N0, D))
    y = np.zeros(N0)
    print("")
    print("###################################################")
    print(" Making initial training data")
    print("###################################################")
    for i in range(N0):
        print(f"Generating training data set #{i}.")
        # Randomly determine input vector consisting of (0 or 1)
    
        x = input_generator()
        # If the same input vector has already exist, regenerate the input values
        is_identical = True
        while is_identical:
            is_identical = False
            for j in range(i + 1):
                if np.all(x == X[j, :]):
                    x = input_generator()
                    is_identical = True
                    break
        # Evaluate input value x with objective function f(x) and copy input-output pair (x,f(x)) to teacher data
        X[i, :] = x
        y[i] = true_func(x)
        print(f"------------------------")
    return X, y

from amplify import VariableGenerator, solve, one_hot, sum
import matplotlib.pyplot as plt
import sys
    
# FMQA for non-binary inputs
    
    
class FMQA_NB:
    def __init__(self, D, N, N0, k, true_func, client, encoder_decoder) -> None:
        assert N0 < N
        self.D = D
        self.N = N
        self.N0 = N0
        self.k = k
        self.true_func = true_func
        self.client = client
        self.nbins_array = encoder_decoder.nbins_array
        self.encoder = encoder_decoder.encoder
        self.decoder = encoder_decoder.decoder
        self.y = None
    
    # A member function that repeatedly performs (N-N0)x FMQA based on the training data with adding new training data
    def cycle(self, X, y, log=False) -> np.ndarray:
        # Weight for the one-hot constraints
        constraint_weight = max([1, int(np.abs(y).max() + 0.5)])
        print("")
        print("###################################################")
        print(f" Starting FMQA cycles... {constraint_weight=}")
        print("###################################################")
    
        pred_x = X[0]
        pred_y = 1e18
        for i in range(self.N - self.N0):
            print(f"FMQA Cycle #{i} ")
            try:
                x_hat = self.step(X, y, constraint_weight)
            except RuntimeError:
                sys.exit(f"Unknown error, i = {i}")
            # If an input value identical to the found x_hat already exists in the current training data set, a neighboring value is used as the new x_hat.
            is_identical = True
            while is_identical:
                is_identical = False
                for j in range(i + self.N0):
                    if np.all(x_hat == X[j, :]):
                        # Decode binary vector to integer vector and copy to inputss
                        inputs = self.decoder(x_hat)
                        # Get the surrounding values (inputs is an integer vector)
                        inputs += np.random.randint(-1, 2, len(inputs))
                        for i_inp in range(len(inputs)):
                            if inputs[i_inp] < 0:
                                inputs[i_inp] = 0
                            elif inputs[i_inp] > self.nbins_array[i_inp] - 1:
                                inputs[i_inp] = self.nbins_array[i_inp] - 1
                        # Encode from integer vector to binary vector and copy to x_hat
                        x_hat = self.encoder(inputs)
                        if log:
                            print(f"{i=}, Identical x is found, {x_hat=}")
                        is_identical = True
                        break
            # Evaluate objective function f() with x_hat
            y_hat = self.true_func(x_hat)
            # Add the input-output pair [x_hat, y_hat] to the training data set
            X = np.vstack((X, x_hat))
            y = np.append(y, y_hat)
            # Copy the input-output pair to [pred_x, pred_y] when the evaluated value of the objective function updates the minimum value
            if pred_y > y_hat:
                pred_y = y_hat
                pred_x = x_hat
                print(f"variable updated, {pred_y=}")
            else:
                print("")
            print(f"------------------------")
        self.y = y
        return pred_x
    
    # Member function to perform one FMQA cycle
    def step(self, X, y, constraint_weight) -> np.ndarray:
        # Train FM
        model = train(
            X,
            y,
            model_class=TorchFM,
            model_params={"d": self.D, "k": self.k},
            batch_size=8,
            epochs=2000,
            criterion=nn.MSELoss(),
            optimizer_class=torch.optim.AdamW,
            # Set scheduler to reduce learning rate with number of epochs
            opt_params={"lr": 0.5},
            lr_sche_class=torch.optim.lr_scheduler.StepLR,
            lr_sche_params={"step_size": 50, "gamma": 0.9},
        )
        # Extract FM parameters from the trained FM model
        v, w, w0 = list(model.parameters())
        v = v.detach().numpy()
        w = w.detach().numpy()[0]
        w0 = w0.detach().numpy()[0]
        # Solve a QUBO problem using a quantum annealing or Ising machine
        gen = VariableGenerator()  # Declare a variable generator
        q = gen.array("Binary", self.D)  # Generate binary decision variables
        cost = self.__FM_as_QUBO(q, w0, w, v)  # Define FM as a QUBO equation
        # Construct one-hot constraints for non-binary variables
        constraint_list = []
        ista = 0
        iend = 0
        for i in range(len(self.nbins_array)):
            iend += self.nbins_array[i]
            constraint_list.append(one_hot(q[ista:iend]))
            ista = iend
        constraints = sum(constraint_list)
        # Add up the objective function and constraints, and pass it to the Fixstars Amplify solver
        model = cost + constraint_weight * constraints
        result = solve(model, self.client)  # Pass the objective function to Amplify
        if len(result.solutions) == 0:
            raise RuntimeError("No solution was found.")
        q_values = q.evaluate(result.best.values)
        return q_values
    
    # A function that defines FM as a QUBO equation from FM parameters. As with the previously defined TorchFM class, the formula is written as per the acquisition function form of g(x).
    def __FM_as_QUBO(self, x, w0, w, v):
        lin = w0 + (x.T @ w)
        out_1 = np.array([(x * v[:, i]).sum() ** 2 for i in range(self.k)]).sum()
        # Note that x[j] = x[j]^2 since x[j] is a binary variable in the following equation.
        out_2 = np.array([(x * v[:, i] * v[:, i]).sum() for i in range(self.k)]).sum()
        return lin + (out_1 - out_2) / 2
    
    """The sum_poly used in __FM_as_QUBO above is inefficient in terms of computation speed and
    memory. In the case of FM, where the interaction terms of the decision variables are generally
    nonzero, the following implementation using BinaryMatrix is more efficient. Here, the quadratic
    terms in BinaryMatrix correspond to the non-diagonal terms represented by the upper triangular
    matrix, so x(1/2) for the quadratic terms in the FM formula is unnecessary. Also, although x is
    taken as an argument just to match the function signature with __FM_as_QUBO above (implementation
    using sum_poly), it is not needed in this implementation using BinaryMatrix.
    def __FM_as_QUBO(self, x, w0, w, v):
        out_1_matrix = v @ v.T
        out_2_matrix = np.diag((v * v).sum(axis=1))
        matrix = BinaryMatrix(out_1_matrix - out_2_matrix + np.diag(w))
        # Do not forget to put the constant term w0 in the second argument of BinaryQuadraticModel.
        model = BinaryQuadraticModel(matrix, w0)
        return model
    """
    
    # A function to plot the history of i-th objective function evaluations performed within the initial training data construction (blue) and during FMQA cycles (red).
    def plot_history(self):
        assert self.y is not None
        fig = plt.figure(figsize=(6, 4))
        plt.plot(
            [i for i in range(self.N0)],
            self.y[: self.N0],
            marker="o",
            linestyle="-",
            color="b",
        )  # Objective function evaluation values at the time of initial training data generation (random process)
        plt.plot(
            [i for i in range(self.N0, self.N)],
            self.y[self.N0 :],
            marker="o",
            linestyle="-",
            color="r",
        )  # Objective function evaluation values during the FMQA cycles (FMQA cycle process)
        plt.xlabel("i-th evaluation of f(x)", fontsize=18)
        plt.ylabel("f(x)", fontsize=18)
        plt.tick_params(labelsize=18)
        return fig

class EncDecOneHot:
    def __init__(self, D: int, nbins_array):
        self.D = D
        self.nbins_array = nbins_array
    
    def gen_random_input(self):
        ista = 0
        iend = 0
        x = np.zeros(self.D, int)
        for i in range(len(self.nbins_array)):
            iend += self.nbins_array[i]
            idx = np.random.randint(ista, iend, 1)[0]
            x[idx] = 1
            ista = iend
        return x
    
    def encoder(self, inputs):
        if len(inputs) != len(self.nbins_array):
            raise RuntimeError("inputs should be the same length as nbins_array!")
        x = np.zeros(self.D, int)
        i_offset = 0
        for i in range(len(inputs)):
            x[inputs[i] + i_offset] = 1
            i_offset += self.nbins_array[i]
        return x
    
    def decoder(self, x):
        bit_locs = np.where(x == 1)[0]
        if len(bit_locs) != len(self.nbins_array):
            raise RuntimeError("bit_locs should be the same length as nbins_array!")
        i_offset = 0
        for i in range(1, len(bit_locs)):
            i_offset += self.nbins_array[i - 1]
            bit_locs[i] -= i_offset
        return bit_locs

# Initialize random seed value
seed_everything()
    
# Instantiate airfoil generator class, joukowski.WingGenerator
wing_generator = joukowski.WingGenerator()
    
# Display integer index table of airfoil generator parameters
wing_generator.show_index_table()
    
# One-hot encoding bin size for each variable
nbins_array = np.array(
    [wing_generator.nbins_xi0, wing_generator.nbins_eta0, wing_generator.nbins_alpha]
)
    
# Overall size of binary variables (= size of binary decision variables handled by QUBO)
D = wing_generator.nbins_xi0 + wing_generator.nbins_eta0 + wing_generator.nbins_alpha
    
# Instantiate one-hot encoder-decoder class (EncDecOneHot)
enc_dec_one_hot = EncDecOneHot(D, nbins_array)
    
# Counter for the number of objective function evaluations
n_eval = 0
    
# Objective function (construct airfoil → run simulation → compute objective function value)
    
    
def obj_func(x):
    global n_eval
    # Decode binary variables to integers (mapping parameter indices)
    i_xi0, i_eta0, i_alpha = enc_dec_one_hot.decoder(x)
    # Perform airfoil generation (mapping transformation). xi0, eta0, alpha are integer indices.
    # See the output of wing_generator.show_index_table() for table of indices
    wing = wing_generator.generate_wing(i_xi0, i_eta0, i_alpha)
    # Show the mapping parameters for the transformation
    wing_generator.show_parameters("My wing")
    # Simulation grid size (nx: horizontal, ny: vertical)
    nx, ny = 400, 160
    # Instantiate the lbm.Solver class.
    solver = lbm.Solver(nx, ny)
    # Add the airfoil, wing, constructed above as a boundary condition (airfoil model)
    solver.add_model(wing)
    # Run 5000-time-step simulation
    solver.integrate(steps=5000)
    # Negative value of the lift-drag ratio
    cost = -solver.fy_average / solver.fx_average
    # If the fluid simulation diverges
    if np.isnan(cost):
        cost = 0
    print(f"#{n_eval} evaluation finished. cost:{cost:.2f}")
    n_eval += 1
    return cost
    
    
N = 5  # 40 is used for Sec. 3.3.  # Number of times the objective function can be evaluated
N0 = 3  # 20 is used for Sec. 3.3.  # Number of samples of initial training data
k = 20  # Dimension of vector in FM (hyperparameters)
    
# Generate initial training data
X, y = gen_training_data_flow_onehot(D, N0, obj_func, enc_dec_one_hot.gen_random_input)
    
# Instantiate FMQA
fmqa_solver = FMQA_NB(D, N, N0, k, obj_func, client, enc_dec_one_hot)
    
# Perform FMQA cycle
pred_x = fmqa_solver.cycle(X, y)
    
# Output optimization results
print("###################################################")
print(" Optimal input values and corresponding output")
print("###################################################")
print("pred x:", pred_x)
print("pred value:", obj_func(pred_x))

fig = fmqa_solver.plot_history()

N = 40  # Number of times the objective function can be evaluated
N0 = 20  # Number of samples of initial training data

Black-Box Optimization of the Airfoil Geometry with Fluid Flow Simulation¶

1. Airfoil definition and fluid flow simulation¶

1.1. Airfoil generator class¶

1.2. Fluid flow simulator class¶

2. FMQA program implementation (non-binary decision variables)¶

2.1. Random seed initialization¶

2.2. Configuration of Amplify client¶

2.3. Implementing FM with PyTorch¶

2.4. Construction of initial training data¶

2.5. Execution class for FMQA cycle¶

2.6. One-hot encoder/decoder¶

3. Search for the optimal airfoil profile¶

3.1. FMQA execution example¶

3.2. Transition of objective function values during the optimization process¶

3.3. Example output from this FMQA sample program¶

3.4. Commentary¶