ограничение линейного решателя только неотрицательными решениями python - PullRequest
Новая
       23

ограничение линейного решателя только неотрицательными решениями python

0 голосов
/ 03 июля 2018

Я разрабатываю программу, которая преобразует куб данных изображений, снятых солнцем, в температурный профиль для каждого изображения. Работает путем линейной оптимизации матриц.

Итак, я попробовал два разных линейных решателя. Scipy.optimize.linprog симплекс-метод. Этот метод работает, но очень медленный. Я обнаружил, что cvxopt также имеет функцию, которая, как я думал, будет работать - solvers.lp, но должна была внести некоторые изменения в настройку. Я добавил в матрицу ограничений матрицу формы mxm с -1 по диагонали.

НО! этот метод дал решения, которые были отрицательными. Мне нужен решатель, который ограничен неотрицательными решениями. Я смотрел на cvxpy, потому что есть возможность сказать, что nonneg = True, но я не могу понять, как перевести их примеры в то, что я делаю. Вот фрагмент кода соскученного, который работал. Res - это место, где функция линейной алгебры фактически выполняется.

после этого я добавил часть программы, которую мне пришлось изменить, чтобы попытаться использовать cvxopt. 

def DEM_solve(img, Dict, lgtaxis, tolfac):
global basis_funcs
t0 = time.time()
eps = 1e-3
relax = 1
symmbuff = 1.0
adaptive_tolfac = 1
dim = np.array(img).shape
NOCOUNTS = np.where(np.array(img).sum(axis=0) < 10*eps)

ntemp = len(Dict[:,0])
coeffs = np.zeros((dim[2], dim[1], ntemp))
zmax = np.zeros((dim[2], dim[1]))
status = np.zeros((dim[2], dim[1]))
tolmap = np.zeros((dim[2], dim[1]))

zequation = np.zeros(ntemp) #63 elements
zequation[:] = 1.0 #simplex in python is a MINIMIZATION not maximization so the zequation is all +1.0 not -1.0

constraint = np.zeros((ntemp, 8)) # [8 63]
constraint[0:ntemp, 0:4] = Dict # row 1-64 and columns 0-4
constraint[0:ntemp, 4:8] = -Dict # row 1-64 and columns 5-8
constraint = np.transpose(constraint)
B = np.zeros(8)

X, Y = np.meshgrid(np.arange(-int(len(img[0,:,0]))/2, int(len(img[0,0,:])/2)), np.arange(-int(len(img[0,:,0]))/2, int(len(img[0,0,:])/2))) 
xsi, ysi = np.where(np.hypot(X,Y) < 20) #limits number of pixels looked at for now since the solver is so slow
xyzip = list(zip_longest(xsi, ysi))
for i, j in xyzip:
        temp = []
        tol = []
        for k in range(4):
           if img[k][i][j] < 0:
               temp.append(0)
           else:
               temp.append(img[k][i][j])
        for k in temp:
           if k > 1.0:
               tol.append(tolfac*0.8*math.sqrt(k))
           else:
               tol.append(tolfac*0.8)

        B[0:4] = [a+b for a, b in zip(tol, temp)]
        B[4:8] = list(map(lambda x: -x if x > 0 else 0,  map(sub, temp, tol)))


        res = scipy.optimize.linprog(zequation,
                                     A_ub = constraint,
                                     b_ub = B,
                                     options = {'tol': eps*np.max(temp)*8e-4})
        coeff = res.x
        soln = -res.fun
        t1p = time.time()
        if isinstance(coeff, float): #if nan 
            new = [0]*ntemp
            s = 10
        else: # successfully provided a list
            new = [soln] + coeff.tolist()
            s = res.status

        # if one of the new solutions is less than zero
        # make all zero, otherwise keep them
        # if nan befor or < 0 flags as unsuccessful
        if min(new[1:]) < 0:
            coeffs[i][j][0:ntemp] = 0.0
            s = 10
        else:
            coeffs[i][j][0:ntemp-1] = new[1:ntemp]

        zmax[i][j] = new[0]
        status[i][j] = s

if len(NOCOUNTS[0]) != 0:
    status[NOCOUNTS] = 11.0
    for i in range(len(lgtaxis)):
        coeffsnew = np.squeeze(coeffs[:,:,i])
        coeffsnew[NOCOUNTS] = 0.0
        coeffs[:,:,i] = coeffsnew
oem = np.zeros((dim[1], dim[2], len(lgtaxis)))
for i, j in xyzip:
      oem[i, j, :] = np.squeeze(np.matmul(np.squeeze(coeffs[i, j, :]), np.transpose(basis_funcs)))

return(oem)


zequation = np.zeros(ntemp) #63 elements
zequation[:] = 1.0 #simplex in python is a MINIMIZATION not maximization so the zequation is all +1.0 not -1.0
zequation = matrix(zequation)
constraint = np.zeros((ntemp, 8)) # 63 rows by 8 columns
constraint[0:ntemp, 0:4] = Dict # row 1-64 and columns 0-4
constraint[0:ntemp, 4:8] = -Dict # row 1-64 and columns 5-8
identity = np.zeros((ntemp, ntemp))
np.fill_diagonal(identity, -1)
constraint = np.concatenate((constraint, identity), axis = 1)
constraint = matrix(np.transpose(constraint)) # 71 rows by 63 columns
B = np.zeros(8+63)
X, Y = np.meshgrid(np.arange(-int(len(img[0,:,0]))/2, int(len(img[0,0,:])/2)), np.arange(-int(len(img[0,:,0]))/2, int(len(img[0,0,:])/2))) 
xsi, ysi = np.where(np.hypot(X,Y) < 5)
xyzip = list(zip_longest(xsi, ysi))
for i, j in xyzip:
        temp = []
        tol = []
        for k in range(4):
           if img[k][i][j] < 0:
               temp.append(0)
           else:
               temp.append(img[k][i][j])
        for k in temp:
           if k > 1.0:
               tol.append(tolfac*0.8*math.sqrt(k))
           else:
               tol.append(tolfac*0.8)

        B[0:4] = [a+b for a, b in zip(tol, temp)]
        B[4:8] = list(map(lambda x: -x if x > 0 else 0,  map(sub, temp, tol)))
        B[9:] = 0
        B = matrix(B)
        sol = solvers.lp(zequation, constraint, B, solver = 'glpk')
        coeff = list(sol['x'])
        print(coeff)
        soln = sol['primal objective']
...