1% BEGIN LICENSE BLOCK
2% Version: CMPL 1.1
3%
4% The contents of this file are subject to the Cisco-style Mozilla Public
5% License Version 1.1 (the "License"); you may not use this file except
6% in compliance with the License.  You may obtain a copy of the License
7% at www.eclipse-clp.org/license.
8% 
9% Software distributed under the License is distributed on an "AS IS"
10% basis, WITHOUT WARRANTY OF ANY KIND, either express or implied.  See
11% the License for the specific language governing rights and limitations
12% under the License. 
13% 
14% The Original Code is  The ECLiPSe Constraint Logic Programming System. 
15% The Initial Developer of the Original Code is  Cisco Systems, Inc. 
16% Portions created by the Initial Developer are
17% Copyright (C) 1995 - 2009 Cisco Systems, Inc.  All Rights Reserved.
18% 
19% Contributor(s): 	Helmut Simonis, 4C, Univerity College Cork, Cork
20%			Kish Shen
21% 
22% END LICENSE BLOCK
23% ----------------------------------------------------------------------
24
25:-export(alldifferent/1).
26:-export(matching/2).
27:-export(gcc/2).
28:-export(same/2).
29:-export(inverse/2).
30:-export(lex_le/2).
31:-export(lex_lt/2).
32
33%----------------------------------------------------------------------
34% Output transformations
35%----------------------------------------------------------------------
36
37:-pragma(nodebug).
38
39:-comment(tr_global_gac_out/2,hidden).
40:- export tr_global_gac_out/2.
41
42:- export portray(update_alldifferent/4, tr_global_gac_out/2, [goal]).
43
44tr_global_gac_out(update_alldifferent(_, _,List,_), alldifferent(List)).
45
46
47:- import flow_constraints_support.
48
49/***************************************************************
50
51Structure definitions
52
53***************************************************************/
54        
55% common structure for all constraints
56%:-local struct(low_high(low,high,mapping)).
57:-local struct(mapping(n,
58                       value2node,
59                       var_array)).
60
61:-local struct(remember(matching,
62                        solved)).
63
64/*******************************************************************************
65
66alldifferent
67
68*******************************************************************************/
69
70:-comment(alldifferent/1,[summary:"GAC version of alldifferent",
71                          amode:alldifferent(+),
72                          args:["L":"List of integers or domain"
73                                    " variables, or a collection a la collection_to_list/2"],
74                          kind:[constraint:[root:[ic,fd]]],
75                          desc:html("This predicate implements a GAC"
76                                    " (generalized arc consistency)"
77                                    " version of the alldifferent"
78                                    " constraint. It uses the"
79                                    " classical bitpartite matching"
80                                    " implementation using the"
81                                    " graph_algorithms library. This"
82                                    " version often removes more"
83                                    " values than the bound consistent"
84                                    " alldifferent in the ic_global"
85                                    " library, or the forward checking"
86                                    " variant in the ic library, but"
87                                    " may spend much more time doing"
88                                    " this."),
89                          fail_if:"fails if there is no bipartite matching"
90                                  " between all variables and the"
91                                  " possible values",
92                          see_also:[matching/2,
93                                    ic:alldifferent/1,
94                                    ic_global:alldifferent/1]
95                         ]).
96
97:-comment(matching/2,[summary:"Get a matching between a list of"
98                             " domain variables and their possible"
99                              " values",
100                      amode:matching(+,-),
101                      args:["L":"A list of integers or domain"
102                                " variables, or a collection a la collection_to_list/2",
103                            "K":"A free variable, will be bound to a"
104                                " list of integers"
105                           ],
106                      fail_if:"The predicate fails if no matching"
107                              " exists",
108                      desc:html("This predicate can be used to get the"
109                                " matching into the user program."
110                                " Sometimes it is a good starting"
111                                " point for heuristics. It only gets"
112                                " the current matching and does not do"
113                                " any propagation, the matching is not"
114                                " updated when values are removed, the"
115                                " predicate must be called again in"
116                                " the user program if this is"
117                                " required"),
118                      see_also:[alldifferent/1]
119                     ]).
120
121% this is the top-level entry point
122% it sets up a delay on each variable, and then calls the checker
123alldifferent(Vars):-
124        collection_to_list(Vars,L),
125        init_remember(Remember),
126        call_priority((
127            check_alldifferent(L,Remember),
128            shrink(L,L1),
129            (L1 = [] ->
130                true
131            ; L1 = [_] ->
132                true
133            ;
134                suspend(update_alldifferent(L1,Remember,L,Susp),4,[L1->any],Susp)
135            )
136        ), 2).
137
138
139shrink(L,L1):-
140        (foreach(X,L),
141         fromto(L1,A1,A,[]) do
142            (var(X) ->
143                A1 = [X|A]
144            ;
145                A1 = A
146            )
147        ).
148
149% this is called whenever a variable changes
150% it removes the constraint once it is solved
151% the third argument is passed to allow output of original argument
152:-demon(update_alldifferent/4).
153update_alldifferent(L,Remember,_OrigArg,Susp):-  
154        check_alldifferent(L,Remember),
155        (is_solved(Remember) ->
156            kill_suspension(Susp)
157        ;
158            unschedule_suspension(Susp) 
159        ).
160
161% this is the actual constraint check code
162check_alldifferent(_L,Remember):-
163        is_solved(Remember),
164        !.
165check_alldifferent(L,Remember):-
166        length(L,N),
167%        writeln(check),
168        % create a list of the variable nodes
169        % create a list of edges
170        create_nodes_and_edges(L,N,Mapping,VarNodes,Edges,LastNode,Ground),
171        (Ground = 1 ->
172            mark_solved(Remember),
173%            writeln(ground),
174            true
175        ;
176            N1 is N+1,
177            % build a list of the value nodes in the graph
178            value_nodes(N1,LastNode,ValueNodes),
179            % test if there are at least as many values as variables
180            % if not the constraint must fail here
181            LastNode - N >= N,
182            % create graph data structure
183            make_graph(LastNode,Edges,Graph),
184            remembered_matching(Mapping,L,Remember,OldMatching,
185                                SizeOld),
186%            writeln(old(SizeOld,OldMatching)),
187            % call the matching algorithm, if old matching does not work
188            (SizeOld = N ->
189                MaximalM = OldMatching
190            ;
191                maximum_matching_hopcroft_karp(Graph, VarNodes,
192                                               ValueNodes, 
193                                               MaximalM),
194           
195                length(MaximalM,SizeM),
196%                writeln(matching(SizeM,MaximalM)),
197                % check that every variable is matched, otherwise fail
198                SizeM = N,
199                remember_matching(N,MaximalM,Remember,Mapping)
200            ),
201            % invert the edges not in the matching
202            invert_edges(Edges,MaximalM,InvertedEdges,Edges1),
203            % build a new graph with matching and inverted edges
204            make_graph(LastNode,Edges1,Graph1),
205            % search for strongly connected components in new graph
206            strong_components(Graph1, StrongComponents),
207%            writeln(strong(StrongComponents)),
208            % extract edges between nodes in different components
209            mark_components(StrongComponents,LastNode,InvertedEdges,
210                        NotMarked),
211%            writeln(not_marked(NotMarked)),
212            % find value nodes not in matching, they are possible starts
213            % of alternating paths; for permutations there should be none
214            unmatched_value_nodes(N1,LastNode,MaximalM,MFree),
215            % MFree is a list of unmatched value nodes
216            % find edges on alternate paths starting in unmatched values
217            % and mark them
218%            writeln(mfree(MFree)),
219            alternate_paths(MFree,Graph1,NotMarked,FinalNotMarked),
220            % remove the values which correspond to unmarked edges in the
221            % graph; note that values are the sources of edges
222            remove_unmarked_edges(FinalNotMarked,Mapping)
223        ).
224
225
226
227create_nodes_and_edges(L,N,Mapping,VarNodes,Edges,LastNode,Ground):-
228        mapping_create(N,Mapping),
229        hash_create(UniqueValues),
230        (foreach(X,L),
231         foreach(VarNode,VarNodes),
232         fromto([],A,A1,Edges),
233         fromto(1,B,B1,Ground),
234         fromto(N,NextV,NextV1,LastNode),
235         count(VarNode,1,_),
236         param(Mapping,UniqueValues) do
237            % create a node for each variable
238            node_map(Mapping,VarNode,X),
239            (integer(X) ->
240                B1 = B,
241                (hash_find(UniqueValues,X,_) ->
242                    fail
243                ;
244                    hash_add(UniqueValues,X,X)
245                )
246            ;
247                B1  = 0
248            ),
249            get_full_domain_as_list(X,Dom),
250            % create an edge for every value in the domain of a variable
251            (foreach(V,Dom),
252             fromto(NextV,NextValue,NextValue1,NextV1),
253             fromto(A,B,[e(VarNode,ValueNode,V)|B],A1),
254             param(Mapping,VarNode) do
255                % create a node for the value, if not already there
256                value_insert(Mapping,V,ValueNode,NextValue,NextValue1)
257            )
258        ).
259
260value_nodes(N1,LastNode,ValueNodes):-
261        (for(J,N1,LastNode),
262         foreach(J,ValueNodes) do
263            true
264        ).
265
266unmatched_value_nodes(N1,LastNode,MaximalM,MFree):-
267        dim(MNodes,[LastNode]),
268        (foreach(e(_From,To,_),MaximalM),
269         param(MNodes) do
270            subscript(MNodes,[To],1)
271        ),
272        (for(J,N1,LastNode),
273         fromto([],A,A1,MFree),
274         param(MNodes) do
275            subscript(MNodes,[J],V),
276            (var(V) ->
277                A1 = [J|A]
278            ;
279                A1 = A
280            )
281        ).
282
283% special case for no unmatched value nodes = permutation
284alternate_paths([],_Graph,NotMarked,NotMarked):-
285        !.
286alternate_paths(MFree,Graph,NotMarked,FinalNotMarked):-
287        % hash table with nodes as key, dummy values
288        hash_create(CheckedNodes),
289        % hash table with Edges as key, dummy values 
290        hash_create(MarkedEdges),
291        scan_edges(MFree,CheckedNodes,Graph,MarkedEdges),
292        reduce_unmarked_edges(NotMarked,MarkedEdges,FinalNotMarked).
293
294scan_edges([],_CheckedNodes,_Graph,_MarkedEdges).
295scan_edges([H|T],CheckedNodes,Graph,MarkedEdges):-
296        hash_create(NewNodes),
297        scan([H|T],NewNodes,CheckedNodes,Graph,MarkedEdges),
298        hash_list(NewNodes,NewNodesList,_),
299        scan_edges(NewNodesList,CheckedNodes,Graph,MarkedEdges).
300
301scan([],_NewNodes,_CheckedNodes,_Graph,_MarkedEdges).
302scan([H|T],NewNodes,CheckedNodes,Graph,MarkedEdges):-
303        graph_get_adjacent_edges(Graph,H,Edges),
304        scan_outbound(Edges,NewNodes,CheckedNodes,MarkedEdges),
305        hash_add(CheckedNodes,H,H), 
306        scan(T,NewNodes,CheckedNodes,Graph,MarkedEdges).
307
308scan_outbound([],_NewNodes,_CheckedNodes,_MarkedEdges).
309scan_outbound([Edge|Edge1],NewNodes,CheckedNodes,MarkedEdges):-
310        consider_edge(Edge,NewNodes,CheckedNodes,MarkedEdges),
311        scan_outbound(Edge1,NewNodes,CheckedNodes,MarkedEdges).
312
313consider_edge(Edge,NewNodes,CheckedNodes,MarkedEdges):-
314        Edge = e(_From,To,_Info),
315        % if node has been checked before or
316        % node is already on ToDo list, do not add
317        ((hash_find(CheckedNodes,To,_);hash_find(NewNodes,To,_)) ->
318            true
319        ;
320            hash_add(NewNodes,To,To)
321        ),
322        hash_add(MarkedEdges,Edge,0).
323
324% if edge is in hash table, then it is marked, and should not be removed
325reduce_unmarked_edges(L,Hash,K):-
326        (foreach(Edge,L),
327         fromto(K,A1,A,[]),
328         param(Hash) do
329            (hash_find(Hash,Edge,_) ->
330                A1 = A
331            ;
332                A1 = [Edge|A]
333            )
334        ).
335
336% remove values from variables corresponding to unmarked edges
337remove_unmarked_edges(FinalNotMarked,Mapping):-
338        (foreach(e(_ValueN,VarN,Value),FinalNotMarked),
339         param(Mapping) do
340            node_map(Mapping,VarN,Var),
341%           writeln(rem(VarN,Var,ValueN,Value)),
342            excl(Var,Value)
343        ).
344/*
345
346Dealing with var/value <-> node mapping
347vars are numbered 1..n
348values are numbered n+1..lastnode
349
350this allows for arbitrary values in the domains
351we can not use nodenames for this as we deal with variables
352*/
353
354mapping_create(N,mapping{n:N,
355                         value2node:Value2Node,
356                         var_array:VarArray}):-
357        hash_create(Value2Node),
358        N1 is N+1,
359        dim(VarArray,[N1]).
360
361node_map(mapping{var_array:VarArray},Node,Var):-
362        subscript(VarArray,[Node],Var).
363
364% find the node for a given value
365value_insert(mapping{value2node:Hash},V,ValueNode,Next,Next):-
366        hash_find(Hash,V,ValueNode),
367        !.
368value_insert(mapping{value2node:V2N},V,Next1,Next,Next1):-
369        Next1 is Next+1,
370        hash_add(V2N,V,Next1).
371
372value_lookup(mapping{value2node:Hash},V,ValueNode):-
373        hash_find(Hash,V,ValueNode),
374        !.
375
376
377matching(L,K):-
378        length(L,N),
379        % create a list of the variable nodes
380        % create a list of edges
381        create_nodes_and_edges(L,N,_Mapping,VarNodes,Edges,LastNode,_),
382        N1 is N+1,
383        % build a list of the value nodes in the graph
384        value_nodes(N1,LastNode,ValueNodes),
385        % test if there are at least as many values as variables
386        % if not the constraint must fail here
387        LastNode - N >= N,
388        % create graph data structure
389        make_graph(LastNode,Edges,Graph),
390        % call the matching algorithm
391        maximum_matching_hopcroft_karp(Graph, VarNodes,ValueNodes, MaximalM),
392        length(MaximalM,SizeM),
393        % check that every variable is matched, otherwise fail
394        SizeM = N,
395        sort(1,=<,MaximalM,Sorted),
396        (foreach(e(_From,_To,Value),Sorted),
397         foreach(Value,K) do
398            true
399        ).
400
401/*
402invert the edges in a list which are not in matching
403uses hash table to avoid list lookup
404also returns all edges in new graph
405*/
406invert_edges(Edges,MaximalM,InvertedEdges,AllEdges):-
407        hash_create(Match),
408        (foreach(Edge,MaximalM),
409         param(Match) do
410            hash_add(Match,Edge,1)
411        ),
412        (foreach(Edge,Edges),
413         fromto([],A,A1,InvertedEdges),
414         fromto(MaximalM,B,B1,AllEdges),
415         param(Match) do
416            direct_edge(Edge,Match,A,A1,B,B1)
417        ).
418
419% ignore edge if it is in the matching
420% otherwise invert its direction and put into accumulator
421direct_edge(Edge,Match,A,A,B,B):-
422        hash_find(Match,Edge,_),
423        !.
424direct_edge(e(From,To,W),_Match,A,[e(To,From,W)|A],
425            B,[e(To,From,W)|B]).
426
427
428init_remember(remember{}).
429
430remembered_matching(_Mapping,_L,remember{matching:Array},[],0):-
431        var(Array),
432        !.
433remembered_matching(Mapping,L,remember{matching:Array},
434                    OldMatching,
435                    SizeOld):-
436        
437        (foreach(X,L),
438         count(J,1,_),
439         fromto(OldMatching,A1,A,[]),
440         fromto(0,B,B1,SizeOld),
441         param(Mapping,Array) do
442            subscript(Array,[J],Old),
443            (check_in(Old,X) ->
444                value_lookup(Mapping,Old,To),
445                A1 = [e(J,To,Old)|A],
446                B1 is B+1
447            ;
448                A1 = A,
449                B1 = B
450            )
451        ).
452
453remember_matching(N,MaximalM,Remember,_Mapping):-
454        N1 is N+1,
455        dim(OldMatch,[N1]),
456        (foreach(e(Node,_,Value),MaximalM),
457         param(OldMatch) do
458            subscript(OldMatch,[Node],Value)
459        ),
460        setarg(1,Remember,OldMatch).
461
462mark_solved(remember{solved:1}).
463
464is_solved(remember{solved:Flag}):-
465        integer(Flag).
466
467
468
469/***************************************************************
470
471GCC with fixed bounds
472
473***************************************************************/
474
475:- comment(gcc/2, [
476        amode: gcc(++,+),
477        args: ["Bounds":"A list of elements of the form gcc(Low,High,Value),"
478                        " where Low, High and Value are integers, and High and"
479                        " Low are non-negative (High >= Low), and Value must"
480                        " be different from other Values",
481               "Vars":"A collection of different variables or integers"
482              ],
483        summary:"Constrain the cardinality of each Value specified in Bound's"
484                " gcc(Low,High,Value) to be between Low and High in Vars",
485        kind:[constraint:[root:[ic,fd]]],
486        desc:html("\
487<P>
488    This constraint ensures that the cardinality (the number of occurrences)
489    of values in Vars conforms to the specifications in Bounds. Bounds is a
490    list of triples in the form gcc(Low,High,Value) where Value is an integer,
491    a value that Vars is to be assigned to, and must occur only once as a
492    Value in Bounds, and whose cardinality |Value| is specified by 
493    Low =< |Value| =< High, where Low and High are non-negative integers.
494    Note that all values that Vars can take must be specified in Bounds.
495</P><P>
496    This is currently a prototype -- the constraint has not been tested
497    very extensively and little effort has been spent to optimise performance.
498    We welcome any feedback on using this constraint.
499</P><P>
500    This constraint is known as global_cardinality_low_up in the global
501    constraint catalog. The algorithm implemented is described in 
502    J.-C. Regin's paper 'Generalized Arc Consistency for Global Cardinality
503    Constraint', published in AAAI-1996. 
504")
505                  ]).
506
507gcc(Bounds,Vars):-
508        collection_to_list(Vars,Variables),
509        length(Variables,N),
510        length(Bounds,M),
511        call_priority((
512            check_gcc(Bounds,Variables,N,M),
513            (ground(Variables) ->
514                true
515            ;
516                suspend(update_gcc(Bounds,Variables,N,M,Susp),9,
517                        [Variables->[inst,any]],Susp)
518
519            )
520         ), 2).
521
522:-demon(update_gcc/5).
523update_gcc(Bounds,Variables,N,M,Susp):-
524        check_gcc(Bounds,Variables,N,M),
525        (ground(Variables) ->
526            kill_suspension(Susp)
527        ;
528            unschedule_suspension(Susp)
529        ).
530
531check_gcc(Bounds,Variables,N,M):-
532        create_gcc_edges(Variables,Bounds,N,M,
533                          Mapping,
534                          SourceNode,SinkNode,Edges),
535%        writeln(Mapping),
536        make_graph(SinkNode,Edges,Graph),
537        feas_flow_with_lb(Graph, low of low_high, high of low_high,
538                          SourceNode, SinkNode, 
539                          MaxFlowValue, MaxFlowEdges, _),
540%        writeln(N-MaxFlowValue-MaxFlowEdges),
541        MaxFlowValue >= N, % may fail
542%        is_feasible_flow(MaxFlowEdges,Edges), % testing only, remove later
543        residual_graph(MaxFlowEdges,Edges,ResidualEdges,MidEdges,_),
544        make_graph(SinkNode,ResidualEdges,ResidualGraph),
545        strong_components(ResidualGraph,StrongComponents),
546%        writeln(strong(StrongComponents)),
547        mark_components(StrongComponents,SinkNode,MidEdges,NotMarked),
548%        writeln(unmarked(NotMarked)),
549        gcc_remove_unmarked_edges(NotMarked,Mapping).
550
551create_gcc_edges(Variables,Bounds,N,M,Mapping,
552                  SourceNode,SinkNode,Edges):-
553        SourceNode is N+M+1,
554        SinkNode is N+M+2,
555        dim(Mapping,[SinkNode]),
556        arg(SourceNode,Mapping,source),
557        arg(SinkNode,Mapping,sink),
558        hash_create(Hash),
559        hash_add(Hash,node(SourceNode),source),
560        hash_add(Hash,node(SinkNode),sink),
561        (foreach(gcc(Low,High,Value),Bounds),
562         fromto([],E,[e(J,SinkNode,
563                        low_high{low:Low,high:High,mapping:0})|E],
564                EdgesSink),
565         count(J,N+1,_),
566         param(Mapping,Hash,SinkNode) do
567            (hash_find(Hash,Value,_) ->
568                writeln(no_value(Value)),
569                abort
570            ;
571                true
572            ),
573            hash_add(Hash,Value,J),
574            arg(J,Mapping,Value)
575        ),
576                 
577        (foreach(X,Variables),
578         count(J,1,_),
579         fromto([],E,E1,EdgesMain),
580         fromto([],F,[e(SourceNode,J,
581                        low_high{low:1,high:1,mapping:0})|F],EdgesSource),
582         param(Mapping,Hash,SourceNode) do
583            arg(J,Mapping,X),
584            get_full_domain_as_list(X,Domain),
585            (foreach(V,Domain),
586             fromto(E,Edges,Edges1,E1),
587             param(Hash,J,X) do
588                (hash_find(Hash,V,Node) ->
589                    Edges1 = [e(J,Node,
590                                low_high{low:0,high:1,mapping:0})|Edges]
591                ;
592                    excl(X,V),
593                    Edges1 = Edges
594                )
595            )
596        ),
597%        writeln(EdgesSource),
598%        writeln(EdgesMain),
599%        writeln(EdgesSink),
600        append(EdgesSource,EdgesMain,EE),
601        append(EdgesSink,EE,Edges).
602
603% remove values from variables corresponding to unmarked edges
604gcc_remove_unmarked_edges(NotMarked,Mapping):-
605        (foreach(e(VarN,ValueN,_),NotMarked),
606         param(Mapping) do
607            gcc_node_map(Mapping,VarN,Var),
608            gcc_node_map(Mapping,ValueN,Value),
609            ((atom(Var);atom(Value)) ->
610                true
611            ;
612                
613%            writeln(rem(VarN,Var,ValueN,Value)),
614                excl(Var,Value)
615            )
616        ).
617
618gcc_node_map(Mapping,Node,Value):-
619        arg(Node,Mapping,Value).
620
621
622
623/*********************************************************************
624
625SAME
626
627
628*********************************************************************/
629
630:- comment(same/2, [
631        amode: same(+,+),
632        args: ["Vars1":"A collection of N different variables or integers",
633               "Vars2":"A collection of N different variables or integers"
634              ],
635        summary: "Vars1 and Vars2 are constrained to be a permutation of each"
636                 " other in the values taken by the variables.",
637        kind:[constraint:[root:[ic,fd],extra:[gccat:same]]],
638        desc: html("\
639<P>
640    This constraint ensures that the values taken by the variables in Vars1
641    and Vars2 are permutations of each other. Vars1 and Vars must be the same
642    length.
643</P><P>
644    This is currently a prototype -- the constraint has not been tested
645    very extensively and little effort has been spent to optimise performance.
646    We welcome any feedback on using this constraint.
647</P><P>
648    This constraint is also known as same in the global constraint catalog. 
649    The implementation is the generalised arc-consistent version described
650    in the catalog. 
651")]).
652
653same(L1,L2):-
654        collection_to_list(L1,Vars1),
655        collection_to_list(L2,Vars2),
656        length(Vars1,N),
657        length(Vars2,N),
658        append(Vars1,Vars2,Variables),
659        call_priority((
660            check_same(Vars1,Vars2,N),
661            (ground(Variables) ->
662                true
663            ;
664                suspend(update_same(Variables,Vars1,Vars2,N,Susp),
665                        0,[Variables->any],Susp)
666            )
667        ), 2).
668
669:-demon(update_same/5).
670update_same(Variables,Vars1,Vars2,N,Susp):-
671        check_same(Vars1,Vars2,N),
672        (ground(Variables) ->
673            kill_suspension(Susp)
674        ;
675            unschedule_suspension(Susp)
676        ).
677
678check_same(Vars1,Vars2,N):-
679        same_create_edges(Vars1,Vars2,N,NrNodes,
680                          CenterEdges,SourceNode,SinkNode,Edges,
681                          Mapping,ValueHash),
682        make_graph(NrNodes,Edges,Graph),
683        feas_flow_with_lb(Graph, low of low_high,high of low_high, 
684                          SourceNode, SinkNode, 
685                          MaxFlowValue, MaxFlowEdges, _),
686        MaxFlowValue >= N,
687        same_residual_graph(MaxFlowEdges,Edges,ResidualEdges),
688        make_graph(NrNodes,ResidualEdges,ResidualGraph),
689        strong_components(ResidualGraph,StrongComponents),
690        mark_components(StrongComponents,NrNodes,CenterEdges,NotMarked),
691        remove_same_inconsistent(NotMarked,Mapping,ValueHash).
692
693same_create_edges(Vars1,Vars2,N,NrNodes,
694                  CenterEdges,
695                  SourceNode,SinkNode,Edges,Mapping,Hash):-
696        N2 is 2*N,
697        SourceNode is N2+1,
698        SinkNode is N2+2, 
699        FreeNode is SinkNode+1,
700        dim(Mapping,[N2]),
701        hash_create(Hash),
702        (foreach(X,Vars1),
703         count(J,1,_),
704         fromto([],A,[e(SourceNode,J,low_high{low:1,high:1,mapping:0})|A],LeftEdges),
705         fromto([],B,B1,MidLeftEdges),
706         fromto(FreeNode,C,C1,LastNode1),
707         param(SourceNode,Mapping,Hash) do
708            arg(J,Mapping,X),
709            get_full_domain_as_list(X,List),
710            (foreach(V,List),
711             fromto(B,BB,[e(J,Node,low_high{low:0,high:1,mapping:0})|BB],B1),
712             fromto(C,CC,CC1,C1),
713             param(J,Hash) do
714                new_node(Hash,V,Node,CC,CC1)
715            )
716        ),
717        (foreach(X,Vars2),
718         count(J,N+1,_),
719         fromto([],A,[e(J,SinkNode,low_high{low:1,high:1,mapping:0})|A],RightEdges),
720         fromto([],B,B1,MidRightEdges),
721         fromto(LastNode1,C,C1,NrNodes),
722         param(SinkNode,Mapping,Hash) do
723            arg(J,Mapping,X),
724            get_full_domain_as_list(X,List),
725            (foreach(V,List),
726             fromto(B,BB,[e(Node,J,low_high{low:0,high:1,mapping:0})|BB],B1),
727             fromto(C,CC,CC1,C1),
728             param(J,Hash) do
729                new_node(Hash,V,Node,CC,CC1)
730            )
731        ),
732        append(MidLeftEdges,MidRightEdges,CenterEdges),
733        append(LeftEdges,CenterEdges,Edges1),
734        append(RightEdges,Edges1,Edges).
735
736new_node(Hash,V,Node,CC,CC):-
737        hash_find(Hash,V,Node),
738        !.
739new_node(Hash,V,CC,CC,CC1):-
740        hash_add(Hash,V,CC),
741        hash_add(Hash,node(CC),V),
742        CC1 is CC+1.
743
744
745                
746same_residual_graph(MaxFlowEdges,Edges,ResidualEdges):-
747        (foreach(_-X,MaxFlowEdges),
748         fromto(Edges,A,A1,ResidualEdges) do
749            invert_edge(X,A,A1)
750        ).
751
752invert_edge(e(From,To,Cap),A,[e(To,From,Cap)|A]).
753
754
755remove_same_inconsistent(NotMarked,Mapping,Hash):-
756        dim(Mapping,[N2]),
757        (foreach(e(From,To,_),NotMarked),
758         param(Mapping,Hash,N2) do
759            ((From =< N2,hash_find(Hash,node(To),Value)) ->
760                arg(From,Mapping,X),
761                excl(X,Value)
762            ;(To =< N2,hash_find(Hash,node(From),Value)) ->
763                arg(To,Mapping,X),
764                excl(X,Value)
765            ;
766                writeln(inconsistent(From,To)),
767                abort
768            )
769        ).
770
771/************************************************************************
772
773inverse
774
775************************************************************************/
776
777:- comment(inverse/2, [
778        amode: inverse(+,+),
779        args: ["Succ":"A collection of N different variables or integers",
780               "Pred":"A collection  of N different variables or integers"
781              ],
782        summary: "Constrains elements of Succ to be the successors and"
783                 " Pred to be the predecessors of nodes in a digraph",
784        kind:[constraint:[root:[ic,fd]]],
785        desc: html("\
786<P>
787     Succ and Pred are list of N elements, representing a digraph of N nodes,
788     where the i'th element of Succ and Pred represents the successor and
789     predecessor of the node i respectively. The constraint enforces each
790     node in the digraph to have one successor and one predessor node, and
791     that if node y is the successor of node x, then node x is the
792     predecessor of node y.
793</P><P>
794    This is currently a prototype -- the constraint has not been tested
795    very extensively and little effort has been spent to optimise performance.
796    We welcome any feedback on using this constraint.
797</P><P>
798     This constraint is known as inverse in the global constraint catalog,
799     but with implicit node index based on the position in the list.  
800")]).
801
802inverse(XL,YL):-
803        collection_to_list(XL,Vars1),
804        collection_to_list(YL,Vars2),
805        length(Vars1,N),
806        length(Vars2,N),
807        append(Vars1,Vars2,Variables),
808        Variables :: 1..N,
809        call_priority((
810            check_inverse(Vars1,Vars2,N),
811            (ground(Variables) ->
812                true
813            ;
814                suspend(update_inverse(Variables,Vars1,Vars2,N,Susp),
815                        10,[Variables->any],Susp)
816            )
817        ), 2).
818
819
820:-demon(update_inverse/5).
821update_inverse(Variables,Vars1,Vars2,N,Susp):-
822        check_inverse(Vars1,Vars2,N),
823        (ground(Variables) ->
824            kill_suspension(Susp)
825        ;
826            unschedule_suspension(Susp)
827        ).
828
829check_inverse(Vars1,Vars2,N):-
830        dim(Matrix,[N,N]),
831        (for(I,1,N),
832         foreach(X,Vars1),
833         param(Matrix) do
834            get_full_domain_as_list(X,Dom),
835            (foreach(V,Dom),
836             param(I,Matrix) do
837                subscript(Matrix,[I,V],1-_)
838            )
839        ),
840        (for(J,1,N),
841         foreach(Y,Vars2),
842         param(Matrix) do
843            get_full_domain_as_list(Y,Dom),
844            (foreach(V,Dom),
845             param(Y,J,Matrix) do
846                subscript(Matrix,[V,J],Z-1),
847                (var(Z) ->
848                    excl(Y,V)
849                ;
850                    true
851                )
852            )
853        ),
854        (for(I,1,N),
855         foreach(X,Vars1),
856         param(Matrix) do
857            get_full_domain_as_list(X,Dom),
858            (foreach(V,Dom),
859             param(X,I,Matrix) do
860                subscript(Matrix,[I,V],1-Z),
861                (var(Z) ->
862                    excl(X,V)
863                ;
864                    true
865                )
866            )
867        ).
868
869/************************************************************************
870
871lex_le/2, lex_lt/2
872
873************************************************************************/
874
875:-local struct(store(alpha,beta,n)).
876
877:- comment(lex_le/2, [
878    summary:"List1 is lexicographically less or equal to List2",
879    amode:lex_le(+,+),
880    args:[
881	"List1":"List of integers or domain variables",
882	"List2":"List of integers or domain variables"
883    ],
884    kind:[constraint:[root:[ic,fd]]],
885    desc:html("\
886    	Imposes a lexicographic ordering between the two lists. 
887	I.e. either is the first element of List1 strictly smaller
888	than the first element of List2, or the first elements are
889	equal and the lexicographic order holds between the two list
890	tails. A non-existing element (i.e. when the end of list is 
891        reached)is strictly smaller than any existing element.
892</P><P>
893    This is currently a prototype -- the constraint has not been tested
894    very extensively and little effort has been spent to optimise performance.
895    We welcome any feedback on using this constraint.
896</P><P>
897        This constraint is known as lex_lesseq in the global constraint
898        catalog. The implementation here maintains generalised arc
899        consistency and uses the algorithm described in:
900        Z. Kiziltan, 'Symmetry Breaking Ordering Constraints, Ph.D thesis,
901        Uppsala University, 2004.
902")
903]).
904
905lex_le(XVector,YVector):-
906        collection_to_list(XVector,XList),
907        collection_to_list(YVector,YList),
908        setup_lex_gac(XList,YList,'=<').
909
910:- comment(lex_lt/2, [
911    summary:"List1 is lexicographically less than  List2",
912    amode:lex_lt(+,+),
913    args:[
914	"List1":"List of integers or domain variables",
915	"List2":"List of integers or domain variables"
916    ],
917    kind:[constraint:[root:[ic,fd]]],
918    desc:html("\
919    	Imposes a lexicographic ordering between the two lists. 
920	I.e.  either is the first element of List1 strictly smaller
921	than the first element of List2, or the first elements are
922	equal and the lexicographic order holds between the two list
923	tails. A non-existing element (i.e. when the end of list is 
924        reached)is strictly smaller than any existing element.
925</P><P>
926    This is currently a prototype -- the constraint has not been tested
927    very extensively and little effort has been spent to optimise performance.
928    We welcome any feedback on using this constraint.
929</P><P>
930        This constraint is known as lex_less in the global constraint
931        catalog. The implementation here maintains generalised arc
932        consistency and uses the algorithm described in:
933        Z. Kiziltan, 'Symmetry Breaking Ordering Constraints, Ph.D thesis,
934        Uppsala University, 2004.
935")
936]).
937
938lex_lt(XVector,YVector):-
939        collection_to_list(XVector,XList),
940        collection_to_list(YVector,YList),
941        call_priority(setup_lex_gac(XList,YList,'<'), 2).
942
943setup_lex_gac(XList,YList,Variant):-
944        length(XList,N),
945        length(YList,N),
946        find_alpha(XList,YList,1,A,XAlpha,YAlpha),
947        (A =< N ->
948            find_beta(XAlpha,YAlpha,A,I,-1,BB),
949            ((Variant = '=<',I =:= N+1) ->
950                B = 1.0Inf
951            ; BB = -1 ->
952                B = I
953            ;
954                B = BB
955            ),
956            A < B, % may fail
957            copy_array(XList,YList,N,XArray,YArray),
958            Store = store{alpha:A,beta:B,n:N},
959            lex_gac(XArray,YArray,Store,A,Variant),
960%            writeln(Store-XArray-YArray),
961            Upper is fix(min(B-1,N)), % convert to integer required
962%            writeln(Upper is min(B,N)),
963            (for(J,A,Upper),
964             param(XArray,YArray,Store,Variant) do
965                arg(J,XArray,Xj),
966                arg(J,YArray,Yj),
967%                writeln(j(J,Xj,Yj)),
968                (var(Xj) ->
969                    suspend(update_lex_gac(Xj,XArray,YArray,
970                                           Store,J,Variant,SuspX),4,
971                            [Xj->[inst,min]],
972                            SuspX)
973                ;
974                    true
975                ),
976                (var(Yj) ->
977                    suspend(update_lex_gac(Yj,XArray,YArray,
978                                           Store,J,Variant,SuspY),4,
979                            [Yj->[inst,max]],
980                            SuspY)
981                ;
982                    true
983                )
984            )
985        ;
986            % all elements equal
987            (Variant = '<' ->
988                fail
989            ;
990                true
991            )
992        ).
993
994:-demon update_lex_gac/7.
995update_lex_gac(Var,XArray,YArray,Store,J,Variant,Susp):-
996        lex_gac(XArray,YArray,Store,J,Variant),
997        (ground(Var) ->
998            kill_suspension(Susp)
999        ;
1000            unschedule_suspension(Susp)
1001        ).
1002
1003lex_gac(XArray,YArray,Store,I,Variant):-
1004        Store = store{alpha:A,beta:B},
1005        ((A =< I,I< B) ->
1006            arg(I,XArray,Xi),
1007            arg(I,YArray,Yi),
1008            ((I = A, I+1 =:= B) ->
1009                Xi #< Yi
1010            ;(I = A, I+1 < B) ->
1011                Xi #=< Yi,
1012                (Xi == Yi -> % ????
1013                    I1 is I+1,
1014                    update_alpha(XArray,YArray,I1,B,Store,Variant)
1015                ;
1016                    true
1017                )
1018            ; (A < I,I<B) ->
1019                get_lwb(Xi,Xmin),
1020                get_upb(Yi,Ymax),
1021                (((I =:= B-1,Xmin = Ymax);Xmin>Ymax) ->
1022                    II is I-1,
1023                    update_beta(XArray,YArray,A,II,Store,Variant)
1024                ;
1025                    true
1026                )
1027            ;
1028                writeln(not_reached(A,I,B)),
1029                abort
1030            )
1031        ;
1032            % the update index is outside the current [a,b) range, ignore
1033%            writeln(outside),
1034            true
1035        ).
1036
1037update_alpha(XArray,YArray,A,B,Store,Variant):-
1038        Store = store{n:N},
1039        A < B,
1040        (A =:= N+1 ->
1041            (Variant = '<' ->
1042                fail
1043            ;
1044                true
1045            )
1046        ;
1047            (A = B ->
1048                true
1049            ;
1050                arg(A,XArray,Xa),
1051                arg(A,YArray,Ya),
1052                (Xa == Ya ->
1053                    A1 is A+1,
1054                    update_alpha(XArray,YArray,A1,B,Store,Variant)
1055                ;
1056                    setarg(alpha of store,Store,A),
1057                    lex_gac(XArray,YArray,Store,A,Variant)
1058                )
1059            )
1060        ).
1061
1062update_beta(XArray,YArray,A,B,Store,Variant):-
1063        B+1 =\= A,
1064        arg(B,XArray,Xb),
1065        arg(B,YArray,Yb),
1066        get_lwb(Xb,Xmin),
1067        get_upb(Yb,Ymax),
1068        (Xmin < Ymax ->
1069            B1 is B+1,
1070            setarg(beta of store,Store,B1),
1071            lex_gac(XArray,YArray,Store,B,Variant)
1072        ;
1073            (Xmin = Ymax ->
1074                BB is B-1,
1075                update_beta(XArray,YArray,A,BB,Store,Variant)
1076            ;
1077                true
1078            )
1079        ).
1080
1081             
1082find_alpha([X|X1],[Y|Y1],I,A,Xr,Yr):-
1083        X == Y,
1084        !,
1085        I1 is I+1,
1086        find_alpha(X1,Y1,I1,A,Xr,Yr).
1087find_alpha(Xr,Yr,A,A,Xr,Yr).
1088
1089find_beta([X|X1],[Y|Y1],I,Iend,B,Bend):-
1090        get_lwb(X,Xmin),
1091        get_upb(Y,Ymax),
1092        Xmin =< Ymax,
1093        !,
1094        (Xmin = Ymax ->
1095            (B = -1 ->
1096                B1 = I
1097            ;
1098                B1 = B
1099            )
1100        ;
1101            B1 = -1
1102        ),
1103        I1 is I+1,
1104        find_beta(X1,Y1,I1,Iend,B1,Bend).
1105find_beta(_,_,I,I,B,B).
1106
1107              
1108copy_array(XList,YList,N,XArray,YArray):-
1109            dim(XArray,[N]),
1110            dim(YArray,[N]),
1111            (foreach(X,XList),
1112             foreach(Y,YList),
1113             foreacharg(X,XArray),
1114             foreacharg(Y,YArray) do
1115                true
1116            ).
1117
1118
1119