索引错误:在scipi . optimization中数组的索引太多

I'm trying to debbug some code with Scipy.Optimize. The bug comes from the constante: the optimisation works fine without it. The constante itself seems to works fine outside scipy.optimize (the variable testconst is computed normally). The code is the following:

我想用scipi . optimization去掉一些代码。缺陷来自于constante:优化没有它就可以正常工作。君士坦丁似乎在scipy之外运行得很好。优化(通常计算变量testconst)。守则如下:

from scipy.optimize import minimize
import numpy as np


def totaldist(dy):
    n = np.shape(dy)[0]
    temp = 0
    for i in range(n):
        temp += dy[i] ** 2

    return -0.5 * temp


def create_bond(dy_max):

    n = np.shape(dy_max)[0]
    bond = np.zeros((n, 2))
    for i in range(n):
        bond[i, :] = [0, dy_max[i]]

    tot = tuple([tuple(row) for row in bond])
    return tot


# def create_const(type_x, dx, gamma, P):
def create_const(dy, *args):
    arg = np.asarray(args)
    n = np.shape(dy)[0]
    dx = np.zeros((n, 2))
    bnd = np.zeros((n, 2))

    # from args to numpy array
    type_x = np.zeros(n)
    dP = 0
    delta1 = np.zeros(n)
    delta2 = np.zeros(n)
    gamma = np.zeros((n, n))

    for i in range(n):
        a, b = bndr(arg[0, i])
        delta1[i] = arg[0, i + n + 1]
        delta2[i] = arg[0, i + 2*n + 1]
        dx[i, 0] = (b - a) * dy[i]

    gamma = GammaApprox(delta1, delta2, dx[:, 1], dx[:, 0])

    d = np.dot(delta2, dx[:, 0])
    g = np.dot(dx[:, 0], gamma)
    g = np.dot(g, dx[:, 0])
    dP = float(arg[0, n])
    return d + 0.5 * g - dP


def GammaApprox(delta1, delta2, x1, x2):
    n = np.shape(delta1)[0]
    gamma = np.zeros((n, n))
    for i in range(n):
        if x2[i] == x1[i]:
            gamma[i, i] = 0
        else:
            gamma[i, i] = (delta2[i] - delta1[i]) / (x2[i] - x1[i])

    return gamma


def GetNewPoint(x1, x2, delta1, delta2, type_x, P):
    n = np.shape(delta1)[0]
    dmax = np.zeros(n)
    dy0 = np.zeros(n)

    # create the inequality data and initial points
    for i in range(n):
        a, b = bndr(type_x[i])
        if x2[i] > x1[i]:
            dmax[i] = (x2[i] - x1[i])/(b - a)
            dy0[i] = 1 / (b - a) * (x2[i] - x1[i]) / 2
        else:
            dmax[i] = (x1[i] - x2[i])/(b - a)
            dy0[i] = 1 / (b - a) * (x1[i] - x2[i]) / 2

    bond = create_bond(dmax)

    # create the args tuple
    arg = ()
    # type x
    for i in range(n):
        arg = arg + (type_x[i],)
    # dP
    arg = arg + (abs(P[0] - P[1]), )

    # delta1
    for i in range(n):
        arg = arg + (delta1[i], )
    # delta1
    for i in range(n):
        arg = arg + (delta2[i], )

    testconst = create_const(dy0, arg)

    # create the equality constraint
    con1 = {'type': 'eq', 'fun': create_const}
    cons = ([con1, ])
    solution = minimize(totaldist, dy0, args=arg, method='SLSQP', bounds=bond, constraints=cons, options={'disp': True})
    x = solution.x
    print(x)
    return x


def bndr(type_x):
    if type_x == 'normal':
        x_0 = -5
        x_f = 1.5
    if type_x == 'lognorm':
        x_0 = 0.0001
        x_f = 5
    if type_x == 'chisquare':
        x_0 = 0.0001
        x_f = (0.8 * (10 ** .5))

    return x_0, x_f


def test():
    x1 = np.array([0.0001, 0.0001, -5])
    x2 = np.array([1.6673, 0.84334, -5])
    delta1 = np.array([0, 0, 0])
    delta2 = np.array([2.44E-7, 2.41E-6, 4.07E-7])
    type_x = np.array(['lognorm', 'chisquare', 'normal'])
    P = (0, 6.54E-8)
    f = GetNewPoint(x1, x2, delta1, delta2, type_x, P)
    return f


test()

the error message is the following:

错误信息如下:

Traceback (most recent call last):
  File "D:/Anaconda Project/TestQP - Simplified/QP.py", line 134, in <module>
    test()
  File "D:/Anaconda Project/TestQP - Simplified/QP.py", line 130, in test
    f = GetNewPoint(x1, x2, delta1, delta2, type_x, P)
  File "D:/Anaconda Project/TestQP - Simplified/QP.py", line 103, in GetNewPoint
    solution = minimize(totaldist, dy0, args=arg, method='SLSQP', bounds=bond, constraints=cons, options={'disp': True})
  File "C:\Program Files\Anaconda\lib\site-packages\scipy\optimize\_minimize.py", line 458, in minimize
    constraints, callback=callback, **options)
  File "C:\Program Files\Anaconda\lib\site-packages\scipy\optimize\slsqp.py", line 311, in _minimize_slsqp
    meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['eq']]))
  File "C:\Program Files\Anaconda\lib\site-packages\scipy\optimize\slsqp.py", line 311, in <listcomp>
    meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['eq']]))
  File "D:/Anaconda Project/TestQP - Simplified/QP.py", line 40, in create_const
    a, b = bndr(arg[0, i])
IndexError: too many indices for array

I find roughly similar error in the website like: IndexError: index 1 is out of bounds for axis 0 with size 1/ForwardEuler

我在网站上发现了类似的错误:索引错误:索引1超出了axis 0的范围，大小为1/转发。

...but I failed to see it's really the same problem.

…但我没发现这其实是同一个问题。

1 个解决方案

#1

args is not passed to constraint-functions (automatically)!

args不会被传递给约束函数(自动)!

This is indicated in the docs:

这在文件中有说明:

args : tuple, optional

参数:tuple,可选的

Extra arguments passed to the objective function and its derivatives (Jacobian, Hessian).

额外的参数传递给目标函数及其导数(雅可比矩阵，黑森)。

You can see the problem easily by adding a print:

你可以很容易地看到问题，添加一个打印:

def create_const(dy, *args):
    print('args:')
    print(args)
    arg = np.asarray(args)
    ...

which will output something like:

会输出如下内容:

args:
(('lognorm', 'chisquare', 'normal', 6.54e-08, 0, 0, 0, 2.4400000000000001e-07, 2.4099999999999998e-06, 4.0699999999999998e-07),)
args:
()
ERROR...

If you remove your test (which is manually passing args; which works) testconst = create_const(dy0, arg), you will see only the non-working output:

如果您删除您的测试(手动传递args;testconst = create_const(dy0, arg)，您只会看到非工作输出:

args:
()
ERROR...

Constraints have their own mechanism of passing args as described in the docs:

约束有它们自己的传递args的机制，如文档所述:

constraints : dict or sequence of dict, optional

约束:命令或命令序列，可选

Constraints definition (only for COBYLA and SLSQP). Each constraint is defined in a dictionary with fields:

约束定义(仅适用于COBYLA和SLSQP)。每个约束都定义在带有字段的字典中:

type : str Constraint type: ‘eq’ for equality, ‘ineq’ for inequality.

类型:str约束类型:' eq ' for equality， ' ineq ' for inequality。

fun : callable The function defining the constraint.

有趣:可调用定义约束的函数。

jac : callable, optional The Jacobian of fun (only for SLSQP).

jac:可调用的，可选的雅可比矩阵(仅用于SLSQP)。

args : sequence, optional Extra arguments to be passed to the function and Jacobian.

args:序列，可选的额外参数传递给函数和雅可比矩阵。

Equality constraint means that the constraint function result is to be zero whereas inequality means that it is to be non-negative. Note that COBYLA only supports inequality constraints.

等式约束意味着约束函数的结果为零，而不等式则意味着它是非负的。注意，COBYLA只支持不等式约束。

In your case:

在你的例子:

con1 = {'type': 'eq', 'fun': create_const}                  # incomplete!
con1 = {'type': 'eq', 'fun': create_const, 'args': (arg,)}  # (,)
                                                            # to make it behave as needed
                                                            # for your code!

This will make it run until some other problem occurs!

这将使它运行，直到出现其他问题!