Make default optimizer for acv use differentialevolution to get initial guess

jdjakem · jdjakem · commit f0c55a972a96 · 2025-10-13T12:37:01.000-06:00
diff --git a/environment.yml b/environment.yml
@@ -19,6 +19,7 @@ dependencies:
   - networkx
   - pytorch
   - numba
+  - cvxpy
   # make html in docs folder
   - sphinx
   - sphinx-gallery
diff --git a/environment_w_extras.yml b/environment_w_extras.yml
@@ -19,6 +19,7 @@ dependencies:
   - networkx
   - pytorch
   - numba
+  - cvxpy
   # gpu support
   - torchvision
   - torchaudio
diff --git a/examples/expdesign/plot_bayesoed4params.py b/examples/expdesign/plot_bayesoed4params.py
@@ -138,7 +138,7 @@
 
 # %%
 # Plot the OED
-# -----------
+# ------------
 axs = plt.subplots(1, 2, figsize=(8, 6), sharey=True)[1]
 # Reshape weights so that each row correspond to one of the OED states
 design_weights = design_weights.reshape((2, -1))
diff --git a/examples/expdesign/plot_bayesoed4pred.py b/examples/expdesign/plot_bayesoed4pred.py
@@ -170,7 +170,7 @@
 
 # %%
 # Plot the OED
-# -----------
+# ------------
 fig, axs = plt.subplots(1, 2, figsize=(8, 6), sharey=True)
 # Reshape weights so that each row correspond to one of the OED states
 design_weights = design_weights.reshape((2, -1))
diff --git a/examples/interface/plot_model_timing.py b/examples/interface/plot_model_timing.py
@@ -1,7 +1,7 @@
 r"""
 Model Timing
 ------------
- It is often useful to be able to track the time needed to evaluate a function. 
+It is often useful to be able to track the time needed to evaluate a function.
 This tutorial demonstrats how to use model to time function calls.
 All model class can track time.
 
diff --git a/examples/surrogates/plot_barycentric_interpolation.py b/examples/surrogates/plot_barycentric_interpolation.py
@@ -27,7 +27,7 @@
 marginals = [stats.uniform(-1, 2), stats.uniform(-1, 2)]
 quad_rules_1d = [GaussQuadratureRule(marginal) for marginal in marginals]
 bases_1d = [
-    UnivariateLagrangeBasis(quad_rule, 1) for quad_rule in quad_rules_1d
+    UnivariateLagrangeBasis(quad_rule, 2) for quad_rule in quad_rules_1d
 ]
 
 
diff --git a/pyapprox/multifidelity/acv.py b/pyapprox/multifidelity/acv.py
@@ -1364,12 +1364,16 @@ def plot_recursion_dag(self, ax):
     #     partition_ratios.grad.zero_()
     #     return val.item(), grad
 
-    def get_npartition_bounds(self) -> Array:
+    def get_npartition_bounds(self, total_cost: float) -> Array:
         nunknowns = self._npartitions - 1
+        # tighter bounds helps global optimizer so insted of using inf,
+        # set upperbounds to be maximum number of samples if all samples
+        # were allocated to a single model, for each model
         bounds = self._bkd.stack(
             (
                 self._bkd.full((nunknowns,), self._npartitions_lower_bound),
-                self._bkd.full((nunknowns,), np.inf),
+                # self._bkd.full((nunknowns,), np.inf),
+                total_cost / self._costs[1:],
             ),
             axis=1,
         )
@@ -1384,9 +1388,17 @@ def get_default_optimizer(self) -> Optimizer:
         # init_gen.set_bounds(self.get_npartition_bounds())
         # optimizer = ConstrainedMultiStartOptimizer(local_optimizer)
         # optimizer.set_initial_iterate_generator(init_gen)
-        global_optimizer = ScipyConstrainedNelderMeadOptimizer(
-            opts={"maxiter": 500}
+        # global_optimizer = ScipyConstrainedNelderMeadOptimizer(
+        #     opts={"maxiter": 500}
+        # )
+        from pyapprox.optimization.scipy import (
+            ScipyConstrainedDifferentialEvolutionOptimizer,
+        )
+
+        global_optimizer = ScipyConstrainedDifferentialEvolutionOptimizer(
+            opts={"maxiter": 3}
         )
+
         local_optimizer = ScipyConstrainedOptimizer()
         optimizer = ChainedOptimizer(global_optimizer, local_optimizer)
         optimizer.set_verbosity(0)
@@ -1415,7 +1427,7 @@ def _allocate_samples_minimize(
         self._optimizer.set_objective_function(objective)
         constraints = [ACVPartitionConstraint(self, target_cost)]
         self._optimizer.set_constraints(constraints)
-        self._optimizer.set_bounds(self.get_npartition_bounds())
+        self._optimizer.set_bounds(self.get_npartition_bounds(target_cost))
         init_iterate = self._bkd.full((self._nmodels - 1, 1), 1.0)
         result = self._optimizer.minimize(init_iterate)
         return result
diff --git a/pyapprox/multifidelity/tests/test_acv.py b/pyapprox/multifidelity/tests/test_acv.py
@@ -724,17 +724,6 @@ def test_best_model_subset_estimator(self):
         target_cost = 10
         est._save_candidate_estimators = True
         est.allocate_samples(target_cost)
-        # {
-        #     "verbosity": 1,
-        #     "nprocs": 1,
-        #     "scaling": 1,
-        #     "init_guess": {
-        #         "disp": True,
-        #         "maxiter": 300,
-        #         "lower_bound": 1e-10,
-        #     },
-        # },
-
         criteria = bkd.array(
             [e[0]._optimized_criteria for e in est._candidate_estimators]
         )
diff --git a/pyapprox/surrogates/affine/linearsystemsolvers.py b/pyapprox/surrogates/affine/linearsystemsolvers.py
@@ -1358,6 +1358,8 @@ class BasisPursuitDensoisingCVXRegressionSolver(
     CVXOPTOptionsMixin, LinearSystemSolver
 ):
     def __init__(self, penalty: float, backend: BackendMixin = NumpyMixin):
+        if not package_available("cvxopt"):
+            raise ImportError("cvxopt is not installed")
         self._penalty = penalty
         super().__init__(backend)
 
diff --git a/pyapprox/surrogates/affine/tests/test_linsolvers.py b/pyapprox/surrogates/affine/tests/test_linsolvers.py
@@ -148,6 +148,7 @@ def test_linearly_constrained_lstsq_solver(self):
         assert bkd.allclose(C @ coef, d, atol=1e-15)
         assert bkd.allclose(A @ coef, y, atol=1e-1)
 
+    @unittest.skipIf(not package_available("cvxpy"), "cvxpy not installed")
     def test_quantile_regression(self):
         bkd = self.get_backend()
         quantile = 0.5
diff --git a/pyapprox/surrogates/kernels/kernels.py b/pyapprox/surrogates/kernels/kernels.py
@@ -284,6 +284,11 @@ def __init__(
         super().__init__(backend)
         if transform is None:
             transform = IdentityHyperParameterTransform(backend=self._bkd)
+        if not self._bkd.bkd_equal(backend, transform._bkd):
+            raise TypeError(
+                f"Transform backend {transform._bkd.__name__} "
+                f"does not macth backend {backend.__name__}"
+            )
         if constant_bounds is None:
             constant_bounds = [-self._bkd.inf(), self._bkd.inf()]
         self._const = HyperParameter(
diff --git a/pyapprox/surrogates/univariate/orthopoly.py b/pyapprox/surrogates/univariate/orthopoly.py
@@ -89,7 +89,7 @@ def _ncoefs(self) -> int:
     def _get_recursion_coefficients(self, ncoefs: int) -> Array:
         raise NotImplementedError
 
-    def set_nterms(self, nterms):
+    def set_nterms(self, nterms: int):
         """Compute and set the recursion coefficients of the polynomial."""
         if self._rcoefs is None or self._ncoefs() < nterms:
             # TODO implement increment of recursion coefficients
diff --git a/pyproject.toml b/pyproject.toml
@@ -35,6 +35,7 @@ dependencies = [
     'numba',
     'scikit-fem',
     'umbridge',
+    'cvxpy',
 ]
 
 
diff --git a/setup.py b/setup.py
@@ -58,6 +58,7 @@
         "numba",
         "scikit-fem",
         "umbridge",
+        "cvxpy",
     ],
     extras_require={
         "docs": docs_extras,
@@ -134,6 +135,12 @@
 
 # pip install "cvxpy[CVXOPT]"
 
+# when running python script.py and recieve error
+# OMP: Error #15: Initializing libomp.dylib, but found libomp.dylib already initialized.
+# Then run command with
+# KMP_DUPLICATE_LIB_OK=TRUE  python script.py
+# warning above is a hack and may not always work, it is better to fix install of packages, which can require recreating environment
+
 # add the following to start of script to help debug
 # torch grad errors
 # import torch
diff --git a/tutorials/multifidelity_surrogates/plot_multifidelity_gp.py b/tutorials/multifidelity_surrogates/plot_multifidelity_gp.py
@@ -116,7 +116,6 @@
 
 from pyapprox.variables.joint import IndependentMarginalsVariable
 from pyapprox.surrogates.gaussianprocess.exactgp import (
-    # GreedyMultifidelityIntegratedVarianceSampler,
     SequentialMultiLevelGaussianProcess,
     MOExactGaussianProcess,
     ExactGaussianProcess,
@@ -128,29 +127,31 @@
 from pyapprox.surrogates.univariate.base import Monomial1D
 from pyapprox.surrogates.affine.basis import MultiIndexBasis
 from pyapprox.expdesign.sequences import HaltonSequence
+from pyapprox.util.backends.torch import TorchMixin as bkd
 
 np.random.seed(1)
 nvars = 1
 variable = IndependentMarginalsVariable(
-    [stats.uniform(0, 1) for ii in range(nvars)]
+    [stats.uniform(0, 1) for ii in range(nvars)], backend=bkd
 )
 
 nmodels = 2
-sigmas = [1, 0.1]
+sigmas = [1.0, 0.1]
 length_scale = np.hstack([np.full(nvars, 0.3)] * nmodels)
 kernels = [
     ConstantKernel(
         sigmas[nn],
         (1e-3, 1e1),
-        transform=LogHyperParameterTransform(),
+        transform=LogHyperParameterTransform(backend=bkd),
+        backend=bkd,
     )
-    * MaternKernel(np.inf, length_scale[nn], (0.1, 1.0), nvars)
+    * MaternKernel(np.inf, length_scale[nn], (0.1, 1.0), nvars, backend=bkd)
     for nn in range(nmodels)
 ]
 
-rho = np.array([[1.0]])
+rho = bkd.array([[1.0]])
 nterms_1d = 1
-basis = MultiIndexBasis([Monomial1D() for ii in range(nvars)])
+basis = MultiIndexBasis([Monomial1D(backend=bkd) for ii in range(nvars)])
 basis.set_tensor_product_indices([nterms_1d] * nvars)
 kernel_scalings = [MonomialExpansion(basis) for nn in range(nmodels - 1)]
 [scaling.set_coefficients(rho) for scaling in kernel_scalings]
@@ -159,7 +160,7 @@
 # sort for plotting covariance matrices
 nsamples_per_model = [100, 50]
 train_samples_per_model = [
-    np.sort(
+    bkd.sort(
         HaltonSequence(nvars, start_idx=1000 * nn + 1, variable=variable).rvs(
             nsamples_per_model[nn]
         )
@@ -175,7 +176,7 @@
 # %%
 # Now lets plot the structure of the multi-fidelity covariance matrix when using two models a linear polynomial scaling :math:`\rho(x)=a+bx`. In comparison to a scalar scaling the correlation between the low and high-fidelity models changes as a function of the input :math:`\rv`.
 ax = plt.subplots(1, 1, figsize=(1 * 8, 6))[1]
-rho_2 = np.array([0.0, 1.0])[:, None]  # [a, b]
+rho_2 = bkd.array([0.0, 1.0])[:, None]  # [a, b]
 basis_2 = MultiIndexBasis([Monomial1D() for ii in range(nvars)])
 basis_2.set_tensor_product_indices([2] * nvars)
 kernel_scalings_2 = [MonomialExpansion(basis_2) for nn in range(nmodels - 1)]
@@ -198,7 +199,7 @@ def scale(degree, x, rho, kk):
 def f0(x):
     # y = x.sum(axis=0)[:, None]
     y = x[0:1].T
-    return ((y * 3 - 1) ** 2 + 1) * np.sin((y * 10 - 2)) / 5
+    return ((y * 3 - 1) ** 2 + 1) * bkd.sin((y * 10 - 2)) / 5
 
 
 def f1(degree, rho, x):
@@ -211,7 +212,7 @@ def f1(degree, rho, x):
 # Now build a GP with 8 samples from the low-fidelity model and four samples from the high-fidelity model
 nsamples_per_model = [8, 4]
 train_samples_per_model = [
-    np.sort(
+    bkd.sort(
         HaltonSequence(nvars, start_idx=1000 * nn + 1, variable=variable).rvs(
             nsamples_per_model[nn]
         ),
@@ -221,7 +222,7 @@ def f1(degree, rho, x):
 ]
 degree = 0
 
-true_rho = np.full((nmodels - 1) * degree + 1, 0.9)
+true_rho = bkd.full(((nmodels - 1) * degree + 1,), 0.9)
 models = [f0, partial(f1, degree, true_rho)]
 train_values_per_model = [
     f(x) for f, x in zip(models, train_samples_per_model)
@@ -234,7 +235,7 @@ def f1(degree, rho, x):
 ax = plt.subplots(1, 1, figsize=(8, 6))[1]
 ax.plot(train_samples_per_model[0][0], train_values_per_model[0], "ko")
 ax.plot(train_samples_per_model[1][0], train_values_per_model[1], "rs")
-xx = np.linspace(0, 1, 101)[None, :]
+xx = bkd.linspace(0, 1, 101)[None, :]
 ax.plot(xx[0], models[0](xx), "k-", label=r"$f_0$")
 ax.plot(xx[0], models[1](xx), "r-", label=r"$f_1$")
 gp.plot_1d(
@@ -261,8 +262,9 @@ def f1(degree, rho, x):
 sf_kernel = ConstantKernel(
     sigmas[0],
     (1e-3, 1e1),
-    transform=LogHyperParameterTransform(),
-) * MaternKernel(np.inf, length_scale[0], (1e-1, 1.0), nvars)
+    transform=LogHyperParameterTransform(backend=bkd),
+    backend=bkd,
+) * MaternKernel(bkd.inf(), length_scale[0], (1e-1, 1.0), nvars, backend=bkd)
 sf_gp = ExactGaussianProcess(nvars, sf_kernel)
 sf_gp.set_optimizer(ncandidates=3)
 sf_gp.fit(train_samples_per_model[1], train_values_per_model[1])
@@ -316,14 +318,17 @@ def f1(degree, rho, x):
 # The following plot compares a sequential multi-level GP with the exact co-kriging on the toy problem introduced earlier.
 
 sml_kernels = [
-    MaternKernel(np.inf, 0.1, (1e-10, 1), nvars) for ii in range(nmodels)
+    MaternKernel(bkd.inf(), 0.1, (1e-10, 1), nvars, backend=bkd)
+    for ii in range(nmodels)
 ]
-sml_scaling_basis = MultiIndexBasis([Monomial1D() for ii in range(nvars)])
+sml_scaling_basis = MultiIndexBasis(
+    [Monomial1D(backend=bkd) for ii in range(nvars)]
+)
 sml_scaling_basis.set_tensor_product_indices([1] * nvars)
 sml_scalings = [
     MonomialExpansion(sml_scaling_basis) for nn in range(nmodels - 1)
 ]
-[scaling.set_coefficients(np.array([[1.0]])) for scaling in sml_scalings]
+[scaling.set_coefficients(bkd.array([[1.0]])) for scaling in sml_scalings]
 sml_gp = SequentialMultiLevelGaussianProcess(sml_kernels, sml_scalings)
 sml_gp.fit(train_samples_per_model, train_values_per_model)
 # gps = sml_gp.gaussian_processes_per_level()
@@ -402,8 +407,8 @@ def f1(degree, rho, x):
 # The following plots the correlations between the models used previously in this tutorial with the non-linear ensemble used here
 axs = plt.subplots(1, 2, figsize=(2 * 8, 6))[1]
 nonlinear_models = [
-    lambda xx: np.sin(np.pi * 8 * xx.T),
-    lambda xx: (xx.T - np.sqrt(2)) * np.sin(np.pi * 8 * xx.T) ** 2,
+    lambda xx: bkd.sin(np.pi * 8 * xx.T),
+    lambda xx: (xx.T - np.sqrt(2)) * bkd.sin(np.pi * 8 * xx.T) ** 2,
 ]
 axs[0].plot(models[0](xx), models[1](xx))
 axs[1].plot(nonlinear_models[0](xx), nonlinear_models[1](xx))
diff --git a/tutorials/sciml/plot_learning_greens_functions.py b/tutorials/sciml/plot_learning_greens_functions.py
@@ -142,7 +142,10 @@ def greens_solution(kernel, forc, xx):
     basis, basis, 1 / bkd.arange(1, nterms + 1, dtype=float), [1e-2, 1]
 )
 const_kernel = ConstantKernel(
-    10.0, [1e-2, 1e4], transform=LogHyperParameterTransform(), backend=bkd
+    10.0,
+    [1e-2, 1e4],
+    transform=LogHyperParameterTransform(backend=bkd),
+    backend=bkd,
 )
 final_kernel = const_kernel * hs_kernel
 green_sol_hs = greens_solution(final_kernel, forc_fun, plot_xx)

Original file line number	Diff line number	Diff line change
`@@ -27,7 +27,7 @@`
`27`	`27`	`marginals = [stats.uniform(-1, 2), stats.uniform(-1, 2)]`
`28`	`28`	`quad_rules_1d = [GaussQuadratureRule(marginal) for marginal in marginals]`
`29`	`29`	`bases_1d = [`
`30`		`- UnivariateLagrangeBasis(quad_rule, 1) for quad_rule in quad_rules_1d`
	`30`	`+ UnivariateLagrangeBasis(quad_rule, 2) for quad_rule in quad_rules_1d`
`31`	`31`	`]`
`32`	`32`
`33`	`33`