1/*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 *
4 * Use of this software is governed by the MIT license
5 *
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
9
10#include <assert.h>
11#include <string.h>
12#include <isl_map_private.h>
13#include <isl/aff.h>
14#include <isl/set.h>
15#include "isl_sample.h"
16#include "isl_scan.h"
17#include <isl_seq.h>
18#include <isl_ilp_private.h>
19#include <isl/printer.h>
20#include <isl_point_private.h>
21#include <isl_vec_private.h>
22#include <isl/options.h>
23#include <isl_config.h>
24
25/* The input of this program is the same as that of the "example" program
26 * from the PipLib distribution, except that the "big parameter column"
27 * should always be -1.
28 *
29 * Context constraints in PolyLib format
30 * -1
31 * Problem constraints in PolyLib format
32 * Optional list of options
33 *
34 * The options are
35 *	Maximize	compute maximum instead of minimum
36 *	Rational	compute rational optimum instead of integer optimum
37 *	Urs_parms	don't assume parameters are non-negative
38 *	Urs_unknowns	don't assume unknowns are non-negative
39 */
40
41struct options {
42	struct isl_options	*isl;
43	unsigned		 verify;
44	unsigned		 format;
45};
46
47#define FORMAT_SET	0
48#define FORMAT_AFF	1
49
50struct isl_arg_choice pip_format[] = {
51	{"set",		FORMAT_SET},
52	{"affine",	FORMAT_AFF},
53	{0}
54};
55
56ISL_ARGS_START(struct options, options_args)
57ISL_ARG_CHILD(struct options, isl, "isl", &isl_options_args, "isl options")
58ISL_ARG_BOOL(struct options, verify, 'T', "verify", 0, NULL)
59ISL_ARG_CHOICE(struct options, format, 0, "format",
60	pip_format, FORMAT_SET, "output format")
61ISL_ARGS_END
62
63ISL_ARG_DEF(options, struct options, options_args)
64
65static __isl_give isl_basic_set *set_bounds(__isl_take isl_basic_set *bset)
66{
67	isl_size nparam;
68	int i, r;
69	isl_point *pt, *pt2;
70	isl_basic_set *box;
71
72	nparam = isl_basic_set_dim(bset, isl_dim_param);
73	if (nparam < 0)
74		return isl_basic_set_free(bset);
75	r = nparam >= 8 ? 4 : nparam >= 5 ? 6 : 30;
76
77	pt = isl_basic_set_sample_point(isl_basic_set_copy(bset));
78	pt2 = isl_point_copy(pt);
79
80	for (i = 0; i < nparam; ++i) {
81		pt = isl_point_add_ui(pt, isl_dim_param, i, r);
82		pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, r);
83	}
84
85	box = isl_basic_set_box_from_points(pt, pt2);
86
87	return isl_basic_set_intersect(bset, box);
88}
89
90static __isl_give isl_basic_set *to_parameter_domain(
91	__isl_take isl_basic_set *context)
92{
93	isl_size dim;
94
95	dim = isl_basic_set_dim(context, isl_dim_set);
96	if (dim < 0)
97		return isl_basic_set_free(context);
98	context = isl_basic_set_move_dims(context, isl_dim_param, 0,
99		    isl_dim_set, 0, dim);
100	context = isl_basic_set_params(context);
101	return context;
102}
103
104/* If "context" has more parameters than "bset", then reinterpret
105 * the last dimensions of "bset" as parameters.
106 */
107static __isl_give isl_basic_set *move_parameters(__isl_take isl_basic_set *bset,
108	__isl_keep isl_basic_set *context)
109{
110	isl_size nparam, nparam_bset, dim;
111
112	nparam = isl_basic_set_dim(context, isl_dim_param);
113	nparam_bset = isl_basic_set_dim(bset, isl_dim_param);
114	if (nparam < 0 || nparam_bset < 0)
115		return isl_basic_set_free(bset);
116	if (nparam == nparam_bset)
117		return bset;
118	dim = isl_basic_set_dim(bset, isl_dim_set);
119	if (dim < 0)
120		return isl_basic_set_free(bset);
121	bset = isl_basic_set_move_dims(bset, isl_dim_param, 0,
122					    isl_dim_set, dim - nparam, nparam);
123	return bset;
124}
125
126/* Plug in the initial values of "params" for the parameters in "bset" and
127 * return the result.  The remaining entries in "params", if any,
128 * correspond to the existentially quantified variables in the description
129 * of the original context and can be ignored.
130 */
131static __isl_give isl_basic_set *plug_in_parameters(
132	__isl_take isl_basic_set *bset, __isl_take isl_vec *params)
133{
134	int i;
135	isl_size n;
136
137	n = isl_basic_set_dim(bset, isl_dim_param);
138	if (n < 0)
139		bset = isl_basic_set_free(bset);
140	for (i = 0; i < n; ++i)
141		bset = isl_basic_set_fix(bset,
142					 isl_dim_param, i, params->el[1 + i]);
143
144	bset = isl_basic_set_remove_dims(bset, isl_dim_param, 0, n);
145
146	isl_vec_free(params);
147
148	return bset;
149}
150
151/* Plug in the initial values of "params" for the parameters in "set" and
152 * return the result.  The remaining entries in "params", if any,
153 * correspond to the existentially quantified variables in the description
154 * of the original context and can be ignored.
155 */
156static __isl_give isl_set *set_plug_in_parameters(__isl_take isl_set *set,
157	__isl_take isl_vec *params)
158{
159	int i;
160	isl_size n;
161
162	n = isl_set_dim(set, isl_dim_param);
163	if (n < 0)
164		set = isl_set_free(set);
165	for (i = 0; i < n; ++i)
166		set = isl_set_fix(set, isl_dim_param, i, params->el[1 + i]);
167
168	set = isl_set_remove_dims(set, isl_dim_param, 0, n);
169
170	isl_vec_free(params);
171
172	return set;
173}
174
175/* Compute the lexicographically minimal (or maximal if max is set)
176 * element of bset for the given values of the parameters, by
177 * successively solving an ilp problem in each direction.
178 */
179static __isl_give isl_vec *opt_at(__isl_take isl_basic_set *bset,
180	__isl_take isl_vec *params, int max)
181{
182	isl_size dim;
183	isl_ctx *ctx;
184	struct isl_vec *opt;
185	struct isl_vec *obj;
186	int i;
187
188	dim = isl_basic_set_dim(bset, isl_dim_set);
189	if (dim < 0)
190		goto error;
191
192	bset = plug_in_parameters(bset, params);
193
194	ctx = isl_basic_set_get_ctx(bset);
195	if (isl_basic_set_plain_is_empty(bset)) {
196		opt = isl_vec_alloc(ctx, 0);
197		isl_basic_set_free(bset);
198		return opt;
199	}
200
201	opt = isl_vec_alloc(ctx, 1 + dim);
202	assert(opt);
203
204	obj = isl_vec_alloc(ctx, 1 + dim);
205	assert(obj);
206
207	isl_int_set_si(opt->el[0], 1);
208	isl_int_set_si(obj->el[0], 0);
209
210	for (i = 0; i < dim; ++i) {
211		enum isl_lp_result res;
212
213		isl_seq_clr(obj->el + 1, dim);
214		isl_int_set_si(obj->el[1 + i], 1);
215		res = isl_basic_set_solve_ilp(bset, max, obj->el,
216						&opt->el[1 + i], NULL);
217		if (res == isl_lp_empty)
218			goto empty;
219		assert(res == isl_lp_ok);
220		bset = isl_basic_set_fix(bset, isl_dim_set, i, opt->el[1 + i]);
221	}
222
223	isl_basic_set_free(bset);
224	isl_vec_free(obj);
225
226	return opt;
227error:
228	isl_basic_set_free(bset);
229	isl_vec_free(params);
230	return NULL;
231empty:
232	isl_vec_free(opt);
233	opt = isl_vec_alloc(ctx, 0);
234	isl_basic_set_free(bset);
235	isl_vec_free(obj);
236
237	return opt;
238}
239
240struct isl_scan_pip {
241	struct isl_scan_callback callback;
242	isl_basic_set *bset;
243	isl_set *sol;
244	isl_set *empty;
245	int stride;
246	int n;
247	int max;
248};
249
250/* Check if the "manually" computed optimum of bset at the "sample"
251 * values of the parameters agrees with the solution of pilp problem
252 * represented by the pair (sol, empty).
253 * In particular, if there is no solution for this value of the parameters,
254 * then it should be an element of the parameter domain "empty".
255 * Otherwise, the optimal solution, should be equal to the result of
256 * plugging in the value of the parameters in "sol".
257 */
258static isl_stat scan_one(struct isl_scan_callback *callback,
259	__isl_take isl_vec *sample)
260{
261	struct isl_scan_pip *sp = (struct isl_scan_pip *)callback;
262	struct isl_vec *opt;
263
264	sp->n--;
265
266	opt = opt_at(isl_basic_set_copy(sp->bset), isl_vec_copy(sample), sp->max);
267	assert(opt);
268
269	if (opt->size == 0) {
270		isl_point *sample_pnt;
271		sample_pnt = isl_point_alloc(isl_set_get_space(sp->empty), sample);
272		assert(isl_set_contains_point(sp->empty, sample_pnt));
273		isl_point_free(sample_pnt);
274		isl_vec_free(opt);
275	} else {
276		isl_set *sol;
277		isl_set *opt_set;
278		opt_set = isl_set_from_basic_set(isl_basic_set_from_vec(opt));
279		sol = set_plug_in_parameters(isl_set_copy(sp->sol), sample);
280		assert(isl_set_is_equal(opt_set, sol));
281		isl_set_free(sol);
282		isl_set_free(opt_set);
283	}
284
285	if (!(sp->n % sp->stride)) {
286		printf("o");
287		fflush(stdout);
288	}
289
290	return sp->n >= 1 ? isl_stat_ok : isl_stat_error;
291}
292
293static void check_solution(isl_basic_set *bset, isl_basic_set *context,
294	isl_set *sol, isl_set *empty, int max)
295{
296	struct isl_scan_pip sp;
297	isl_int count, count_max;
298	int i, n;
299	int r;
300
301	context = set_bounds(context);
302	context = isl_basic_set_underlying_set(context);
303
304	isl_int_init(count);
305	isl_int_init(count_max);
306
307	isl_int_set_si(count_max, 2000);
308	r = isl_basic_set_count_upto(context, count_max, &count);
309	assert(r >= 0);
310	n = isl_int_get_si(count);
311
312	isl_int_clear(count_max);
313	isl_int_clear(count);
314
315	sp.callback.add = scan_one;
316	sp.bset = bset;
317	sp.sol = sol;
318	sp.empty = empty;
319	sp.n = n;
320	sp.stride = n > 70 ? 1 + (n + 1)/70 : 1;
321	sp.max = max;
322
323	for (i = 0; i < n; i += sp.stride)
324		printf(".");
325	printf("\r");
326	fflush(stdout);
327
328	isl_basic_set_scan(context, &sp.callback);
329
330	printf("\n");
331
332	isl_basic_set_free(bset);
333}
334
335int main(int argc, char **argv)
336{
337	struct isl_ctx *ctx;
338	struct isl_basic_set *context, *bset, *copy, *context_copy;
339	struct isl_set *set = NULL;
340	struct isl_set *empty;
341	isl_pw_multi_aff *pma = NULL;
342	int neg_one;
343	char s[1024];
344	int urs_parms = 0;
345	int urs_unknowns = 0;
346	int max = 0;
347	int rational = 0;
348	int n;
349	struct options *options;
350
351	options = options_new_with_defaults();
352	assert(options);
353	argc = options_parse(options, argc, argv, ISL_ARG_ALL);
354
355	ctx = isl_ctx_alloc_with_options(&options_args, options);
356
357	context = isl_basic_set_read_from_file(ctx, stdin);
358	assert(context);
359	n = fscanf(stdin, "%d", &neg_one);
360	assert(n == 1);
361	assert(neg_one == -1);
362	bset = isl_basic_set_read_from_file(ctx, stdin);
363
364	while (fgets(s, sizeof(s), stdin)) {
365		if (strncasecmp(s, "Maximize", 8) == 0)
366			max = 1;
367		if (strncasecmp(s, "Rational", 8) == 0) {
368			rational = 1;
369			bset = isl_basic_set_set_rational(bset);
370		}
371		if (strncasecmp(s, "Urs_parms", 9) == 0)
372			urs_parms = 1;
373		if (strncasecmp(s, "Urs_unknowns", 12) == 0)
374			urs_unknowns = 1;
375	}
376	if (!urs_parms)
377		context = isl_basic_set_intersect(context,
378		isl_basic_set_positive_orthant(isl_basic_set_get_space(context)));
379	context = to_parameter_domain(context);
380	bset = move_parameters(bset, context);
381	if (!urs_unknowns)
382		bset = isl_basic_set_intersect(bset,
383		isl_basic_set_positive_orthant(isl_basic_set_get_space(bset)));
384
385	if (options->verify) {
386		copy = isl_basic_set_copy(bset);
387		context_copy = isl_basic_set_copy(context);
388	}
389
390	if (options->format == FORMAT_AFF) {
391		if (max)
392			pma = isl_basic_set_partial_lexmax_pw_multi_aff(bset,
393								context, &empty);
394		else
395			pma = isl_basic_set_partial_lexmin_pw_multi_aff(bset,
396								context, &empty);
397	} else {
398		if (max)
399			set = isl_basic_set_partial_lexmax(bset,
400								context, &empty);
401		else
402			set = isl_basic_set_partial_lexmin(bset,
403								context, &empty);
404	}
405
406	if (options->verify) {
407		assert(!rational);
408		if (options->format == FORMAT_AFF)
409			set = isl_set_from_pw_multi_aff(pma);
410		check_solution(copy, context_copy, set, empty, max);
411		isl_set_free(set);
412	} else {
413		isl_printer *p;
414		p = isl_printer_to_file(ctx, stdout);
415		if (options->format == FORMAT_AFF)
416			p = isl_printer_print_pw_multi_aff(p, pma);
417		else
418			p = isl_printer_print_set(p, set);
419		p = isl_printer_end_line(p);
420		p = isl_printer_print_str(p, "no solution: ");
421		p = isl_printer_print_set(p, empty);
422		p = isl_printer_end_line(p);
423		isl_printer_free(p);
424		isl_set_free(set);
425		isl_pw_multi_aff_free(pma);
426	}
427
428	isl_set_free(empty);
429	isl_ctx_free(ctx);
430
431	return 0;
432}
433