image-match/python/py/Lib/site-packages/scipy/_lib/pyprima/common/selectx.py

'''
This module provides subroutines that ensure the returned X is optimal among all the calculated
points in the sense that no other point achieves both lower function value and lower constraint
violation at the same time. This module is needed only in the constrained case.

Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA.

Dedicated to late Professor M. J. D. Powell FRS (1936--2015).

Python translation by Nickolai Belakovski.
'''

import numpy as np
import numpy.typing as npt
from .consts import DEBUGGING, EPS, CONSTRMAX, REALMAX, FUNCMAX
from .present import present

def isbetter(f1: float, c1: float, f2: float, c2: float, ctol: float) -> bool:
    '''
    This function compares whether FC1 = (F1, C1) is (strictly) better than FC2 = (F2, C2), which
    basically means that (F1 < F2 and C1 <= C2) or (F1 <= F2 and C1 < C2).
    It takes care of the cases where some of these values are NaN or Inf, even though some cases
    should never happen due to the moderated extreme barrier.
    At return, BETTER = TRUE if and only if (F1, C1) is better than (F2, C2).
    Here, C means constraint violation, which is a nonnegative number.
    '''

    # Preconditions
    if DEBUGGING:
        assert not any(np.isnan([f1, c1]) | np.isposinf([f2, c2]))
        assert not any(np.isnan([f2, c2]) | np.isposinf([f2, c2]))
        assert c1 >= 0 and c2 >= 0
        assert ctol >= 0

    #====================#
    # Calculation starts #
    #====================#

    is_better = False
    # Even though NaN/+Inf should not occur in FC1 or FC2 due to the moderated extreme barrier, for
    # security and robustness, the code below does not make this assumption.
    is_better = is_better or (any(np.isnan([f1, c1]) | np.isposinf([f1, c1])) and not any(np.isnan([f2, c2]) | np.isposinf([f2, c2])))

    is_better = is_better or (f1 < f2 and c1 <= c2)
    is_better = is_better or (f1 <= f2 and c1 < c2)

    # If C1 <= CTOL and C2 is significantly larger/worse than CTOL, i.e., C2 > MAX(CTOL,CREF),
    # then FC1 is better than FC2 as long as F1 < REALMAX. Normally CREF >= CTOL so MAX(CTOL, CREF)
    # is indeed CREF. However, this may not be true if CTOL > 1E-1*CONSTRMAX.
    cref = 10 * max(EPS, min(ctol, 1.0E-2 * CONSTRMAX))  # The MIN avoids overflow.
    is_better = is_better or (f1 < REALMAX and c1 <= ctol and (c2 > max(ctol, cref) or np.isnan(c2)))

    #==================#
    # Calculation ends #
    #==================#

    # Postconditions
    if DEBUGGING:
        assert not (is_better and f1 >= f2 and c1 >= c2)
        assert is_better or not (f1 <= f2 and c1 < c2)
        assert is_better or not (f1 < f2 and c1 <= c2)

    return is_better


def savefilt(cstrv, ctol, cweight, f, x, nfilt, cfilt, ffilt, xfilt, constr=None, confilt=None):
    '''
    This subroutine saves X, F, and CSTRV in XFILT, FFILT, and CFILT (and CONSTR in CONFILT
    if they are present), unless a vector in XFILT[:, :NFILT] is better than X.
    If X is better than some vectors in XFILT[:, :NFILT] then these vectors will be
    removed. If X is not better than any of XFILT[:, :NFILT], but NFILT == MAXFILT,
    then we remove a column from XFILT according to the merit function
    PHI = FFILT + CWEIGHT * max(CFILT - CTOL, 0)
    N.B.:
    1. Only XFILT[:, :NFILT] and FFILT[:, :NFILT] etc contains valid information,
    while XFILT[:, NFILT+1:MAXFILT] and FFILT[:, NFILT+1:MAXFILT] etc are not
    initialized yet.
    2. We decide whether and X is better than another by the ISBETTER function
    '''

    # Sizes
    if present(constr):
        num_constraints = len(constr)
    else:
        num_constraints = 0
    num_vars = len(x)
    maxfilt = len(ffilt)

    # Preconditions
    if DEBUGGING:
        # Check the size of X.
        assert num_vars >= 1
        # Check CWEIGHT and CTOL
        assert cweight >= 0
        assert ctol >= 0
        # Check NFILT
        assert nfilt >= 0 and nfilt <= maxfilt
        # Check the sizes of XFILT, FFILT, CFILT.
        assert maxfilt >= 1
        assert np.size(xfilt, 0) == num_vars and np.size(xfilt, 1) == maxfilt
        assert np.size(cfilt) == maxfilt
        # Check the values of XFILT, FFILT, CFILT.
        assert not (np.isnan(xfilt[:, :nfilt])).any()
        assert not any(np.isnan(ffilt[:nfilt]) | np.isposinf(ffilt[:nfilt]))
        assert not any(cfilt[:nfilt] < 0 | np.isnan(cfilt[:nfilt]) | np.isposinf(cfilt[:nfilt]))
        # Check the values of X, F, CSTRV.
        # X does not contain NaN if X0 does not and the trust-region/geometry steps are proper.
        assert not any(np.isnan(x))
        # F cannot be NaN/+Inf due to the moderated extreme barrier.
        assert not (np.isnan(f) | np.isposinf(f))
        # CSTRV cannot be NaN/+Inf due to the moderated extreme barrier.
        assert not (cstrv < 0 | np.isnan(cstrv) | np.isposinf(cstrv))
        # Check CONSTR and CONFILT.
        assert present(constr) == present(confilt)
        if present(constr):
            # CONSTR cannot contain NaN/-Inf due to the moderated extreme barrier.
            assert not any(np.isnan(constr) | np.isneginf(constr))
            assert np.size(confilt, 0) == num_constraints and np.size(confilt, 1) == maxfilt
            assert not (np.isnan(confilt[:, :nfilt]) | np.isneginf(confilt[:, :nfilt])).any()

    #====================#
    # Calculation starts #
    #====================#

    # Return immediately if any column of XFILT is better than X.
    if any((isbetter(ffilt_i, cfilt_i, f, cstrv, ctol) for ffilt_i, cfilt_i in zip(ffilt[:nfilt], cfilt[:nfilt]))) or \
        any(np.logical_and(ffilt[:nfilt] <= f, cfilt[:nfilt] <= cstrv)):
        return nfilt, cfilt, ffilt, xfilt, confilt

    # Decide which columns of XFILT to keep.
    keep = np.logical_not([isbetter(f, cstrv, ffilt_i, cfilt_i, ctol) for ffilt_i, cfilt_i in zip(ffilt[:nfilt], cfilt[:nfilt])])

    # If NFILT == MAXFILT and X is not better than any column of XFILT, then we remove the worst column
    # of XFILT according to the merit function PHI = FFILT + CWEIGHT * MAX(CFILT - CTOL, ZERO).
    if sum(keep) == maxfilt:  # In this case, NFILT = SIZE(KEEP) = COUNT(KEEP) = MAXFILT > 0.
        cfilt_shifted = np.maximum(cfilt - ctol, 0)
        if cweight <= 0:
            phi = ffilt
        elif np.isposinf(cweight):
            phi = cfilt_shifted
            # We should not use CFILT here; if MAX(CFILT_SHIFTED) is attained at multiple indices, then
            # we will check FFILT to exhaust the remaining degree of freedom.
        else:
            phi = np.maximum(ffilt, -REALMAX)
            phi = np.nan_to_num(phi, nan=-REALMAX)  # Replace NaN with -REALMAX and +/- inf with large numbers
            phi += cweight * cfilt_shifted
        # We select X to maximize PHI. In case there are multiple maximizers, we take the one with the
        # largest CSTRV_SHIFTED; if there are more than one choices, we take the one with the largest F;
        # if there are several candidates, we take the one with the largest CSTRV; if the last comparison
        # still leads to more than one possibilities, then they are equally bad and we choose the first.
        # N.B.:
        # 1. This process is the opposite of selecting KOPT in SELECTX.
        # 2. In finite-precision arithmetic, PHI_1 == PHI_2 and CSTRV_SHIFTED_1 == CSTRV_SHIFTED_2 do
        # not ensure that F_1 == F_2!
        phimax = max(phi)
        cref = max(cfilt_shifted[phi >= phimax])
        fref = max(ffilt[cfilt_shifted >= cref])
        kworst = np.ma.array(cfilt, mask=(ffilt > fref)).argmax()
        if kworst < 0 or kworst >= len(keep):  #  For security. Should not happen.
            kworst = 0
        keep[kworst] = False

    # Keep the good xfilt values and remove all the ones that are strictly worse than the new x.
    nfilt = sum(keep)
    index_to_keep = np.where(keep)[0]
    xfilt[:, :nfilt] = xfilt[:, index_to_keep]
    ffilt[:nfilt] = ffilt[index_to_keep]
    cfilt[:nfilt] = cfilt[index_to_keep]
    if confilt is not None and constr is not None:
        confilt[:, :nfilt] = confilt[:, index_to_keep]

    # Once we have removed all the vectors that are strictly worse than x,
    # we add x to the filter.
    xfilt[:, nfilt] = x
    ffilt[nfilt] = f
    cfilt[nfilt] = cstrv
    if confilt is not None and constr is not None:
        confilt[:, nfilt] = constr
    nfilt += 1  # In Python we need to increment the index afterwards

    #==================#
    # Calculation ends #
    #==================#

    # Postconditions
    if DEBUGGING:
        # Check NFILT and the sizes of XFILT, FFILT, CFILT.
        assert nfilt >= 1 and nfilt <= maxfilt
        assert np.size(xfilt, 0) == num_vars and np.size(xfilt, 1) == maxfilt
        assert np.size(ffilt) == maxfilt
        assert np.size(cfilt) == maxfilt
        # Check the values of XFILT, FFILT, CFILT.
        assert not (np.isnan(xfilt[:, :nfilt])).any()
        assert not any(np.isnan(ffilt[:nfilt]) | np.isposinf(ffilt[:nfilt]))
        assert not any(cfilt[:nfilt] < 0 | np.isnan(cfilt[:nfilt]) | np.isposinf(cfilt[:nfilt]))
        # Check that no point in the filter is better than X, and X is better than no point.
        assert not any([isbetter(ffilt_i, cfilt_i, f, cstrv, ctol) for ffilt_i, cfilt_i in zip(ffilt[:nfilt], cfilt[:nfilt])])
        assert not any([isbetter(f, cstrv, ffilt_i, cfilt_i, ctol) for ffilt_i, cfilt_i in zip(ffilt[:nfilt], cfilt[:nfilt])])
        # Check CONFILT.
        if present(confilt):
            assert np.size(confilt, 0) == num_constraints and np.size(confilt, 1) == maxfilt
            assert not (np.isnan(confilt[:, :nfilt]) | np.isneginf(confilt[:, :nfilt])).any()


    return nfilt, cfilt, ffilt, xfilt, confilt


def selectx(fhist: npt.NDArray, chist: npt.NDArray, cweight: float, ctol: float):
    '''
    This subroutine selects X according to FHIST and CHIST, which represents (a part of) history
    of F and CSTRV. Normally, FHIST and CHIST are not the full history but only a filter, e.g. ffilt
    and CFILT generated by SAVEFILT. However, we name them as FHIST and CHIST because the [F, CSTRV]
    in a filter should not dominate each other, but this subroutine does NOT assume such a property.
    N.B.: CTOL is the tolerance of the constraint violation (CSTRV). A point is considered feasible if
    its constraint violation is at most CTOL. Not that CTOL is absolute, not relative.
    '''

    # Sizes
    nhist = len(fhist)

    # Preconditions
    if DEBUGGING:
        assert nhist >= 1
        assert np.size(chist) == nhist
        assert not any(np.isnan(fhist) | np.isposinf(fhist))
        assert not any(chist < 0 | np.isnan(chist) | np.isposinf(chist))
        assert cweight >= 0
        assert ctol >= 0

    #====================#
    # Calculation starts #
    #====================#

    # We select X among the points with F < FREF and CSTRV < CREF.
    # Do NOT use F <= FREF, because F == FREF (FUNCMAX or REALMAX) may mean F == INF in practice!
    if any(np.logical_and(fhist < FUNCMAX, chist < CONSTRMAX)):
        fref = FUNCMAX
        cref = CONSTRMAX
    elif any(np.logical_and(fhist < REALMAX, chist < CONSTRMAX)):
        fref = REALMAX
        cref = CONSTRMAX
    elif any(np.logical_and(fhist < FUNCMAX, chist < REALMAX)):
        fref = FUNCMAX
        cref = REALMAX
    else:
        fref = REALMAX
        cref = REALMAX

    if not any(np.logical_and(fhist < fref, chist < cref)):
        kopt = nhist - 1
    else:
        # Shift the constraint violations by ctol, so that cstrv <= ctol is regarded as no violation.
        chist_shifted = np.maximum(chist - ctol, 0)
        # cmin is the minimal shift constraint violation attained in the history.
        cmin = np.min(chist_shifted[fhist < fref])
        # We consider only the points whose shifted constraint violations are at most the cref below.
        # N.B.: Without taking np.maximum(EPS, .), cref would be 0 if cmin = 0. In that case, asking for
        # cstrv_shift < cref would be WRONG!
        cref = np.maximum(EPS, 2*cmin)
        # We use the following phi as our merit function to select X.
        if cweight <= 0:
            phi = fhist
        elif np.isposinf(cweight):
            phi = chist_shifted
            # We should not use chist here; if np.minimum(chist_shifted) is attained at multiple indices, then
            # we will check fhist to exhaust the remaining degree of freedom.
        else:
            phi = np.maximum(fhist, -REALMAX) + cweight * chist_shifted
            # np.maximum(fhist, -REALMAX) makes sure that phi will not contain NaN (unless there is a bug).

        # We select X to minimize phi subject to f < fref and cstrv_shift <= cref (see the comments
        # above for the reason of taking "<" and "<=" in these two constraints). In case there are
        # multiple minimizers, we take the one with the least cstrv_shift; if there is more than one
        # choice, we take the one with the least f; if there are several candidates, we take the one
        # with the least cstrv; if the last comparison still leads to more than one possibility, then
        # they are equally good and we choose the first.
        # N.B.:
        # 1. This process is the opposite of selecting kworst in savefilt
        # 2. In finite-precision arithmetic, phi_2 == phi_2 and cstrv_shift_1 == cstrv_shifted_2 do
        # not ensure thatn f_1 == f_2!
        phimin = np.min(phi[np.logical_and(fhist < fref, chist_shifted <= cref)])
        cref = np.min(chist_shifted[np.logical_and(fhist < fref, phi <= phimin)])
        fref = np.min(fhist[chist_shifted <= cref])
        # Can't use argmin here because using it with a mask throws off the index
        kopt = np.ma.array(chist, mask=(fhist > fref)).argmin()

    #==================#
    # Calculation ends #
    #==================#

    # Postconditions
    if DEBUGGING:
        assert kopt >= 0 and kopt < nhist
        assert not any([isbetter(fhisti, chisti, fhist[kopt], chist[kopt], ctol) for fhisti, chisti in zip(fhist[:nhist], chist[:nhist])])

    return kopt