Merge branch 'merging_carrier_signal_sim' into 'master'

Merging carrier and signal simulation

See merge request !49

playground/micro_resonator.py

@@ -46,25 +46,26 @@ class MicroResonator(Connector):
         self.mech.F_z.frequencies = self.frequencies
 
     def _get_workspace(self, sim):
-        if sim.is_audio:
+        if sim.carrier.any_frequencies_changing:
+            # ddb - This case probably needs some more thought on what beatings change
+            raise NotImplementedError(
+                "Changing carrier frequencies whilst using a MicroResonant not supported yet."
+            )
+
+        if sim.signal:
             refill = (
                 sim.model.fsig.f.is_changing
                 or self.mass.is_changing
                 or self.f.is_changing
                 or self.Q.is_changing
             )
-            ws = MicroResonatorWorkspace(self, sim, refill)
-            ws.set_fill_fn(self.fill)
-            ws.Fz_frequencies = sim.mechanical_frequencies[self.mech.F_z]
-            ws.z_frequencies = sim.mechanical_frequencies[self.mech.z]
+            ws = MicroResonatorWorkspace(self, sim, refill, refill)
+            ws.signal.set_fill_function(self.fill)
+            ws.Fz_frequencies = sim.signal.mechanical_frequencies[self.mech.F_z]
+            ws.z_frequencies = sim.signal.mechanical_frequencies[self.mech.z]
             return ws
-        else:
-            if sim.any_frequencies_changing:
-                # ddb - This case probably needs some more thought on what beatings change
-                raise NotImplementedError(
-                    "Changing carrier frequencies whilst using a MicroResonant not supported yet."
-                )
-            return None
+        
+        return None
 
     def _couples_frequency(self, ws, connection, frequency_in, frequency_out):
         return True
@@ -88,7 +89,7 @@ class MicroResonator(Connector):
         for fi in ws.Fz_frequencies.frequencies:
             for fo in ws.z_frequencies.frequencies:
                 if fi.f == fo.f:
-                    with ws.sim.component_edge_fill3(
+                    with ws.sim.signal.component_edge_fill3(
                         ws.owner_id, ws.connections.F_2_Z_idx, fi.index, fo.index
                     ) as mat:
                         mat[:] = Hz(2 * np.pi * fo.f)

src/finesse/analysis/actions.py

@@ -40,9 +40,10 @@ def get_sweep_array(start: float, stop: float, steps: int, mode="lin"):
 
 
 class ActionWorkspace:
-    def __init__(self, s_prev, model):
+    def __init__(self, s_prev, sim):
         self.s_prev = s_prev
-        self.model = model
+        self.sim = sim
+        self.model = sim.model
         self.fn_do = None
 
 
@@ -89,9 +90,9 @@ class Action:
     def fill_info(self, p_info):
         p_info.add(self.copy_info())
 
-    def setup(self, s_prev, model):
+    def setup(self, s_prev, sim):
         # By default simple actions should just run their 'do' method
-        ws = ActionWorkspace(s_prev, model)
+        ws = ActionWorkspace(s_prev, sim)
         ws.fn_do = self.do
         return ws
 
@@ -112,10 +113,10 @@ class Action:
             p.is_tunable = True
             initial_param_values[p] = p.value
 
-        with model.built():
+        with model.built() as sim:
             try:
                 s = BaseSolution(None, None)
-                ws = self.setup(s, model)
+                ws = self.setup(s, sim)
                 ws.fn_do(ws)
             except StopIteration:
                 raise Exception("Should we reach this point? Probably unexpected")
@@ -185,8 +186,7 @@ class ABCD(Action):
 
 
 class StepParamNDWorkspace(ActionWorkspace):
-    def __init__(self):
-        pass
+    pass
 
 
 class StepParamND(Action):
@@ -227,47 +227,37 @@ class StepParamND(Action):
         if self.on_complete:
             Folder("on_complete", self.on_complete).fill_info(info)
 
-    def setup(self, s_prev, model: finesse.model.Model):
-        ws = StepParamNDWorkspace()
-        ws.s_prev = s_prev
-        ws.model = model
+    def setup(self, s_prev, sim):
+        ws = StepParamNDWorkspace(s_prev, sim)
         ws.info = self.copy_info()
         ws.fn_do = self.do
 
-        ws.params = tuple(get_param(model, p) for p in self._info.parameters_changing)
+        ws.params = tuple(get_param(ws.model, p) for p in self._info.parameters_changing)
         for p in ws.params:
             if not p.is_tunable:
                 raise ParameterLocked(
                     f"{repr(p)} must set as tunable " "before building the simulation"
                 )
-
-        ws.carrier = model.carrier_simulation
-        try:
-            ws.signal = model.signal_simulation
-        except AttributeError:
-            ws.signal = None
-
         return ws
 
     def do(self, ws: StepParamNDWorkspace):
         ws.sol = ArraySolution(
-            self.name, ws.s_prev, ws.model, self.out_shape, self.axes, ws.params
+            self.name, ws.s_prev, ws.sim.detector_workspaces, self.out_shape, self.axes, ws.params
         )
 
         if self.pre_step:
-            ws.pre_step = Folder("pre_step", self.pre_step).setup(ws.sol, ws.model)
+            ws.pre_step = Folder("pre_step", self.pre_step).setup(ws.sol, ws.sim)
         else:
             ws.pre_step = None
 
         if self.post_step:
-            ws.post_step = Folder("post_step", self.post_step).setup(ws.sol, ws.model)
+            ws.post_step = Folder("post_step", self.post_step).setup(ws.sol, ws.sim)
         else:
             ws.post_step = None
 
         # Now we loop over the actual simulation and run each point
         run_axes_scan(
-            ws.carrier,
-            ws.signal,
+            ws.sim,
             self.axes,
             ws.params,
             self.offsets,
@@ -278,7 +268,7 @@ class StepParamND(Action):
         )
 
         if self.on_complete:
-            ws = Folder("on_complete", self.on_complete).setup(ws.sol, ws.model)
+            ws = Folder("on_complete", self.on_complete).setup(ws.sol, ws.sim)
             ws.fn_do(ws)
 
 
@@ -319,29 +309,32 @@ class XNaxis(StepParamND):
             on_complete=on_complete,
         )
 
-    def setup(self, s_prev, model):
+    def setup(self, s_prev, sim):
         # If the model has locks, set them up to happen on the
         # pre-step of StepParamND action
-        if len(model.locks) > 0:
-            self.pre_step = RunLocks(*model.locks)
+        if len(sim.model.locks) > 0:
+            self.pre_step = RunLocks(*sim.model.locks)
 
-        return super().setup(s_prev, model)
+        return super().setup(s_prev, sim)
 
     def do(self, ws: StepParamNDWorkspace):
-        model = ws.model.deepcopy()
+        model = ws.sim.model.deepcopy()
 
         params = (
             get_param(model, pstr) for pstr in ws.info.get_all_parameters_changing()
         )
 
-        for p in params:
-            p.is_tunable = True
+        try:
+            for p in params:
+                p.is_tunable = True
 
-        with model.built():
-            super().do(ws)
+            with model.built():
+                super().do(ws)
+                
+        finally:
+            for p in params:
+                p.is_tunable = False
 
-        for p in params:
-            p.is_tunable = False
 
     def __getattr__(self, key):
         res = re.match("(parameter|mode|start|stop|steps|offset)([0-9]*)", key)
@@ -524,15 +517,15 @@ class Serial(Action):
             if len(p_info.children) == 0:
                 raise Exception("Analysis information object should have been made")
 
-    def setup(self, s_prev, model):
-        ws = ActionWorkspace(s_prev, model)
+    def setup(self, s_prev, sim):
+        ws = ActionWorkspace(s_prev, sim)
         ws.wss = []
         ws.fn_do = self.do
         curr_children = len(ws.s_prev.children)
         # Here we get workspaces for each of the
         # actions we need to run
         for arg in self.args:
-            ws.wss.append(arg.setup(s_prev, model))
+            ws.wss.append(arg.setup(s_prev, sim))
             if len(ws.s_prev.children) > curr_children:
                 # If a solution was made and added as a child in the previous coroutine
                 # then that becomes the next solution in the serial chain
@@ -556,10 +549,10 @@ class Folder(Action):
         self.action.fill_info(info)
         self.folder = None
 
-    def setup(self, s_prev, model):
-        ws = ActionWorkspace(s_prev, model)
+    def setup(self, s_prev, sim):
+        ws = ActionWorkspace(s_prev, sim)
         ws.folder = BaseSolution(self.name, ws.s_prev)
-        ws.action_ws = self.action.setup(ws.folder, model)
+        ws.action_ws = self.action.setup(ws.folder, sim)
         ws.fn_do = self.do
         return ws
 
@@ -704,24 +697,23 @@ class RunLocks(Action):
     def fill_info(self, p_info):
         p_info.add(self.copy_info())
 
-    def setup(self, s_prev, model: finesse.model.Model):
-        ws = RunLocksWorkspace(s_prev, model)
-        ws.locks = tuple(model.elements[l] for l in self.locks)
+    def setup(self, s_prev, sim):
+        ws = RunLocksWorkspace(s_prev, sim)
+        ws.locks = tuple(sim.model.elements[l] for l in self.locks)
         ws.det_ws = [None,] * self.num_locks
 
         for i, l in enumerate(ws.locks):
-            for j, d in enumerate(model.detectors):
+            for j, d in enumerate(sim.model.detectors):
                 if l.error_signal.name == d.name:
-                    ws.det_ws[i] = model._detector_workspaces[j]
+                    ws.det_ws[i] = sim.detector_workspaces[j]
 
         if any(_ is None for _ in ws.det_ws):
             raise Exception("Could not find detector workspaces for all locks")
 
         ws.info = self.copy_info()
         ws.s_prev = s_prev
-        ws.sim = model.carrier_simulation
         ws.fn_do = do_lock
-        ws.params = tuple(get_param(model, p) for p in self._info.parameters_changing)
+        ws.params = tuple(get_param(sim.model, p) for p in self._info.parameters_changing)
         ws.max_iterations = self.max_iterations
 
         for p in ws.params:
@@ -737,7 +729,7 @@ def do_lock(ws: RunLocksWorkspace):
     iters = 0
     while recompute and iters < ws.max_iterations:
         iters += 1
-        ws.sim.run(True)
+        ws.sim.run_carrier()
         recompute = False
         for i, dws in enumerate(ws.det_ws):
             acc = ws.locks[i].accuracy

src/finesse/analysis/runners.pyx

@@ -8,7 +8,7 @@ from cpython.ref cimport PyObject
 from finesse.cymath cimport complex_t
 
 from finesse.parameter cimport Parameter
-from finesse.simulations cimport BaseSimulation
+from finesse.simulations.basematrix cimport CarrierSignalMatrixSimulation
 from finesse.solutions.array import ArraySolution
 from finesse.solutions.array cimport ArraySolution
 
@@ -18,8 +18,7 @@ LOGGER = logging.getLogger(__name__)
 @cython.wraparound(False)
 @cython.initializedcheck(False)
 cpdef run_fsig_sweep(
-        BaseSimulation carrier,
-        BaseSimulation signal,
+        CarrierSignalMatrixSimulation sim,
         double[::1] axis,
         long[::1] input_rhs_indices,
         long[::1] output_rhs_indices,
@@ -39,7 +38,7 @@ cpdef run_fsig_sweep(
         int No = len(output_rhs_indices)
         int i, o, j
         complex_t denom
-        Parameter f = signal.model.fsig.f
+        Parameter f = sim.model.fsig.f
 
     if out is None:
         out = np.zeros((Na, Ni, No), dtype=np.complex_128)
@@ -48,36 +47,39 @@ cpdef run_fsig_sweep(
         assert(out.shape[1] == Na)
         assert(out.shape[2] == No)
 
+    # We'll be making our own RHS inputs for this simulation
+    sim.signal.manual_rhs = True
+
     for i in range(Ni):
         for j in range(Na):
             f.set_double_value(axis[j])
-            signal.model_data.fsig = axis[j]
-            signal._clear_rhs()
-            signal.set_source_fast_2(
+            sim.signal.model_data.fsig = axis[j]
+            sim.signal.clear_rhs()
+            sim.signal.set_source_fast_2(
                 input_rhs_indices[i], 1
             )
-            signal.run(False)
+            sim.signal.run()
             if not compute_open_loop:
                 for o in range(No):
-                    out[i][j][o] = signal.out_view[output_rhs_indices[o]]
+                    out[i][j][o] = sim.signal.out_view[output_rhs_indices[o]]
             else:
                 for o in range(No):
                     # We can divide out the 1/(1-H) closed loop behaviours by
                     # using the coupling computed back into the same input node
-                    denom = signal.out_view[input_rhs_indices[i]]
+                    denom = sim.signal.out_view[input_rhs_indices[i]]
                     if denom.real == denom.imag == 0:
                         out[i][j][o] = 0
                     else:
-                        out[i][j][o] = signal.out_view[output_rhs_indices[o]] / denom
+                        out[i][j][o] = sim.signal.out_view[output_rhs_indices[o]] / denom
 
+    sim.signal.manual_rhs = False
     return out
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.initializedcheck(False)
 def run_axes_scan(
-        BaseSimulation carrier,
-        BaseSimulation signal,
+        CarrierSignalMatrixSimulation sim,
         tuple axes,
         tuple params,
         double[:] offsets,
@@ -151,11 +153,10 @@ def run_axes_scan(
             # ------------------------------------------------------
             # DO STEP
             # ------------------------------------------------------
-            mask_this_point = carrier.run(True)
+            mask_this_point = sim.run_carrier()
 
-            if not mask_this_point and signal is not None:
-                signal.model_data.fsig = signal.model.fsig.f.value
-                mask_this_point = signal.run(True)
+            if not mask_this_point and sim.signal:
+                sim.run_signal()
             # ------------------------------------------------------
             # POST STEP
             # ------------------------------------------------------

src/finesse/components/beamsplitter.py

@@ -586,16 +586,17 @@ class Beamsplitter(Surface):
 
     def _get_workspace(self, sim):
         from finesse.components.modal.beamsplitter import (
-            beamsplitter_fill,
+            beamsplitter_carrier_fill,
+            beamsplitter_signal_fill,
             BeamsplitterWorkspace,
         )
 
         _, is_changing = self._eval_parameters()
 
         refill = (
-            (sim.is_audio and sim.model.fsig.f.is_changing)
+            (sim.signal and sim.model.fsig.f.is_changing)
             or self in sim.trace_forest.changing_components
-            or sim.any_frequencies_changing
+            or sim.carrier.any_frequencies_changing
             or (len(is_changing) and is_changing.issubset(self.__changing_check))
         )
 
@@ -604,7 +605,8 @@ class Beamsplitter(Surface):
         ws.nr1 = refractive_index(self.p1) or refractive_index(self.p2) or 1
         ws.nr2 = refractive_index(self.p3) or refractive_index(self.p4) or 1
 
-        ws.set_fill_fn(beamsplitter_fill)
+        ws.carrier.set_fill_function(beamsplitter_carrier_fill)
+        ws.signal.set_fill_function(beamsplitter_signal_fill)
 
         # Initialise the ABCD matrix memory-views
         if sim.is_modal:
@@ -650,24 +652,4 @@ class Beamsplitter(Surface):
             ws.abcd_p4p2_x = self.ABCD(self.p4.i, self.p2.o, "x", copy=False)
             ws.abcd_p4p2_y = self.ABCD(self.p4.i, self.p2.o, "y", copy=False)
 
-        return ws
-
-    def _fill_qnoise_rhs(self, sim):
-        ws = sim.ws[self]
-        if sim.is_modal:
-            for freq in sim.frequencies:
-                for k in range(sim.nhoms):
-                    L1 = self.R * ws.K21_loss[k] + self.T * ws.K31_loss[k]
-                    L2 = self.R * ws.K12_loss[k] + self.T * ws.K42_loss[k]
-                    L3 = self.R * ws.K43_loss[k] + self.T * ws.K13_loss[k]
-                    L4 = self.R * ws.K34_loss[k] + self.T * ws.K24_loss[k]
-                    # TODO, can't access matrix like this!
-                    sim.Mq[sim.field(self.p1.o, freq.index, k)] = (L1 + self.L) / 2
-                    sim.Mq[sim.field(self.p2.o, freq.index, k)] = (L2 + self.L) / 2
-                    sim.Mq[sim.field(self.p3.o, freq.index, k)] = (L3 + self.L) / 2
-                    sim.Mq[sim.field(self.p4.o, freq.index, k)] = (L4 + self.L) / 2
-        else:
-            for freq in sim.frequencies:
-                # TODO, can't access matrix like this!
-                sim.Mq[sim.field(self.p1.o, freq.index, 0)] = self.L / 2
-                sim.Mq[sim.field(self.p2.o, freq.index, 0)] = self.L / 2
+        return ws
\ No newline at end of file

src/finesse/components/directional_beamsplitter.py

@@ -75,8 +75,8 @@ class DirectionalBeamsplitter(Connector):
 
     def _get_workspace(self, sim):
         ws = DBSWorkspace(self, sim, False)
-        ws.I = np.eye(sim.nhoms, dtype=np.complex128)
-        ws.set_fill_fn(self._fill_matrix)
+        ws.I = np.eye(sim.model_data.num_HOMs, dtype=np.complex128)
+        ws.set_fill_function(self._fill_matrix)
         return ws
 
     def _fill_matrix(self, ws):

src/finesse/components/dof.py

@@ -72,9 +72,9 @@ class DegreeOfFreedom(Connector):
         return self.__amplitudes
 
     def _get_workspace(self, sim):
-        if sim.is_audio:
+        if sim.signal:
             ws = DOFWorkspace(self, sim, False)
-            ws.set_fill_fn(self.__fill)
+            ws.set_fill_function(self.__fill)
             ws.drives = self.drives
             ws.amplitudes = self.amplitudes
             return ws

src/finesse/components/electronics.py

@@ -34,12 +34,12 @@ class ZPKNodeActuator(Connector):
         self._register_node_coupling("P1_ACT", self.p1.i, mechanical_node)
 
     def _get_workspace(self, sim):
-        if sim.is_audio:
+        if sim.signal:
             refill = sim.model.fsig.f.is_changing or any(
                 p.is_changing for p in self.parameters
             )
-            ws = FilterWorkspace(self, sim, refill)
-            ws.set_fill_fn(self.fill)
+            ws = FilterWorkspace(self, sim, refill, refill)
+            ws.signal.set_fill_function(self.fill)
             ws.frequencies = sim.electrical_frequencies[self.p1.i].frequencies
             return ws
         else:
@@ -48,7 +48,7 @@ class ZPKNodeActuator(Connector):
     def fill(self, ws):
         for _ in ws.frequencies:
             # Assumes only couples to first mech frequency
-            with ws.sim.component_edge_fill3(
+            with ws.sim.signal.component_edge_fill3(
                 ws.owner_id, ws.connections.P1_ACT_idx, _.index, 0,
             ) as mat:
                 mat[:] = ws.values.gain
@@ -69,12 +69,12 @@ class Amplifier(Connector):
         self._register_node_coupling("P1_P2", self.p1.i, self.p2.o)
 
     def _get_workspace(self, sim):
-        if sim.is_audio:
+        if sim.signal:
             refill = sim.model.fsig.f.is_changing or any(
                 p.is_changing for p in self.parameters
             )
-            ws = FilterWorkspace(self, sim, refill)
-            ws.set_fill_fn(self.fill)
+            ws = FilterWorkspace(self, sim, refill, refill)
+            ws.signal.set_fill_function(self.fill)
             ws.frequencies = sim.electrical_frequencies[self.p1.i].frequencies
             return ws
         else:
@@ -82,7 +82,7 @@ class Amplifier(Connector):
 
     def fill(self, ws):
         for _ in ws.frequencies:
-            with ws.sim.component_edge_fill3(
+            with ws.sim.signal.component_edge_fill3(
                 ws.owner_id, ws.connections.P1_P2_idx, 0, 0,
             ) as mat:
                 mat[:] = ws.values.gain
@@ -106,13 +106,13 @@ class Filter(Connector):
         self._register_node_coupling("P1_P2", self.p1.i, self.p2.o)
 
     def _get_workspace(self, sim):
-        if sim.is_audio:
+        if sim.signal:
             refill = sim.model.fsig.f.is_changing or any(
                 p.is_changing for p in self.parameters
             )
-            ws = FilterWorkspace(self, sim, refill)
-            ws.set_fill_fn(self.fill)
-            ws.frequencies = sim.electrical_frequencies[self.p1.i].frequencies
+            ws = FilterWorkspace(self, sim, refill, refill)
+            ws.signal.set_fill_function(self.fill)
+            ws.frequencies = sim.signal.electrical_frequencies[self.p1.i].frequencies
             return ws
         else:
             return None
@@ -121,7 +121,7 @@ class Filter(Connector):
         Hz = self.eval(ws.sim.model_data.fsig)
 
         for _ in ws.frequencies:
-            with ws.sim.component_edge_fill3(
+            with ws.sim.signal.component_edge_fill3(
                 ws.owner_id, ws.connections.P1_P2_idx, 0, 0,
             ) as mat:
                 mat[:] = Hz

src/finesse/components/isolator.py

@@ -66,7 +66,7 @@ class Isolator(Connector):
         ws.nr1 = refractive_index(self.p1) or 1
         ws.nr2 = refractive_index(self.p2) or 1
 
-        ws.set_fill_fn(self._fill_matrix)
+        ws.set_fill_function(self._fill_matrix)
         return ws
 
     def _fill_matrix(self, ws):

src/finesse/components/laser.py

@@ -9,11 +9,6 @@ from finesse.cymath.complex import crotate
 from finesse.parameter import model_parameter, Rebuild
 
 from finesse.components.general import Connector, FrequencyGenerator
-from finesse.components.modal.laser import (
-    laser_fill_rhs,
-    laser_fill_signal,
-    LaserWorkspace,
-)
 from finesse.components.node import NodeType, NodeDirection
 
 
@@ -84,7 +79,7 @@ class Laser(Connector, FrequencyGenerator):
     def __find_src_freq(self, sim):
         # if it's tunable we want to look for the symbol that is just this
         # lasers frequency, as it will be changing
-        for f in sim.frequencies:
+        for f in sim.optical_frequencies.frequencies:
             if not self.f.is_changing:
                 # Don't match changing frequency bins if ours won't match
                 if not f.symbol.is_changing and (
@@ -101,48 +96,47 @@ class Laser(Connector, FrequencyGenerator):
         return None
 
     def _get_workspace(self, sim):
+        from finesse.components.modal.laser import (
+            laser_carrier_fill_rhs,
+            laser_fill_signal,
+            LaserWorkspace,
+            laser_set_gouy
+        )
         ws = LaserWorkspace(self, sim, True)
-        ws.node_id = sim.node_id(self.p1.o)
+        ws.node_car_id = sim.carrier.node_id(self.p1.o)
         ws.fsrc_car_idx = -1
-
-        if not sim.is_audio: