/*

  Rectangle placement with rotation in Picat.

  This is model that place rectangular shapes in a grid of certain
  width and height, including rotation of the rectangle.

  Compare with the CSPLib #9 for perfect square placement problem.
  https://www.csplib.org/Problems/prob009/
  This problem is implemented in 
  http://hakank.org/picat/perfect_square_placement.pi
  and is the basis of this model.

  The biggest difference between these two problems is that
  the perfect square problem ensures that there are no unoccupied
  places (holes) while the rectangular problem do not has this
  guarantee. Another difference is that widhth and height might
  differ.



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

*/

import sat.
% import cp.

main => go.

%
% Show first solution.
%
go ?=>

   Map = get_global_map(),
   Map.put(count,0),

   % problem(0,X,Y,Width,Height),
   problem(3,X,Y,Width,Height),
   % problem(5,X,Y,Width,Height),

   time2(solve_and_print(X,Y,Width,Height)),
   % rectangle_placement(X,Y,Width,Height, X2,Y2,StartX,StartY,Rotation),
   % print_solution(X2,Y2,Width,Height, StartX,StartY,Rotation),

   flush(stdout),
   Map.put(count,Map.get(count)+1),
   println(count=Map.get(count)),

   fail,
   
   nl.

go => println(count=get_global_map().get(count)).


%
% Solve all instances (first solutions) for the problems defined below:
%   problem/5.
%
% Standard version                       Shapes
% Problem  Time to
%          first solution  
%          (SAT solver)    
% ----------------------------------------
%    1     0.007s          X = [2, 2, 2, 2]             Y = [2, 2, 2, 2],
%    2     0.032s          X = [2, 4, 2, 3, 3, 2, 2],   Y = [4, 2, 4, 6, 3, 2, 2],
%    3     0.952s          X = [1,2,3,4,5,6,7,8,9,10],  Y = [5,3,3,2,2,2,2,2,1,1]
%    4     0.094s          X = [1,2,3,4,5,6,7,8,9,10],  Y = [2,2,2,2,2,2,2,2,2,2], 
%    5     5.809s          Problem #18 from CSPLib problem #9 
%    6     0.107s          X = 1..10                    Y = 1..10
%    7     0.166s          X = 1..11                    Y = 1..11
%    8     0.342s          X = 1..12                    Y = 1..12
%    9     1.381s          X = 1..15                    Y = 1..15
%   10     6.083s          X = 1..20                    Y = 1..20         
%   11     4.944s          X = 1..21                    Y = 1..21
%   12     4.921s          X = 1..22                    Y = 1..22
%   13     10.759s         X = 1..23                    Y = 1..23
%   14     2578.15         X = 1..24                    Y = 1..24
% 
%
go2 ?=>

   % Problems = find_all(P,problem(P,_X,_Y,_Height,_Width)),
   % foreach(Problem in Problems, Problem != 5)
   %    println(problem=Problem),
   %    problem(Problem,X,Y,Width,Height),
   %    time2(solve_and_print(X,Y,Width,Height,MustTouch)),
   %    nl
   % end,
   member(Problem,1..14),
   println(problem=Problem),
   problem(Problem,X,Y,Width,Height),
   
   time2(once(solve_and_print(X,Y,Width,Height))),

   flush(stdout),   
   fail,
   
   nl.

go2 => true.


%
% Solve all Korf instances (first instance)
%
% Perhaps it's silly to test a rectangle placement model
% with rotation on these square packing problems.
% But they have also the right to be tested. :-)

% Problem  Time to        
%          first solution  
%          (SAT solver)    
% ----------------------------------------
%    1     0.002s
%    2     0.002s
%    3     0.005s
%    4     0.011s
%    5     0.018s
%    6     0.034s
%    7a    0.043s
%    7b    0.046s
%    8     0.077s
%    9     0.085s
%   10     0.111s   
%   11     0.332s
%   12     0.872s
%   13     0.536s
%   14     0.577s
%   15     1.1s
%   16a    2.134s
%   16b    4.999s
%   17     3.098s
%   18     2.718s
%   19     11.047s
%   20     15.985s
%   21     "Impossible" CHECK THIS
%   22     57.446s
%   23     1704.12s
%   24     2640.5s
%   25     -
% 
%
go3 ?=>

   korf(N,Width,Height),
   println([n=N,width=Width,height=Height]),
   X = 1..N,
   Y = 1..N,   
   time2(once(solve_and_print(X,Y,Width,Height))),

   flush(stdout),   
   fail,
   
   nl.

go3 => true.


%
% Testing
%
go4 ?=>

   Map = get_global_map(),
   Map.put(count,0),

   % _ = random2(),

   N = 6, % number of shapes
   MX = 10,
   MY = 10, % max length
   X = [random(1,MX) : _ in 1..N],
   Y = [random(1,MY) : _ in 1..N],
   println(x=X),
   println(y=Y),
   RectangleSum = sum([X[I]*Y[I] : I in 1..N]),
   println(rectangleSum=RectangleSum),
   All = [[W*H,W,H] : W in 1..MX*MY,H in W..MX*MY, W*H >= RectangleSum].sort,
   % println(all=All),
   member([_Size,Width,Height],All),
   % println([width=Width,height=Height]),
   time2(once(solve_and_print(X,Y,Width,Height))),
   
   flush(stdout),
   Map.put(count,Map.get(count)+1),
   println(count=Map.get(count)),

   fail,
   
   nl.

go4 => println(count=get_global_map().get(count)).


 
%
% Solve the problem:
% Input:
%  - X,Y: the rectangle size in x and y positions
%  - Width, Height: the size of the main grid
%    must touch (left, right, up, down, or by diagonals)
% Output:
% - StartX, StartY: the coordinates of the start positions
%
rectangle_placement(X,Y,Width,Height, X2,Y2,StartX2,StartY2,Rotation) =>

   N = length(X),
   println([n=N,width=Width,height=Height]),

   SumXY = sum([X[I]*Y[I] : I in 1..N]) ,
   println([sumXY=SumXY,'Width*Height'=(Width*Height)]),
   if SumXY > Width*Height then
     printf("Impossible problem: sum X and Y (%d) >  width*height (%d)\n",SumXY,Width*Height),
     fail
   end,

   MaxWidthHeight = max([Width,Height]),

   StartX = new_list(N),
   StartX :: 1..Width,

   StartY = new_list(N),
   StartY :: 1..Height,

   Rotation = new_list(N),
   Rotation :: 0..1, % Is the I'th shape rotated?

   StartX2 = new_list(N),
   StartX2 :: 1..MaxWidthHeight,

   StartY2 = new_list(N),
   StartX2 :: 1..MaxWidthHeight,

   XY = (X++Y).sort_remove_dups,
   X2 = new_list(N),
   X2 :: XY,
   Y2 = new_list(N),   
   Y2 :: XY,
   
   foreach(I in 1..N)
     % If it's a square we don't have to rotate.
     if X[I] == Y[I] then
       Rotation[I] #= 0
     end,
     (Rotation[I] #= 0) #=> (StartX2[I] #= StartX[I] #/\ StartY2[I] #= StartY[I] #/\ X2[I] #= X[I] #/\ Y2[I] #= Y[I]),
     (Rotation[I] #= 1) #=> (StartX2[I] #= StartY[I] #/\ StartY2[I] #= StartX[I] #/\ X2[I] #= Y[I] #/\ Y2[I] #= X[I])
   end,


   % CP seems to be faster with this, but SAT is slower
   % This don't work 
   % if membchk(cp, loaded_modules()) then
   %   cumulative(StartX2,X2,Y2,Height),
   %   cumulative(StartY2,Y2,X2,Width)
   % end,

   diffn4(StartX2,StartY2,X2,Y2),

   %% Tests (see below for implementations)
   % diffn_me(StartX,StartY,X,Y),
   % diffn_nonstrict(StartX,StartY,X,Y),     

   foreach(I in 1..N)
     StartX2[I] + X2[I] #=< Width+1
   end,
   foreach(J in 1..N)
     StartY2[J] + Y2[J] #=< Height+1
   end,

   Vars = StartX ++ StartY ++ Rotation ++ StartX2 ++ StartY2 ++ X2 ++ Y2,

   println(solve),
   % println(vars=Vars),
   solve($[degree,split], Vars),
   println(rotation=Rotation).
   % solve($[ff,split], Vars).   
   
%
% Print solution
%
print_solution(X,Y,Width,Height,StartX,StartY,Rotation) =>

   N = length(X), % number of rectangles
   println([n=N,width=Width,height=Height]),
   println(x=X),
   println(y=Y),
   println(startX=StartX),
   println(startY=StartY),
   println(rotationFromOrig=Rotation),
   
   % Print the grid
   M = new_array(Width,Height), % bind_vars(M,0),
   foreach(S in 1..N)
     % println([s=S,x=X[S]..X[S]+A[S]-1,y=Y[S]..Y[S]+A[S]-1]),
     foreach(I in StartX[S]..min(StartX[S]+X[S]-1,Width),J in StartY[S]..min(StartY[S]+Y[S]-1,Height))
       % Debugging: Check that we got a proper solution, i.e. that there is no
       % overlapping of values. 
       if nonvar(M[I,J]) then
          println([s=S,[I,J],already_defined,M[I,J]])
          % , halt
       end,
       M[I,J] = S
     end
   end,

   NumEmpty = 0,
   % Numeric solution
   foreach(I in 1..Width)
     foreach(J in 1..Height)
       V = M[I,J],
       if var(V) then
         print(" _ "),
         NumEmpty := NumEmpty + 1
       else
         printf("%2d ",M[I,J])
       end
     end,
     nl
   end,

   if N <= 26 then
     % _ is the unoccupied slots
     Alpha = "ABCDEFGHIJKLMNOPQRSTUVXYZ",
     foreach(I in 1..Width)
       foreach(J in 1..Height)
         V = M[I,J],
         if var(V) then
           print("_")
         else
           printf("%w",Alpha[M[I,J]])
         end
       end,
       nl
     end
   end,
   nl,
   println(num_empty=NumEmpty),
   nl.


solve_and_print(X,Y,Width,Height) =>
   rectangle_placement(X,Y,Width,Height, X2,Y2,StartX,StartY,Rotation),
   print_solution(X2,Y2,Width,Height, StartX,StartY,Rotation).


% This is from fzn_picat_cp.pi
diffn4(VecX,VecY,VecDX,VecDY) =>
    Rects = [[VecX[I],VecY[I],VecDX[I],VecDY[I]] : I in 1 .. length(VecX)],
    diffn(Rects).


%
% From MiniZinc's diffn_nonstrict/4
%
diffn_nonstrict(X,Y,DX,DY) =>
  N = X.len,
  foreach(I in 1..N, J in I+1..N)
    DX[I] #= 0 #\/
    DX[J] #= 0 #\/
    DY[I] #= 0 #\/
    DY[J] #= 0 #\/
    X[I] + DX[I] #=< X[J] #\/
    Y[I] + DY[I] #=< Y[J] #\/
    X[J] + DX[J] #=< X[I] #\/
    Y[J] + DY[J] #=< Y[I]
  end.

/*
From MiniZinc's diffn/4
*/
diffn_me(X,Y,DX,DY) =>
  N = X.len,
  foreach(I in 1..N, J in I+1..N)
    X[I] + DX[I] #<= X[J] #\/
    Y[I] + DY[I] #<= Y[J] #\/
    X[J] + DX[J] #<= X[I] #\/
    Y[J] + DY[J] #<= Y[I]
  end.



%
% time2 + time_out as a function.
%
time2f(Goal,Timeout) = [End,Backtracks,Status] =>
    statistics(runtime,_),
    statistics(backtracks, Backtracks1),
    time_out(Goal,Timeout,Status),
    statistics(backtracks, Backtracks2),
    statistics(runtime, [_,End]),
    Backtracks = Backtracks2 - Backtracks1.


%
% Problem instances
%

% Test problems
problem(0,X,Y,Width,Height) :- 
   X = [1,3,2,2,1],
   Y = [4,2,1,2,5],
   Width = 5,   
   Height = 5.


% This is a variant of the problem from
% https://stackoverflow.com/questions/70483069/n-rectangles-inside-the-square-minizinc
problem(1,X,Y,Width,Height) :- 
   X = [2, 2, 2, 2],
   Y = [2, 2, 2, 2],
   % Width = 8,   
   % Height = 8.
   Width = 4,   
   Height = 5.


% A slightly larger problem: 
problem(2,X,Y,Width,Height) :- 
   X = [2, 4, 2, 3, 3, 2, 2],
   Y = [4, 2, 4, 6, 3, 2, 2],
   Width = 8,
   Height = 8.


% 
problem(3,X,Y,Width,Height) :- 
  X = [1,2,3,4,5,6,7,8,9,10],
  Y = [5,3,3,2,2,2,2,2,1,1], 
  Width = 10,
  Height = 10.


% 
problem(4,X,Y,Width,Height) :- 
   X = [1,2,3,4,5,6,7,8,9,10],
   Y = [2,2,2,2,2,2,2,2,2,2], 
   Width = 11,
   Height = 11.

% Problem #18 from CSPLib problem #9 (perfect square placement)
% 8 solutions.
% Use the SAT solver for this.
problem(5,X,Y,Width,Height) :-
  X = [1,3,4,9,10,11,12,16,17,18,22,23,24,40,41,60,66,68,70,71,74,76,79],
  Y = [1,3,4,9,10,11,12,16,17,18,22,23,24,40,41,60,66,68,70,71,74,76,79],   
  Width = 215,
  Height = 215.


%
% Here are some instances with has squares
% but with (random) holes in the grid.
%

% Time to first solution: 0.289s
problem(6,X,Y,Width,Height) :- 
   X = 1..10,
   Y = 1..10,
   Width = 21,
   Height = 21.

problem(7,X,Y,Width,Height) :- 
   X = 1..11,
   Y = 1..11,
   Width = 21,
   Height = 26.

problem(8,X,Y,Width,Height) :- 
   X = 1..12,
   Y = 1..12,
   Width = 24,
   Height = 29.


problem(9,X,Y,Width,Height) :- 
   X = 1..15,
   Y = 1..15,
   Width = 33,
   Height = 39.


problem(10,X,Y,Width,Height) :-
   X = 1..20,
   Y = 1..20,
   Width = 42,
   Height = 70.

% 2.704s
problem(11,X,Y,Width,Height) :-
   X = 1..21,
   Y = 1..21,
   Width = 44,
   Height = 77.

% 9.207s
problem(12,X,Y,Width,Height) :-
   X = 1..22,
   Y = 1..22,
   Width = 56,
   Height = 69.


% 
problem(13,X,Y,Width,Height) :-
   X = 1..23,
   Y = 1..23,
   Width = 48,
   Height = 92.

% 
problem(14,X,Y,Width,Height) :-
   X = 1..24,
   Y = 1..24,
   Width = 56,
   Height = 88.


% See https://sofdem.github.io/gccat/gccat/Ksmallest_rectangle_area.html
% 225.123s
problem(20,X,Y,Width,Height) :-
   X = 1..35,
   Y = 1..35,
   Width = 123,
   Height = 123.



% These (optimal) widths and heights for 1..N packing  are from
% Richard E. Korf "Optimal Rectangle Packing: New Results"
% https://www.aaai.org/Papers/ICAPS/2004/ICAPS04-019.pdf

% korf(N,Width,Height).
korf(1,1,1).
korf(2,2,3).
korf(3,3,5).
korf(4,5,7).
korf(5,5,12).
korf(6,9,11).
korf(7,7,22). % yes, there are two solutions for 7x7!
korf(7,11,14). 
korf(8,14,15).
korf(9,15,20).
korf(10,15,27).
korf(11,19,27).
korf(12,23,29).
korf(13,22,38).
korf(14,23,45).
korf(15,23,55).
korf(16,27,56). % Two of 16x16
korf(16,28,54).
korf(17,39,46).
korf(18,31,69).
korf(19,47,53).
korf(20,34,85).
korf(21,38,85). % This is impossible! sum xs ys=3311 > width*height=3230
korf(22,39,98).
korf(23,64,68).
korf(24,56,88).
korf(25,43,129).