289 lines
14 KiB
Python
289 lines
14 KiB
Python
'''
|
|
This module contains subroutines concerning the update of the interpolation set.
|
|
|
|
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.
|
|
'''
|
|
|
|
from ..common.consts import DEBUGGING
|
|
from ..common.infos import DAMAGING_ROUNDING, INFO_DEFAULT
|
|
from ..common.linalg import isinv, matprod, outprod, inprod, inv, primasum
|
|
import numpy as np
|
|
|
|
|
|
def updatexfc(jdrop, constr, cpen, cstrv, d, f, conmat, cval, fval, sim, simi):
|
|
'''
|
|
This function revises the simplex by updating the elements of SIM, SIMI, FVAL, CONMAT, and CVAL
|
|
'''
|
|
|
|
# Local variables
|
|
itol = 1
|
|
|
|
# Sizes
|
|
num_constraints = np.size(constr)
|
|
num_vars = np.size(sim, 0)
|
|
|
|
# Preconditions
|
|
if DEBUGGING:
|
|
assert num_constraints >= 0
|
|
assert num_vars >= 1
|
|
assert jdrop >= 0 and jdrop <= num_vars + 1
|
|
assert not any(np.isnan(constr) | np.isneginf(constr))
|
|
assert not (np.isnan(cstrv) | np.isposinf(cstrv))
|
|
assert np.size(d) == num_vars and all(np.isfinite(d))
|
|
assert not (np.isnan(f) | np.isposinf(f))
|
|
assert np.size(conmat, 0) == num_constraints and np.size(conmat, 1) == num_vars + 1
|
|
assert not (np.isnan(conmat) | np.isneginf(conmat)).any()
|
|
assert np.size(cval) == num_vars + 1 and not any(cval < 0 | np.isnan(cval) | np.isposinf(cval))
|
|
assert np.size(fval) == num_vars + 1 and not any(np.isnan(fval) | np.isposinf(fval))
|
|
assert np.size(sim, 0) == num_vars and np.size(sim, 1) == num_vars + 1
|
|
assert np.isfinite(sim).all()
|
|
assert all(primasum(abs(sim[:, :num_vars]), axis=0) > 0)
|
|
assert np.size(simi, 0) == num_vars and np.size(simi, 1) == num_vars
|
|
assert np.isfinite(simi).all()
|
|
assert isinv(sim[:, :num_vars], simi, itol)
|
|
|
|
#====================#
|
|
# Calculation starts #
|
|
#====================#
|
|
|
|
|
|
# Do nothing when JDROP is None. This can only happen after a trust-region step.
|
|
if jdrop is None: # JDROP is None is impossible if the input is correct.
|
|
return conmat, cval, fval, sim, simi, INFO_DEFAULT
|
|
|
|
sim_old = sim
|
|
simi_old = simi
|
|
if jdrop < num_vars:
|
|
sim[:, jdrop] = d
|
|
simi_jdrop = simi[jdrop, :] / inprod(simi[jdrop, :], d)
|
|
simi -= outprod(matprod(simi, d), simi_jdrop)
|
|
simi[jdrop, :] = simi_jdrop
|
|
else: # jdrop == num_vars
|
|
sim[:, num_vars] += d
|
|
sim[:, :num_vars] -= np.tile(d, (num_vars, 1)).T
|
|
simid = matprod(simi, d)
|
|
sum_simi = primasum(simi, axis=0)
|
|
simi += outprod(simid, sum_simi / (1 - sum(simid)))
|
|
|
|
# Check whether SIMI is a poor approximation to the inverse of SIM[:, :NUM_VARS]
|
|
# Calculate SIMI from scratch if the current one is damaged by rounding errors.
|
|
itol = 1
|
|
erri = np.max(abs(matprod(simi, sim[:, :num_vars]) - np.eye(num_vars))) # np.max returns NaN if any input is NaN
|
|
if erri > 0.1 * itol or np.isnan(erri):
|
|
simi_test = inv(sim[:, :num_vars])
|
|
erri_test = np.max(abs(matprod(simi_test, sim[:, :num_vars]) - np.eye(num_vars)))
|
|
if erri_test < erri or (np.isnan(erri) and not np.isnan(erri_test)):
|
|
simi = simi_test
|
|
erri = erri_test
|
|
|
|
# If SIMI is satisfactory, then update FVAL, CONMAT, CVAL, and the pole position. Otherwise restore
|
|
# SIM and SIMI, and return with INFO = DAMAGING_ROUNDING.
|
|
if erri <= itol:
|
|
fval[jdrop] = f
|
|
conmat[:, jdrop] = constr
|
|
cval[jdrop] = cstrv
|
|
# Switch the best vertex to the pole position SIM[:, NUM_VARS] if it is not there already
|
|
conmat, cval, fval, sim, simi, info = updatepole(cpen, conmat, cval, fval, sim, simi)
|
|
else:
|
|
info = DAMAGING_ROUNDING
|
|
sim = sim_old
|
|
simi = simi_old
|
|
|
|
#==================#
|
|
# Calculation ends #
|
|
#==================#
|
|
|
|
# Postconditions
|
|
if DEBUGGING:
|
|
assert np.size(conmat, 0) == num_constraints and np.size(conmat, 1) == num_vars + 1
|
|
assert not (np.isnan(conmat) | np.isneginf(conmat)).any()
|
|
assert np.size(cval) == num_vars + 1 and not any(cval < 0 | np.isnan(cval) | np.isposinf(cval))
|
|
assert np.size(fval) == num_vars + 1 and not any(np.isnan(fval) | np.isposinf(fval))
|
|
assert np.size(sim, 0) == num_vars and np.size(sim, 1) == num_vars + 1
|
|
assert np.isfinite(sim).all()
|
|
assert all(primasum(abs(sim[:, :num_vars]), axis=0) > 0)
|
|
assert np.size(simi, 0) == num_vars and np.size(simi, 1) == num_vars
|
|
assert np.isfinite(simi).all()
|
|
assert isinv(sim[:, :num_vars], simi, itol) or info == DAMAGING_ROUNDING
|
|
|
|
return sim, simi, fval, conmat, cval, info
|
|
|
|
def findpole(cpen, cval, fval):
|
|
'''
|
|
This subroutine identifies the best vertex of the current simplex with respect to the merit
|
|
function PHI = F + CPEN * CSTRV.
|
|
'''
|
|
|
|
# Size
|
|
num_vars = np.size(fval) - 1
|
|
|
|
# Preconditions
|
|
if DEBUGGING:
|
|
assert cpen > 0
|
|
assert np.size(cval) == num_vars + 1 and not any(cval < 0 | np.isnan(cval) | np.isposinf(cval))
|
|
assert np.size(fval) == num_vars + 1 and not any(np.isnan(fval) | np.isposinf(fval))
|
|
|
|
#====================#
|
|
# Calculation starts #
|
|
#====================#
|
|
|
|
# Identify the optimal vertex of the current simplex
|
|
jopt = np.size(fval) - 1
|
|
phi = fval + cpen * cval
|
|
phimin = min(phi)
|
|
# Essentially jopt = np.argmin(phi). However, we keep jopt = num_vars unless there
|
|
# is a strictly better choice. When there are multiple choices, we choose the jopt
|
|
# with the smallest value of cval.
|
|
if phimin < phi[jopt] or any((cval < cval[jopt]) & (phi <= phi[jopt])):
|
|
# While we could use argmin(phi), there may be two places where phi achieves
|
|
# phimin, and in that case we should choose the one with the smallest cval.
|
|
jopt = np.ma.array(cval, mask=(phi > phimin)).argmin()
|
|
|
|
#==================#
|
|
# Calculation ends #
|
|
#==================#
|
|
|
|
# Postconditions
|
|
if DEBUGGING:
|
|
assert jopt >= 0 and jopt < num_vars + 1
|
|
assert jopt == num_vars or phi[jopt] < phi[num_vars] or (phi[jopt] <= phi[num_vars] and cval[jopt] < cval[num_vars])
|
|
return jopt
|
|
|
|
|
|
def updatepole(cpen, conmat, cval, fval, sim, simi):
|
|
#--------------------------------------------------------------------------------------------------!
|
|
# This subroutine identifies the best vertex of the current simplex with respect to the merit
|
|
# function PHI = F + CPEN * CSTRV, and then switch this vertex to SIM[:, NUM_VARS], which Powell called
|
|
# the "pole position" in his comments. CONMAT, CVAL, FVAL, and SIMI are updated accordingly.
|
|
#
|
|
# N.B. 1: In precise arithmetic, the following two procedures produce the same results:
|
|
# 1) apply UPDATEPOLE to SIM twice, first with CPEN = CPEN1 and then with CPEN = CPEN2;
|
|
# 2) apply UPDATEPOLE to SIM with CPEN = CPEN2.
|
|
# In finite-precision arithmetic, however, they may produce different results unless CPEN1 = CPEN2.
|
|
#
|
|
# N.B. 2: When JOPT == N+1, the best vertex is already at the pole position, so there is nothing to
|
|
# switch. However, as in Powell's code, the code below will check whether SIMI is good enough to
|
|
# work as the inverse of SIM(:, 1:N) or not. If not, Powell's code would invoke an error return of
|
|
# COBYLB; our implementation, however, will try calculating SIMI from scratch; if the recalculated
|
|
# SIMI is still of poor quality, then UPDATEPOLE will return with INFO = DAMAGING_ROUNDING,
|
|
# informing COBYLB that SIMI is poor due to damaging rounding errors.
|
|
#
|
|
# N.B. 3: UPDATEPOLE should be called when and only when FINDPOLE can potentially returns a value
|
|
# other than N+1. The value of FINDPOLE is determined by CPEN, CVAL, and FVAL, the latter two being
|
|
# decided by SIM. Thus UPDATEPOLE should be called after CPEN or SIM changes. COBYLA updates CPEN at
|
|
# only two places: the beginning of each trust-region iteration, and when REDRHO is called;
|
|
# SIM is updated only by UPDATEXFC, which itself calls UPDATEPOLE internally. Therefore, we only
|
|
# need to call UPDATEPOLE after updating CPEN at the beginning of each trust-region iteration and
|
|
# after each invocation of REDRHO.
|
|
|
|
# Local variables
|
|
itol = 1
|
|
|
|
# Sizes
|
|
num_constraints = conmat.shape[0]
|
|
num_vars = sim.shape[0]
|
|
|
|
# Preconditions
|
|
if DEBUGGING:
|
|
assert num_constraints >= 0
|
|
assert num_vars >= 1
|
|
assert cpen > 0
|
|
assert np.size(conmat, 0) == num_constraints and np.size(conmat, 1) == num_vars + 1
|
|
assert not (np.isnan(conmat) | np.isneginf(conmat)).any()
|
|
assert np.size(cval) == num_vars + 1 and not any(cval < 0 | np.isnan(cval) | np.isposinf(cval))
|
|
assert np.size(fval) == num_vars + 1 and not any(np.isnan(fval) | np.isposinf(fval))
|
|
assert np.size(sim, 0) == num_vars and np.size(sim, 1) == num_vars + 1
|
|
assert np.isfinite(sim).all()
|
|
assert all(primasum(abs(sim[:, :num_vars]), axis=0) > 0)
|
|
assert np.size(simi, 0) == num_vars and np.size(simi, 1) == num_vars
|
|
assert np.isfinite(simi).all()
|
|
assert isinv(sim[:, :num_vars], simi, itol)
|
|
|
|
#====================#
|
|
# Calculation starts #
|
|
#====================#
|
|
|
|
# INFO must be set, as it is an output.
|
|
info = INFO_DEFAULT
|
|
|
|
# Identify the optimal vertex of the current simplex.
|
|
jopt = findpole(cpen, cval, fval)
|
|
|
|
# Switch the best vertex to the pole position SIM[:, NUM_VARS] if it is not there already and update
|
|
# SIMI. Before the update, save a copy of SIM and SIMI. If the update is unsuccessful due to
|
|
# damaging rounding errors, we restore them and return with INFO = DAMAGING_ROUNDING.
|
|
sim_old = sim.copy()
|
|
simi_old = simi.copy()
|
|
if 0 <= jopt < num_vars:
|
|
# Unless there is a bug in FINDPOLE it is guaranteed that JOPT >= 0
|
|
# When JOPT == NUM_VARS, there is nothing to switch; in addition SIMI[JOPT, :] will be illegal.
|
|
# fval[[jopt, -1]] = fval[[-1, jopt]]
|
|
# conmat[:, [jopt, -1]] = conmat[:, [-1, jopt]] # Exchange CONMAT[:, JOPT] AND CONMAT[:, -1]
|
|
# cval[[jopt, -1]] = cval[[-1, jopt]]
|
|
sim[:, num_vars] += sim[:, jopt]
|
|
sim_jopt = sim[:, jopt].copy()
|
|
sim[:, jopt] = 0 # np.zeros(num_constraints)?
|
|
sim[:, :num_vars] -= np.tile(sim_jopt, (num_vars, 1)).T
|
|
# The above update is equivalent to multiplying SIM[:, :NUM_VARS] from the right side by a matrix whose
|
|
# JOPT-th row is [-1, -1, ..., -1], while all the other rows are the same as those of the
|
|
# identity matrix. It is easy to check that the inverse of this matrix is itself. Therefore,
|
|
# SIMI should be updated by a multiplication with this matrix (i.e. its inverse) from the left
|
|
# side, as is done in the following line. The JOPT-th row of the updated SIMI is minus the sum
|
|
# of all rows of the original SIMI, whereas all the other rows remain unchanged.
|
|
# NDB 20250114: In testing the cutest problem 'SYNTHES2' between the Python implementation and
|
|
# the Fortran bindings, I saw a difference between the following for loop and the
|
|
# np.sum command. The differences were small, on the order of 1e-16, i.e. epsilon.
|
|
# According to numpy documentation, np.sum sometimes uses partial pairwise summation,
|
|
# depending on the memory layout of the array and the axis specified.
|
|
# for i in range(simi.shape[1]):
|
|
# simi[jopt, i] = -sum(simi[:, i])
|
|
simi[jopt, :] = -primasum(simi, axis=0)
|
|
|
|
# Check whether SIMI is a poor approximation to the inverse of SIM[:, :NUM_VARS]
|
|
# Calculate SIMI from scratch if the current one is damaged by rounding errors.
|
|
erri = np.max(abs(matprod(simi, sim[:, :num_vars]) - np.eye(num_vars))) # np.max returns NaN if any input is NaN
|
|
itol = 1
|
|
if erri > 0.1 * itol or np.isnan(erri):
|
|
simi_test = inv(sim[:, :num_vars])
|
|
erri_test = np.max(abs(matprod(simi_test, sim[:, :num_vars]) - np.eye(num_vars)))
|
|
if erri_test < erri or (np.isnan(erri) and not np.isnan(erri_test)):
|
|
simi = simi_test
|
|
erri = erri_test
|
|
|
|
|
|
# If SIMI is satisfactory, then update FVAL, CONMAT, and CVAL. Otherwise restore SIM and SIMI, and
|
|
# return with INFO = DAMAGING_ROUNDING.
|
|
if erri <= itol:
|
|
if 0 <= jopt < num_vars:
|
|
fval[[jopt, num_vars]] = fval[[num_vars, jopt]]
|
|
conmat[:, [jopt, num_vars]] = conmat[:, [num_vars, jopt]]
|
|
cval[[jopt, num_vars]] = cval[[num_vars, jopt]]
|
|
else: # erri > itol or erri is NaN
|
|
info = DAMAGING_ROUNDING
|
|
sim = sim_old
|
|
simi = simi_old
|
|
|
|
#==================#
|
|
# Calculation ends #
|
|
#==================#
|
|
|
|
# Postconditions
|
|
if DEBUGGING:
|
|
assert findpole(cpen, cval, fval) == num_vars or info == DAMAGING_ROUNDING
|
|
assert np.size(conmat, 0) == num_constraints and np.size(conmat, 1) == num_vars + 1
|
|
assert not (np.isnan(conmat) | np.isneginf(conmat)).any()
|
|
assert np.size(cval) == num_vars + 1 and not any(cval < 0 | np.isnan(cval) | np.isposinf(cval))
|
|
assert np.size(fval) == num_vars + 1 and not any(np.isnan(fval) | np.isposinf(fval))
|
|
assert np.size(sim, 0) == num_vars and np.size(sim, 1) == num_vars + 1
|
|
assert np.isfinite(sim).all()
|
|
assert all(primasum(abs(sim[:, :num_vars]), axis=0) > 0)
|
|
assert np.size(simi, 0) == num_vars and np.size(simi, 1) == num_vars
|
|
assert np.isfinite(simi).all()
|
|
# Do not check SIMI = SIM[:, :num_vars]^{-1}, as it may not be true due to damaging rounding.
|
|
assert isinv(sim[:, :num_vars], simi, itol) or info == DAMAGING_ROUNDING
|
|
|
|
return conmat, cval, fval, sim, simi, info
|