/* 

  Magic cube in Picat.

  Here are two models:
  - From the F1 model
    http://www.f1compiler.com/samples/Magic%20Cube%203x3x3.f1.html

  - From Pitrat's MALICE system
    and Pitrat's blog post
    http://bootstrappingartificialintelligence.fr/WordPress3/2015/01/caia-surprised-me-part-i/

  The first one is very fast.

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

*/


% import util.
import cp.
% import sat.



main => go.



%
% From the F1 model
% http://www.f1compiler.com/samples/Magic%20Cube%203x3x3.f1.html
%
% Representation:
% """
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%
% 1-----10-----19
% | 2     11    |20
% |   3-----12-----21
% 4   | 13     22   |
% | 5 |   14    |23 |
% |   6     15  |  24
% 7---|-16-----25   |
%   8 |   17     26 |
%     9-----18-----27
%
% 
% """
%
% All solutions have C = 42, hence we can fix it 
% (though it might be considered cheating).
%
% CP : fixed C: 0.02s  (unfixed C: 0.4s)
% SAT: fixed C: 3.7s   (unfixed C: 14.7s)
%
% There are 4 solutions with symmetry breaking:
%  [1,23,18,17,3,22,24,16,2,15,7,20,19,14,9,8,21,13,26,12,4,6,25,11,10,5,27]
%  [2,24,16,18,1,23,22,17,3,15,7,20,19,14,9,8,21,13,25,11,6,5,27,10,12,4,26]
%  [4,26,12,18,1,23,20,15,7,17,3,22,19,14,9,6,25,11,21,13,8,5,27,10,16,2,24]
%  [6,26,10,17,1,24,19,15,8,16,3,23,21,14,7,5,25,12,20,13,9,4,27,11,18,2,22]
% 
% Without symmetry breaking there are 192 solutions (found in fixed: 0.8s / unfixed: 10.7s)
% 
go ?=>
  nolog,
  F = new_list(27),
  F:: 1..27,

  F = [F01,F02,F03,F04,F05,F06,F07,F08,F09,F10,
       F11,F12,F13,F14,F15,F16,F17,F18,F19,F20,
       F21,F22,F23,F24,F25,F26,F27],

  C #= 42,
  % C :: 6..78,
  % C #!= 42, % "prove" that C must be 42, i.e.this give no solution

  all_different(F),
  % all_distinct(F),

  % Remove symmetries...
  F01 #< F03, F01 #< F07, F01 #< F09, F01 #< F19,  F01 #< F21, F01 #< F25, F01 #< F27,
  F02 #> F04, F04 #> F10,


  % Setup constraints between elements
  F01 + F02 + F03 #= C, F01 + F10 + F19 #= C, F01 + F04 + F07 #= C,
  F04 + F05 + F06 #= C, F02 + F11 + F20 #= C, F02 + F05 + F08 #= C,
  F07 + F08 + F09 #= C, F03 + F12 + F21 #= C, F03 + F06 + F09 #= C,

  F10 + F11 + F12 #= C, F04 + F13 + F22 #= C, F10 + F13 + F16 #= C,
  F13 + F14 + F15 #= C, F05 + F14 + F23 #= C, F11 + F14 + F17 #= C,
  F16 + F17 + F18 #= C, F06 + F15 + F24 #= C, F12 + F15 + F18 #= C,
  F19 + F20 + F21 #= C, F07 + F16 + F25 #= C, F19 + F22 + F25 #= C,

  F22 + F23 + F24 #= C, F08 + F17 + F26 #= C, F20 + F23 + F26 #= C,
  F25 + F26 + F27 #= C, F09 + F18 + F27 #= C, F21 + F24 + F27 #= C,
  F01 + F14 + F27 #= C, F03 + F14 + F25 #= C,
  F07 + F14 + F21 #= C, F09 + F14 + F19 #= C,


  % stats(F,Size,Degree,Min,Max),
  % println(F),

  Vars = F,

  solve($[split], Vars),

  println(c=C),
  % println(f=F),

  printMC3(F,0),
  nl,

  fail,

  nl.

go => true.

printMC3(MC,R) =>
   if R < 3 then 
        print("\n"),
        printDigit2(MC[3*R+1])  , printDigit2(MC[3*R+2])  , printDigit2(MC[3*R+3]) ,
        print("      ") ,
        printDigit2(MC[3*R+10])  , printDigit2(MC[3*R+11]) , printDigit2(MC[3*R+12]) ,
        print("      ") ,
        printDigit2(MC[3*R+19]) , printDigit2(MC[3*R+20]) , printDigit2(MC[3*R+21]) ,
        printMC3(MC,R+1)
    else
        print("\n") 
    end.

printDigit2(D) =>
  printf("%2d ", D).


%
% Problem formulation from Pitrat's MALICE.
%
% From Pitrat's blog post
% http://bootstrappingartificialintelligence.fr/WordPress3/2015/01/caia-surprised-me-part-i/
%
%
% SAT finds these solutions after about 3 minute (and then interrupted) with
% TheSum #= 126
%
%   [8,9,20,7,26,18,12,16,10,14,23,6,25,15,2,24,4,13,22,5,19,11,1,21,3,27,17]
%   [16,3,24,12,23,1,22,20,5,11,17,18,7,10,25,4,15,19,27,8,2,13,9,26,14,21,6]
%   [14,17,21,10,4,15,12,25,8,24,5,6,23,18,1,3,19,27,2,11,26,22,20,13,16,7,9]
%   [22,18,23,21,1,4,14,8,15,3,2,11,7,16,19,17,24,27,12,26,9,20,25,13,10,6,5]
%   [22,8,25,5,2,26,19,18,1,16,6,3,20,13,9,12,23,24,11,14,21,17,27,7,4,15,10]
%   [16,22,8,12,5,18,7,13,25,26,9,6,11,14,17,21,19,3,4,20,15,2,23,24,27,1,10]
%  
%
go2 =>
  nolog,
  F = new_list(27),
  F:: 1..27,

  % TheSum :: 45..207,
  TheSum #= 126, % hard coded...

  % all_different(F),
  all_distinct(F),

  % symmetry breaking.
  % Note :I'm not sure that all these are correct...
  % F[1] #< F[3], F[1] #< F[7], F[1] #< F[9], F[1] #< F[19], F[1] #< F[21],  
  % F[1] #< F[25], F[1]#< F[27], F[2] #> F[4], F[4] #> F[10],

  sum([F[1],F[2],F[3],F[4],F[5],F[6],F[7],F[8],F[9]]) #= TheSum,
  sum([F[10],F[11],F[12],F[13],F[14],F[15],F[16],F[17],F[18]]) #= TheSum,
  sum([F[19],F[20],F[21],F[22],F[23],F[24],F[25],F[26],F[27]]) #= TheSum,
  sum([F[1],F[4],F[7],F[10],F[13],F[16],F[19],F[22],F[25]]) #= TheSum,
  sum([F[2],F[5],F[8],F[11],F[14],F[17],F[20],F[23],F[26]]) #= TheSum,
  sum([F[3],F[6],F[9],F[12],F[15],F[18],F[21],F[24],F[27]]) #= TheSum,
  sum([F[1],F[2],F[3],F[10],F[11],F[12],F[19],F[20],F[21]]) #= TheSum,
  sum([F[4],F[5],F[6],F[13],F[14],F[15],F[22],F[23],F[24]]) #= TheSum,
  sum([F[7],F[8],F[9],F[16],F[17],F[18],F[25],F[26],F[27]]) #= TheSum,
  sum([F[1],F[5],F[9],F[10],F[14],F[18],F[19],F[23],F[27]]) #= TheSum,
  sum([F[3],F[5],F[7],F[12],F[14],F[16],F[21],F[23],F[25]]) #= TheSum,
  sum([F[1],F[4],F[7],F[11],F[14],F[17],F[21],F[24],F[27]]) #= TheSum,
  sum([F[3],F[6],F[9],F[11],F[14],F[17],F[19],F[22],F[25]]) #= TheSum,
  sum([F[7],F[8],F[9],F[13],F[14],F[15],F[19],F[20],F[21]]) #= TheSum,
  sum([F[1],F[2],F[3],F[13],F[14],F[15],F[25],F[26],F[27]]) #= TheSum,


  % From Pitrat's blog post
  % http://bootstrappingartificialintelligence.fr/WordPress3/2015/01/caia-surprised-me-part-i/

  % % The three constraints for the horizontal planes are:
  % F[1]+F[2]+F[3]+F[4]+F[5]+F[6]+F[7]+F[8]+F[9] #= TheSum,
  % F[10]+F[11]+F[12]+F[13]+F[14]+F[15]+F[16]+F[17]+F[18] #= TheSums,
  % F[19]+F[20]+F[21]+F[22]+F[23]+F[24]+F[25]+F[26]+F[27] #= TheSums,

  % % The six constraints for the vertical planes:
  % F[1]+F[2]+F[3]+F[10]+F[11]+F[12]+F[19]+F[20]+F[21] #= TheSums,
  % F[4]+F[5]+F[6]+F[13]+F[14]+F[15]+F[22]+F[23]+F[24] #= TheSums,
  % F[7]+F[8]+F[9]+F[16]+F[17]+F[18]+F[25]+F[26]+F[27] #= TheSums,
  % F[1]+F[4]+F[7]+F[10]+F[13]+F[16]+F[19]+F[22]+F[25] #= TheSums,
  % F[2]+F[5]+F[8]+F[11]+F[14]+F[17]+F[20]+F[23]+F[26] #= TheSums,
  % F[3]+F[6]+F[9]+F[12]+F[15]+F[18]+F[21]+F[24]+F[27] #= TheSums,

  % % The last six constraints are for the diagonal planes:
  % F[1]+F[2]+F[3]+F[13]+F[14]+F[15]+F[25]+F[26]+F[27] #= TheSums,
  % F[7]+F[8]+F[9]+F[13]+F[14]+F[15]+F[19]+F[20]+F[21] #= TheSums,
  % F[1]+F[4]+F[7]+F[11]+F[14]+F[17]+F[21]+F[24]+F[27] #= TheSums,
  % F[3]+F[6]+F[9]+F[11]+F[14]+F[17]+F[19]+F[22]+F[25] #= TheSums,
  % F[1]+F[5]+F[9]+F[10]+F[14]+F[18]+F[19]+F[23]+F[27] #= TheSums,
  % F[3]+F[5]+F[7]+F[12]+F[14]+F[16]+F[21]+F[23]+F[25] #= TheSums, 

  % stats(F,Size,Degree,Min,Max),

  println(F),

  Vars = F,

  % solve($[ff,updown], Vars),
  solve($[ff,split], Vars),

  println(theSum=TheSum),
  println(F),

  fail,

  nl.



% only for the cp solver
% stats(Vars, Size,Degree,Min,Max) =>
%     Size   = [fd_size(V) : V in Vars],
%     Degree = [fd_degree(V) : V in Vars],
%     Min    = [fd_min(V) : V in Vars],
%     Max    = [fd_max(V) : V in Vars],
%     % println(vars=Vars),
%     println(size=Size),
%     println(degree=Degree),
%     println(min=Min),
%     println(max=Max),
%     nl.