Merge pull request #6 from esa/khan-boundaries

Khan boundaries gradient fix

Author: waldemar-martens

.gitlab-ci.yml

@@ -163,12 +163,30 @@ test-wheels:
   needs:
     - python-wheels:manylinux2014
 
-publish-test:
+publish-gitlab:
   stage: deploy
   image: python:3.9
   script:
     - pip install twine
-    # upload to the project pypi (whl) (change testpypi to pypi to upload to PyPi)
+    # upload to repo (ESA Gitlab)
+    - twine upload --non-interactive --disable-progress-bar --verbose -u gitlab-ci-token -p $CI_JOB_TOKEN --repository-url https://gitlab.esa.int/api/v4/projects/$CI_PROJECT_ID/packages/pypi pyoptgra-*.whl
+    - twine upload --non-interactive --disable-progress-bar --verbose -u gitlab-ci-token -p $CI_JOB_TOKEN --repository-url https://gitlab.esa.int/api/v4/projects/$CI_PROJECT_ID/packages/pypi pyoptgra-*.tar.gz
+
+  needs:
+    - python-sdist
+    - test-wheels
+    - python-wheels:manylinux2014
+  only:
+    refs:
+      - master
+  when: manual
+
+publish-testpypi:
+  stage: deploy
+  image: python:3.9
+  script:
+    - pip install twine
+    # upload to the project pypi (whl)
     - twine upload --non-interactive --disable-progress-bar --verbose --repository testpypi -u __token__ -p $TEST_PYPI_TOKEN pyoptgra-*.whl
     # upload to the project pypi (sdist)
     - twine upload --non-interactive --disable-progress-bar --verbose --repository testpypi -u __token__ -p $TEST_PYPI_TOKEN pyoptgra-*.tar.gz

pyoptgra/optgra.py

@@ -13,7 +13,6 @@
 # and https://essr.esa.int/license/european-space-agency-community-license-v2-4-weak-copyleft
 
 from math import isfinite
-from collections import deque
 from typing import Callable, List, Optional, Tuple, Union
 
 import numpy as np
@@ -29,7 +28,7 @@
 
 
 class khan_function:
-    """Function to smothly enforce optimisation parameter bounds as Michal Khan used to do:
+    r"""Function to smothly enforce optimisation parameter bounds as Michal Khan used to do:
 
     .. math::
 
@@ -59,9 +58,14 @@ def __init__(self, lb: List[float], ub: List[float], unity_gradient: bool = True
         self._ub = np.asarray(ub)
         self._nx = len(lb)
 
-        # determine finite lower/upper bounds
-        finite_lb = np.isfinite(self._lb)
-        finite_ub = np.isfinite(self._ub)
+        # determine finite lower/upper bounds\
+        def _isfinite(a: np.ndarray):
+            """Custom _ function"""
+            almost_infinity = 1e300
+            return np.logical_and(np.isfinite(a), np.abs(a) < almost_infinity)
+
+        finite_lb = _isfinite(self._lb)
+        finite_ub = _isfinite(self._ub)
 
         # we only support cases where both lower and upper bounds are finite if given
         check = np.where(finite_lb != finite_ub)[0]
@@ -71,8 +75,14 @@ def __init__(self, lb: List[float], ub: List[float], unity_gradient: bool = True
                 "must be finite."
                 f"Detected mismatch at decision vector indices: {check}"
             )
+
+        # also exclude parameters where lower and upper bounds are identical
+        with np.errstate(invalid="ignore"):
+            # we ignore RuntimeWarning: invalid value encountered in subtract
+            nonzero_diff = abs(self._lb - self._ub) > 1e-9
+
         # store the mask of finite bounds
-        self.mask = finite_ub
+        self.mask = np.logical_and(finite_ub, nonzero_diff)
         self._lb_masked = self._lb[self.mask]
         self._ub_masked = self._ub[self.mask]
 
@@ -93,12 +103,14 @@ def _eval_inv(self, x_masked: np.ndarray) -> np.ndarray:
         arg = (2 * x_masked - self._ub_masked - self._lb_masked) / (
             self._ub_masked - self._lb_masked
         )
-        if np.any((arg < -1.0) | (arg > 1.0)):
+
+        clip_value = 1.0 - 1e-8  # avoid boundaries
+        if np.any((arg < -clip_value) | (arg > clip_value)):
             print(
                 "WARNING: Numerical inaccuracies encountered during khan_function inverse.",
                 "Clipping parameters to valid range.",
             )
-            arg = np.clip(arg, -1.0, 1.0)
+            arg = np.clip(arg, -clip_value, clip_value)
         return (np.arcsin(arg) - self._b) / self._a
 
     def _eval_grad(self, x_khan_masked: np.ndarray) -> np.ndarray:
@@ -249,31 +261,26 @@ def wrapped_fitness(x):
                 # if Khan function is used, we first need to convert to pagmo parameters
                 x = khanf.eval(x_khan=x)
 
-            fixed_x = x
-            lb, ub = problem.get_bounds()
-
             if force_bounds:
-                for i in range(problem.get_nx()):
-                    if x[i] < lb[i]:
-                        fixed_x[i] = lb[i]
-                    if x[i] > ub[i]:
-                        fixed_x[i] = ub[i]
+                fixed_x = np.clip(x, lb, ub)
+            else:
+                fixed_x = x
 
-            result = deque(problem.fitness(fixed_x))
+            # call pagmo fitness function
+            result = problem.fitness(fixed_x)
 
             # add constraints derived from box bounds
             if bounds_to_constraints:
-                for i in range(len(lb)):
-                    if isfinite(lb[i]):
-                        result.append(x[i] - lb[i])
+                # Add (x[i] - lb[i]) for finite lb[i] and (ub[i] - x[i]) for finite ub[i]
+                result = np.concatenate(
+                    [result, (x - lb)[np.isfinite(lb)], (ub - x)[np.isfinite(ub)]]
+                )
 
-                for i in range(len(ub)):
-                    if isfinite(ub[i]):
-                        result.append(ub[i] - x[i])
+            # reorder constraint order, optgra expects the merit function last, pagmo has it first
+            # equivalent to rotating in a dequeue
+            result = np.concatenate([result[1:], result[0:1]])
 
-            # optgra expects the fitness last, pagmo has the fitness first
-            result.rotate(-1)
-            return list(result)
+            return result.tolist()  # return a list
 
         return wrapped_fitness
 
@@ -284,17 +291,14 @@ def _wrap_gradient_func(
         force_bounds=False,
         khanf: Optional[khan_function] = None,
     ):
+
         # get the sparsity pattern to index the sparse gradients
         sparsity_pattern = problem.gradient_sparsity()
         f_indices, x_indices = sparsity_pattern.T  # Unpack indices
 
         # expected shape of the non-sparse gradient matrix
         shape = (problem.get_nf(), problem.get_nx())
 
-        # get problem parameters
-        lb, ub = problem.get_bounds()
-        nx = problem.get_nx()
-
         def wrapped_gradient(x):
 
             # we are using vectorisation internally -> convert to ndarray
@@ -304,6 +308,10 @@ def wrapped_gradient(x):
                 # if Khan function is used, we first need to convert to pagmo parameters
                 x = khanf.eval(x_khan=x)
 
+            # get problem parameters
+            lb, ub = problem.get_bounds()
+            nx = problem.get_nx()
+
             # force parameters to lower and upper bounds if needed
             if force_bounds:
                 fixed_x = np.clip(x, lb, ub)
@@ -321,15 +329,15 @@ def wrapped_gradient(x):
             result[f_indices, x_indices] = sparse_values
 
             # add box-derived constraints
+            result = result.tolist()
             if bounds_to_constraints:
-
                 # lower bound gradients
                 finite_indices = np.isfinite(lb)  # Boolean mask for valid indices
                 box_lb_grads = np.eye(nx)[finite_indices]
 
                 # upper bound gradients
                 finite_indices = np.isfinite(ub)  # Boolean mask for valid indices
-                box_ub_grads = np.eye(nx)[finite_indices]
+                box_ub_grads = -1.0 * np.eye(nx)[finite_indices]
 
                 # append box bounds to gradient matrix
                 result = np.concatenate([result, box_lb_grads, box_ub_grads])

pyproject.toml

@@ -12,7 +12,7 @@ license = {text = "GPL-3.0 or ESCL-2.4"}
 name = "pyoptgra"
 readme = "README.rst"
 requires-python = ">=3.9"
-version = "1.0.0"
+version = "1.0.1"
 
 [build-system]
 build-backend = "scikit_build_core.build"

tests/python/test.py

@@ -308,7 +308,7 @@ def gradient_no_constraints_test(self):
     def gradient_with_constraints_test(self):
         # 1. Run Luksan-Vlcek problem with optgra
         prob = pygmo.problem(luksan_vlcek())
-        prob.c_tol = 1e-6
+        prob.c_tol = 1e-7
         og = pyoptgra.optgra(
             optimization_method=1,
             max_iterations=100,
@@ -325,6 +325,8 @@ def gradient_with_constraints_test(self):
 
         # objective function
         self.assertLess(pop.champion_f[0], 2.26)
+        # checking exact value as regression test
+        self.assertEqual(pop.champion_f[0], 0.82929210248477)
 
         # equality constraints
         for i in [1, 2, 3, 4]:
@@ -351,6 +353,8 @@ def gradient_with_constraints_test(self):
 
         # objective function
         self.assertLess(pop2.champion_f[0], 2.26)
+        # checking exact value as regression test
+        self.assertEqual(pop2.champion_f[0], 0.8292921025820391)
 
         # equality constraints
         for i in [1, 2, 3, 4]: