Merge pull request #515 from sandialabs/simplify-TimeIndependentMDCObjectiveFunction

Simplify the TimeIndependentMDCObjectiveFunction class

Author: sserita

pygsti/objectivefns/objectivefns.py

@@ -1124,13 +1124,22 @@ def dclip(ar): return ar
                                 printer.log("      Accepted%s! gain ratio=%g  mu * %g => %g"
                                             % (" UPHILL" if uphill_ok else "", dF / dL, mu_factor, mu), 2)
                                 last_accepted_dx = dx.copy()
-                                if new_x_is_known_inbounds and norm_f < min_norm_f:
-                                    min_norm_f = norm_f
-                                    best_x[:] = x[:]
-                                    best_x_state = (mu, nu, norm_f, f.copy(), spow, None)
-                                    #Note: we use rawJTJ=None above because the current `JTJ` was evaluated
-                                    # at the *last* x-value -- we need to wait for the next outer loop
-                                    # to compute the JTJ for this best_x_state
+                                if norm_f < min_norm_f:
+                                    if not new_x_is_known_inbounds:
+                                        try:
+                                            _ = obj_fn(global_x, oob_check=True)
+                                            # ^ Dead-store the return value.
+                                            new_x_is_known_inbounds = True
+                                        except ValueError:
+                                            # Then we keep new_x_is_known_inbounds==False.
+                                            pass
+                                    if new_x_is_known_inbounds:
+                                        min_norm_f = norm_f
+                                        best_x[:] = x[:]
+                                        best_x_state = (mu, nu, norm_f, f.copy(), spow, None)
+                                        # ^ Note: we use rawJTJ=None above because the current `JTJ` was evaluated
+                                        #   at the *last* x-value -- we need to wait for the next outer loop
+                                        #   to compute the JTJ for this best_x_state
 
                                 #assert(_np.isfinite(x).all()), "Non-finite x!" # NaNs tracking
                                 #assert(_np.isfinite(f).all()), "Non-finite f!" # NaNs tracking

pygsti/optimize/customlm.py

@@ -600,7 +600,6 @@ def simplish_leastsq(
 
             #determing increment using adaptive damping
             while True:  # inner loop
-
                 if profiler: profiler.memory_check("simplish_leastsq: begin inner iter")
 
                 # ok if assume fine-param-proc.size == 1 (otherwise need to sync setting local JTJ)
@@ -783,17 +782,30 @@ def simplish_leastsq(
                 norm_f = norm_new_f
                 global_x[:] = global_new_x[:]
                 printer.log("      Accepted%s! gain ratio=%g  mu * %g => %g" % ("", dF / dL, mu_factor, mu), 2)
-                if new_x_is_known_inbounds and norm_f < min_norm_f:
-                    min_norm_f = norm_f
-                    best_x[:] = x[:]
-                    best_x_state = (mu, nu, norm_f, f.copy())
+                if norm_f < min_norm_f:
+                    if not new_x_is_known_inbounds:
+                        try:
+                            _ = obj_fn(global_x, oob_check=True)
+                            # ^ Dead-store the return value.
+                            new_x_is_known_inbounds = True
+                        except ValueError:
+                            # Then we keep new_x_is_known_inbounds==False.
+                            pass
+                    if new_x_is_known_inbounds:
+                        min_norm_f = norm_f
+                        best_x[:] = x[:]
+                        best_x_state = (mu, nu, norm_f, f.copy())
 
                 #assert(_np.isfinite(x).all()), "Non-finite x!" # NaNs tracking
                 #assert(_np.isfinite(f).all()), "Non-finite f!" # NaNs tracking
 
                 break 
                 # ^ exit inner loop normally ...
             # end of inner loop
+            #
+            # x[:] = best_x[:]
+            # mu, nu, norm_f, f[:] = best_x_state
+            #
         # end of outer loop
         else:
             #if no break stmt hit, then we've exceeded max_iter

pygsti/optimize/simplerlm.py

@@ -29,4 +29,4 @@ export MKL_NUM_THREADS=1
 # Note: This flags are useful on Kahuna to avoid error messages
 # But the --mca flags are not necessary for performance
 mpirun -np ${NUM_PROCS} --mca pml ucx --mca btl '^openib' \
-  python ./mpi_test.py &> ${PREFIX}.out
+  python ./run_me_with_mpirun.py &> ${PREFIX}.out

test/performance/mpi_2D_scaling/run.sh

@@ -1,4 +1,4 @@
-# This file is designed to be run via: mpiexec -np 4 python -W ignore testMPI.py
+# This file is designed to be run via: mpiexec -np 4 python -W ignore run_me_with_mpiexec.py
 # This does not use nosetests because I want to set verbosity differently based on rank (quiet if not rank 0)
 # By wrapping asserts in comm.rank == 0, only rank 0 should fail (should help with output)
 # Can run with different number of procs, but 4 is minimum to test all modes (pure MPI, pure shared mem, and mixed)
@@ -226,7 +226,7 @@ def run_fills(self, sim, natoms, nparams):
         else:
             raise RuntimeError("Improper sim type passed by test_fills_generator")
 
-        serial_layout = mdl.sim.create_layout(circuits, array_types=('E','EP','EPP'), derivative_dimension=nP)
+        serial_layout = mdl.sim.create_layout(circuits, array_types=('E','EP','EPP'), derivative_dimensions=(nP,))
 
         nE = serial_layout.num_elements
         nC = len(circuits)
@@ -246,7 +246,7 @@ def run_fills(self, sim, natoms, nparams):
         global_serial_layout = serial_layout.global_layout
 
         #Use a parallel layout to compute the same probabilities & their derivatives
-        local_layout = mdl.sim.create_layout(circuits, array_types=('E','EP','EPP'), derivative_dimension=nP,
+        local_layout = mdl.sim.create_layout(circuits, array_types=('E','EP','EPP'), derivative_dimensions=(nP,),
                                              resource_alloc=self.ralloc)
 
         vp_local = local_layout.allocate_local_array('e', 'd')
@@ -334,7 +334,7 @@ def setup_class(cls):
     tester = PureMPIParallel_Test()
     tester.setup_class()
     tester.ralloc = pygsti.baseobjs.ResourceAllocation(wcomm)
-    #tester.run_objfn_values('matrix','logl',4)
+    tester.run_objfn_values('matrix','logl',4)
     tester.run_fills('map', 1, None)
     tester.run_fills('map', 4, None)
     tester.run_fills('matrix', 4, 15)

test/performance/mpi_2D_scaling/mpi_test.py renamed to test/performance/mpi_2D_scaling/run_me_with_mpirun.py

@@ -0,0 +1,26 @@
+import subprocess
+import pytest
+import os
+from pathlib import Path
+
+try:
+    from mpi4py import MPI
+except (ImportError, RuntimeError):
+    MPI = None
+
+
+class MPITester:
+
+    @pytest.mark.skipif(MPI is None, reason="mpi4py could not be imported")
+    def test_all(self, capfd: pytest.LogCaptureFixture):
+        current_filepath = Path(os.path.abspath(__file__))
+        to_run = current_filepath.parents[0] / Path('run_me_with_mpiexec.py')
+        subprocess_args = (f"mpiexec -np 4 python -W ignore {str(to_run)}").split(' ')
+        
+        result = subprocess.run(subprocess_args, capture_output=False, text=True)
+        out, err = capfd.readouterr()
+        if len(out) + len(err) > 0:
+            msg = out + '\n'+ 80*'-' + err
+            raise RuntimeError(msg)
+        return
+

test/test_packages/mpi/test_mpi.py renamed to test/unit/mpi/run_me_with_mpiexec.py

@@ -300,10 +300,12 @@ def test_derivative(self):
                                          places=3)  # each *element* should match to 3 places
 
             if self.computes_lsvec:
-                lsvec = objfn.lsvec().copy()
-                dlsvec = objfn.dlsvec().copy()
-                self.assertArraysAlmostEqual(dterms / nEls, 2 * lsvec[:, None] * dlsvec / nEls,
-                                             places=4)  # each *element* should match to 4 places
+                arg1   = dterms / nEls
+                lsvec  = objfn.lsvec(v0).copy()
+                dlsvec = objfn.dlsvec(v0).copy()
+                arg2   = 2 * lsvec[:, None] * dlsvec / nEls
+                self.assertArraysAlmostEqual(arg1, arg2, places=4)  # each *element* should match to 4 places
+                return
 
     def test_approximate_hessian(self):
         if not self.enable_hessian_tests: