https://github.com/baharev/ManiSolve
Tip revision: 4df5c5f8c67154d8faf93dafcf5df85ce9899e1a authored by Ali Baharev on 18 November 2019, 10:48:51 UTC
Adding a CITATION file
Adding a CITATION file
Tip revision: 4df5c5f
nl_interpreter.py
# Copyright (C) 2016, 2017 University of Vienna
# All rights reserved.
# BSD license.
# Author: Ali Baharev <ali.baharev@gmail.com>
from __future__ import print_function, division
from collections import namedtuple, OrderedDict
from itertools import chain
from toolz import partitionby
from auxfiles import get_col_names, get_row_names, get_def_var_names, nl_fname,\
get_solutions
from opcodes import ARITY, NAME
from py3compat import imap, irange, izip
from utils import warning, as_pairs, get_lines
Problem = namedtuple('Problem', '''name nl_header segments col_names
def_var_names row_names solutions''')
Header = namedtuple('Header', 'n_vars n_cons n_nzeros n_objs')
#-------------------------------------------------------------------------------
def get_nl_file_content(full_nl_file_name):
lines = get_lines(full_nl_file_name)
lines = remove_comments(lines)
header, content = split_to_header_and_content(lines)
nl_header = parse_nl_header(header)
return nl_header, get_segments(nl_header, content)
def parse_nl_file(probname):
return get_nl_file_content(nl_fname(probname))
def remove_comments(lines):
# 'v1 #x[6]' -> 'v1'
return list(imap(lambda l: l.rsplit('#', 1)[0].rstrip(), lines))
#-------------------------------------------------------------------------------
# Stuff that mainly concerns the header
#
# 0: g3 1 1 0 # problem silly
# 1: 9 4 3 1 1 2 # vars, constraints, objectives, ranges, eqns, lcons
# 2: 3 3 # nonlinear constraints, objectives
# 3: 0 0 # network constraints: nonlinear, linear
# 4: 7 8 5 # nonlinear vars in constraints, objectives, both
# 5: 0 1 0 1 # linear network variables; functions; arith, flags
# 6: 0 0 0 0 0 # discrete variables: binary, integer, nonlinear (b,c,o)
# 7: 17 5 # nonzeros in Jacobian, gradients
# 8: 5 4 # max name lengths: constraints, variables
# 9: 2 0 0 3 1 # common exprs: b,c,o,c1,o1
# The lines below the header have a non-blank first character
def split_to_header_and_content(lines):
if lines[0][0] != 'g':
raise IOError('ASCII .nl file is expected.\nPlease generate the'
' .nl file with "write gname;" in the AMPL .mod file\n'
'or as "ampl -ogname name.mod" in the command line.')
parts = list(partitionby(lambda l: l[0]==' ', lines))
header, content = list(chain(parts[0], parts[1])), parts[2]
return header, content
def parse_nl_header(header):
assert len(header) == 10, 'See expected format above in comment'
header = list(imap(to_ints, header[1:])) # <- discarding the first line!
n_vars, n_cons, n_objs = header[0][:3]
nzeros = header[6][0]
if any(header[5]):
warning('The model has discrete variables but we treat them as continuous!')
return Header(n_vars=n_vars, n_cons=n_cons, n_nzeros=nzeros, n_objs=n_objs)
def to_ints(line):
# '1 2' -> [1, 2]
entries = line.split()
return list(imap(int, entries))
#-------------------------------------------------------------------------------
# F imported function description
# S suffix values
# V defined variable definition (must precede V,C,L,O segments where used)
# C algebraic constraint body
# L logical constraint expression
# O objective function
# d dual initial guess
# x primal initial guess
# r bounds on algebraic constraint bodies ("ranges")
# b bounds on variable
# k Jacobian column counts (must precede all J segments)
# J Jacobian sparsity, linear terms
# G Gradient sparsity, linear terms
SEGMENTS = {'F', 'S', 'V', 'C', 'L', 'O', 'd', 'x', 'r', 'b', 'k', 'J', 'G'}
# FIXME Start documenting the shape of the members!
Segments = namedtuple('Segments', '''def_vars cons objs var_bnds
initial_guess con_jacobian eval_order
con_blocks var_blocks''')
def get_segments(nl_header, content):
# Returns Segments.
n_cons = nl_header.n_cons
n_objs = nl_header.n_objs
allsegs = Segments(def_vars=OrderedDict(), cons=[None]*n_cons,
objs=[None]*n_objs, var_bnds=[], initial_guess={},
con_jacobian=[], eval_order=[], con_blocks=[],
var_blocks=[])
# Creates a sequence of tuples. Two adjacent tuples are the header and the
# body of the segment:
# V68 0 0
# o3 -> ('V68 0 0',), ('o3', 'v2', 'v53')
# v2
# v53
segments_w_header = list(partitionby(lambda l: l[0] in SEGMENTS, content))
segments_w_header = fix_empty_segments(segments_w_header)
# Assumption: Constraints precede r segments. Enforce it:
segments_w_header = move_r_beyond_last_C(segments_w_header)
# Populate allsegs by dispatching to the corresponding segment handler.
for header, segment in segments_w_header:
ADD_SEGMENT.get(header[0], ignored)(header, segment, allsegs)
return allsegs
def move_r_beyond_last_C(segments_w_header):
r_index = first_index_of(segments_w_header, 'r')
if r_index is None:
warning('Apparently the problem does not have r segments.')
return segments_w_header
last_C_index = last_index_of(segments_w_header, 'C')
assert last_C_index is not None
if r_index > last_C_index:
return segments_w_header
pos = last_C_index + 1
new_lst = segments_w_header[:r_index] # elems before r segment
new_lst.extend(segments_w_header[r_index+1:pos]) # elems after r till last C
new_lst.append(segments_w_header[r_index]) # move r segment here
new_lst.extend(segments_w_header[pos:]) # copy the rest after last C
assert len(new_lst) == len(segments_w_header)
return new_lst
def first_index_of(segments_w_header, header):
# returns None if not found
return next((i for i, e in enumerate(segments_w_header) if e[0][0]==header), None)
def last_index_of(segments_w_header, header):
# returns None if not found
idx = first_index_of(segments_w_header[::-1], header)
return len(segments_w_header) - idx - 1 if idx is not None else None
# TODO Write an appropriate algorithm to replace partitionby, which then
# eliminates the need for this post-processing.
def fix_empty_segments(segments_w_header):
# The assumption is that each header is followed by a segment. However,
# if the segment is empty, then the header is followed by an another header.
new_segments_with_header = [ ]
for headers, segment in as_pairs(segments_w_header):
n_headers = len(headers)
assert n_headers
if n_headers == 1:
# Everything is fine, just one header with a segment
new_segments_with_header.append((headers[0], segment))
else:
# Insert fake empty segments after each header
segments = [tuple() for _ in irange(n_headers-1)]
segments.append(segment)
for h, s in izip(headers, segments):
new_segments_with_header.append((h, s))
return new_segments_with_header
def V_segment(header, segment, allsegs):
index, length, _k = to_ints(header[1:])
linear_part, nl_part = segment[:length], segment[length:]
linear_part = parse_linear_part(linear_part)
V_segs = allsegs.def_vars
assert index not in V_segs, index
V_segs[index] = (nl_part, linear_part)
allsegs.eval_order.append(('V', index))
def C_segment(header, segment, allsegs):
i = int(header[1:])
cons = allsegs.cons
assert cons[i] is None
cons[i] = (segment, None, None)
allsegs.eval_order.append(('C', i))
def J_segment(header, segment, allsegs):
i, j = to_ints(header[1:])
indices, coefficients = parse_linear_part(segment)
assert len(indices) == j, i
nl_part, lin_part, rng = allsegs.cons[i]
assert lin_part is None, (i, lin_part)
allsegs.cons[i] = (nl_part, (indices, coefficients), rng)
def O_segment(header, segment, allsegs):
# similar to the C segment (nl_part of the constraints)
i, j = to_ints(header[1:])
O_segs = allsegs.objs
assert O_segs[i] is None
direction = 'minimize' if j == 0 else 'maximize'
O_segs[i] = (direction, segment, None)
allsegs.eval_order.append(('O', i))
def G_segment(header, segment, allsegs):
# similar to the J segment (lin_part of the constraints)
i, j = to_ints(header[1:])
indices, coefficients = parse_linear_part(segment)
assert len(indices) == j, i
direction, nl_part, lin_part = allsegs.objs[i]
assert lin_part is None, (i, lin_part)
allsegs.objs[i] = (direction, nl_part, (indices, coefficients))
def r_segment(_header, segment, allsegs):
cons = allsegs.cons
for i, rng in enumerate(imap(to_range, segment)):
nl_part, lin_part, _None = cons[i]
cons[i] = (nl_part, lin_part, rng)
def b_segment(_header, segment, allsegs):
allsegs.var_bnds.extend(imap(to_range, segment))
def x_segment(_header, segment, allsegs):
allsegs.initial_guess.update(imap(to_index_value, segment))
def k_segment(_header, segment, allsegs):
allsegs.con_jacobian.extend(imap(int, segment))
def S_segment(header, segment, allsegs):
kind, _length, name = header.split()
if name != 'blockid':
warning('S segment "%s" ignored' % header)
return
kind = int(kind[1:])
assert kind == 0 or kind == 1, (kind, header)
block_ids = allsegs.con_blocks if kind == 1 else allsegs.var_blocks
block_ids.extend(parse_blockids(segment))
def ignored(header, *_args, **_kwargs):
warning('%s segment ignored' % header[0])
def to_index_value(line):
parts = line.split()
assert len(parts) == 2
return int(parts[0]), parts[1]
def parse_linear_part(segment):
indices, coeffs = [ ], [ ]
for i, c in imap(lambda l: l.split(), segment):
indices.append(int(i)) # make the indices integers
coeffs.append(c) # but keep the coefficients as strings
return indices, coeffs
def parse_blockids(segment):
return [(i, j) for i, j in imap(to_ints, segment)]
ADD_SEGMENT = {
'C': C_segment,
'V': V_segment,
'J': J_segment,
'O': O_segment,
'G': G_segment,
'r': r_segment,
'b': b_segment,
'x': x_segment,
'k': k_segment,
'S': S_segment,
}
class Range:
lb_ub = 0 # lb <= body <= ub
ub = 1 # body <= ub
lb = 2 # lb <= body
no_bnd = 3 # body is unbounded
eq = 4 # body = c
compl = 5 # "5 k i" in the .nl file, complementary constraint, ignoring it
TO_RANGE_KIND = {
'0': Range.lb_ub,
'1': Range.ub,
'2': Range.lb,
'3': Range.no_bnd,
'4': Range.eq,
'5': Range.compl
}
def to_range(line):
parts = line.split()
return TO_RANGE_KIND[parts[0]], tuple(parts[1:])
#-------------------------------------------------------------------------------
# Interpreting nonlinear constraint segments
# https://en.wikipedia.org/wiki/Reverse_Polish_notation#Postfix_algorithm
#
# While there are input tokens left
# Read the next token from input.
# If the token is a value
# Push it onto the stack.
# Otherwise, the token is an operator (operator here includes both operators and functions).
# It is already known that the operator takes n arguments.
# If there are fewer than n values on the stack
# (Error) The user has not input sufficient values in the expression.
# Else, Pop the top n values from the stack.
# Evaluate the operator, with the values as arguments.
# Push the returned results, if any, back onto the stack.
# If there is only one value in the stack
# That value is the result of the calculation.
# Otherwise, there are more values in the stack
# (Error) The user input has too many values.
def interpret(body, counter):
segment = list(body)
stack = []
evaluated = ('n', 'v', 't')
while segment:
token = segment.pop()
kind = token[0]
if kind in evaluated or token.isdigit():
stack.append(token)
else:
assert kind == 'o', kind
assert token in NAME, token
arity = ARITY[token]
if not arity: # n-ary operators have their arity as the first argument
arity = int(stack.pop())
assert len(stack) >= arity
args = [stack.pop() for _ in irange(arity)]
tmp = 't%d' % next(counter)
#print(tmp, '=', NAME[token], ' '.join(args))
yield tmp, NAME[token], args
stack.append(tmp)
assert len(stack) == 1, stack
assert stack[0][0] in evaluated, stack
yield stack[0]
#-------------------------------------------------------------------------------
def get_problem(name):
nl_header, segments = parse_nl_file(name)
col_names = get_col_names(name)
def_var_names = get_def_var_names(name)
row_names = get_row_names(name)
solutions = get_solutions(name, col_names)
return Problem(name=name, nl_header=nl_header, segments=segments,
col_names=col_names, def_var_names=def_var_names,
row_names=row_names, solutions=solutions)
def get_J_rowwise(problem):
cons = problem.segments.cons
# Get the sparsity pattern as a sparse matrix rowwise
return [indices for _nl_part, (indices, _coefficients), _rng in cons]