Feature/vispi final

fourier_generator.py

@@ -0,0 +1,73 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+class FourierSeriesGenerator:
+
+    def __init__(self, total_time=400, time_step=0.002):
+        self.default_total_time = total_time
+        self.default_time_step = time_step
+
+    @staticmethod
+    def _normalize(x):
+        return 2 * (x - np.min(x)) / (np.max(x) - np.min(x)) - 1
+
+    @staticmethod
+    def _rescale(x, min_value, max_value):
+        return x * (max_value - min_value) + min_value
+
+    def fourier_series(self, x, params, rescale_mag=600, rescale_amplitude=50):
+        x = self._normalize(x)
+
+        n, freq, amplitude, phase, trend, seasonality, = params
+        n = int(self._rescale(n, 0, 10))
+        freq = self._rescale(freq, 0, 10)
+        amplitude = self._rescale(amplitude, 0, 10)
+        phase = self._rescale(phase, 0, 10000)
+        trend = self._rescale(trend, -500, 500)
+        seasonality = self._rescale(seasonality, 0, 200)
+
+        sum = np.zeros_like(x)
+        for i in range(1, n + 1, 2):
+            term = (1 / i) * np.sin(2 * np.pi * freq * i * x + phase)
+            sum += term
+
+        y = amplitude * (2 / np.pi) * sum
+        if np.sum(y) == 0:
+            return np.zeros_like(x) + 600
+        else:
+            y = (y - np.min(y)) / (np.max(y) - np.min(y))
+            y = (y * rescale_amplitude) + rescale_mag
+
+        y += trend * x
+        y += seasonality * np.sin(2 * np.pi * x)
+        return y
+
+    def plot_and_save(self, params, base_path, iteration, total_time=None, time_step=None):
+        if total_time is None:
+            total_time = self.default_total_time
+        if time_step is None:
+            time_step = self.default_time_step
+
+        folder_name = f"Iteration_{iteration}"
+        save_directory = os.path.join(base_path, folder_name)
+        if not os.path.exists(save_directory):
+            os.makedirs(save_directory)  # Create directory if it doesn't exist
+
+        x = np.linspace(0, total_time, int(total_time / time_step))
+        y = self.fourier_series(x, params)
+
+        plt.plot(x, y, label=f'Curve')
+        plt.xlabel("Time")
+        plt.ylabel("Laser Power")
+        plt.legend()
+        image_path = os.path.join(save_directory, "plot.png")
+        plt.savefig(image_path)
+        plt.show()
+
+        output_string = "laser_power,time_elapsed\n"
+        for i in range(len(x)):
+            output_string += f"{y[i]:.15f},{x[i]:.2f}\n"
+        csv_path = os.path.join(save_directory, "data.csv")
+        with open(csv_path, "w") as f:
+            f.write(output_string)

gamma_model_simulator.py

@@ -0,0 +1,37 @@
+# gamma_model_simulator.py
+
+import os
+import numpy as np
+import cupy as cp
+import gamma.interface as rs
+from multiprocessing import Process
+import time
+
+class GammaModelSimulator:
+    def __init__(self, input_data_dir, sim_dir_name, laser_file, VtkOutputStep=1., ZarrOutputStep=0.02, outputVtkFiles=True, verbose=True):
+        self.input_data_dir = input_data_dir
+        self.sim_dir_name = sim_dir_name
+        self.laser_file = laser_file
+        self.VtkOutputStep = VtkOutputStep
+        self.ZarrOutputStep = ZarrOutputStep
+        self.outputVtkFiles = outputVtkFiles
+        self.verbose = verbose
+
+        self.sim_itr = None
+
+    def setup_simulation(self):
+        self.sim_itr = rs.FeaModel(
+                    input_data_dir=self.input_data_dir,
+                    geom_dir=self.sim_dir_name,
+                    laserpowerfile=self.laser_file,
+                    VtkOutputStep=self.VtkOutputStep,
+                    ZarrOutputStep=self.ZarrOutputStep,
+                    outputVtkFiles=self.outputVtkFiles,
+                    verbose=self.verbose)
+
+    def run_simulation(self):
+        if self.sim_itr:
+            self.sim_itr.run()
+        else:
+            raise ValueError("Simulation is not setup yet. Call setup_simulation() first.")
+

gamma_simulator_runner.py

@@ -0,0 +1,25 @@
+from gamma_model_simulator import GammaModelSimulator
+
+class GammaSimulatorRunner:
+    def __init__(self, input_data_dir, sim_dir_name, laser_file):
+        self.simulator = GammaModelSimulator(
+            input_data_dir=input_data_dir,
+            sim_dir_name=sim_dir_name,
+            laser_file=laser_file
+        )
+
+    def run(self):
+        # Set up the simulation
+        self.simulator.setup_simulation()
+
+        # Execute the simulation
+        self.simulator.run_simulation()
+
+if __name__ == "__main__":
+    # Define your parameters
+    INPUT_DATA_DIR = "/home/vnk3019/ded_dt_thermomechanical_solver/examples/data"
+    SIM_DIR_NAME = "thin_wall"
+    LASER_FILE = "/home/vnk3019/ded_dt_thermomechanical_solver/examples/data/laser_inputs/thin_wall/LP_2"
+
+    runner = GammaSimulatorRunner(INPUT_DATA_DIR, SIM_DIR_NAME, LASER_FILE)
+    runner.run()

run_analysis.py

@@ -0,0 +1,127 @@
+import numpy as np
+import pandas as pd
+import zarr
+from botorch.fit import fit_gpytorch_model
+import torch
+from botorch.models import SingleTaskGP
+from gpytorch.mlls import ExactMarginalLogLikelihood
+from botorch.optim import optimize_acqf
+from botorch.acquisition import UpperConfidenceBound
+from tqdm import tqdm
+import pandas as pd
+import numpy as np
+import os
+import warnings
+warnings.filterwarnings("ignore")
+
+# Define the neural network parameters in the following function:
+def objective(params, iteration_number, min_temp=893, max_temp=993):
+    from fourier_generator import FourierSeriesGenerator
+    from temperature_analyzer import TemperatureAnalyzer
+
+
+    # Generate and save the Fourier series
+    generator = FourierSeriesGenerator()
+    base_path = "/home/vnk3019/ded_dt_thermomechanical_solver/examples/data/laser_inputs/thin_wall"
+    generator.plot_and_save(params, base_path, iteration_number, total_time=300, time_step=0.002)
+
+    # Paths for the simulator
+    INPUT_DATA_DIR = "/home/vnk3019/ded_dt_thermomechanical_solver/examples/data"
+    SIM_DIR_NAME = "thin_wall"
+
+    # Modify LASER_FILE to reflect the correct iteration and CSV filename
+    LASER_FILE = os.path.join(base_path, f"Iteration_{iteration_number}", "data")
+
+    # Modify ZARR_LOCATION to reflect the correct iteration
+    ZARR_LOCATION = os.path.join(base_path, f"Iteration_{iteration_number}", "data.zarr", "ff_dt_temperature")
+
+    selected_nodes_list = ["45003"]  # As an example
+
+    analyzer = TemperatureAnalyzer(
+        INPUT_DATA_DIR,
+        SIM_DIR_NAME,
+        LASER_FILE
+    )
+
+    # Call function
+    avg_time = analyzer.run_simulation_and_analyze(ZARR_LOCATION, min_temp, max_temp, selected_nodes_list, plot_graph=True)
+    return torch.tensor(avg_time)  # Now returns a 1D tensor
+
+
+
+# Bayesian Optimization functions
+# -----------------------------------------------------------
+def initialize_model():
+    train_X = pd.read_excel("optimized_params.xlsx")
+    train_X_np = train_X[["n", "freq", "amplitude", "phase", "trend", "seasonality"]].values
+    train_X_torch = torch.tensor(train_X_np, dtype=torch.float32)
+    train_Y = pd.read_excel("avg_heat_treatment_times.xlsx")
+    train_Y_np = train_Y[["Average Heat Treatment Time"]].values
+    train_Y_torch = torch.tensor(train_Y_np, dtype=torch.float32)
+    gp = SingleTaskGP(train_X_torch.float(), train_Y_torch.float())
+    mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
+    fit_gpytorch_model(mll)
+    return gp
+
+def optimize(bounds, n_steps=50):
+    gp = initialize_model()
+    best_value = -float('inf')
+    best_params = None
+
+    param_history = []
+    value_history = []
+    uncertainty_history = []
+
+    # Create an empty dataframe with the required columns
+    df = pd.DataFrame(columns=['Parameters', 'Objective Value', 'Uncertainty'])
+
+    for i in tqdm(range(n_steps)):
+        UCB = UpperConfidenceBound(gp, beta=0.5)
+        candidate, _ = optimize_acqf(UCB, bounds=bounds, q=1, num_restarts=5, raw_samples=20)
+        candidate_numpy = candidate.detach().numpy().flatten()
+        new_Y = objective(candidate_numpy, iteration_number=i).unsqueeze(0).unsqueeze(-1)
+
+        variance = gp.posterior(candidate).variance
+        uncertainty_history.append(variance.item())
+
+        print(new_Y)
+        print(candidate_numpy)
+
+
+        if new_Y.item() > best_value:
+            best_value = new_Y.item()
+            best_params = candidate_numpy
+
+        param_history.append(candidate_numpy)
+        value_history.append(new_Y.item())
+
+        # Update the Gaussian Process model
+        gp = SingleTaskGP(
+            torch.cat([gp.train_inputs[0], torch.tensor(candidate_numpy.astype(np.float32)).unsqueeze(0)]),
+            torch.cat([gp.train_targets.unsqueeze(-1).float(), new_Y.float()], dim=0)
+        )
+        mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
+        fit_gpytorch_model(mll)
+
+        # Save the current iteration results to the dataframe
+        df = df.append({
+            'Current Best Parameters': best_params,
+            'Current Best Value': best_value, 
+            'Parameters': candidate_numpy.tolist(),
+            'Objective Value': new_Y.item(),
+            'Uncertainty': variance.item()
+        }, ignore_index=True)
+
+        # Save the dataframe to an Excel file
+        df.to_csv('bayesian_optimization_results.csv', index=False)
+
+    return gp, best_params, best_value, param_history, value_history, uncertainty_history
+# -----------------------------------------------------------
+
+if __name__ == "__main__":
+    input_size = 6  # Assuming 6 parameters
+    bounds = torch.tensor([[0]*input_size, [1]*input_size], dtype=torch.float32)
+    optimized_model, best_params, best_value, param_history, value_history, uncertainty_history = optimize(bounds)
+
+    print(f'Best parameters: {best_params}, Best value: {best_value}')
+    print('Uncertainty at each step:', uncertainty_history)