/*

  Some CP utilities in Picat.

  Note: Some of the predicates/global constraints are now implemented
  in Picat (and thus are commented out).

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

module cp_utils.

import cp.
% import util.
import random_v3.


% Different implementations of matrix_element(X,I,J,Val)
% to handle
%    Val = X[I,J]
% which is not available when I or J are CP-variables.
% 

matrix_element1(X, I, J, Val) =>
  element(I, X, Row),
  element(J, Row, Val).

matrix_element2(X, I, J, Val) =>
  nth(I, X, Row),
  element(J, Row, Val).

matrix_element3(X, I, J, Val) =>
     freeze(I, (nth(I, X, Row),freeze(J,nth(J,Row,Val)))).

matrix_element4(X, I, J, Val) =>
   freeze(I, (element(I, X, Row),freeze(J,element(J,Row,Val)))).

matrix_element5(X, I, J, Val) =>
  nth(I, X, Row),
  nth(J, Row, Val).

matrix_element6(X, I, J, Val) =>
   freeze(I, (nth(I, X, Row),freeze(J,element(J,Row,Val)))).

matrix_element7(X, I, J, Val) =>
 foreach (Row in 1..len(X), Col in 1..len(X[1]))
    (I #= Row #/\ J #= Col) #=> (Val #= X[Row,Col])
 end.


%
% converts a number Num to/from a list of integer List given a base Base
%
to_num(List, Base, Num) =>
        Len = length(List),
        Num #= sum([List[I]*Base**(Len-I) : I in 1..Len]).

to_num(List, Num) =>
       to_num(List, 10, Num).

%
% Ensure a Latin square, 
% i.e. all rows and all columns are different
%
latin_square(Board) =>
   foreach(Row in Board) all_different(Row) end,
   foreach(Column in transpose(Board)) all_different(Column) end.

%
% Ensure that List is increasing
%
% increasing(List) =>
%    foreach(I in 2..List.length) List[I-1] #=< List[I] end.

% increasing_strict(List) =>
%    foreach(I in 2..List.length) List[I-1] #< List[I] end.

% 
% Reified version of increasing
%
increasing_reif(L,B) =>
  (B #= 1) #<=> (sum([L[I-1] #> L[I] : I in 2..L.len]) #= 0).


%
% Ensure that values != 0 is increasing
%
increasing_except_0(List) =>
   Len = List.length,
   foreach(I in 1..Len, J in 1..Len, I < J) 
      (List[I] #!= 0 #/\ List[J] #!= 0) #=> List[I] #=< List[J] 
   end.


%
% Ensure that List is decreasing
%
% decreasing(List) =>
%    foreach(I in 2..List.length) List[I-1] #>= List[I] end.

% decreasing_strict(List) =>
%    foreach(I in 2..List.length) List[I-1] #> List[I] end.


%
% Ensure that values != 0 is decreasing
%
decreasing_except_0(List) =>
   Len = List.length,
   foreach(I in 1..Len, J in 1..Len, I < J) 
      (List[I] #!= 0 #/\ List[J] #!= 0) #=> List[I] #>= List[J] 
   end.

% 
% reified version of decreasing
%
decreasing_reif(L,B) =>
  (B #= 1) #<=> (sum([L[I-1] #< L[I] : I in 2..L.len]) #= 0).


%
% increasing_values(X, Values)
% Ensure that values that are in Values should be increasing in X.
%
% It's a generalization (and the oppsite of) the
% global constraint increasing_except_0(X).
%
% See increasing_values.pi for some examples
% 
increasing_values(X,Values) =>
  Len = X.length,
  foreach(I in 1..Len, J in 1..Len, I < J)
    ( sum( [X[I] #= V : V in Values] ) #> 0 #/\ sum([X[J] #= V : V in Values]) #> 0) #=> X[I] #=< X[J] 
  end.


% 
% All values in List are equal
%
all_equal(List) =>
  foreach(I in 2..List.len) List[I-1] #= List[I] end.


%
% reverse an array
%
reverse(X, Rev) =>
  Len = X.len,
  foreach(I in 1..Len)
    Rev[I] #= X[Len-I+1]
  end.

%
% Ensure that X is an square
%
% square(X) =>
%    Max = fd_max(X),
%    Y :: 0..ceiling(sqrt(Max)),
%    Y*Y #= X.

%
% Ensure that X is an square. Better Min limit.
%
square(X) =>
   fd_min_max(X,Min,Max),
   Y :: ceiling(sqrt(Min))..ceiling(sqrt(Max)),
   Y*Y #= X.


%
% Scalar product of the list A and X
%
% scalar_product(A, X, Product) => 
%    Product #= sum([A[I]*X[I] : I in 1..A.length]).

%
% scalar product of List A and X, with one of the relations: 
% 
%   #=, #<, #>, #<=, #>=, #!=
%
% scalar_product(A, X, Rel, Product) => 
%    scalar_product(A, X, P),
%    call(Rel,P,Product).


%%
%% Product = prod(List)
%%
%% returns the product of the numbers in List
%%
%% Note: This is now a built-in in Picat.
%
% prod(List, Product) =>
%    Product1 = 1, 
%    foreach(L in List)
%       % We must "taint" the variable to be a CP variable
%       Product1 := $Product1 * L
%    end,
%    Product #= Product1.


%
% Requires that all values in Xs != 0 must be distinct.
%
% Note: This is now a built-in predicate.
%
% alldifferent_except_0(Xs) =>
%   foreach(I in 1..Xs.length, J in 1..I-1)
%     (Xs[I] #!= 0 #/\ Xs[J] #!= 0) #=> (Xs[I] #!= Xs[J])
%    end.

% all_different_except_0(Xs) =>
%   foreach(I in 1..Xs.length, J in 1..I-1)
%      (Xs[I] #!= 0 #/\ Xs[J] #!= 0) #=> (Xs[I] #!= Xs[J])
%    end.

%
% Requires that all values in Xs != I must be distinct.
%
all_different_except_c(Xs,C) =>
  foreach(I in 1..Xs.length, J in 1..I-1)
     (Xs[I] #!= C #/\ Xs[J] #!= C) #=> (Xs[I] #!= Xs[J])
   end.


%
% reified version of all_different/1
%
all_different_reif(L,B) =>
  (B #= 1) #<=> (sum([L[I] #= L[J] : I in 2..L.len, J in 1..I-1]) #= 0).


%
% Ensure that the minumum value (> 0) is MinVal.
%
min_except_0(X,MinVal) =>
  Len = X.length,
  between(1,Len,I),
  MinVal #= X[I],
  foreach(J in 1..Len)
     MinVal #=< X[J] #\/ X[J] #= 0
  end,
  MinVal #> 0.


%
% nvalue(?N,?X)
%
% Requires that there are N distinct values in X.
%
% nvalue(N, X) =>
%   Len = length(X),
%   N #= sum([ (sum([ (X[J] #= I) : J in 1..Len]) #> 0) : I in 1..Len]).

% nvalue(N, X) =>
%    [Min, Max] = fd_min_max_array(X),
%    N #= sum([ (sum([ (X[J] #= I) : J in 1..X.length]) #> 0) : I in Min..Max]).


%
% Maximum/minimum domain value in array X
%
fd_max_array(X) = max([fd_max(V) : V in X]).
fd_min_array(X) = min([fd_min(V) : V in X]).

% combine them
fd_min_max_array(X) = [fd_min_array(X), fd_max_array(X)].

% Get the full domain of X
fd_min_max_array_dom(X) = fd_min_array(X)..fd_max_array(X).

%
% Get Min and Max for an array/list
%
/*
fd_min_max_array(X) = [Min,Max] =>
   Max = fd_max(X[1]),
   Min = fd_min(X[1]),
   foreach(Y in X) 
      if fd_min(Y) < Min then Min = fd_min(Y) end,
      if fd_max(Y) > Max then Max = fd_max(Y) end
   end.
*/


%
% nvalues(X,Op,N)
%
% Requires that the number of distinct values in the array X is 
%    Op N 
% where
% Op is either one of 
%   #=, #<m, #=<, #>=, #>
% (this is not checked though)    
%
nvalues(X, Op, N) =>
   nvalue(M,X),
   call(Op, M, N).

%
% global_cardinality(A, Gcc)
%
% This version is bidirectional but limited:
%
% Both A and Gcc are (plain) lists.
%  
% The list A can contain only values 1..Max (i.e. the length of Gcc).
% This means that the caller must know the max values of A.
% Or rather: if A contains another values they will not be counted.
% 
global_cardinality2(A, Gcc) =>
   Len = length(A),
   Max = length(Gcc),
   Gcc :: 0..Len,
   foreach(I in 1..Max) count(I,A,#=,Gcc[I]) end.

%
% global_cardinality_low_up(V,C,Low,Up)
%
% ensure that the occurrences of C[I] in the list V
% are between Low[I] and Up[I]
%
global_cardinality_low_up(V,C,Low,Up) =>
  T = new_list(C.len),
  foreach(I in 1..C.len)
    T[I] :: Low[I]..Up[I]
  end,
  Gcc = $[C[I]-T[I] : I in 1..C.len],
  global_cardinality(V,Gcc).


global_cardinality_table(Variables,Values) =>
  % sanity/symmetry
  increasing_strict([Values[J,1] :J in 1..Values.len]),
  
  foreach(I in 1..Values.len) 
    Values[I,2] #= sum([Variables[J] #= Values[I,1] : J in 1..Variables.len])
  end.



% From MiniZinc's lex2.mzn 
% """
%-----------------------------------------------------------------------------%
% Require adjacent rows and adjacent columns in the array 'x' to be
% lexicographically ordered.  Adjacent rows and adjacent columns may be equal.
%-----------------------------------------------------------------------------%
% """
% This use lex_le/1
lex2(X) =>
   Len1 = X.len,
   Len2 = X[1].length,
   foreach(I in 2..Len1) 
      lex_le([X[I-1, J] : J in 1..Len2], [X[I, J] : J in 1..Len2])
   end,
   foreach(J in 2..Len2)
      lex_le([X[I ,J-1] : I in 1..Len1], [X[I, J] : I in 1..Len1])      
   end.

% Note: This use lex_lt/1.
lex2lt(X) =>
   Len1 = X.len,
   Len2 = X[1].length,
   foreach(I in 2..Len1) 
      lex_lt([X[I-1, J] : J in 1..Len2], [X[I, J] : J in 1..Len2])
   end,
   foreach(J in 2..Len2)
      lex_lt([X[I ,J-1] : I in 1..Len1], [X[I, J] : I in 1..Len1])      
   end.





%% lex_lt/2 and lex_le/2 are now built-in constraints.
/*

% Port of MiniZinc's lex_less_int.mzn
% """
%-----------------------------------------------------------------------------%
% Requires that the array 'x' is strictly lexicographically less than array 'y'.
% Compares them from first to last element, regardless of indices
%-----------------------------------------------------------------------------%
% """
lex_less(X,Y) =>
   LX = 1,
   UX = X.length,
   LY = 1,
   UY = Y.length,
   Size = max(UX-LX,UY-LY),
   B = new_list(Size+2), % (Note: The MiniZinc version is 0-based.)
   B :: 0..1,
   B[1] #= 1,
   foreach(I in 1..Size+1) 
      B[I] #= ( X[I] #<= Y[I] #/\ (X[I] #< Y[I] #\/ B[I+1] #= 1) )
   end,
   B[Size + 2] #= (Size - 1 #< Size - 1).


% Port of MiniZinc's lex_lesseq_in.mzn
% """
%-----------------------------------------------------------------------------%
% Requires that the array 'x' is lexicographically less than or equal to
% array 'y'.  Compares them from first to last element, regardless of indices
%-----------------------------------------------------------------------------%
% """
lex_lesseq(X,Y) =>
   LX = 1,
   UX = X.length,
   LY = 1,
   UY = Y.length,
   Size = max(UX-LX,UY-LY),
   B = new_list(Size+2), % (Note: The MiniZinc version is 0-based.)
   B :: 0..1,
   B[1] #= 1,
   foreach(I in 1..Size+1) 
      B[I] #= ( X[I] #<= Y[I] #/\ (
                ((I #= Size) #=> 1)
                #/\ 
                ((I #< Size) #=> X[I] #< Y[I] #\/ B[I+1] #= 1) 
                )
              )
   end.
*/

lex_greater(X,Y) => lex_less(Y,X).
lex_greatereq(X,Y) => lex_lesseq(Y,X).



%
% Alternative approach where we convert to two
% numbers and ensure that the first number is < second number.
%
lex_less2(X,Y,Base) =>
  to_num(X,Base,NumX),
  to_num(Y,Base,NumY),
  NumX #<= NumY.


%
% exactly(?N,?X,?V)
%
% Requires that exactly N variables in X takes the value V.
%
% exactly(N, X, V) =>
%   count(V,X,#=,N).

%
% atmost(?N,?X,?V)
%
% Requires that atmost N variables in X takes the value V.
%
% atmost(N,X,V) => 
%   count(V,X,#=<,N).

%
% atleast(?N,?X,?V)
%
% Requires that atleast N variables in X takes the value V.
%
% atleast(N,X,V) => 
%   count(V,X,#>=,N).


% The permutation from A <-> B using the permutation P
permutation3(A,P,B) =>
   foreach(I in 1..A.length)
       %  B[I] #= A[P[I]]
       PI #= P[I],
       BI #= B[I],
       element(PI, A, BI)
   end.


%
% knapsack: 
%  given Weights, Values and decision variable list Take
%  - ensure that the total weights <= WeightMax
%  - calculate the profit (which will be maximized)
%
knapsack(Weights, Values,Take, WeightMax,Profit) =>
   sum([W*T : {W,T} in zip(Weights, Take)]) #=< WeightMax,
   Profit #= sum([V*T : {V,T} in zip(Values,Take)]).


value_precede_chain(C, X) =>
  foreach(I in 2..C.length)
    value_precede(C[I-1], C[I], X)
  end.

% This definition is inspired by 
% MiniZinc definition value_precede_int.mzn
value_precede(S,T,X) =>
   XLen = X.length,
   B = new_list(XLen+1),
   B :: 0..1,
   foreach(I in 1..XLen)
     % Xis :: 0..1,
     Xis #= (X[I] #= S),
     (Xis #=> (B[I+1] #= 1))
     #/\ ((#~ Xis #= 1) #=> (B[I] #= B[I+1]))
     #/\ ((#~ B[I] #= 1) #=> (X[I] #!= T))
   end,
   B[1] #= 0.



%
% M1 x M2 = M3
%
% assume square matrices
%
% (not very efficient...)
%
matrix_multi_cp(M1,M2,M3) => 
   N = M1.length,
   Dim = 1..N,
   foreach(I in Dim, J in Dim)
      Zip = [T : {P,Q} in zip([M1[I,K] : K in Dim],[M2[L,J] : L in Dim]), T #= P*Q],
      M3[I,J] #= sum(Zip)
   end.


% assume square matrix
matrix_power_cp(M,M2) =>
   matrix_multi_cp(M,M,M2).


% chain(+List, +Relation)
% 
% Chain establishes the Relation constraint between each pair of adjacent 
% elements in List.
% 
% Example:
%   
%   L = new_list(3),
%   L :: 1..3,
%   chain(L,#>=),
%   All = solve_all(L).
% 
% ->
%  All = [[1,1,1],[2,1,1],[2,2,1],[2,2,2],[3,1,1],[3,2,1],[3,2,2],[3,3,1],[3,3,2],[3,3,3]]
%%
chain(L,Rel) => 
  foreach(I in 2..L.len)
    call(Rel,L[I-1],L[I])
  end.


%
% The sum of Seq numbers must be between Low and Up.
%
% Note: Seq must be instantiated, but neither Low or Up has
% to be (the result may be weird unless they are, though).
%
sliding_sum(Low, Up, Seq, Variables) =>
   foreach(I in 1..Variables.length-Seq+1)
      Sum #= sum([Variables[J] : J in I..I+Seq-1]),
      Sum #>= Low,
      Sum #=< Up
   end.

%
% The number of X's in a sequence of Seq numbers must be between Low and Up.
% Note: Seq must be instantiated, but neither Low or Up has
% to be (the result may be weird unless they are, though).
%
sliding_count(X,Low, Up, Seq, Variables) =>
   foreach(I in 1..Variables.length-Seq+1)
      Sum #= sum([Variables[J] #= X : J in I..I+Seq-1]),
      Sum #>= Low,
      Sum #=< Up
   end.


%
% A is the number of fix_points in the list X
%
fix_points(A, X) =>
    A #= sum([X[I] #= I : I in 1..X.len]).



%
% gcd(X,Y,G)
% G is the gcd of X and Y.
%
% Not very efficient.
%
gcd(X, Y, GCD) =>
  fd_min_max(X,_MinX,MaxX),
  fd_min_max(Y,_MinY,MaxY),
  PMax = min(MaxX,MaxY),
  % PMin = max(MinX,MinY),
  P :: 1..PMax,

  foreach (I in 1..PMax)
    I #= P #=> (X mod I #= 0 #/\ Y mod I #= 0 #/\ GCD #= I),
    foreach(J in I+1..PMax) 
      I#>= P #=> (X mod J #!= 0 #\/ Y mod J #!= 0)
    end
  end.


%
% lcm(X,Y,LCM)
% cm is the lcm of x and y
%
lcm(X, Y, LCM) =>
  fd_min_max(X,_MinX,MaxX),
  fd_min_max(Y,_MinY,MaxY),
  PMax = min(MaxX,MaxY),
  GCD :: 1..PMax,
  gcd(X,Y,GCD),
  LCM #= X*Y div GCD.


%
% All 0s must be placed last
%
collect_0s_last(X) =>
  N = X.len,
  foreach(I in 1..N)
    X[I] #= 0 #=> sum([X[J] #> 0 : J in I+1..N]) #= 0
  end.



%
% Find the first position of value Val in the list Y.
% Note that we assume that Val is in the list.
%
first_pos(Val, Y, Ix) =>
  N = Y.len,
  Ix :: 1..N,
  sum([
       Y[I] #= Val #/\ Ix #= I #/\
       sum([Y[J] #= Val : J in 1..I-1]) #= 0
      : I in 1..N]) #= 1.

%
% Find the first position (FirstPos[I]) of value Vals[I] in the list Y.
% Note that we assume that Val is in the list.
%
first_pos_all(Vals,Y,FirstPos) =>
  foreach({I,V} in zip(1..Vals.len,Vals))
    first_pos(V, Y, FirstPos[I])
  end.



% From what_are_the_five_numbers.pi
% Assume integer values
mean(X,Mean) =>
  Mean*X.len #= sum(X).

median(X,Median) =>
  Median #= X[ceiling(X.len / 2)].

mode(X,Mode) =>
  element(_,X,Mode), % Mode must be in the list
  % Count the number of occurrences of each number
  MaxVal = max([fd_max(X[I]) : I in 1..X.len]),
  GCC = new_list(MaxVal+1), GCC :: 0..X.len,
  global_cardinality(X,$[I-GCC[I] : I in 1..MaxVal]),
  % the max value is unique
  element(Mode,GCC,ModeNum), 
  ModeNum #> 1,  
  sum([V #= ModeNum : V in GCC ]) #= 1.


%
% product(L,X,Product)
% better handling than prod/1 for 0 in X and especially L
% See the_weekly_challenge_218.pi for an example
%
product(L,X,Product) =>
  product(L,X,1,Product).
product([],[],Product,Product).
product([H|T],[X|Xs],Product0,Product) :-
  (X #= 1 #/\ H #!= 0) #=> Product1 #= (X*H)*Product0,
  (X #= 0 #\/ H #= 0) #=> Product1 #= Product0,
  product(T,Xs,Product1,Product).



%
% We assume integer values.
%
mean(X,Mean) =>
  Mean*X.len #= sum(X).

median(X,Median) =>
  Median #= X[ceiling(X.len / 2)].

mode(X,Mode) =>
  element(_,X,Mode), % Mode must be in the list
  % Count the number of occurrences of each number
  MaxVal = max([fd_max(X[I]) : I in 1..X.len]),
  GCC = new_list(MaxVal+1), GCC :: 0..X.len,
  global_cardinality(X,$[I-GCC[I] : I in 1..MaxVal]),
  % the max value is unique
  element(Mode,GCC,ModeNum), 
  ModeNum #> 1,  
  sum([V #= ModeNum : V in GCC ]) #= 1.