rewrites code generation to use gathered params, static sparse jacobians

pyCamSet/optimisation/optimisation_handling.py

@@ -94,12 +94,10 @@ def run_bundle_adjustment(param_handler: TemplateBundleHandler,
         jac= bundle_jac if bundle_jac is not None else "2-point", #pass the function for the jacobian if it exists
         max_nfev=param_handler.problem_opts["max_nfev"],
         # loss = "cauchy"
-        # x_scale='jac',
+        x_scale='jac',
     )
     end = time.time()
 
-
-
     final_euclid = np.mean(np.linalg.norm(np.reshape(optimisation.fun, (-1, 2)), axis=1))
     logging.info(f'Final Euclidean error: {final_euclid:.2f} px')
     logging.info(f'Optimisation took {end - start: .2f} seconds.')

pyCamSet/optimisation/standard_bundle_handler.py

@@ -8,6 +8,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.optimize import least_squares
+from scipy.sparse import csr_array
 import pyvista as pv
 from itertools import combinations
 
@@ -167,6 +168,13 @@ def __init__(self,
         self.op_fun: fb.optimisation_function = fb.projection() + fb.extrinsic3D() + fb.rigidTform3d() +  fb.free_point()
 
     def make_loss_fun(self, threads):
+        """
+        Describes and writes the loss function of the loss function represented by self.
+        Wraps the loss function to account for the fixed parameters of the optimisation.
+
+        :params threads: the number of threads to use.
+
+        """
         target_shape = self.target.point_data.shape
         dd = self.detection.return_flattened_keys(target_shape[:-1]).get_data()
         temp_loss = self.op_fun.make_full_loss_fn(dd, threads)
@@ -177,9 +185,15 @@ def loss_fun(params):
         return loss_fun
 
     def make_loss_jac(self, threads): 
+        """
+        Describes and writes the jacobian of the loss function described in self.
+        Wraps the jacobian of the loss function to deal with the fixed parameters of the optimisation.
+
+        :params threads: the number of threads to use for the optimisation.
+        :returns jac_fn: a callable jacobian function that returns the jacobian of the given paramaters.
+        """
         target_shape = self.target.point_data.shape
         dd = self.detection.return_flattened_keys(target_shape[:-1]).get_data()
-        temp_loss = self.op_fun.make_jacobean(dd, threads)
         mask = np.concatenate(
             ( 
                 np.repeat(self.intr_unfixed, 9),
@@ -189,18 +203,13 @@ def make_loss_jac(self, threads):
             ), axis=0
         )
 
-        if np.all(mask):
-            def jac_fn(params):
-                inps = self.get_bundle_adjustment_inputs(params) #return proj, extr, poses
-                param_str = self.op_fun.build_param_list(*inps)
-                return temp_loss(param_str, )
-            return jac_fn
-        else:
-            def jac_fn(params):
-                inps = self.get_bundle_adjustment_inputs(params) #return proj, extr, poses
-                param_str = self.op_fun.build_param_list(*inps)
-                return temp_loss(param_str)[:, mask]
-            return jac_fn
+        temp_loss = self.op_fun.make_jacobean(dd, threads, unfixed_params=mask)
+        def jac_fn(params):
+            inps = self.get_bundle_adjustment_inputs(params) #return proj, extr, poses
+            param_str = self.op_fun.build_param_list(*inps)
+            d, c, rp = temp_loss(param_str)
+            return csr_array((d,c,rp), shape=(2*dd.shape[0], params.shape[0]))
+        return jac_fn
 
     def get_bundle_adjustment_inputs(self, x) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
         """

pyCamSet/optimisation/template_handler.py

@@ -7,6 +7,7 @@
 import numpy as np
 
 from typing import TYPE_CHECKING
+from scipy.sparse import csr_array
 
 import pyCamSet.utils.general_utils as gu
 import pyCamSet.optimisation.compiled_helpers as ch
@@ -170,31 +171,23 @@ def loss_fun(params):
     def make_loss_jac(self, threads): 
         #TODO implement proper culling
         obj_data = self.target.point_data.reshape((-1, 3))
-
         target_shape = self.target.point_data.shape
         dd = self.detection.return_flattened_keys(target_shape[:-1]).get_data()
-        temp_loss = self.op_fun.make_jacobean(dd, threads)
         mask = np.concatenate(
             ( 
                 np.repeat(self.intr_unfixed, 9),
                 np.repeat(self.extr_unfixed, 6),
                 np.repeat(self.pose_unfixed, 6),
             ), axis=0
         )
-        if np.all(mask):
-
-            def jac_fn(params):
-                inps = self.get_bundle_adjustment_inputs(params) #return proj, extr, poses
-                param_str = self.op_fun.build_param_list(*inps)
-                return temp_loss(param_str, obj_data)
-            return jac_fn
-        else:
-
-            def jac_fn(params):
-                inps = self.get_bundle_adjustment_inputs(params) #return proj, extr, poses
-                param_str = self.op_fun.build_param_list(*inps)
-                return temp_loss(param_str, obj_data)[:, mask]
-            return jac_fn
+
+        temp_loss = self.op_fun.make_jacobean(dd, threads, unfixed_params=mask)
+        def jac_fn(params):
+            inps = self.get_bundle_adjustment_inputs(params) #return proj, extr, poses
+            param_str = self.op_fun.build_param_list(*inps)
+            d, c, rp = temp_loss(param_str, obj_data)
+            return csr_array((d,c,rp), shape=(2*dd.shape[0], params.shape[0]))
+        return jac_fn
 
     def special_plots(self, params):
         """

pyCamSet/utils/general_utils.py

@@ -88,9 +88,9 @@ def run_benchmark():
         median = np.median(times) * ranges[mode]
         max_t = min(mean + 3*stdev,np.amax(times) * ranges[mode])
         print(f"Mean: {mean:.2f} {mode}, median: {median:.2f} {mode}, stdev: {stdev:.2f} {mode}")
-        histogram(times*ranges[mode], bins=50,
+        histogram(times*ranges[mode], bins=20,
                   bins_min=max(mean- 3*stdev, 0),
-                  x_max= max_t,
+                  x_max = min(mean + 5*stdev, max_t),
                   height = 3,
                   color = True,
                   y_unit=" freq",

setup_scripts/graph_runtime_calib.py

@@ -0,0 +1,40 @@
+from pathlib import Path
+from cv2 import aruco
+from multiprocessing import cpu_count
+
+from pyCamSet import Ccube, load_CameraSet
+from pyCamSet.utils.general_utils import benchmark
+
+from pyCamSet.optimisation.template_handler import TemplateBundleHandler
+from pyCamSet.optimisation.standard_bundle_handler import SelfBundleHandler
+from pyCamSet.optimisation.optimisation_handling import run_bundle_adjustment, make_optimisation_function
+
+# load scripts that represent all of the calibration examples
+# then, drag race the functions!
+
+target = Ccube(n_points=10, length=40,  aruco_dict=aruco.DICT_6X6_1000, border_fraction=0.2)
+loc=Path("tests/test_data/calibration_ccube")
+test_cams = loc/'self_calib_test.camset'
+debug = True
+
+cams = load_CameraSet(test_cams)
+
+param_handler = SelfBundleHandler(
+    detection=cams.calibration_handler.detection, target=target, camset=cams,
+    options={'max_nfev':100}
+)
+
+param_handler.set_from_templated_camset(cams)
+
+loss, jac, init = make_optimisation_function(
+    param_handler=param_handler,
+    threads = cpu_count(),
+)
+print("initialising")
+loss(init)
+jac(init)
+print("Testing self calib loss")
+benchmark(lambda :loss(init), repeats=100, mode='ms')
+
+print("Testing self calib jac")
+benchmark(lambda :jac(init), repeats=100, mode='ms')