Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove old _optimize method and the test #1066

Merged
merged 1 commit into from
Jun 19, 2024
Merged

Remove old _optimize method and the test #1066

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
158 changes: 0 additions & 158 deletions desc/equilibrium/equilibrium.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@
from desc.profiles import PowerSeriesProfile, SplineProfile
from desc.transform import Transform
from desc.utils import (
Timer,
check_nonnegint,
check_posint,
copy_coeffs,
Expand Down Expand Up @@ -1990,163 +1989,6 @@ def optimize(

return things[0], result

def _optimize( # noqa: C901
self,
objective,
constraint=None,
ftol=1e-6,
xtol=1e-6,
maxiter=50,
verbose=1,
copy=False,
solve_options=None,
perturb_options=None,
):
"""Optimize an equilibrium for an objective.

Parameters
----------
objective : ObjectiveFunction
Objective function to optimize.
constraint : ObjectiveFunction
Objective function to satisfy. Default = fixed-boundary force balance.
ftol : float
Relative stopping tolerance on objective function value.
xtol : float
Stopping tolerance on optimization step size.
maxiter : int
Maximum number of optimization steps.
verbose : int
Level of output.
copy : bool, optional
Whether to update the existing equilibrium or make a copy (Default).
solve_options : dict
Dictionary of additional options used in Equilibrium.solve().
perturb_options : dict
Dictionary of additional options used in Equilibrium.perturb().

Returns
-------
eq_new : Equilibrium
Optimized equilibrium.

"""
import inspect
from copy import deepcopy

from desc.optimize.tr_subproblems import update_tr_radius
from desc.optimize.utils import check_termination
from desc.perturbations import optimal_perturb

solve_options = {} if solve_options is None else solve_options
perturb_options = {} if perturb_options is None else perturb_options

if constraint is None:
constraint = get_equilibrium_objective(eq=self)

timer = Timer()
timer.start("Total time")

if not objective.built:
objective.build()
if not constraint.built:
constraint.build()

cost = objective.compute_scalar(objective.x(self))
perturb_options = deepcopy(perturb_options)
tr_ratio = perturb_options.get(
"tr_ratio",
inspect.signature(optimal_perturb).parameters["tr_ratio"].default,
)

if verbose > 0:
objective.print_value(objective.x(self))

params = orig_params = self.params_dict.copy()

iteration = 1
success = None
while success is None:
timer.start("Step {} time".format(iteration))
if verbose > 0:
print("====================")
print("Optimization Step {}".format(iteration))
print("====================")
print("Trust-Region ratio = {:9.3e}".format(tr_ratio[0]))

# perturb + solve
(_, predicted_reduction, dc_opt, dc, c_norm, bound_hit) = optimal_perturb(
self,
constraint,
objective,
copy=False,
**perturb_options,
)
self.solve(objective=constraint, **solve_options)

# update trust region radius
cost_new = objective.compute_scalar(objective.x(self))
actual_reduction = cost - cost_new
trust_radius, ratio = update_tr_radius(
tr_ratio[0] * c_norm,
actual_reduction,
predicted_reduction,
np.linalg.norm(dc_opt),
bound_hit,
)
tr_ratio[0] = trust_radius / c_norm
perturb_options["tr_ratio"] = tr_ratio

timer.stop("Step {} time".format(iteration))
if verbose > 0:
objective.print_value(objective.x(self))
print("Predicted Reduction = {:10.3e}".format(predicted_reduction))
print("Reduction Ratio = {:+.3f}".format(ratio))
if verbose > 1:
timer.disp("Step {} time".format(iteration))

# stopping criteria
success, message = check_termination(
actual_reduction,
cost,
np.linalg.norm(dc),
c_norm,
np.inf, # TODO: add g_norm
ratio,
ftol,
xtol,
0, # TODO: add gtol
iteration,
maxiter,
0,
np.inf,
)
if actual_reduction > 0:
params = self.params_dict.copy()
cost = cost_new
else:
# reset equilibrium to last good params
self.params_dict = params
if success is not None:
break

iteration += 1

timer.stop("Total time")
print("====================")
print("Done")
if verbose > 0:
print(message)
if verbose > 1:
timer.disp("Total time")

if copy:
eq = self.copy()
self.params = orig_params
else:
eq = self
return eq

def perturb(
self,
deltas,
Expand Down
20 changes: 0 additions & 20 deletions tests/test_examples.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,6 @@
)
from desc.optimize import Optimizer
from desc.profiles import FourierZernikeProfile, PowerSeriesProfile
from desc.vmec_utils import vmec_boundary_subspace

from .utils import area_difference_desc, area_difference_vmec

Expand Down Expand Up @@ -179,25 +178,6 @@ def test_1d_optimization():
np.testing.assert_allclose(eq.compute("R0/a")["R0/a"], 2.5, rtol=2e-4)


@pytest.mark.regression
@pytest.mark.optimize
def test_1d_optimization_old():
"""Tests 1D optimization for target aspect ratio."""
eq = get("SOLOVEV")
objective = ObjectiveFunction(AspectRatio(eq=eq, target=2.5))
eq._optimize(
objective,
copy=False,
solve_options={"verbose": 0},
perturb_options={
"dZb": True,
"subspace": vmec_boundary_subspace(eq, ZBS=[0, 1]),
},
)

np.testing.assert_allclose(eq.compute("R0/a")["R0/a"], 2.5, rtol=1e-3)


def run_qh_step(n, eq):
"""Run 1 step of the precise QH optimization example from Landreman & Paul."""
print(f"==========QH step {n+1}==========")
Expand Down
Loading