1% ----------------------------------------------------------------------
2% BEGIN LICENSE BLOCK
3% Version: CMPL 1.1
4%
5% The contents of this file are subject to the Cisco-style Mozilla Public
6% License Version 1.1 (the "License"); you may not use this file except
7% in compliance with the License.  You may obtain a copy of the License
8% at www.eclipse-clp.org/license.
9%
10% Software distributed under the License is distributed on an "AS IS"
11% basis, WITHOUT WARRANTY OF ANY KIND, either express or implied.  See
12% the License for the specific language governing rights and limitations
13% under the License.
14%
15% The Original Code is  The ECLiPSe Constraint Logic Programming System.
16% The Initial Developer of the Original Code is  Cisco Systems, Inc.
17% Portions created by the Initial Developer are
18% Copyright (C) 1999-2006 Cisco Systems, Inc.  All Rights Reserved.
19%
20% Contributor(s): IC-Parc, Imperal College London
21%
22% END LICENSE BLOCK
23%
24% System:	ECLiPSe Constraint Logic Programming System
25% Version:	$Id: branch_and_bound.pl,v 1.7 2015/04/04 22:05:50 jschimpf Exp $
26% ----------------------------------------------------------------------
27
28:- module(branch_and_bound).
29
30:- comment(categories, ["Algorithms"]).
31:- comment(summary, "Generic branch-and-bound optimization").
32:- comment(copyright, "Cisco Systems, Inc").
33:- comment(author, "Joachim Schimpf, Vassilis Liatsos, IC-Parc, Imperial College, London").
34:- comment(date, "$Date: 2015/04/04 22:05:50 $").
35:- comment(index, ["branch-and-bound","dichotomic search","optimization"]).
36:- comment(desc, html("<P>
37    This is a solver-independent library implementing branch-and-bound
38    optimization. It can be used with any nondeterministic search routine
39    that instantiates a cost variable when a solution is found. The cost
40    variable can be an arbitrary numerical domain variable or even a
41    simple domain-less Prolog variable.
42</P><P>
43    The main primitive is bb_min/3.  Assume we have the following
44    collection of facts:
45<PRE>
46        % item(Food, Calories, Price)
47        item(bread,  500, 1.50).
48        item(eggs,   600, 1.99).
49        item(milk,   400, 0.99).
50        item(apples, 200, 1.39).
51        item(butter, 800, 1.89).
52</PRE>
53   Then we can find a minimum-calorie solution as follows:
54<PRE>
55        ?- bb_min(item(Food,Cal,Price), Cal, _).
56        Found a solution with cost 500
57        Found a solution with cost 400
58        Found a solution with cost 200
59        Found no solution with cost -1.0Inf .. 199
60
61        Food = apples
62        Cal = 200
63        Price = 1.39
64        Yes (0.00s cpu)
65</PRE>
66    In this example, the item/3 predicate serves as a nondeterministic
67    generator of solutions with different values for the variable Cal,
68    which we have chosen as our cost variable.  As can be seen from the
69    progress messages, the optimization procedure registers increasingly
70    good solutions (i.e. solutions with smaller cost), and finally delivers
71    the minimum-cost solution with Cal=200.
72</P><P>
73    Alternatively, we can minimize the item price:
74<PRE>
75        ?- bb_min(item(Food,Cal,Price), Price, bb_options{delta:0.05}).
76        Found a solution with cost 1.5
77        Found a solution with cost 0.99
78        Found no solution with cost -1.0Inf .. 0.94
79
80        Food = milk
81        Cal = 400
82        Price = 0.99
83        Yes (0.00s cpu)
84</PRE>
85    Because the price is non-integral, we had to adjust the step-width
86    of the optimization procedure using the delta-option.
87</P>
88<H3>Optimization with Constraints</H3>
89<P>
90    This library is designed to work together with arbitrary constraint
91    solvers, for instance library(ic).  The principle there is to wrap
92    the solver's nondeterministic search procedure into a bb_min/3 call.
93    This turns a program that finds all solutions into one that finds
94    the best solution.  For example:
95<PRE>
96        ?- [X,Y,Z] #:: 1..5,                   % constraints (model)
97           X+Z #>= Y,
98
99           C #= 3*X - 5*Y + 7*Z,               % objective function
100
101           bb_min(labeling([X,Y,Z]), C, _).    % nondet search + b&b
102
103        Found a solution with cost 5
104        Found a solution with cost 0
105        Found a solution with cost -2
106        Found a solution with cost -4
107        Found a solution with cost -6
108        Found no solution with cost -15.0 .. -7.0
109        X = 4
110        Y = 5
111        Z = 1
112        C = -6
113        Yes (0.00s cpu)
114</PRE>
115    The code shows the general template for such an optimization solver:
116    All constraints should be set up BEFORE the call to bb_min/3,
117    while the nondeterministic search procedure (here labeling/1)
118    must be invoked WITHIN bb_min/3.  The branch-and-bound procedure
119    only works if it envelops all nondeterminism.
120</P><P>
121    The cost variable (here C) must be defined in such a way that it is
122    instantiated (possibly by propagation) whenever the search procedure
123    succeeds with a solution.  Moreover, good, early bounds on the cost
124    variable are important for efficiency, as they help the branch-and-bound
125    procedure to prune the search.  Redundant constraints on the cost
126    variable can sometimes help.
127</P>
128
129<H3>Note on the treatment of bounded reals</H3>
130<P>
131    The library allows the cost to be instantiated to a number of type
132    breal.  This is useful e.g. when using lib(ic) to solve problems
133    with continuous variables.  When the variable domains have been
134    narrowed sufficiently, the problem variables (in particular the
135    cost variable) should be instantiated to a bounded real, e.g.
136    using the following idiom:
137    <PRE>
138	    X is breal_from_bounds(get_min(X),get_max(X))
139    </PRE>
140    Bounded reals contain some uncertainty about their true value. If
141    this uncertainty is too large, the branch-and-bound procedure may
142    not be able to compare the quality of two solutions. In this case,
143    a warning is issued and the search terminated prematurely.  The
144    problem can be solved by increasing the delta-parameter, or by
145    locating the cost value more precisely.
146</P>")).
147
148:- export
149	bb_min/3,		% main predicate
150	bb_min_cost/4,
151	bb_min/6,
152	minimize/2,
153
154	bb_init/2,		% underlying primitives
155	bb_cost/2,
156	bb_solution/2,
157	bb_finish/1,
158	bb_probe/7.
159
160:- meta_predicate((
161	bb_min(0,*,*),
162	bb_min(0,*,*,*,*,*),
163	minimize(0,*),
164	bb_probe(*,*,0,*,*,*,*)
165    )).
166
167:- export struct(bb_options(
168	    strategy,		% atom
169	    from,		% number
170	    to,			% number
171	    delta,		% number
172	    factor,		% number
173	    timeout,		% number
174	    probe_timeout,	% number
175	    solutions,		% atom
176	    report_success,	% GoalPrefix/A
177	    report_failure,	% GoalPrefix/A
178	    report_timeout	% GoalPrefix/A
179	)).
180
181:- lib(timeout).
182
183:- import
184	% general-purpose predicates
185	call_local/1,
186	par_true/0,
187	worker_boundary/0
188    from sepia_kernel.
189
190
191%----------------------------------------------------------------------
192% Toplevel
193%----------------------------------------------------------------------
194
195:- tool(minimize/2, minimize/3).
196minimize(Goal, Cost, Module) :-
197	bb_min(Goal, Cost, _DefaultOptions, Module).
198
199:- tool(bb_min/3, bb_min/4).
200bb_min(Goal, Cost, Options, Module) :-
201	% we used to use the Goal directly as the Template, but that can
202	% cause problems when Goal contains strange stuff, like handles
203	term_variables(Goal, Template),
204	bb_min(Goal, Cost, Template, Solution, Opt, Options, Module),
205	( Options = bb_options{solutions:one} ->
206	    % instantiate to the optimum solution (atomically)
207	    Template-Cost = Solution-Opt
208	;
209	    % impose Opt as a lower bound, because there could be
210	    % better solutions within the delta-range
211	    set_var_bounds(Cost, -1.0Inf, Opt),
212	    get_var_bounds(Cost, _, Max),
213	    % only if bound can't be imposed, instantiate instead
214	    ( Max =< Opt -> true ; Cost = Opt ),
215	    call(Goal)@Module
216	).
217
218:- tool(bb_min_cost/4, bb_min_cost_/5).
219bb_min_cost_(Goal, Cost, Opt, Options, Module) :-
220	bb_min(Goal, Cost, [], [], Opt, Options, Module).
221
222:- tool(bb_min/6, bb_min/7).
223bb_min(Goal, Cost, Template, Solution, Opt, Options, Module) :-
224	( var(Cost) ; number(Cost) ),
225	callable(Goal),
226	default_options(Options),
227	check_options(Options),
228	!,
229	initial_bounds(Cost, Options, From, To),
230	bb_init(To, Handle),
231	Options = bb_options{timeout:Timeout},
232        ( Timeout > 0, Timeout < 1.0Inf ->
233	    timeout(
234		bb_any(From, To, Goal, Template, Cost, Handle, Module, Options),
235		Timeout,
236		report_timeout(Options, From..To, Handle, Module)
237	    )
238	;
239	    bb_any(From, To, Goal, Template, Cost, Handle, Module, Options)
240	),
241	bb_solution(Handle, Solution),	% may fail
242	bb_cost(Handle, Opt),
243	bb_finish(Handle).
244bb_min(Goal, Cost, _Template, _Solution, _Opt, Options, Module) :-
245	error(5, bb_min(Goal, Cost, Options), Module).
246
247    bb_any(From, To, Goal, Template, Cost, Handle, Module, Options) :-
248	% The following line is for parallelism (scheduling and copying)
249	%worker_boundary, par_true,
250	Options = bb_options{strategy:Strategy},
251	( Strategy == continue ->
252	    bb_continue(From, To, Goal, Template, Cost, Handle, Module, Options)
253	; Strategy == restart ->
254	    bb_delta(From, To, Goal, Template, Cost, Handle, Module, Options)
255	; Strategy == step ->
256	    bb_delta(From, To, Goal, Template, Cost, Handle, Module, Options)
257	; Strategy == dichotomic ->
258	    bb_dichotomic(From, To, Goal, Template, Cost, Handle, Module, Options)
259	;
260	    error(6, bb_min(Goal, Cost, Options), Module)
261	).
262    bb_any(_From,_To,_Goal,_Template,_Cost,_Handle,_Module,_Options).
263
264
265default_options(bb_options{
266	    from:From,to:To,factor:Factor,strategy:Strategy,solutions:Solutions,
267	    delta:Delta,timeout:Timeout,probe_timeout:ProbeTimeout}) :-
268	set_default(Strategy, continue),
269	set_default(Solutions, one),
270	set_default(From, -1.0Inf),
271	set_default(To, 1.0Inf),
272	( Strategy==dichotomic -> set_default(Factor, 0.5)
273	; set_default(Factor, 1) ),
274	set_default(Delta, 1),
275	set_default(Timeout, 0),
276	set_default(ProbeTimeout, 0).
277
278    set_default(X, X) :- !.
279    set_default(_, _).
280
281check_options(bb_options{from:From,to:To,factor:Factor,solutions:Solutions,
282	     delta:Delta,timeout:Timeout,probe_timeout:ProbeTimeout}) :-
283	precise_number(From),
284	precise_number(To),
285	precise_number(Delta),
286	precise_number(Factor),
287	(Solutions==one;Solutions==all),
288	From =< To,
289	0 < Factor, Factor =< 1,
290	0 < Delta,
291	0 =< Timeout,
292	0 =< ProbeTimeout.
293
294precise_number(X) :- breal(X), !,
295	printf(error, "branch_and_bound: parameter must not be a bounded real: %w%n", [X]),
296	fail.
297precise_number(X) :- number(X).
298
299%----------------------------------------------------------------------
300% Stepwise improvement
301% bb_delta always fails or aborts!
302%----------------------------------------------------------------------
303
304bb_delta(From, To, Goal, Template, Cost, Handle, Module, Options) :-
305	( bb_probe(From, To, Goal, Template, Cost, Handle, Module,Options) ->	% may fail
306	    Best is bb_cost(Handle),
307	    step(From, Best, Options, NewTo),		% may fail
308	    bb_delta(From, NewTo, Goal, Template, Cost, Handle, Module, Options)
309	;
310	    report_failure(Options, From..To, Handle, Module),
311	    fail
312	).
313
314    % fails if there is no range left to explore
315    % PRE: no breals among the inputs, except Best
316    step(From, Best, bb_options{factor:Factor,delta:Delta}, NewTo) :-
317	To is number_max(Best),	% we are trying to improve the guaranteed best
318	Gap is To - From,
319	Gap > 0,		% termination condition
320	( Gap < 1.0Inf ->
321	    NewTo is min(From + Gap*Factor, To-Delta),
322	    NewTo < To		% can only be violated if precision problems
323	;
324	    NewTo is To - Delta
325	),
326	( NewTo < number_min(Best) ->
327	    true
328	;
329	    % this prevents looping with overlapping costs
330	    writeln(warning_output, "WARNING: bb_min: search terminated prematurely - cost uncertainty too large."),
331	    writeln(warning_output, "Either increase bb_min delta parameter, or compute cost more precisely."),
332	    fail
333	).
334
335    % like breal_max, but preserve the type for non-breals
336    number_max(X, Max) :- breal(X), !, breal_max(X, Max).
337    number_max(X, X).
338
339    % like breal_min, but preserve the type for non-breals
340    number_min(X, Min) :- breal(X), !, breal_min(X, Min).
341    number_min(X, X).
342
343%----------------------------------------------------------------------
344% Stepwise improvement without restart
345% bb_continue always fails or aborts!
346%
347% The search tree is only traversed once. Every time a solution is
348% found, the new cost bound gets imposed dynamically for the rest
349% of the search. To do that, we use the fail-event mechanism to trigger
350% a demon that imposes the current global bound on the cost variable.
351% This is necessary in two situations:
352% 1. If we fail across a change of the global bound (after a new solution)
353% 2. If we fail across the imposition of a new bound on the cost variable
354%    (we have to do it again because the effect is lost)
355% For the sake of completeness, we also take care of the case where
356% the cost variable is not a domain variable and cannot have a bound
357% imposed. In that case, the demon must additionally be woken on
358% instantiation of the cost variable to effect the check.
359%----------------------------------------------------------------------
360
361:- import
362	request_fail_event/3,
363	timestamp_init/2,
364	timestamp_update/2
365    from sepia_kernel.
366
367:- set_event_handler('branch_and_bound:apply_bound', trigger/1).
368
369
370bb_continue(From, To, Goal, Template, Cost, Handle, Module, Options) :-
371	(
372	    Stamp = stamp(_),
373	    timestamp_init(Stamp, 1),
374	    timestamp_update(Stamp, 1),
375	    suspend(bb_impose_bound(Cost, Handle, Stamp, From),
376		2,	% this should be between 1 (debug) and 2 (constraints)
377		trigger('branch_and_bound:apply_bound'), Susp),
378	    ( cannot_impose_bound(Cost) ->
379		insert_suspension(Cost, Susp, inst of suspend, suspend)
380	    ;
381		true
382	    ),
383	    call_local((
384		    % impose initial bounds, if possible
385		    set_var_bounds(Cost, From, To),
386		    Goal,
387		    kill_suspension(Susp)
388		))@Module,
389	    ( var(Cost) ->
390		writeln(error, "bb_min: search did not instantiate cost variable"),
391		abort
392	    ;
393		true
394	    ),
395	    % In case the solution still contains variables,
396	    % we want to strip most of their attributes.
397	    % Otherwise we might copy the whole constraint store!
398	    copy_term(Template, StrippedSolution),
399	    % Compute the new bound
400	    ( step(From, Cost, Options, NewTo) ->
401		% Set all shelf fields atomically (timeout!)
402		shelf_set(Handle, 0, sol(s(StrippedSolution),Cost,NewTo)),
403		request_fail_event(Stamp, 1, 'branch_and_bound:apply_bound')
404	    ;
405		% Set all shelf fields atomically (timeout!)
406		shelf_set(Handle, 0, sol(s(StrippedSolution),Cost,From)),
407	    	% We found an optimal solution, no need to continue search
408		!
409	    ),
410	    report_success(Options, Cost, Handle, Module),
411	    fail	% continue search, fail into Goal
412	;
413	    shelf_get(Handle, 3, NewTo),
414	    report_failure(Options, From..NewTo, Handle, Module),
415	    fail
416	).
417
418
419% This demon will be woken after failures when a global bound needs
420% to be re-imposed. It will itself schedule another event to be raised
421% in case we fail across the current demon invocation.
422% The bound is imposed using set_var_bounds/3. In case the Cost
423% variable does not have a domain which can represent that,
424% we need to wake also when Cost gets instantiated.
425
426:- demon bb_impose_bound/4.
427bb_impose_bound(Cost, _Handle, _Stamp, _From) :-
428	% Optimization: if Cost is a free variable, we can't impose bound
429	% neither now nor further up the search tree. No point scheduling
430	% further fail-events, need to wait for instantiation of Cost.
431	free(Cost), !.
432bb_impose_bound(Cost, Handle, Stamp, From) :-
433	% (number(Cost) ; meta(Cost))
434	% when we fail across this point, we will need to re-apply the bound!
435	request_fail_event(Stamp, 1, 'branch_and_bound:apply_bound'),
436	% this will apply the bound, if Cost is a domain variable or constat
437	NewTo is bb_bound(Handle),
438%	writeln(bb_impose_bound(Cost::From..NewTo)),
439	set_var_bounds(Cost, From, NewTo).
440
441
442% succeed if the argument cannot have a bound imposed
443% This fails for breals, although the bound cannot really be imposed if
444% it overlaps, but we take care of that elsewhere.
445
446cannot_impose_bound(X) :- free(X).
447cannot_impose_bound(X) :- meta(X),
448	get_var_bounds(X, L, H),
449	H =:= 1.0Inf,
450	L =:= -1.0Inf,
451	% try imposing a bound and check if it worked
452	call_priority((
453		set_var_bounds(X, L, 0),
454		get_var_bounds(X, _, H1),
455		H1 =:= 1.0Inf
456	    ), 1).
457
458
459%----------------------------------------------------------------------
460% Dichotomic search
461% bb_dichotomic always fails or aborts!
462%----------------------------------------------------------------------
463
464bb_dichotomic(From, To, Goal, Template, Cost, Handle, Module, Options) :-
465	% look for an initial solution
466	( bb_probe(From, To, Goal, Template, Cost, Handle, Module, Options) ->	% may fail
467	    Best is bb_cost(Handle),
468	    bb_dichotomic1(From, false, Best, Goal, Template, Cost, Handle, Module, Options)
469	;
470	    report_failure(Options, From..To, Handle, Module),
471	    fail
472	).
473
474    % We assume that there is a solution with cost Best.
475    % If FromNoSol==true then we know there is no solution with cost From
476    bb_dichotomic1(From, FromNoSol, Best, Goal, Template, Cost, Handle, Module, Options) :-
477	split(From, FromNoSol, Best, Options, Split), 	% may fail
478	( bb_probe(From, Split, Goal, Template, Cost, Handle, Module, Options) ->
479	    NewBest is bb_cost(Handle),
480	    bb_dichotomic1(From, FromNoSol, NewBest, Goal, Template, Cost, Handle, Module, Options)
481	;
482	    report_failure(Options, From..Split, Handle, Module),
483	    bb_dichotomic1(Split, true, Best, Goal, Template, Cost, Handle, Module, Options)
484	).
485
486    % PRE: no breals among the inputs, except Best
487    split(From, FromNoSol, Best, bb_options{factor:Factor,delta:Delta}, Split) :-
488	To is number_max(Best),	% we are trying to improve the guaranteed best
489	Gap is To - From,
490	( FromNoSol == true ->
491	    Gap > Delta		% termination condition
492	;
493	    Gap >= Delta	% termination condition
494	),
495	( Gap < 1.0Inf ->
496	    % Normally, Split is From + Gap*Factor, but we
497	    % have to consider all sorts of special cases...
498	    FromPlusDelta is From + Delta,
499	    ToMinusDelta is To - Delta,
500	    NormalSplit is From + Gap*Factor,
501	    ( FromPlusDelta > ToMinusDelta ->
502		% Gap size between Delta and 2*Delta: split in the middle,
503		% so both probe-range and improvement are at least Delta/2
504		Split is From + Gap*0.5
505	    ; NormalSplit < FromPlusDelta ->
506		% make sure the probe covers at least a Delta-range,
507		% otherwise proving optimality can take a long time
508	    	Split = FromPlusDelta
509	    ; NormalSplit > ToMinusDelta ->
510		% make sure the improvement is at least Delta,
511		% otherwise we might iterate through too many solutions
512	    	Split = ToMinusDelta
513	    ;
514		% normal case: split according to factor
515	    	Split = NormalSplit
516	    ),
517	    % Following line is normally redundant, but can be violated in
518	    % case of precision/rounding problems. It then prevents looping.
519	    Split < To
520	; From >= 0 ->
521	    Split is From + 10000000
522	; To =< 0 ->
523	    Split is To - 10000000
524	;
525	    Split = 0.0
526	),
527	( Split < number_min(Best) ->
528	    true
529	;
530	    % this prevents looping with overlapping costs
531	    writeln(warning_output, "WARNING: bb_min: search terminated prematurely - cost uncertainty too large."),
532	    writeln(warning_output, "Either increase bb_min delta parameter, or compute cost more precisely."),
533	    fail
534	).
535
536%----------------------------------------------------------------------
537% Primitives
538%----------------------------------------------------------------------
539
540bb_init(ExtremeCost, Handle) :-
541	shelf_create(sol(no_solution,ExtremeCost,ExtremeCost), Handle).
542
543bb_cost(Handle, Cost) :-
544	shelf_get(Handle, 2, Cost).
545
546bb_solution(Handle, Solution) :-
547	shelf_get(Handle, 1, s(Solution)).	% fail here if no solution
548
549bb_bound(Handle, Bound) :-
550	shelf_get(Handle, 3, Bound).
551
552bb_finish(Handle) :-
553	shelf_abolish(Handle).
554
555
556% bb_probe tries to find a solution for Goal in the range From..To.
557% If there is a solution, its Template and Cost are stored in Handle,
558% the computation is undone, and bb_probe succeeds.
559% If there is no solution, Handle is not changed and bb_probe fails.
560
561bb_probe(From, To, Goal, Template, Cost, Handle, Module) :-
562	bb_probe(From, To, Goal, Template, Cost, Handle, Module,
563		bb_options{report_success:true/0}).
564
565bb_probe(From, To, Goal, Template, Cost, Handle, Module, Options) :-
566	(
567	    call_local((
568		    % impose bounds early if possible
569		    set_var_bounds(Cost, From, To),
570		    Goal,
571		    % in case there was no set_bounds handler:
572		    set_var_bounds(Cost, From, To)
573		))@Module
574	->
575	    ( var(Cost) ->
576	    	writeln(error, "bb_min: search did not instantiate cost variable"),
577		abort
578	    ;
579		true
580	    ),
581	    % In case the solution still contains variables,
582	    % we want to strip most of their attributes.
583	    % Otherwise we might copy the whole constraint store!
584	    copy_term(Template, StrippedSolution),
585	    % Set all shelf fields atomically (timeout!)
586	    shelf_set(Handle, 0, sol(s(StrippedSolution),Cost,Cost)),
587	    report_success(Options, Cost, Handle, Module),
588	    fail	% to undo the call-effect and succeed with 2nd clause
589	;
590	    !,
591	    fail
592	).
593bb_probe(_From, _To, _Goal, _Template, _Cost, _Handle, _Module, _Options).
594
595
596% Get an initial lower and upper bound for the search by intersecting
597% the bounds from propagation with the user option input (if any)
598
599initial_bounds(Cost, Options, Lower, Upper) :-
600	Options = bb_options{from:UserFrom, to:UserTo},
601	get_var_bounds(Cost, CostL, CostU),
602	Lower is max(UserFrom, CostL),
603	Upper is min(UserTo, CostU).
604
605
606report_success(bb_options{report_success:Spec}, Cost, Handle, Module) :-
607	report_result(Spec, "Found a solution with cost %q%n", Cost, Handle, Module).
608
609report_failure(bb_options{report_failure:Spec}, Range, Handle, Module) :-
610	report_result(Spec, "Found no solution with cost %q%n", Range, Handle, Module).
611
612report_timeout(bb_options{report_timeout:Spec}, _InitialRange, Handle, Module) :-
613	bb_cost(Handle, Cost),
614	report_result(Spec, "Branch-and-bound timeout while searching for solution better than %q%n%b", Cost, Handle, Module).
615
616    report_result(Prefix/A, Msg, Cost, Handle, Module) ?-
617    	callable(Prefix),
618	!,
619	( A==0 -> call(Prefix)@Module
620	; A==1 -> call(Prefix,Cost)@Module
621	; A==2 -> call(Prefix,Cost,Handle)@Module
622	; A==3 -> call(Prefix,Cost,Handle,Module)@Module
623	; printf(log_output, Msg, Cost)
624	).
625    report_result(Prefix, _Msg, Range, Handle, Module) :-
626	callable(Prefix),
627	!,
628	call(Prefix, Range, Handle, Module)@Module.
629    report_result(_, Msg, Cost, _, _) :-
630	printf(log_output, Msg, Cost).
631
632
633%----------------------------------------------------------------------
634% Documentation
635%----------------------------------------------------------------------
636
637:- comment(minimize/2, [
638    summary:"Find a minimal solution using the branch-and-bound method",
639    desc:html("This is a shorthand for
640    	<PRE>
641	bb_min(+Goal, ?Cost, _DefaultOptions)
642	</PRE>
643	See bb_min/3 for details."),
644    template:"minimize(+Goal, ?Cost)",
645    see_also:[bb_min/3]]).
646
647
648:- comment(bb_min/3, [
649    summary:"Find one or all minimal solutions using the branch-and-bound method",
650    see_also:[bb_min/6],
651    desc:html("\
652	A solution of the goal <EM>Goal</EM> is found that minimizes
653	the value of <EM>Cost</EM>.  <EM>Cost</EM> should be a
654	variable that is affected, and eventually instantiated, by
655	<EM>Goal</EM>.  Usually, <EM>Goal</EM> is the search procedure
656	of a constraint problem and <EM>Cost</EM> is the variable
657	representing the cost.  The solution is found using the branch
658	and bound method:  as soon as a solution is found, it gets
659	remembered and the search is continued or restarted with an
660	additional constraint on the <EM>Cost</EM> variable which
661	requires the next solution to be better than the previous one.
662	Iterating this process finally yields an optimal solution.
663</P><P>
664	The possible options are
665	<DL>
666	<DT><STRONG>strategy:</STRONG></DT><DD>
667	    <DL>
668	    <DT>continue (default)</DT>
669	    	<DD>after finding a solution, continue search with the newly
670		found bound imposed on Cost</DD>
671	    <DT>restart</DT>
672	    	<DD>after finding a solution, restart the whole search with
673		the newly found bound imposed on Cost</DD>
674	    <DT>step</DT>
675	    	<DD>a synonym for 'restart'</DD>
676	    <DT>dichotomic</DT>
677	    	<DD>after finding a solution, split the remaining cost range
678		and restart search to find a solution in the lower sub-range.
679		If that fails, assume the upper sub-range as the remaining
680		cost range and split again.</DD>
681	    </DL>
682	    The new bound (or the split point, respectively), is computed
683	    from the current best solution, taking into account the
684	    parameters delta and factor.
685	    </DD>
686	<DT><STRONG>from:</STRONG></DT>
687	    <DD>number - an initial lower bound for the cost (default -1.0Inf).
688	    Only useful if Cost is not a domain variable.</DD>
689
690	<DT><STRONG>to:</STRONG></DT>
691	    <DD>number - an initial upper bound for the cost (default +1.0Inf).
692	    Only useful if Cost is not a domain variable.</DD>
693
694	<DT><STRONG>delta:</STRONG></DT>
695	    <DD>number - minimal absolute improvement required for each step
696	    (applies to all strategies). The default value of 1.0 is
697	    appropriate for integral costs.  Any solution that improves on
698	    the best solution by less than this value will be missed.</DD>
699
700	<DT><STRONG>factor:</STRONG></DT>
701	    <DD>number - minimal improvement ratio (with respect to the lower
702	    cost bound) for strategies 'continue' and 'restart' (default 1.0),
703	    or split factor for strategy 'dichotomic' (default 0.5)</DD>
704
705	<DT><STRONG>solutions:</STRONG></DT><DD>
706	    <DL>
707	    <DT>one (default)</DT><DD>
708		Compute one (of possibly multiple) optimal solutions.</DD>
709	    <DT>all</DT><DD>
710		Nondeterministically compute all optimal solutions.
711		This has a performance penalty, as the search is restarted
712		one more time after the optimum has been determined.</DD>
713	    </DL>
714	    Note the dependence on the delta-parameter: the costs of these
715	    solutions may deviate by less than delta from the true optimum.</DD>
716
717	<DT><STRONG>timeout:</STRONG></DT>
718	    <DD>number - maximum seconds of cpu time to spend (default: no limit)</DD>
719
720	<DT><STRONG>report_success:</STRONG></DT>
721	    <DD>GoalPrefix/N - this specifies a goal to be invoked whenever
722	    the branch-and-bound process finds a better solution.  GoalPrefix
723	    is a callable term (atom or compound) and N is an integer between
724	    0 and 3.  The invoked goal is constructed by adding N optional
725	    arguments to GoalPrefix: Cost, Handle and Module.  Cost is
726	    a float number representing the cost of the solution found,
727	    Handle is a handle as accepted by bb_cost/2 or bb_solution/2,
728	    and Module is the context module of the minimisation.
729	    To disable any reporting, choose report_success:true/0.
730	    The default handler prints a message to log_output.</DD>
731
732	<DT><STRONG>report_failure:</STRONG></STRONG></DT>
733	    <DD>GoalPrefix/N - this specifies a goal to be invoked whenever
734	    the branch-and-bound process cannot find a solution in a cost
735	    range.  GoalPrefix is a callable term (atom or compound) and
736	    N is an integer between 0 and 3.  The invoked goal is
737	    constructed by adding N optional arguments to GoalPrefix:
738	    Cost, Handle and Module.   Cost is a From..To structure
739	    representing the range of cost in which no solution could be found,
740	    Handle is a handle as accepted by bb_cost/2 or bb_solution/2,
741	    and Module is the context module of the minimisation.
742	    To disable any reporting, choose report_failure:true/0.
743	    The default handler prints a message to log_output.</DD>
744
745	<DT><STRONG>report_timeout:</STRONG></DT>
746	    <DD>GoalPrefix/N - this specifies a goal to be invoked when the
747	    branch-and-bound process times out.  GoalPrefix is a callable
748	    term (atom or compound) and N is an integer between 0 and 3.
749	    The invoked goal is constructed by adding N optional arguments
750	    to GoalPrefix: Cost, Handle and Module.  Cost is a float number
751	    representing the cost of the best solution found, Handle
752	    is a handle as accepted by bb_cost/2 or bb_solution/2,
753	    and Module is the context module of the minimisation.
754	    To disable any reporting, choose report_timeout:true/0.
755	    The default handler prints a message to log_output.</DD>
756	</DL>
757	The default options can be selected by passing a free variable as
758	the Options-argument. To specify other options, pass a bb_options-
759	structure in struct-syntax, e.g.
760	<PRE>
761	    bb_min(..., ..., bb_options{strategy:dichotomic, timeout:60})
762	</PRE>
763</P><P>
764	In order to maximize instead of minimizing, introduce a negated
765	cost variable in your model and minimize that instead, e.g.
766	<PRE>
767	    % maximize Profit
768	    Cost #= -Profit,
769	    bb_min(search(...), Cost, bb_options{}),
770	</PRE>
771</P>"),
772    args:["Goal":"The (nondeterministic) search goal",
773	"Cost":"A (usually numeric domain) variable representing the cost",
774	"Options":"A bb_options structure or variable"],
775    fail_if:"Goal has no solutions",
776    amode:(bb_min(+,?,+) is nondet),
777    eg:"
778% simple minimization with default options
779    ?- bb_min(member(X,[9,6,8,4,7,2,4,7]), X, Options).
780    Found a solution with cost 9
781    Found a solution with cost 6
782    Found a solution with cost 4
783    Found a solution with cost 2
784    Found no solution with cost -1.0Inf .. 1
785    X = 2
786    Options = bb_options(continue, -1.0Inf, 1.0Inf, 1, 1, 0, 0, _, _)
787    yes.
788
789% coarser granularity: faster, but missing the optimum
790    ?- bb_min(member(X,[9,6,8,4,7,2,4,7]), X, bb_options{delta:4}).
791    Found a solution with cost 9
792    Found a solution with cost 4
793    Found no solution with cost -1.0Inf .. 0
794    X = 4
795    yes.
796
797% alternative strategy based on bisecting the cost space
798    ?- bb_min(member(X,[99,60,80,40,70,30,70]), X,
799	    bb_options{strategy:dichotomic, from:0}).
800    Found a solution with cost 99
801    Found a solution with cost 40
802    Found no solution with cost 0.0 .. 20.0
803    Found a solution with cost 30
804    Found no solution with cost 20.0 .. 25.0
805    Found no solution with cost 25.0 .. 27.5
806    Found no solution with cost 27.5 .. 28.75
807    Found no solution with cost 28.75 .. 29.0
808    X = 30
809    yes.
810
811% examples with library(ic) constraints
812    ?- [X,Y,Z] :: 1..5,                    % constraints (model)
813       X+Z #>=Y,
814       C #= 3*X - 5*Y + 7*Z,               % objective function
815       bb_min(labeling([X,Y,Z]), C, _).    % nondet search + b&b
816
817    Found a solution with cost 5
818    Found a solution with cost 0
819    Found a solution with cost -2
820    Found a solution with cost -4
821    Found a solution with cost -6
822    Found no solution with cost -15.0 .. -7.0
823    X = 4
824    Y = 5
825    Z = 1
826    C = -6
827    Yes (0.00s cpu)
828
829
830    ?- [X,Y,Z] :: 1..5,
831       X+Z #>=Y,
832       C #= 3*X - 5*Y + 7*Z,
833       bb_min(search([X,Y,Z],0,input_order,indomain_middle,complete,[]), C, _).
834
835    Found a solution with cost 15
836    Found a solution with cost 8
837    Found a solution with cost 1
838    Found a solution with cost -4
839    Found a solution with cost -6
840    Found no solution with cost -15.0 .. -7.0
841    X = 4
842    Y = 5
843    Z = 1
844    C = -6
845    Yes (0.00s cpu)
846
847
848"]).
849
850
851:- comment(bb_min_cost/4, [
852    summary:"Find the minimal cost using the branch-and-bound method",
853    args:["Goal":"The (nondeterministic) search goal",
854	"Cost":"A (usually numeric domain) variable representing the cost",
855	"Optimum":"A variable which will be set to the optimum value of Cost",
856	"Options":"A bb_options structure or variable"],
857    see_also:[bb_min/3,bb_min/6],
858    desc:html("<P>\
859    Determines the minimum possible value that the variable Cost can
860    attain in any solution of Goal.  This value is returned as Optimum.
861    Neither Cost nor any variable in Goal is instantiated.  The predicate
862    is useful when one is interested in sub-optimal solutions, see example.
863</P><P>
864    This predicate can be defined as
865<PRE>
866	bb_min_cost(Goal, Cost, Optimum, Options) :-
867		bb_min(Goal, Cost, [], _, Optimum, Options).
868</PRE>
869    Options are interpreted in the same way as for bb_min/6
870    (the solutions-parameter is ignored).
871</P>"),
872    fail_if:"Goal has no solutions",
873    amode:(bb_min_cost(+,?,-,+) is semidet),
874    eg:"
875    % A predicate to enumerate solutions in increasing cost order
876    :- lib(ic).
877    :- lib(branch_and_bound).
878
879    ic_increasing_cost(Goal, Cost) :-
880    	bb_min_cost(Goal, Cost, Opt,
881		    bb_options{report_success:true/0,report_failure:true/0}),
882	(
883	    Cost = Opt,
884	    call(Goal)
885	;
886	    Cost #> Opt,
887	    ic_increasing_cost(Goal, Cost)
888	).
889
890    % sample run:
891    ?- ic_increasing_cost(member(C-X,[9-a,4-b,2-c,4-d]), C).
892    C = 2
893    X = c
894    Yes (0.00s cpu, solution 1, maybe more) ? ;
895    C = 4
896    X = b
897    Yes (0.00s cpu, solution 2, maybe more) ? ;
898    C = 4
899    X = d
900    Yes (0.00s cpu, solution 3, maybe more) ? ;
901    C = 9
902    X = a
903    Yes (0.00s cpu, solution 4, maybe more) ? ;
904    No (0.00s cpu)
905    "
906    ]).
907
908:- comment(bb_min/6, [
909    summary:"Find a minimal solution using the branch-and-bound method",
910    see_also:[bb_min/3],
911    desc:html("\
912	A solution of the goal <EM>Goal</EM> is found that minimizes
913	the value of <EM>Cost</EM>.  <EM>Cost</EM> should be a
914	variable that is affected, and eventually instantiated, by
915	<EM>Goal</EM>.  Usually, <EM>Goal</EM> is the search procedure
916	of a constraint problem and <EM>Cost</EM> is the variable
917	representing the cost.  The solution is found using the branch
918	and bound method:  as soon as a solution is found, it gets
919	remembered and the search is continued or restarted with an
920	additional constraint on the <EM>Cost</EM> variable which
921	requires the next solution to be better than the previous one.
922	Iterating this process yields an optimal solution in the end.
923	<P>
924	The possible options are
925	<DL>
926	<DT><STRONG>strategy:</STRONG></DT><DD>
927	    <DL>
928	    <DT>continue (default)</DT>
929	    	<DD>after finding a solution, continue search with the newly
930		found bound imposed on Cost</DD>
931	    <DT>restart</DT>
932	    	<DD>after finding a solution, restart the whole search with
933		the newly found bound imposed on Cost</DD>
934	    <DT>step</DT>
935	    	<DD>a synonym for 'restart'</DD>
936	    <DT>dichotomic</DT>
937	    	<DD>after finding a solution, split the remaining cost range
938		and restart search to find a solution in the lower sub-range.
939		If that fails, assume the upper sub-range as the remaining
940		cost range and split again.</DD>
941	    </DL>
942	    The new bound (or the split point, respectively), are computed
943	    from the current best solution, taking into account the
944	    parameters delta and factor.
945	    </DD>
946	<DT><STRONG>from:</STRONG></DT>
947	    <DD>number - an initial lower bound for the cost (default -1.0Inf)</DD>
948
949	<DT><STRONG>to:</STRONG></DT>
950	    <DD>number - an initial upper bound for the cost (default +1.0Inf)</DD>
951
952	<DT><STRONG>delta:</STRONG></DT>
953	    <DD>number - minimal absolute improvement required for each step
954	    (default 1.0), applies to all strategies</DD>
955
956	<DT><STRONG>factor:</STRONG></DT>
957	    <DD>number - minimal improvement ratio (with respect to the lower
958	    cost bound) for strategies 'continue' and 'restart' (default 1.0),
959	    or split factor for strategy 'dichotomic' (default 0.5)</DD>
960
961	<DT><STRONG>timeout:</STRONG></DT>
962	    <DD>number - maximum seconds of cpu time to spend (default: no limit)</DD>
963
964	<DT><STRONG>report_success:</STRONG></DT>
965	    <DD>GoalPrefix/N - this specifies a goal to be invoked whenever
966	    the branch-and-bound process finds a better solution.  GoalPrefix
967	    is a callable term (atom or compound) and N is an integer between
968	    0 and 3.  The invoked goal is constructed by adding N optional
969	    arguments to GoalPrefix: Cost, Handle and Module.  Cost is
970	    a float number representing the cost of the solution found,
971	    Handle is a handle as accepted by bb_cost/2 or bb_solution/2,
972	    and Module is the context module of the minimisation.
973	    To disable any reporting, choose report_success:true/0.
974	    The default handler prints a message to log_output.</DD>
975
976	<DT><STRONG>report_failure:</STRONG></STRONG></DT>
977	    <DD>GoalPrefix/N - this specifies a goal to be invoked whenever
978	    the branch-and-bound process cannot find a solution in a cost
979	    range.  GoalPrefix is a callable term (atom or compound) and
980	    N is an integer between 0 and 3.  The invoked goal is
981	    constructed by adding N optional arguments to GoalPrefix:
982	    Cost, Handle and Module.   Cost is a From..To structure
983	    representing the range of cost in which no solution could be found,
984	    Handle is a handle as accepted by bb_cost/2 or bb_solution/2,
985	    and Module is the context module of the minimisation.
986	    To disable any reporting, choose report_failure:true/0.
987	    The default handler prints a message to log_output.</DD>
988
989	<DT><STRONG>report_timeout:</STRONG></DT>
990	    <DD>GoalPrefix/N - this specifies a goal to be invoked when the
991	    branch-and-bound process times out.  GoalPrefix is a callable
992	    term (atom or compound) and N is an integer between 0 and 3.
993	    The invoked goal is constructed by adding N optional arguments
994	    to GoalPrefix: Cost, Handle and Module.  Cost is a float number
995	    representing the cost of the best solution found, Handle
996	    is a handle as accepted by bb_cost/2 or bb_solution/2,
997	    and Module is the context module of the minimisation.
998	    To disable any reporting, choose report_timeout:true/0.
999	    The default handler prints a message to log_output.</DD>
1000	</DL>
1001	The default options can be selected by passing a free variable as
1002	the Options-argument. To specify other options, pass a bb_options-
1003	structure in struct-syntax, e.g.
1004	<PRE>
1005	bb_options{strategy:dichotomic, timeout:60}
1006	</PRE>
1007	In order to maximize instead of minimizing, introduce a negated
1008	cost variable in your model and minimize that instead.
1009	<P>
1010	Unlike bb_min/3, bb_min/6 does <STRONG>not</STRONG> affect Goal or Cost after
1011	the optimum has been found. Instead, the optimum cost value is returned
1012	in Optimum, and the Solution argument gets unified with an instance of
1013	Template where the variables have the values that correspond to the
1014	optimal solution. Note that bb_min/3 is actually based on bb_min/6
1015	and can (for the one-solution case) be defined as:
1016	<PRE>
1017	bb_min(Goal, Cost, Options) :-
1018	    bb_min(Goal, Cost, Goal, Goal, Cost, Options).
1019	</PRE>
1020	"),
1021    args:["Goal":"The (nondeterministic) search goal",
1022	"Cost":"A (usually numeric domain) variable representing the cost",
1023	"Template":"A term containing all or some problem variables",
1024	"Solution":"A term which will be unified with the optimized Template",
1025	"Optimum":"A variable which will be set to the optimum value of Cost",
1026	"Options":"A bb_options structure or variable"],
1027    fail_if:"Goal has no solutions",
1028    amode:(bb_min(+,?,?,?,?,?) is semidet)
1029    ]).
1030
1031:- comment(bb_init/2, [
1032    summary:"Low-level primitive for building branch-and-bound-style search procedures",
1033    template:"bb_init(++ExtremeCost, -Handle)",
1034    see_also:[bb_probe/7]
1035    ]).
1036
1037:- comment(bb_cost/2, [
1038    summary:"Low-level primitive for building branch-and-bound-style search procedures",
1039    template:"bb_cost(++Handle, -Cost)",
1040    see_also:[bb_probe/7]
1041    ]).
1042
1043:- comment(bb_solution/2, [
1044    summary:"Low-level primitive for building branch-and-bound-style search procedures",
1045    template:"bb_solution(++Handle, -Solution)",
1046    see_also:[bb_probe/7]
1047    ]).
1048
1049:- comment(bb_finish/1, [
1050    summary:"Low-level primitive for building branch-and-bound-style search procedures",
1051    template:"bb_finish(++Handle)",
1052    see_also:[bb_probe/7]
1053    ]).
1054
1055:- comment(bb_probe/7, [
1056    summary:"Low-level primitive for building branch-and-bound-style search procedures",
1057    template:"bb_probe(++From, ++To, +Goal, ?Template, ?Cost, ++Handle, ++Module)",
1058    desc:html("
1059	bb_probe tries to find a solution for Goal in the range From..To.
1060	If there is a solution, its Template and Cost are stored in Handle,
1061	the computation is undone, and bb_probe succeeds.
1062	If there is no solution, Handle is not changed and bb_probe fails.
1063	The primitive set_var_bounds/3 is used to impose cost bounds
1064	during the search process in a generic way."),
1065    see_also:[bb_init/2,bb_cost/2,bb_solution/2,bb_finish/1,bb_min/3,bb_min/6,
1066    	set_var_bounds/3],
1067    eg:"% a simple branch-and-bound procedure
1068my_minimize(Goal, Cost, Solution, OptCost, Module) :-
1069	bb_init(1000000, Handle),
1070	(
1071	    bb_delta(0, 1000000, Goal, Cost, Handle, Module)
1072	;
1073	    bb_solution(Handle, Solution),
1074	    bb_cost(Handle, OptCost)
1075	),
1076	bb_finish(Handle).
1077
1078bb_delta(From, To, Goal, Cost, Handle, Module) :-
1079	bb_probe(From, To, Goal, Goal, Cost, Handle, Module),
1080	NewTo is bb_cost(Handle) - 1,
1081	bb_delta(From, NewTo, Goal, Cost, Handle, Module).
1082    "]).
1083
1084