How do I create a constraint using a linear symbolic inequality in Mystic?

Question

I am trying to use Mystic to minimize a nonlinear function with linear constraints.

As a simple example, I have the following:

import numpy as np
import mystic.symbolic as ms
from mystic.symbolic import generate_constraint
from mystic.symbolic import generate_solvers
from mystic.symbolic import linear_symbolic
from mystic.monitors import Monitor
from mystic.solvers import LatticeSolver
from mystic.solvers import NelderMeadSimplexSolver
from mystic.termination import CandidateRelativeTolerance as CRT

# diamond-shaped constraint
# same format as output of mystic.linear_symbolic()
basic_constraint = '''
1.0*x0 + 1.0*x1 <= 5
1.0*x0 - 1.0*x1 >= -5
1.0*x0 + 1.0*x1 >= -5
1.0*x0 - 1.0*x1 <= 5
'''[1:]

def basic_objective(x, *args):
    v1 = x[0] * x[1] / (1 + np.abs(x[0] + x[1]))
    v2 = np.min(x)
    return v1 + v2/(1+np.abs(v1))

When trying to run the code, I do the following:

def test_basic():
    stepmon=Monitor()
    nbins = [6,6,]
    solver = LatticeSolver(len(nbins), nbins)
    solver.SetNestedSolver(NelderMeadSimplexSolver)
    print('Generating Solvers')
    constraint_solver = generate_solvers(
        basic_constraint,
        nvars=2
    )
    print(constraint_solver)
    # HERE IS ISSUE, IF COMMENTED ISSUE BELOW
    print(constraint_solver[0](np.ones(2)))
    print('Setting Constraints')
    solver.SetConstraints(
        generate_constraint(constraint_solver)
    )
    solver.SetGenerationMonitor(stepmon)
    solver.SetTermination(CRT())
    print('Solving...')
    # ISSUE APPEARS HERE IF print(constraint_solver[0]...)
    # IS COMMENTED OUT
    solver.Solve(basic_objective)
    solution = solver.Solution()
    print(solution)
    return solution

test_basic()

When I run the above, the error occurs at

print(constraint_solver[0](np.ones(2)))

or, if I comment it out,

solver.Solve(basic_objective)

The only noticeable difference is the size of the call stack.

The error I get is

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 12, in test_basic
  File "<string>", line 4, in solver_139632515562208
  File "<string>", line 1
SyntaxError: cannot assign to operator

This is a result of Mystic trying to compile Python code from a string and encountering a syntax error, but I do not know how to fix this issue.

Answer 1

I'm the mystic author. You are missing one key function, and a second that is not needed in this case but often is.

If you print the doc for your constraint solvers, you'll see that they are not formed well.

>>> constraint_solver = generate_solvers(basic_constraint, nvars=2)
>>> print(constraint_solver[0].__doc__)
1.0*x[0] - 1.0*x[1] = min(5 - (_tol(5,tol,rel) * any(equal(5,[]))), 1.0*x[0] - 1.0*x[1])
>>> 

You need to isolate a single variable on the left-hand side. Hence, we either need to solve or simplify. For inequalities, simplify works better, and for equalities solve generally works. I am not sure the level of documentation that states that. Anyway, I use simplify before building the constraints.

>>> from mystic.symbolic import simplify
>>> constraint_solver = generate_solvers(simplify(basic_constraint), nvars=2)
>>> print(constraint_solver[0].__doc__)
x[0] = max(1.0*x[1] - 5.0 + (_tol(1.0*x[1] - 5.0,tol,rel) * any(equal(1.0*x[1] - 5.0,[]))), x[0])
>>> 
>>> print(constraint_solver[0](np.ones(2)))
[1. 1.]
>>> 

Now, your code works as expected.

However, I'd generally make one other modification.

>>> from mystic.constraints import and_
>>> c = generate_constraint(constraint_solver, join=and_)
>>> c(np.ones(2)*5)
[0.0, 5.0]
>>> print(c.__doc__)
inner: x[0] = max(1.0*x[1] - 5.0 + (_tol(1.0*x[1] - 5.0,tol,rel) * any(equal(1.0*x[1] - 5.0,[]))), x[0])
inner: x[0] = min(1.0*x[1] + 5.0 - (_tol(1.0*x[1] + 5.0,tol,rel) * any(equal(1.0*x[1] + 5.0,[]))), x[0])
inner: x[0] = min(5.0 - 1.0*x[1] - (_tol(5.0 - 1.0*x[1],tol,rel) * any(equal(5.0 - 1.0*x[1],[]))), x[0])
inner: x[0] = max(-1.0*x[1] - 5.0 + (_tol(-1.0*x[1] - 5.0,tol,rel) * any(equal(-1.0*x[1] - 5.0,[]))), x[0])

Without the join=and_, your code still works. The difference is that without an explicit join statement, it's assumed the constraints are independent of each other, and can be solved one at a time. Using join=and_ forces the constraints to be solved simultaneously, which is slower. There's also or_ and other more complex combinations in building constraints, but the default is to assume independence.

Both points are subtle, and, I believe, in the documentation it should state that the constraints solvers need the symbolic equations need to have a single variable isolated on the left-hand side. However, it's probably not obvious as that's often missed.