Make pyprima update constraints when eliminating fixed variables and bring the code to python bindings

[nbelakovski]

This PR covers both #249 and #250. The first commit is for #249 and the second is for #250.

[zaikunzhang]

@nbelakovski Thank you.

@OptHuang Please review this.

The code is copied from COBYQA. Please compare with the corresponding code in PDFO. Any difference is suspicious and hence needs justification.

Thank you.

[OptHuang]

Two typos in comments.

[zaikunzhang]

The code is copied from COBYQA. Please compare with the corresponding code in PDFO. Any difference is suspicious and hence needs justification.

@OptHuang

[nbelakovski]

@OptHuang do you have any further comments?

[OptHuang]

@OptHuang do you have any further comments?

It's good so far. But I haven't found time to compare the corresponding code with that in pdfo.

[nbelakovski]

@OptHuang Do you think you'll get time to do so in the coming week?

[OptHuang]

@OptHuang Do you think you'll get time to do so in the coming week?

Hi, Nickolai @nbelakovski . Sorry for waiting. I checked the code in pdfo about preprocessing fixed variables. I have a few suggestions or questions, which may be wrong since I am not too familiar with the code in prima.

1. When dealing with bounds, is there code in prima that will turn nan to inf or -inf?
2. Is there code in prima that will deal with infeasible bounds?
3. pdfo specially handles the case where all the variables are fixed (directly returns the input x0).
4. There is a dictionary in pdfo called prob_info that will store the information of solving the problem, which also records the information about fixed variables.

That's all my observations and suggestions.

Thanks,
Cunxin

[nbelakovski]

Hi @OptHuang, thank you for taking a look. Below are my responses.

1. When dealing with bounds, is there code in prima that will turn nan to inf or -inf?

Yes, see the relevant code in cobyla.f90 and cobyla.py

2. Is there code in prima that will deal with infeasible bounds?

No. I propose adding the following code to the end of the process_bounds in _bounds.py to address this:

    # Check the infeasibility of the bounds.
    infeasible = np.isposinf(lb) | np.isneginf(ub) | (lb > ub)
    assert not np.any(infeasible), f"Some of the provided bounds are infeasible. {infeasible=} {lb[infeasible]=}, {ub[infeasible]=}"

This code is slightly different from what's in PDFO but I think it's clearer. PDFO uses np.logical_and(lb > 0, np.isinf(lb)) to detect positive infinity in the lower bound but obviously np.isposinf is clearer.

3. pdfo specially handles the case where all the variables are fixed (directly returns the input x0).

At the moment we only handle this in the special case where x0 is a scalar. Here's one proposed way we can handle this for the general case, by inserting the following code just before if any(_fixed_idx): in __init__.py

if all(_fixed_idx):
        result = PRIMAResult()
        result.x = 0.5 * (
            lb[_fixed_idx] + ub[_fixed_idx]
        )
        result.x = np.clip(
            _fixed_values,
            lb[_fixed_idx],
            ub[_fixed_idx],
        )
        result.success = True
        result.status = 0
        result.message = "All variables were fixed by the provided bounds."
        result.fun = np.nan
        result.nfev = 0
        result.maxcv = np.nan
        result.nlconstr = np.nan
        result.method = "N/A"
        return result

None would make more sense for some of the variables, but since PRIMAResult is a structure coming from C++ it has some quirks that make it easier to use np.nan. If we can quickly agree on a way to handle this case within this PR I'd be happy to implement it, but if it starts to become a larger discussion I'd prefer to break it off into a separate issue/PR because this PR is delaying the fixing of a bug which is in turn delaying our ability to create a prima python package.

4. There is a dictionary in pdfo called prob_info that will store the information of solving the problem, which also records the information about fixed variables.

This would be too large of a refactor for right now. We should address the first 3 points first and then work on #252 to consolidate the Python interface code before we think of something like that.

[OptHuang]

Hi @OptHuang, thank you for taking a look. Below are my responses.

1. When dealing with bounds, is there code in prima that will turn nan to inf or -inf?

Yes, see the relevant code in cobyla.f90 and cobyla.py

2. Is there code in prima that will deal with infeasible bounds?

No. I propose adding the following code to the end of the process_bounds in _bounds.py to address this:

    # Check the infeasibility of the bounds.
    infeasible = np.isposinf(lb) | np.isneginf(ub) | (lb > ub)
    assert not np.any(infeasible), f"Some of the provided bounds are infeasible. {infeasible=} {lb[infeasible]=}, {ub[infeasible]=}"

This code is slightly different from what's in PDFO but I think it's clearer. PDFO uses np.logical_and(lb > 0, np.isinf(lb)) to detect positive infinity in the lower bound but obviously np.isposinf is clearer.

3. pdfo specially handles the case where all the variables are fixed (directly returns the input x0).

At the moment we only handle this in the special case where x0 is a scalar. Here's one proposed way we can handle this for the general case, by inserting the following code just before if any(_fixed_idx): in __init__.py

if all(_fixed_idx):
        result = PRIMAResult()
        result.x = 0.5 * (
            lb[_fixed_idx] + ub[_fixed_idx]
        )
        result.x = np.clip(
            _fixed_values,
            lb[_fixed_idx],
            ub[_fixed_idx],
        )
        result.success = True
        result.status = 0
        result.message = "All variables were fixed by the provided bounds."
        result.fun = np.nan
        result.nfev = 0
        result.maxcv = np.nan
        result.nlconstr = np.nan
        result.method = "N/A"
        return result

None would make more sense for some of the variables, but since PRIMAResult is a structure coming from C++ it has some quirks that make it easier to use np.nan. If we can quickly agree on a way to handle this case within this PR I'd be happy to implement it, but if it starts to become a larger discussion I'd prefer to break it off into a separate issue/PR because this PR is delaying the fixing of a bug which is in turn delaying our ability to create a prima python package.

4. There is a dictionary in pdfo called prob_info that will store the information of solving the problem, which also records the information about fixed variables.

This would be too large of a refactor for right now. We should address the first 3 points first and then work on #252 to consolidate the Python interface code before we think of something like that.

Hi Nickolai @nbelakovski ,

Thank you for the detailed explanation and proposed changes. I completely agree with your reasoning on all three points.

I’m happy to proceed with your proposal and will leave the final decision about merging this PR to @zaikunzhang.

Thanks again!

[nbelakovski]

@zaikunzhang Any further comments on this? I believe all the issues raised have been addressed. If you think it would simplify review to split this into two separate PRs (one for pyprima and one for the Python bindings) I could do that.