Feature/tikhonov inverse vc

examples/example09 - virtual_circuits_MASTU.ipynb

@@ -291,6 +291,10 @@
     "\n",
     "$$ V = (S^T S)^{-1} S^T \\in \\Reals^{N_c \\times N_T}.$$\n",
     "\n",
+    "The matrix inversion above is the \"Moore-Penrose\" pseudo inverse which is used by default. We can also use a \"Tikhonov regularised\" inverse by providing an optional array or diagonal matrix of regularisation coeficients. This is computed as \n",
+    "\n",
+    "$$ V = (S^T S + \\Lambda ^T \\Lambda) ^{-1}S^T.$$\n",
+    "\n",
     "This matrix provides a mapping from requested shifts in the targets to the shifts in the coil currents required."
    ]
   },
@@ -328,6 +332,30 @@
     "    )"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# example of tikhonov inverse \n",
+    "lambda_matrix = 1e-12 * np.eye(len(coils)) #simple example proportional to the identity matrix\n",
+    "\n",
+    "VCs.calculate_VC(\n",
+    "    eq=eq,\n",
+    "    profiles=profiles,\n",
+    "    coils=coils,\n",
+    "    target_names=target_names,\n",
+    "    target_calculator=plasma_descriptors,\n",
+    "    starting_dI=None,    \n",
+    "    min_starting_dI=50,\n",
+    "    tikhonov_lambda= lambda_matrix,\n",
+    "    verbose=True,\n",
+    "    \n",
+    "    name=\"VC_for_lower_targets_with_tikhonov\", # name for the VC\n",
+    "    )"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -651,20 +679,6 @@
     "\n",
     "plt.tight_layout()"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

freegsnke/control_loop/pcs.py

@@ -184,7 +184,7 @@ def calculate_ctrl_voltages(
         emulated_VC_targets_calc=None,
         emulator_coils_calc=None,
         emu_inputs=None,
-        vc_update_rate=None,
+        tikhonov_lambda=None,
         verbose=False,
     ):
         """
@@ -253,6 +253,10 @@ def calculate_ctrl_voltages(
         emulator_coils_calc : list of str, optional
             List of coils to use in emulated VC compuation. These are coils to use in computing shape sensitivity matrix.
 
+        tikhonov_lambda : numpy.ndarray , optional
+            Array of regularisation values for Tikhonov regularisation in emualted VC matrix inversion.
+            Must be same length as emulator_coils_calc.
+
         verbose : bool, optional
             If True, prints diagnostic information from subsystem controllers.
 
@@ -331,6 +335,7 @@ def calculate_ctrl_voltages(
                     emulated_VC_targets_calc=emulated_VC_targets_calc,
                     emulator_coils_calc=emulator_coils_calc,
                     emu_inputs=emu_inputs,
+                    tikhonov_lambda=tikhonov_lambda,
                 )
             )
 

freegsnke/control_loop/virtual_circuits_category.py

@@ -156,6 +156,7 @@ def run_control(
         emulated_VC_targets_calc=None,
         emulator_coils_calc=None,
         emu_inputs=None,
+        tikhonov_lambda=None,
         verbose=False,
     ):
         """
@@ -199,6 +200,10 @@ def run_control(
         emu_inputs : np.ndarray , optional
             Array of input values for all input parameters (currents and other plasma parameters) of the Neural Network emulator.
 
+        tikhonov_lambda : numpy.ndarray , optional
+            Array of regularisation values for Tikhonov regularisation in emualted VC matrix inversion.
+            Must be same length as emulator_coils_calc.
+
         verbose : bool
             Print some output if True.
 
@@ -253,6 +258,7 @@ def run_control(
                     coils=self.ctrl_coils,
                     coils_calc=emulator_coils_calc,
                     input_data=emu_inputs,
+                    tikhonov_lambda=tikhonov_lambda,
                 )
                 # update latest vcs/times
                 self.latest_vc_time = 1.0 * t
@@ -272,15 +278,17 @@ def run_control(
                     coils=self.ctrl_coils,
                     coils_calc=emulator_coils_calc,
                     input_data=emu_inputs,
+                    tikhonov_lambda=tikhonov_lambda,
                 )
 
                 # update latest VCs and times
                 self.latest_vc_time = 1.0 * t
                 self.latest_vc = VC_shape_emu
 
                 # store sensitivity matrix (Jacobian)
-                # self.emulated_jacobian_list.append(self.vc_generator.jacobian_matrix)
-                # self.emulated_vc_list.append(self.vc_generator.vc_matrix)
+                # # self.emulated_jacobian_list.append(self.vc_generator.jacobian_matrix)
+                # # self.emulated_vc_list.append(self.vc_generator.vc_matrix)
+                self.emulated_vc_list.append(self.latest_vc)
                 self.emulated_vc_times.append(t)
 
             else:

freegsnke/virtual_circuits.py

@@ -314,6 +314,71 @@ def build_dIydI_j(
 
         return dtargets / final_dI
 
+    @staticmethod
+    def calculate_matrix_inverse(
+        matrix: np.ndarray, tikhonov_lambda: np.ndarray = None
+    ):
+        """
+        Compute inverse of a generically non-square matrix
+        By default Moore Penrose inverse is used (np.pinv).
+        If Tikhonov_lambda is provided then Tikhonov regularisation is applied
+        and inv = [M^T M + diag(lambda)]^-1 M^T
+
+        Inputs:
+        -------
+        matrix : np.ndarray
+            matrix to be inverted
+
+        tikhonov_lambda : np.ndarray, optional
+            1d array of tikhonov coefficients, or 2d diagonal matrix of coefficients.
+            Must have size/shape consistent with matrix.shape[1].
+
+
+        Returns:
+        --------
+        inverse : np.ndarray
+            inverse of matrix
+        """
+        matrix = np.asarray(matrix)  # convert tensorflow to numpy.
+        if tikhonov_lambda is None:
+            # use regular moore-penrose pseudo inverse
+            print("Computing Moore-Penrose pseudoinverse")
+            inverse = np.linalg.pinv(matrix)
+
+        else:
+            print("Computing Tikhonov regularised inverse")
+            tikhonov_lambda = np.asarray(tikhonov_lambda)
+            n_cols = matrix.shape[1]
+
+            if tikhonov_lambda.ndim == 1:
+                if tikhonov_lambda.shape[0] != n_cols:
+                    raise ValueError(
+                        f"tikhonov_lambda length {tikhonov_lambda.shape[0]} "
+                        f"must match matrix column count {n_cols}."
+                    )
+                tikhonov_matrix = np.diag(tikhonov_lambda)
+
+            elif tikhonov_lambda.ndim == 2:
+                if tikhonov_lambda.shape != (n_cols, n_cols):
+                    raise ValueError(
+                        f"tikhonov_lambda shape {tikhonov_lambda.shape} "
+                        f"must be ({n_cols}, {n_cols}) to match matrix column count."
+                    )
+                if not np.allclose(tikhonov_lambda, np.diag(np.diag(tikhonov_lambda))):
+                    raise ValueError(
+                        "tikhonov_lambda 2d array must be a diagonal matrix."
+                    )
+                tikhonov_matrix = tikhonov_lambda
+
+            else:
+                raise ValueError(
+                    f"tikhonov_lambda must be 1d or 2d, got {tikhonov_lambda.ndim}d."
+                )
+
+            inverse = np.linalg.inv(matrix.T @ matrix + tikhonov_matrix) @ matrix.T
+
+        return inverse
+
     def calculate_VC(
         self,
         eq,
@@ -325,6 +390,7 @@ def calculate_VC(
         starting_dI=None,
         min_starting_dI=50,
         verbose=False,
+        tikhonov_lambda=None,
         name=None,
     ):
         """
@@ -454,15 +520,18 @@ def calculate_VC(
         print("Inverting the shape matrix to get the virtual circuit matrix.")
         print(f"VC object stored under name: '{name}'.")
 
+        # vc matrix is pseudo inverse of shape matrix
+        vc_matrix = self.calculate_matrix_inverse(
+            shape_matrix, tikhonov_lambda=tikhonov_lambda
+        )
+
         # store the VC object dynamically
         store_VC = VirtualCircuit(
             name=name,
             eq=eq,
             profiles=profiles,
             shape_matrix=shape_matrix,
-            VCs_matrix=np.linalg.pinv(
-                shape_matrix
-            ),  # "virtual circuits" are the pseudo-inverse of the shape matrix
+            VCs_matrix=vc_matrix,
             target_names=target_names,
             coils=coils,
             target_calculator=target_calculator,