constrain f(x) with g(x) when minimizing f(x) in mystic?

Question

Hi I am new to mystic so I apologize if I ask similar questions that have been answered before.

Say I have an input x (a list of length n), and a function f(x) that will map the input to a m-dimensional output space. I also have a m-dimensional constraint list that also depends on the input x (say it is g(x)), and I want to make sure the output is less than the constraint for each element of the m-dimensional list. How should I specify this in mystic?

simplified example:

x = [1,2,3]
output = [6, 12] # with some f(x)
constraints = [5, 10] # with some g(x)

Answer 1

I'm the mystic author. The easiest way I'd constrain one function's output with another is by using a penalty (soft constraint).

I'll show a simple case, similar to what I think you are looking for:

"""
2-D input x
1-D output y = f(x)
where y > g(x)

with f(x) = x0^(sin(x0)) + x1^4 + 6x1^3 - 5x1^2 - 40x1 + 35
and g(x) = 7x1 - x0 + 5
in x = [0,10]
"""

Let's find the minimum of f, larger than g.

>>> import numpy as np
>>> import mystic as my
>>> 
>>> def f(x):
...     x0,x1 = x
...     return x0**np.sin(x0) + x1**4 + 6*x1**3 - 5*x1**2 - 40*x1 + 35
... 
>>> def g(x):
...     x0,x1 = x
...     return 7*x1 - x0 + 5
... 
>>> def penalty(x):
...     return g(x) - f(x)
... 
>>> @my.penalty.quadratic_inequality(penalty, k=1e12)
... def p(x):
...     return 0.0
... 
>>> mon = my.monitors.Monitor()
>>> my.solvers.fmin(f, [5,5], bounds=[(0,10)]*2, penalty=p, itermon=mon, disp=1)
Optimization terminated successfully.
         Current function value: 11.898373
         Iterations: 89
         Function evaluations: 175
STOP("CandidateRelativeTolerance with {'xtol': 0.0001, 'ftol': 0.0001}")
array([7.81765653, 2.10228969])
>>> 
>>> my.log_reader(mon)

Make a plot to visualize the constrained solver trajectory.

>>> fig = my.model_plotter(f, mon, depth=True, scale=1.0, bounds="0:10, 0:10", out=True)
>>> 
>>> from matplotlib import cm
>>> import matplotlib.pyplot as plt
>>> 
>>> x,y = my.scripts._parse_axes("0:10, 0:10", grid=True)
>>> x, y = np.meshgrid(x, y)
>>> z = 0*x
>>> s,t = x.shape
>>> for i in range(s):
...     for j in range(t):
...         xx,yy = x[i,j], y[i,j]
...         z[i,j] = g([xx,yy])
... 
>>> z = np.log(4*z*1.0+1)+2 # scale=1.0
>>> ax = fig.axes[0]
>>> ax.contour(x, y, z, 50, cmap=cm.cool)
<matplotlib.contour.QuadContourSet object at 0x12bc0d0d0>
>>> plt.show()

From the above plot, you can sort of see that the solver started minimizing f (color=jet), then hit g (color=cool), and traced along the intersection until it hit the minimum.