model_generate.py

#python:
import numpy as np
import random as rd
import csv
import math
import h5py
from scipy.ndimage import rotate, zoom
from user_libs.antennas.GSSI import antenna_like_GSSI_2000

# ----- Helper Functions -----

def get_value(key, file_path):
    with open(file_path, 'r') as file:
        for line in file:
            if line.startswith(key):
                    # Splitting by '=' and stripping to remove any whitespace
                    return line.split("=")[1].strip()
    return None


def clear_csv_file(file_path):
    """
    Clear the contents of a CSV file at the specified file path.
    """
    with open(file_path, 'w') as file:
        pass  # Opening in 'w' mode and closing will clear the file

def append_to_csv(file_path, data):
    """
    Append a row of data to a CSV file.
    """
    with open(file_path, 'a', newline='') as file:
        writer = csv.writer(file)
        writer.writerow(data)

def edit_file_values_in_place(file_path, new_er, new_se):
    """
    Edit the values of er and se in the file.
    """
    try:
        with open(file_path, 'r') as file:
            content = file.read()

        # Replace the values of er and se
        lines = content.split('\n')
        for i, line in enumerate(lines):
            if line.startswith('#material'):
                parts = line.split()
                parts[1] = str(new_er)
                parts[2] = str(new_se)
                lines[i] = ' '.join(parts)
                break

        with open(file_path, 'w') as file:
            file.write('\n'.join(lines))
    except FileNotFoundError:
        print("File not found: ", file_path)


def process_and_save_data(input_file_path, output_file_path, thetas, axes, scale_factors,
                          x_posA, y_posA, pos_z, rock_type):
    """
    Reads a 3D array from an HDF5 file, applies rotation and scaling,
    and saves the processed data back to an HDF5 file.
    """
    with h5py.File(input_file_path, 'r') as file:
        data = file['data'][:]

    # Apply rotation
    for theta, axis in zip(thetas, axes):
        data = rotate(data, angle=np.degrees(theta), axes=axis, reshape=True, order=1, mode='nearest')

    # Apply scaling
    data = zoom(data, scale_factors, order=1, mode='nearest')

    with h5py.File(output_file_path, 'w') as file:
        file.create_dataset('data', data=data)


def GSSI_2(x_posA, y_posA, z_posA):
    """
    Call the GSSI 2 GHz antenna model.
    """
    antenna_like_GSSI_2000(x_posA, y_posA, z_posA, resolution=0.001)


def semi_empirical(T, S, C, rs, rb, fw):
    """
    Semi-empirical dielectric mixing model for soil.
    Based on Peplinski et al. mixing model.
    """
    alpha = 0.65
    ew_real = 80.1
    ew_imag = 0
    es_real = rs
    porosity = 0.5

    eb_real = (1 + (porosity ** alpha) * ((ew_real ** alpha) - 1) + ((1 - porosity) ** alpha) * ((es_real ** alpha) - 1)) ** (1 / alpha)
    eb_imag = 0

    mv = fw
    if mv > 0:
        e_real = (1 + (eb_real ** alpha - 1) + (mv ** alpha) * (ew_real ** alpha) - mv) ** (1 / alpha)
        e_imag = (mv ** alpha) * (ew_imag ** alpha)
    else:
        e_real = eb_real
        e_imag = eb_imag

    conductivity = 2 * np.pi * 2e9 * 8.854e-12 * e_imag

    print("#material: {} {} 1 0 my_sand_semi".format(e_real, conductivity))


# ---- State Variables ----

axes = [(1, 0), (2, 0), (2, 1)]
rock_d1 = 0
rock_d2 = 0
T, S, C, rs, rb = 20, 50, 30, 5, 2.5

# ---- Functions for scan setup ----

def reset_scan_variables():
    global landmine_status, actual_bullets, rock_count
    landmine_status = 0
    actual_bullets = 0
    rock_count = 0
    print("6th scan: Empty model. Resetting landmine_status, rock_count, and actual_bullets to 0")


def configure_landmine_and_sand_properties(ground_depth, landmine_depth_variation, mean_bullet_count, antenna_z_variation):
    global landmine_position, sand_rel_permittivity, sand_rel_conductivity, actual_bullets, landmine_status, z_posA, fw
    landmine_position = (ground_depth + 0.045) + (landmine_depth_variation * rd.uniform(0, 1))
    sand_rel_permittivity = rd.uniform(2.5, 25)
    sand_rel_conductivity = rd.uniform(0.001, 1)
    actual_bullets = mean_bullet_count  # Assuming you want to reset this at the start of each set
    landmine_status = 1  # Ensure landmine is included for the next 5 scans
    z_posA = 0.037 + (antenna_z_variation * rd.uniform(-1, 1))
    fw = rd.uniform(0.0001,0.12) # water fraction

    # Assuming you might want to use these values outside, you can return them

def setup_rock_variables():
    global rock_count, rock1_x_posA, rock1_y_posA, rock_1_pos_z, rock2_x_posA, rock2_y_posA, rock_2_pos_z
    global rock_1_type, rock_2_type, rock_1_permittivity, rock_1_conductivity, rock_2_permittivity, rock_2_conductivity
    global rock_1_thetas, rock_2_thetas, rock_1_scale_factors, rock_2_scale_factors
    global input_file_path_1, output_file_path_1, input_file_path_2, output_file_path_2, file1_path, file2_path

    # rock_count = int(rd.uniform(1, 2))
    rock_count = 2

    rock1_x_posA = Domain_x/2 - 0.0215 + (antenna_x_variation * rd.uniform(-1.1, 1.1))
    rock1_y_posA = Domain_y/2 - 0.0215 + (antenna_y_variation * rd.uniform(-1.1,1.1))
    rock_1_pos_z = ground_depth + rd.uniform(0.03, 0.1)

    rock2_x_posA = Domain_x/2 - 0.0215 + (antenna_x_variation * rd.uniform(-1.1, 1.1))
    rock2_y_posA = Domain_y/2 - 0.0215 + (antenna_y_variation * rd.uniform(-1.1, 1.1))
    rock_2_pos_z = ground_depth + rd.uniform(0.03, 0.1)

    rock_1_type, rock_2_type = int(rd.uniform(1, 3)), int(rd.uniform(1, 3))
    rock_1_permittivity, rock_2_permittivity = rd.uniform(5, 10), rd.uniform(5, 10)
    rock_1_conductivity, rock_2_conductivity = rd.uniform(0.0000001, 1), rd.uniform(0.0000001, 1)

    rock_1_thetas, rock_2_thetas = [rd.uniform(0, np.pi) for _ in range(3)], [rd.uniform(0, np.pi) for _ in range(3)]
    rock_1_scale, rock_2_scale = rd.uniform(0.5, 1.3), rd.uniform(0.5, 1.3)
    rock_1_scale_factors, rock_2_scale_factors = [rock_1_scale] * 3, [rock_2_scale] * 3

    input_file_path_1, output_file_path_1 = f'Rocks/rock{rock_1_type}.h5', f'Rocks/rock_ro_{rock_1_type}.h5'
    input_file_path_2, output_file_path_2 = f'Rocks/rock{rock_2_type}.h5', f'Rocks/rock_ro_{rock_2_type}.h5'
    file1_path, file2_path = f'Rocks/rock{rock_1_type}_new_material.txt', f'Rocks/rock{rock_2_type}_new_material.txt'

    edit_file_values_in_place(file1_path, rock_1_permittivity, rock_1_conductivity)
    edit_file_values_in_place(file2_path, rock_2_permittivity, rock_2_conductivity)

def calculate_casing_positions(Domain_x, Domain_y, ground_depth, bullet_depth, length, b_rad):
    pitch_angle = rd.uniform(-45, 45)  # casing orientation
    yaw_angle = rd.uniform(-180, 180)  # casing orientation
    bpos_x = Domain_x / 2 + rd.uniform(-1, 1) * ((Domain_x / 2) - 0.09)  # casing position
    bpos_y = Domain_y / 2 + rd.uniform(-1, 1) * ((Domain_x / 2) - 0.09)  # casing position
    bpos_z = (ground_depth + bullet_depth) + (rd.uniform(-1, 1) * 0.025)  # casing position

    pitch = math.radians(pitch_angle)
    yaw = math.radians(yaw_angle)

    # Calculate the direction vector based on pitch and yaw angles
    direction_x = math.cos(yaw) * math.cos(pitch)
    direction_y = math.sin(yaw) * math.cos(pitch)
    direction_z = math.sin(pitch)

    # Calculate the end position of the cylinder
    bpos2_x = bpos_x + length * direction_x
    bpos2_y = bpos_y + length * direction_y
    bpos2_z = bpos_z + length * direction_z

    bpos3_x = bpos_x + 0.002 * direction_x
    bpos3_y = bpos_y + 0.002 * direction_y
    bpos3_z = bpos_z + 0.002 * direction_z

    print("#cylinder: {} {} {} {} {} {} {} pec y".format(bpos_x, bpos_y, bpos_z, bpos2_x, bpos2_y, bpos2_z, b_rad))
    print("#cylinder: {} {} {} {} {} {} {} free_space y".format(bpos3_x, bpos3_y, bpos3_z, bpos2_x, bpos2_y, bpos2_z, b_rad - 0.0015))


Dataset_definitions_file_path = "Dataset definitions3.txt"
# Convert these values to floats as they represent numerical data
Training_examples_TGT = int(get_value("Training examples Target", Dataset_definitions_file_path))
Training_examples_MIX = int(get_value("Training examples Mix", Dataset_definitions_file_path))
Training_examples_FA = int(get_value("Training examples False Alarm", Dataset_definitions_file_path))
scans_per_example = int(get_value("Scans per training example", Dataset_definitions_file_path))
mean_bullet_count = int(get_value("Mean bullet count in model", Dataset_definitions_file_path))
Domain_x = float(get_value("Domain_x", Dataset_definitions_file_path))
Domain_y = float(get_value("Domain_y", Dataset_definitions_file_path))
Domain_z = float(get_value("Domain_z", Dataset_definitions_file_path))
Batch_name = get_value("Batch name", Dataset_definitions_file_path)
# Extracting key parameters and assigning variables
Training_scans_TGT = int(scans_per_example * Training_examples_TGT)
Training_scans_MIX = int(scans_per_example * Training_examples_MIX)
Training_scans_FA = int(scans_per_example * Training_examples_FA)
total_scans = int(Training_scans_TGT + Training_scans_MIX + Training_scans_FA)

csv_file_path = f"Current_code/Outputs/Model_data_SAND_MIX_{Batch_name}_.csv"


# gprMax setup parameters

print("#domain: {} {} {}".format(Domain_x, Domain_y, Domain_z))
print("#dx_dy_dz: 0.001 0.001 0.001")
print("#time_window: 6e-9")

# Main logic


if (current_model_run - 1) % (scans_per_example + 1) == scans_per_example:

    reset_scan_variables()

else:
    if (current_model_run - 1) % (scans_per_example + 1) == 0:

        setup_rock_variables()
        # Reset or initialize landmine and bullet properties
        configure_landmine_and_sand_properties(ground_depth, landmine_depth_variation, mean_bullet_count, antenna_z_variation)

# Creating the primary medium for the model

print("#material: {} {} 1 0 my_sand".format(sand_rel_permittivity, sand_rel_conductivity))
semi_empirical(T, S, C, rs, rb, fw)
print("#box: 0 0 {} {} {} {} my_sand".format(ground_depth, Domain_x, Domain_y, Domain_z))


# inserting bullet cases

if current_model_run <= Training_scans_TGT:
    actual_bullets = 0

if Training_scans_TGT < current_model_run <= (Training_scans_TGT + Training_scans_MIX):
    for i in range(actual_bullets):
        if landmine_position >= (ground_depth + 0.045):

            bullet_depth = 0.025

            calculate_casing_positions(Domain_x, Domain_y, ground_depth, bullet_depth, length, b_rad)

    # Log the geometry objects read


if (Training_scans_TGT + Training_scans_MIX) < current_model_run:
    if rock_count >= 1:
        process_and_save_data(input_file_path_2, output_file_path_2, rock_2_thetas, axes, rock_2_scale_factors,
                                            rock1_x_posA, rock1_y_posA, rock_1_pos_z, rock_1_type)
        print("#geometry_objects_read: {} {} {} Rocks/rock_ro_{}.h5 Rocks/rock{}_new_material.txt ".format(rock1_x_posA, rock1_y_posA, rock_1_pos_z, rock_1_type, rock_1_type))
        rock_d1 = rock_1_pos_z

        if rock_count == 2:
            process_and_save_data(input_file_path_2, output_file_path_2, rock_2_thetas, axes, rock_2_scale_factors,
                                                rock2_x_posA, rock2_y_posA, rock_2_pos_z, rock_2_type)
            print("#geometry_objects_read: {} {} {} Rocks/rock_ro_{}.h5 Rocks/rock{}_new_material.txt ".format(rock2_x_posA, rock2_y_posA, rock_2_pos_z, rock_2_type, rock_2_type))
            rock_d2 = rock_2_pos_z
            print("rock two being printed")
            print(rock_d2)


    for i in range(actual_bullets):

        bullet_depth = 0.045
        calculate_casing_positions(Domain_x, Domain_y, ground_depth, bullet_depth, length, b_rad)
# Set landmine target

if current_model_run <= (Training_scans_TGT + Training_scans_MIX):

    if landmine_status == 1:

        print("#geometry_objects_read: {} {} {} PMN.h5 PMN_materials.txt".format(((Domain_x/2) -0.05), ((Domain_y/2) -0.07), landmine_position))

    else:
        landmine_position = 0
else:

    landmine_position = 0

# Load antenna model

GSSI_2(x_posA, y_posA, z_posA)

# print("#geometry_view: 0 0 0 {} {} {} 0.001 0.001 0.001 TRAIN_SAND_ALL_ROCK_ n".format(Domain_x, Domain_y, Domain_z))

# Adding key parameters to a csv file for post analysis

if (current_model_run - 1) == 0:
    clear_csv_file(csv_file_path)
    append_to_csv(csv_file_path, ["current_model_run", "sand_rel_permittivity", "sand_rel_conductivity", "x_posA", "y_posA", "z_posA", "landmine_position", "actual_bullets","rock_d1","rock_d2"])
    new_data = [current_model_run, sand_rel_permittivity, sand_rel_conductivity, x_posA, y_posA, z_posA, landmine_position, actual_bullets, rock_d1, rock_d2]
    append_to_csv(csv_file_path, new_data)

else:
    new_data = [current_model_run, sand_rel_permittivity, sand_rel_conductivity, x_posA, y_posA, z_posA, landmine_position, actual_bullets, rock_d1, rock_d2]
    append_to_csv(csv_file_path, new_data)


#end_python: