Laplacian quality using determinant and more luigi parameters (#18)

.pylintrc

@@ -25,7 +25,7 @@ max-args=10
 # with leading underscore
 ignored-argument-names=_.*
 # Maximum number of locals for function / method body
-max-locals=15
+max-locals=20
 # Maximum number of return / yield for function / method body
 max-returns=6
 # Maximum number of branch for function / method body

netsalt/algorithm.py

@@ -31,6 +31,7 @@ def refine_mode_brownian_ratchet(
     disp=False,
     save_mode_trajectories=False,
     seed=42,
+    quality_method="eigenvalue",
 ):
     """Accurately find a mode from an initial guess, using brownian ratchet algorithm.
 
@@ -54,7 +55,7 @@ def refine_mode_brownian_ratchet(
     if save_mode_trajectories:
         mode_trajectories = [current_mode.copy()]
 
-    initial_quality = mode_quality(current_mode, graph)
+    initial_quality = mode_quality(current_mode, graph, quality_method=quality_method)
     current_quality = initial_quality
 
     search_stepsize = params["search_stepsize"]
@@ -66,7 +67,7 @@ def refine_mode_brownian_ratchet(
             + search_stepsize * current_quality / initial_quality * np.random.uniform(-1, 1, 2)
         )
 
-        new_quality = mode_quality(new_mode, graph)
+        new_quality = mode_quality(new_mode, graph, quality_method=quality_method)
 
         if disp:
             L.debug(
@@ -114,6 +115,7 @@ def refine_mode_brownian_ratchet(
         params,
         disp=disp,
         save_mode_trajectories=save_mode_trajectories,
+        quality_method=quality_method,
     )
 
 

netsalt/modes.py

@@ -34,14 +34,23 @@
 class WorkerModes:
     """Worker to find modes."""
 
-    def __init__(self, estimated_modes, graph, D0s=None, search_radii=None, seed=42):
+    def __init__(
+        self,
+        estimated_modes,
+        graph,
+        D0s=None,
+        search_radii=None,
+        seed=42,
+        quality_method="eigenvalue",
+    ):
         """Init function of the worker."""
         self.graph = graph
         self.params = graph.graph["params"]
         self.estimated_modes = estimated_modes
         self.D0s = D0s
         self.search_radii = search_radii
         self.seed = seed
+        self.quality_method = quality_method
 
     def set_search_radii(self, mode):
         """This fixes a local search region set by search radii."""
@@ -59,26 +68,29 @@ def __call__(self, mode_id):
         mode = self.estimated_modes[mode_id]
         if self.search_radii is not None:
             self.set_search_radii(mode)
-        return refine_mode_brownian_ratchet(mode, self.graph, self.params, seed=self.seed)
+        return refine_mode_brownian_ratchet(
+            mode, self.graph, self.params, seed=self.seed, quality_method=self.quality_method
+        )
 
 
 class WorkerScan:
     """Worker to scan complex frequency."""
 
-    def __init__(self, graph):
+    def __init__(self, graph, quality_method="eigenvalue"):
         self.graph = graph
+        self.quality_method = quality_method
         np.random.seed(42)
 
     def __call__(self, freq):
-        return mode_quality(to_complex(freq), self.graph)
+        return mode_quality(to_complex(freq), self.graph, quality_method=self.quality_method)
 
 
-def scan_frequencies(graph):
+def scan_frequencies(graph, quality_method="eigenvalue"):
     """Scan a range of complex frequencies and return mode qualities."""
     ks, alphas = get_scan_grid(graph)
     freqs = [[k, a] for k in ks for a in alphas]
 
-    worker_scan = WorkerScan(graph)
+    worker_scan = WorkerScan(graph, quality_method=quality_method)
     chunksize = max(1, int(0.1 * len(freqs) / graph.graph["params"]["n_workers"]))
     with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool:
         qualities_list = list(
@@ -104,15 +116,17 @@ def _init_dataframe():
     return pd.DataFrame(columns=indexes)
 
 
-def find_modes(graph, qualities):
+def find_modes(graph, qualities, quality_method="eigenvalue", min_distance=2, threshold_abs=1.0):
     """Find the modes from a scan."""
     ks, alphas = get_scan_grid(graph)
     estimated_modes = find_rough_modes_from_scan(
-        ks, alphas, qualities, min_distance=2, threshold_abs=1.0
+        ks, alphas, qualities, min_distance=min_distance, threshold_abs=threshold_abs
     )
     L.info("Found %s mode candidates.", len(estimated_modes))
     search_radii = [1 * (ks[1] - ks[0]), 1 * (alphas[1] - alphas[0])]
-    worker_modes = WorkerModes(estimated_modes, graph, search_radii=search_radii)
+    worker_modes = WorkerModes(
+        estimated_modes, graph, search_radii=search_radii, quality_method=quality_method
+    )
     with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool:
         refined_modes = list(
             tqdm(
@@ -714,7 +728,7 @@ def compute_modal_intensities(modes_df, max_pump_intensity, mode_competition_mat
     return modes_df
 
 
-def pump_trajectories(modes_df, graph, return_approx=False):
+def pump_trajectories(modes_df, graph, return_approx=False, quality_method="eigenvalue"):
     """For a sequence of D0s, find the mode positions of the modes modes."""
 
     D0s = np.linspace(
@@ -738,7 +752,12 @@ def pump_trajectories(modes_df, graph, return_approx=False):
         for m in range(n_modes):
             pumped_modes_approx[-1][m] = pump_linear(pumped_modes[-1][m], graph, D0s[d], D0s[d + 1])
 
-        worker_modes = WorkerModes(pumped_modes_approx[-1], graph, D0s=n_modes * [D0s[d + 1]])
+        worker_modes = WorkerModes(
+            pumped_modes_approx[-1],
+            graph,
+            D0s=n_modes * [D0s[d + 1]],
+            quality_method=quality_method,
+        )
         with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool:
             pumped_modes.append(list(tqdm(pool.imap(worker_modes, range(n_modes)), total=n_modes)))
         for i, mode in enumerate(pumped_modes[-1]):
@@ -779,7 +798,8 @@ def _get_new_D0(arg, graph=None, D0_steps=0.1):
     return mode_id, new_D0, new_modes_approx
 
 
-def find_threshold_lasing_modes(modes_df, graph):  # pylint:disable=too-many-statements
+def find_threshold_lasing_modes(modes_df, graph, quality_method="eigenvalue"):
+    # pylint:disable=too-many-statements
     """Find the threshold lasing modes and associated lasing thresholds."""
     stepsize = graph.graph["params"]["search_stepsize"]
     D0_steps = graph.graph["params"]["D0_max"] / graph.graph["params"]["D0_steps"]
@@ -818,7 +838,9 @@ def find_threshold_lasing_modes(modes_df, graph):  # pylint:disable=too-many-sta
         )
 
         L.debug("Current search_stepsize: %s", graph.graph["params"]["search_stepsize"])
-        worker_modes = WorkerModes(new_modes_approx, graph, D0s=new_D0s)
+        worker_modes = WorkerModes(
+            new_modes_approx, graph, D0s=new_D0s, quality_method=quality_method
+        )
         new_modes_tmp = np.zeros([len(modes_df), 2])
 
         with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool:

netsalt/quantum_graph.py

@@ -396,7 +396,7 @@ def laplacian_quality(laplacian, method="eigenvalue"):
 
     Args:
         laplacian (sparse matrix): laplacian matrix
-        method (str): either eigenvalue or singular value
+        method (str): method for quality evaluation (eigenvalue, singular value or determinant)
     """
     v0 = np.random.random(laplacian.shape[0])
     if method == "eigenvalue":
@@ -422,6 +422,10 @@ def laplacian_quality(laplacian, method="eigenvalue"):
                 )
             )[0]
 
+    if method == "determinant":
+        logdet = np.linalg.slogdet(laplacian.todense())[1]
+        return np.exp(np.real(logdet / laplacian.shape[0]))
+
     if method == "singularvalue":
         return sc.sparse.linalg.svds(
             laplacian,
@@ -433,12 +437,13 @@ def laplacian_quality(laplacian, method="eigenvalue"):
     return 1.0
 
 
-def mode_quality(mode, graph):
+def mode_quality(mode, graph, quality_method='eigenvalue'):
     """Quality of a mode, small means good quality, thus the mode is close to a correct mode.
 
     Args:
         mode (complex): complex mode
         graph (graph): quantum graph
+        quality_method (str): method for quality evaluation (eig, singular value or det)
     """
     laplacian = construct_laplacian(to_complex(mode), graph)
-    return laplacian_quality(laplacian)
+    return laplacian_quality(laplacian, method=quality_method)

netsalt/tasks/config.py

@@ -19,6 +19,11 @@ class ModeSearchConfig(luigi.Config):
     max_tries_reduction = luigi.IntParameter(default=50)
     reduction_factor = luigi.FloatParameter(default=1.0)
     search_stepsize = luigi.FloatParameter(default=0.001)
+    quality_method = luigi.ChoiceParameter(
+        default="eigenvalue", choices=["eigenvalue", "singularvalue", "determinant"]
+    )
+    threshold_abs = luigi.FloatParameter(default=0.1)
+    min_distance = luigi.FloatParameter(default=2)
 
 
 class PumpConfig(luigi.Config):

netsalt/tasks/lasing.py

@@ -13,6 +13,7 @@
     pump_trajectories,
 )
 from netsalt.pump import make_threshold_pump
+from .config import ModeSearchConfig
 from .netsalt_task import NetSaltTask
 from .passive import CreateQuantumGraph, FindPassiveModes
 from .pump import OptimizePump
@@ -94,7 +95,9 @@ def run(self):
         modes_df = load_modes(self.input()["modes"].path)
         qg = self.get_graph_with_pump(self.input()["graph"].path)
         if not self.skip:
-            modes_df = pump_trajectories(modes_df, qg, return_approx=True)
+            modes_df = pump_trajectories(
+                modes_df, qg, return_approx=True, quality_method=ModeSearchConfig().quality_method
+            )
         save_modes(modes_df, filename=self.output().path)
 
     def output(self):
@@ -120,7 +123,9 @@ def run(self):
         """ """
         qg = self.get_graph_with_pump(self.input()["graph"].path)
         modes_df = load_modes(self.input()["modes"].path)
-        modes_df = find_threshold_lasing_modes(modes_df, qg)
+        modes_df = find_threshold_lasing_modes(
+            modes_df, qg, quality_method=ModeSearchConfig().quality_method
+        )
         save_modes(modes_df, filename=self.output().path)
 
     def output(self):

netsalt/tasks/passive.py

@@ -101,7 +101,7 @@ def requires(self):
     def run(self):
         """ """
         qg = self.get_graph(self.input().path)
-        qualities = scan_frequencies(qg)
+        qualities = scan_frequencies(qg, quality_method=ModeSearchConfig().quality_method)
         save_qualities(qualities, filename=self.output().path)
 
     def output(self):
@@ -122,7 +122,13 @@ def run(self):
         """ """
         qg = self.get_graph(self.input()["graph"].path)
         qualities = load_qualities(filename=self.input()["qualities"].path)
-        modes_df = find_modes(qg, qualities)
+        modes_df = find_modes(
+            qg,
+            qualities,
+            quality_method=ModeSearchConfig().quality_method,
+            min_distance=ModeSearchConfig().min_distance,
+            threshold_abs=ModeSearchConfig().threshold_abs,
+        )
         save_modes(modes_df, filename=self.output().path)
 
     def output(self):