<h3>基于AMPLPY</h3>
<p>正如@denfromufa指出的,有一个<code>AMPL</code>接口到<code>PATH</code>解算器。他建议使用<code>Pyomo</code>,因为它是开源的,能够处理{<cd1>}。然而,<code>Pyomo</code>结果是缓慢和乏味的工作。我最终在cython中编写了自己的接口到<code>PATH</code>解算器,并希望在某个时候发布它,但目前它没有输入卫生设施,速度快而且脏,而且我看不到时间来处理它。在</p>
<p>现在,我可以分享一个使用python扩展<code>AMPL</code>的答案。它不如<code>PATH</code>的直接接口快:对于我们要解决的每个<code>LCP</code>,它都会创建一个(临时)模型文件,运行<code>AMPL</code>,并收集输出。这有点快又脏,但我觉得我至少应该报告几个月来提出这个问题以来的部分结果。在</p>
<pre><code>import os
# PATH license string needs to be updated
os.environ['PATH_LICENSE_STRING'] = "3413119131&Courtesy&&&USR&54784&12_1_2016&1000&PATH&GEN&31_12_2017&0_0_0&5000&0_0"
from amplpy import AMPL, Environment, dataframe
import numpy as np
from scipy import sparse
from tempfile import mkstemp
import os
import sys
import contextlib
class DummyFile(object):
def write(self, x): pass
@contextlib.contextmanager
def nostdout():
save_stdout = sys.stdout
sys.stdout = DummyFile()
yield
sys.stdout = save_stdout
class modFile:
# context manager: create temporary file, insert model data, and supply file name
# apparently, amplpy coders are inable to support StringIO
content = """
set Rn;
param B {Rn,Rn} default 0;
param q {Rn} default 0;
var x {j in Rn};
s.t. f {i in Rn}:
0 <= x[i]
complements
sum {j in Rn} B[i,j]*x[j]
>= -q[i];
"""
def __init__(self):
self.fd = None
self.temp_path = None
def __enter__(self):
fd, temp_path = mkstemp()
file = open(temp_path, 'r+')
file.write(self.content)
file.close()
self.fd = fd
self.temp_path = temp_path
return self
def __exit__(self, exc_type, exc_val, exc_tb):
os.close(self.fd)
os.remove(self.temp_path)
def solveLCP(B, q, x=None, env=None, binaryDirectory=None, pathOptions={'logfile':'logpath.tmp' }, verbose=False):
# x: initial guess
if binaryDirectory is not None:
env = Environment(binaryDirectory='/home/foo/amplide.linux64/')
if verbose:
pathOptions['output'] = 'yes'
ampl = AMPL(environment=env)
# read model
with modFile() as mod:
ampl.read(mod.temp_path)
n = len(q)
# prepare dataframes
dfQ = dataframe.DataFrame('Rn', 'c')
for i in np.arange(n):
# shitty amplpy does not support np.float
dfQ.addRow(int(i)+1, np.float(q[i]))
dfB = dataframe.DataFrame(('RnRow', 'RnCol'), 'val')
if sparse.issparse(B):
if not isinstance(B, sparse.lil_matrix):
B = B.tolil()
dfB.setValues({
(i+1, j+1): B.data[i][jPos]
for i, row in enumerate(B.rows)
for jPos, j in enumerate(row)
})
else:
r = np.arange(n) + 1
Rrow, Rcol = np.meshgrid(r, r, indexing='ij')
#dfB.setColumn('RnRow', [np.float(x) for x in Rrow.reshape((-1), order='C')])
dfB.setColumn('RnRow', list(Rrow.reshape((-1), order='C').astype(float)))
dfB.setColumn('RnCol', list(Rcol.reshape((-1), order='C').astype(float)))
dfB.setColumn('val', list(B.reshape((-1), order='C').astype(float)))
# set values
ampl.getSet('Rn').setValues([int(x) for x in np.arange(n, dtype=int)+1])
if x is not None:
dfX = dataframe.DataFrame('Rn', 'x')
for i in np.arange(n):
# shitty amplpy does not support np.float
dfX.addRow(int(i)+1, np.float(x[i]))
ampl.getVariable('x').setValues(dfX)
ampl.getParameter('q').setValues(dfQ)
ampl.getParameter('B').setValues(dfB)
# solve
ampl.setOption('solver', 'pathampl')
pathOptions = ['{}={}'.format(key, val) for key, val in pathOptions.items()]
ampl.setOption('path_options', ' '.join(pathOptions))
if verbose:
ampl.solve()
else:
with nostdout():
ampl.solve()
if False:
bD = ampl.getParameter('B').getValues().toDict()
qD = ampl.getParameter('q').getValues().toDict()
xD = ampl.getVariable('x').getValues().toDict()
BB = ampl.getParameter('B').getValues().toPandas().values.reshape((n, n,), order='C')
qq = ampl.getParameter('q').getValues().toPandas().values[:, 0]
xx = ampl.getVariable('x').getValues().toPandas().values[:, 0]
ineq2 = BB.dot(xx) + qq
print((xx * ineq2).min(), (xx * ineq2).max() )
return ampl.getVariable('x').getValues().toPandas().values[:, 0]
if __name__ == '__main__':
# solve problem from the Julia port at https://github.com/chkwon/PATHSolver.jl
n = 4
B = np.array([[0, 0, -1, -1], [0, 0, 1, -2], [1, -1, 2, -2], [1, 2, -2, 4]])
q = np.array([2, 2, -2, -6])
BSparse = sparse.lil_matrix(B)
env = Environment(binaryDirectory='/home/foo/amplide.linux64/')
print(solveLCP(B, q, env=env))
print(solveLCP(BSparse, q, env=env))
# to test licensing
from numpy import random
n = 1000
B = np.diag((random.randn(n)))
q = np.ones((n,))
print(solveLCP(B, q, env=env))
</code></pre>