Examples¶
NLP Interface¶
Basic introduction to defining an own NLP and using the NLP interface.
# ---
# jupyter:
# jupytext:
# formats: ipynb,py:percent
# text_representation:
# extension: .py
# format_name: percent
# format_version: '1.3'
# jupytext_version: 1.19.1
# kernelspec:
# display_name: optsam
# language: python
# name: python3
# ---
# %% [markdown]
# # NLP interface & Solver basics
#
# THis tutorial illustrates the generic NLP interface, access to a solver, and plotting the traces.
# %%
import optsam as op
import numpy as np
import matplotlib.pyplot as plt
# %% [markdown]
# Define a problem, here: 2 dimensional, with 2 SOS features, and bounds [-1,-2]x[2,3]
# %%
class MyNLP:
def __init__(self):
self.dimension = 2
self.types = [op.OT.sos] * 2
self.bounds = np.array([[-2,-1], [2,3]])
self.b = 3
def evaluate(self, x):
phi = np.array([ x[0]-1,
self.b*(x[1]-x[0]**2) ])
J = np.array([[ 1, 0 ],
[ -2*self.b*x[0], self.b ]])
return phi, J
nlp = MyNLP()
# %% [markdown]
# Create a solver
# %%
sol = op.NLP_Solver()
sol.setPyProblem(nlp)
sol.setOptions(stepMax=.5, damping=1e-4)
sol.setTracing(trace_x=True, trace_errs=True)
# %% [markdown]
# Call the solver. Here 20 times in a row, each time automatically initialized with uniform in the bounds (default implementation of nlp.getInitializationSample)
# %%
trace_x = []
trace_err = []
for i in range(20):
ret = sol.solve(1)
print(ret)
trace_x.append(sol.getTrace_x())
trace_err.append(sol.getTrace_errs())
sol.clearTracing()
trace_f = [np.sum(E, axis=1) for E in trace_err]
# %% [markdown]
# The following creates a grid X of input points, and evaluates the fct on X
# %%
nlp = sol.getProblem()
B = nlp.bounds
X = [None] * nlp.dimension
for i in range(nlp.dimension):
X[i] = np.linspace(B[0][i], B[1][i], 30)
X = np.stack(np.meshgrid(*X, indexing='ij'), axis=-1)
fX = np.array([nlp.eval_scalar(x)[0] for x in X.reshape(-1, nlp.dimension)])
fX = fX.reshape(X.shape[:-1])
# %% [markdown]
# ... to prepare plotting.
# %%
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(121)
ax1.contour(X[:,:,0], X[:,:,1], fX, 200)
for x in trace_x:
ax1.plot(x[:,0], x[:,1], 'o-r', ms=3)
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_wireframe(X[:,:,0], X[:,:,1], fX)
for x,f in zip(trace_x, trace_f):
ax2.plot(x[:,0], x[:,1], f, 'o-r', ms=3)
plt.show()
# %% [markdown]
# Finally, an example to check the derivatives (Jacobian of all problem features) at random initialization points:
# %%
for i in range(20):
x = sol.getProblem().getInitializationSample()
r = sol.getProblem().checkJacobian(x, 1e-6)
# print(r, x)
# %%
Test Problems¶
Browsing through all test problems.
# ---
# jupyter:
# jupytext:
# cell_metadata_filter: -all
# formats: ipynb,py:percent
# text_representation:
# extension: .py
# format_name: percent
# format_version: '1.3'
# jupytext_version: 1.19.1
# kernelspec:
# display_name: venv (3.12.3)
# language: python
# name: python3
# ---
# %% [markdown]
# # Test Problems
#
# This script loops through all test problems, displays their signature, cost function, and -- if robotics problem -- komo scene.
# %%
import optsam as op
import numpy as np
import matplotlib.pylab as plt
import time
def scalar_objective(nlp: op.NLP, x):
phi, _ = nlp.evaluate(x)
ty = nlp.getTypes()
if len(phi) != len(ty):
return nlp.eval_scalar(x)[0]
cost = 0.0
for value, feature_type in zip(phi, ty):
if feature_type == op.OT.f:
cost += value
elif feature_type == op.OT.sos:
cost += value * value
return cost
# %% [markdown]
#
# The following is a plotting helper problem: It creates a 2D grid; if the problem is 2D it evaluates on that grid; otherwise on a random hyperplane (determined by x0). Then plots.
# %%
def display_2d_unconstrained(nlp: op.NLP, p_name, resolution=30):
B = nlp.bounds
x0 = nlp.getInitializationSample()
dim = 2
X = [None] * dim
for i in range(dim):
X[i] = np.linspace(B[0][i], B[1][i], resolution)
X = np.stack(np.meshgrid(*X, indexing='ij'), axis=-1)
if nlp.dimension==2:
fX = np.array([scalar_objective(nlp, x) for x in X.reshape(-1, dim)])
else:
f = lambda x: scalar_objective(nlp, np.concatenate((x, x0[2:])))
fX = np.array([f(x) for x in X.reshape(-1, dim)])
fX = fX.reshape(X.shape[:-1])
fig, ax = plt.subplots()
ax.contour(X[:,:,0], X[:,:,1], fX, 200)
# ax.plot(x[:,0], x[:,1], 'o-r', ms=3)
ax.set_title(p_name)
plt.show()
# %% [markdown]
#
# The following displays a problem more generically: If it has constraints, it is converted to unconstraint (using the Augmented Lagrangian). If it is a robotics problem, it also displays the komo scene.
# %%
def display_any(nlp: op.NLP, p_name):
ty = nlp.types
komo = nlp.as_KOMO()
if komo is not None:
for _ in range(5):
# x = nlp.getUniformSample()
x = nlp.getInitializationSample()
nlp.evaluate(x)
komo.view(False, 'random init')
time.sleep(.2)
if (op.OT.eq in ty) or (op.OT.ineq in ty):
nlp_org = nlp
nlp = nlp.aug_lag(1e1, -1.)
display_2d_unconstrained(nlp, p_name)
# %% [markdown]
#
# We can now loop through all test problems, print their signatur, and display.
# %%
def main():
problems = op.get_NLP_Problem_names() # some pre-defined benchmark problems
print(problems)
# p = problems[5]
for p in problems:
nlp = op.make_NLP_Problem(p)
print('===\n', p, nlp.report(1), '===')
display_any(nlp, p)
#
if __name__ == "__main__":
main()
Comparing Solvers¶
Looping through solvers and test problems to compare.
# ---
# jupyter:
# jupytext:
# cell_metadata_filter: -all
# formats: ipynb,py:percent
# text_representation:
# extension: .py
# format_name: percent
# format_version: '1.3'
# jupytext_version: 1.19.1
# ---
# %% [markdown]
# # Solvers
#
# This script runs several solvers on several test problems, displaying their performance.
# %%
import optsam as op
import numpy as np
import numpy.typing as npt
import matplotlib.pylab as plt
import time
method_info = {
'none': (False),
'GradientDescent': (False), 'Rprop': (False), 'LBFGS': (False), 'Newton': (False),
'AugmentedLag': (True), 'LogBarrier': (True), 'slackGN_logBarrier': (True), 'SquaredPenalty': (True), 'singleSquaredPenalty': (True),
'slackGN': (True),
'NLopt': (True), 'Ipopt': (True), 'slackGN_Ipopt': (True), 'Ceres': (True),
'LSBO': (False), 'greedy': (False), 'NelderMead': (False),
'CMA': (False), 'LS_CMA': (False), 'ES': (False),
}
def run(nlp: op.NLP, method: op.OptMethod, fixed_x0: npt.NDArray[np.float64]):
sol = op.NLP_Solver()
sol.setProblem(nlp)
sol.setOptions(method=method, stopTolerance=1e-4, damping=1e-3)
sol.setInitialization(fixed_x0)
ret = sol.solve()
print(f'-- method {method}: {ret}')
return sol.getTrace_best()
def main():
problems = op.get_NLP_Problem_names()
print('-- all problems:', problems)
# problems = ['square', 'Rugged', 'Rastrigin', 'Rosenbrock', 'Ackley', 'Himmelblau', 'Box', 'Modes', 'Wedge', 'HalfCircle', 'LinearProgram', 'IK', 'IKobstacle', 'IKtorus', 'PushToReach', 'StableSphere', 'SpherePacking', 'MinimalConvexCore']
problems = ['square', 'Modes', 'LinearProgram', 'IK', 'IKobstacle', 'SpherePacking', 'MinimalConvexCore']
methods = [op.OptMethod.AugmentedLag, op.OptMethod.LBFGS, op.OptMethod.Rprop, op.OptMethod.Newton ]
probs = problems
traces = dict()
seeds = 2
n = len(probs)
fig, axes = plt.subplots(4, (n+3)//4, figsize=(20,15))
for ax, p in zip(axes.reshape(-1), probs):
nlp = op.make_NLP_Problem(p)
ty = nlp.getTypes()
is_constrained = (op.OT.eq in ty) or (op.OT.ineq in ty)
print('===', p, '===')
print(nlp.report(1))
print('=============')
for m in methods:
traces[m.name] = []
for s in range(seeds):
print('--', s)
x0 = nlp.getInitializationSample()
for m in methods:
is_constrained_method = method_info[m.name]
print('--', m)
if is_constrained and not is_constrained_method:
nlp_tmp = nlp.aug_lag(muSquaredPenalty=1e2)
else:
nlp_tmp = nlp
best_trace = run(nlp_tmp, m, x0)
traces[m.name].append(best_trace)
print('-- plot')
ax.set_title(f'{p} ({nlp.dimension}D)')
color = 0
x_min = 1
for m in methods:
for i,t in enumerate(traces[m.name]):
ax.plot(t, color=f'C{color}', label=(str(m.name) if i==0 else None))
x_min = min(x_min, np.min(t))
color += 1
if x_min>1e-10:
ax.set_yscale('log')
ax.legend(loc="upper right")
del nlp
fig.tight_layout()
plt.show()
if __name__ == "__main__":
main()
NLP Sampling¶
This creates the ‘IK’ problem, and shows how to call the default NLP sampler, as well as display and evaluate the samples.
import optsam as op
import numpy as np
import time
def generate_samples(nlp : op.NLP, n):
solver = op.NLP_Sampler(nlp)
#solver.setOptions() #access to many options - let's take default
data = []
while len(data)<n:
ret = solver.sample()
print(ret)
if ret.feasible:
data.append(ret.x)
return np.stack(data)
def display_samples(nlp : op.NLP, data):
komo = nlp.as_KOMO()
assert komo is not None, 'this works only for komo problems'
for x in data:
phi, J = nlp.evaluate(x)
err_eq = np.sum(phi[nlp.types==op.OT.eq])
print('data:', x, '\nphi:', phi, '\nerr:', err_eq)
komo.view(False, f'data {x}')
time.sleep(1.)
def main():
problems = op.get_NLP_Problem_names()
print('== these are all problems:', problems)
nlp = op.make_NLP_Problem('IK')
print('== we pick "IK" as a problem:\n', nlp.report(1), '===')
data = generate_samples(nlp, 10)
display_samples(nlp, data)
if __name__ == "__main__":
main()