esa · gomezzz · May 6, 2022 · Apr 29, 2022 · May 2, 2022 · May 2, 2022
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
@@ -427,7 +427,7 @@ Training Properties
 - model_type (str) : "siren" or "nerf". What type of model architecture to use
 - iterations (int) : number of iterations in training the neural network
 - learning_rate (float) : The learning rate for the training of the neural network, a higher value leads to larger steps along the gradient. Reduce if problems occur.
-- type (str) : "classification" or "regression". For more details see Example_Notebook.ipynb
+- type (str) : "classification" or "regression". For more details see notebooks in the repository.
 - reg_loss_weight (float) :  weighting of the regularization loss, only used in classification to force the network to pick one material.
 - use_regularization_loss (bool) : Only relevant for classification, activates the regularization.
 

diff --git a/nidn/fdtd/grid.py b/nidn/fdtd/grid.py
@@ -134,7 +134,7 @@ def __init__(
         # save electric and magnetic field
         self.E = bd.zeros((self.Nx, self.Ny, self.Nz, 3))
         self.H = bd.zeros((self.Nx, self.Ny, self.Nz, 3))
-        self.absorption_factor = bd.zeros(self.Nx, self.Ny, self.Nz, 3)
+        self.absorption_factor = bd.zeros((self.Nx, self.Ny, self.Nz, 3))
 
         # save the inverse of the relative permittiviy and the relative permeability
         # these tensors can be anisotropic!

diff --git a/nidn/fdtd_integration/compute_spectrum_fdtd.py b/nidn/fdtd_integration/compute_spectrum_fdtd.py
@@ -1,7 +1,7 @@
 from dotmap import DotMap
 from tqdm import tqdm
-import torch
 from loguru import logger
+import torch
 
 from nidn.fdtd_integration.constants import FDTD_GRID_SCALE
 from nidn.utils.global_constants import UNIT_MAGNITUDE
@@ -11,25 +11,15 @@
     calculate_transmission_reflection_coefficients,
 )
 from .init_fdtd import init_fdtd
+from ..fdtd.backend import backend as bd
 
 
-def compute_spectrum_fdtd(permittivity, cfg: DotMap):
-    """Generates a spectrum of transmission and reflection coefficients for the specified wavelengths in the cfg, by using FDTD simulations.
+def _check_grid_scale(cfg):
+    """Checks whether grid scale allows specified sizes.
 
     Args:
-        permitivity (torch.tensor): Array of permittivity for each layer
-        cfg (DotMap): Configurations needed to perform the simulations
-
-    Returns:
-        tuple[array, array]: Transmission spectrum and reflection spectrum
+        cfg (DotMap): Configurations for the simulation
     """
-    transmission_spectrum = []
-    reflection_spectrum = []
-    physical_wavelengths, norm_freq = get_frequency_points(cfg)
-    logger.debug("Wavelenghts in spectrum : ")
-    logger.debug(physical_wavelengths)
-    logger.debug("Number of layers: ")
-    logger.debug(len(permittivity[0, 0, :, 0]))
 
     scaling = max(
         UNIT_MAGNITUDE / (cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE),
@@ -53,6 +43,26 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
                     cfg.PER_LAYER_THICKNESS[i],
                 )
             )
+
+
+def compute_spectrum_fdtd(permittivity, cfg: DotMap):
+    """Generates a spectrum of transmission and reflection coefficients for the specified wavelengths in the cfg, by using FDTD simulations.
+
+    Args:
+        permitivity (torch.tensor): Array of permittivity for each layer
+        cfg (DotMap): Configurations needed to perform the simulations
+
+    Returns:
+        tuple[array, array]: Transmission spectrum and reflection spectrum
+    """
+    transmission_spectrum = []
+    reflection_spectrum = []
+    physical_wavelengths, _ = get_frequency_points(cfg)
+    logger.debug("Wavelenghts in spectrum : ")
+    logger.debug(physical_wavelengths)
+    logger.debug("Number of layers: " + str(len(permittivity[0, 0, :, 0])))
+    _check_grid_scale(cfg)
+
     # For each wavelength, calculate transmission and reflection coefficents
     disable_progress_bar = logger._core.min_level >= 20
     for i in tqdm(range(len(physical_wavelengths)), disable=disable_progress_bar):
@@ -65,7 +75,7 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
             cfg,
             include_object=False,
             wavelength=physical_wavelengths[i],
-            permittivity=permittivity[:, :, :, i],
+            permittivity=None,
         )
         grid.run(cfg.FDTD_niter, progress_bar=False)
         transmission_free_space, reflection_free_space = _get_detector_values(
@@ -142,7 +152,7 @@ def _get_abs_value_from_3D_signal(signal):
     """
     signal = _average_along_detector(signal)
 
-    abs_value = torch.zeros(len(signal))
+    abs_value = bd.zeros(len(signal))
     for i in range(len(signal)):
         # Added 1e-16 to prevent gradient flow from breaking, without significantly changing the result
         squared_value = torch.square(signal[i] + 1e-16)
@@ -163,9 +173,9 @@ def _average_along_detector(signal):
     Returns:
         Array[timesteps, 3]: averaged signal along detector
     """
-    avg = torch.zeros([len(signal), 3])
+    avg = bd.zeros([len(signal), 3])
     for i in range(len(signal)):
-        s = torch.zeros(3)
+        s = bd.zeros(3)
         for p in signal[i]:
             s[0] += p[0] / len(signal[i])
             s[1] += p[1] / len(signal[i])

diff --git a/nidn/fdtd_integration/init_fdtd.py b/nidn/fdtd_integration/init_fdtd.py
@@ -1,5 +1,6 @@
 from dotmap import DotMap
 from loguru import logger
+import torch
 
 from ..fdtd import (
     AbsorbingObject,
@@ -30,81 +31,102 @@ def init_fdtd(cfg: DotMap, include_object, wavelength, permittivity):
     set_backend("torch")
     if (len(cfg.PER_LAYER_THICKNESS) == 1) and (cfg.N_layers > 1):
         cfg.PER_LAYER_THICKNESS = cfg.PER_LAYER_THICKNESS * cfg.N_layers
-    # The scaling is the number of grid points per unit magnitude. This is the maximum of the reation between the unit magnitude and 1/10th of the smallest wavelength,
+    thicknesses = torch.tensor(cfg.PER_LAYER_THICKNESS)
+    # The scaling is the number of grid points per unit magnitude. This is the maximum of the relation between the unit magnitude and 1/10th of the smallest wavelength,
     # and a constant which is defaulted to 10. If this scaling becomes too low, i.e. below 2, there might be some errors in creating the grid,
-    # as there is too feew grid points for ceartian elements to be placed correctly.
-    scaling = max(
-        UNIT_MAGNITUDE / (cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE),
-        cfg.FDTD_min_gridpoints_per_unit_magnitude,
+    # as there are too few grid points for certain elements to be placed correctly.
+    scaling = torch.maximum(
+        torch.tensor(
+            UNIT_MAGNITUDE / (cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE)
+        ),
+        torch.tensor(cfg.FDTD_min_gridpoints_per_unit_magnitude),
     )
 
-    x_grid_size = int(
+    x_grid_size = (
         scaling
         * (
             cfg.FDTD_pml_thickness * 2
             + cfg.FDTD_free_space_distance * 2
-            + sum(cfg.PER_LAYER_THICKNESS)
+            + torch.sum(thicknesses)
         )
-    )
+    ).int()
+
     y_grid_size = 3
     logger.debug(
         "Initializing FDTD grid with size {} by {} grid points, with a scaling factor of {} grid points per um".format(
             x_grid_size, y_grid_size, scaling
         )
     )
     grid = Grid(
-        (x_grid_size, y_grid_size, 1),
+        (x_grid_size.item(), y_grid_size, 1),
         grid_spacing=UNIT_MAGNITUDE / scaling,
         permittivity=1.0,
         permeability=1.0,
     )
     grid = _add_boundaries(grid, int(cfg.FDTD_pml_thickness * scaling))
+
+    # Determine detector positions
+    detector_x = (
+        cfg.FDTD_pml_thickness + cfg.FDTD_free_space_distance + torch.sum(thicknesses)
+    )
+
+    detector_y = torch.tensor(cfg.FDTD_pml_thickness + cfg.FDTD_free_space_distance)
+    detector_x *= scaling
+    detector_y *= scaling
+    detector_x += cfg.FDTD_gridpoints_from_material_to_detector
+    detector_y -= cfg.FDTD_gridpoints_from_material_to_detector
+
+    logger.trace(
+        "Placing detectors at "
+        + str(detector_x.int().item())
+        + ","
+        + str(detector_y.int().item())
+    )
+
     grid, t_detector, r_detector = _add_detectors(
         grid,
-        int(
-            scaling
-            * (
-                cfg.FDTD_pml_thickness
-                + cfg.FDTD_free_space_distance
-                + sum(cfg.PER_LAYER_THICKNESS)
-            )
-            + cfg.FDTD_gridpoints_from_material_to_detector
-        ),
-        int(
-            scaling * (cfg.FDTD_pml_thickness + cfg.FDTD_free_space_distance)
-            - cfg.FDTD_gridpoints_from_material_to_detector
-        ),
+        detector_x.int().item(),
+        detector_y.int().item(),
+    )
+
+    source_x = (cfg.FDTD_source_position[0] * scaling).int()
+    source_y = (cfg.FDTD_source_position[1] * scaling).int()
+
+    logger.trace(
+        "Placing source at "
+        + str(source_x.int().item())
+        + ","
+        + str(source_y.int().item())
     )
 
     grid = _add_source(
         grid,
-        int(cfg.FDTD_source_position[0] * scaling),
-        int(cfg.FDTD_source_position[1] * scaling),
+        source_x.item(),
+        source_y.item(),
         wavelength / SPEED_OF_LIGHT,
         cfg.FDTD_pulse_type,
         cfg.FDTD_source_type,
     )
+
     if include_object:
-        for i in range(cfg.N_layers):
-            x_start = cfg.FDTD_pml_thickness + cfg.FDTD_free_space_distance
-            x_end = x_start
-            if i == 0:
-                x_end += cfg.PER_LAYER_THICKNESS[0]
-            elif i == cfg.N_layers - 1:
-                x_start += sum(cfg.PER_LAYER_THICKNESS[:i])
-                x_end += sum(cfg.PER_LAYER_THICKNESS)
-            else:
-                x_start += sum(cfg.PER_LAYER_THICKNESS[:i])
-                x_end += sum(cfg.PER_LAYER_THICKNESS[: i + 1])
+        x_start = torch.tensor(cfg.FDTD_pml_thickness + cfg.FDTD_free_space_distance)
+        x_end = x_start + thicknesses[0]
+        for i in range(permittivity.shape[2]):
+            # TODO: Implement possibility for patterned grid, currently uniform layer is used
             grid = _add_object(
                 grid,
-                int(scaling * x_start),
-                int(scaling * x_end),
-                permittivity[0][0][
-                    i
-                ],  # TODO: Implement possibility for patterned grid, currently uniform layer is used
+                (scaling * x_start).int(),
+                (scaling * x_end).int(),
+                permittivity[0][0][i],
                 frequency=SPEED_OF_LIGHT / wavelength,
             )
+
+            if i < permittivity.shape[2] - 1:
+                x_start += thicknesses[i]
+                x_end += thicknesses[i + 1]
+            else:
+                break
+
     return grid, t_detector, r_detector
 
 

diff --git a/nidn/materials/data/zirconium.csv → nidn/materials/data/removed/zirconium.csv b/nidn/materials/data/zirconium.csv → nidn/materials/data/removed/zirconium.csv
diff --git a/nidn/materials/find_closest_material.py b/nidn/materials/find_closest_material.py
@@ -3,22 +3,20 @@
 import torch
 
 
-def _find_closest_material(eps, r
A377
un_cfg):
+def _find_closest_material(eps, run_cfg, material_collection):
     """Finds the closest matching material and distance to that (in epsilon)
     from a given epsilon grid.
 
     Args:
         eps (torch.tensor): A tensor of epsilon values (Nx,Ny,N_layers,N_freq).
         run_cfg (DotMap): Run configuration.
+        material_collection (MaterialCollection): The material collection.
 
     Returns:
         tuple: The closest materials and distances to that.
     """
     Nx, Ny, N_layers, N_freq = run_cfg.Nx, run_cfg.Ny, run_cfg.N_layers, run_cfg.N_freq
 
-    # Initiate material collection
-    material_collection = MaterialCollection(run_cfg.target_frequencies)
-
     comparisons = torch.zeros(
         [
             Nx,

diff --git a/nidn/materials/material_collection.py b/nidn/materials/material_collection.py
@@ -13,11 +13,12 @@ class MaterialCollection:
     material_names = []
     N_materials = 0
 
-    def __init__(self, target_frequencies):
+    def __init__(self, target_frequencies, to_skip=None):
         """Initalizes the MaterialCollection. Loads all materials from the data folder.
 
         Args:
             target_frequencies (list): Frequencies we are targeting. Closest ones in the data will be used.
+            to_skip (list): List of materials to skip.
         """
         logger.trace("Initializing material collection")
         self.epsilon_matrix = None
@@ -28,7 +29,7 @@ def __init__(self, target_frequencies):
 
        self.target_frequencies = target_frequencies
         self.materials_folder = os.path.dirname(__file__) + "/data/"
-        self._load_materials_folder()
+        self._load_materials_folder(to_skip)
 
     def __getitem__(self, key):
         """Given the name of a material will return the corresponding epsilon tensor.
@@ -42,11 +43,16 @@ def __getitem__(self, key):
         if key in self.material_names:
             return self.epsilon_matrix[self.material_names.index(key)]
         else:
-            raise KeyError(f"Material '{key}' not found in the material collection.")
+            raise KeyError(
+                f"Material '{key}' not found in the material collection. Available are: {self.material_names}"
+            )
 
-    def _load_materials_folder(self):
+    def _load_materials_folder(self, to_skip=None):
         """Loads all csv files from folder "data"
         and sets up the EPS_MATRIX and MATERIAL_NAMES
+
+        Args:
+            to_skip (list): List of materials to skip.
         """
         logger.trace(f"Loading materials from folder: {self.materials_folder}")
 
@@ -57,12 +63,19 @@ def _load_materials_folder(self):
 
         eps_list = []
         for file in files:
+            file = file.replace("\\", "/")
+            name = file.split("/")[-1].split(".csv")[0]
+
+            # Skip materials we don't want
+            if to_skip is not None and name in to_skip:
+                logger.debug(f"Skipping material {name}")
+                continue
+
             # Load material epsilon (permittivity)
             eps_list.append(self._load_material_data(file))
 
             # Remember file name as material name
-            file = file.replace("\\", "/")
-            self.material_names.append(file.split("/")[-1].split(".csv")[0])
+            self.material_names.append(name)
 
         # Create a single tensor of all materials
         logger.trace("Creating material tensor")
@@ -78,7 +91,7 @@ def _load_material_data(self, name):
         Returns:
             torch.tensor: Epsilon for the material (permittivity)
         """
-        logger.trace(f"Loading material {name}")
+        logger.debug(f"Loading material {name}")
         csv_data = pandas.read_csv(name, delimiter="\t")
 
         eps = []
@@ -94,5 +107,5 @@ def _load_material_data(self, name):
             eps.append([real + imag * 1.0j])
 
         eps = torch.tensor(eps)
-        logger.debug(f"Epsilon for material {name} is: {eps}")
+        logger.trace(f"Epsilon for material {name} is: {eps}")
         return eps