/* 

  Recover origin of pair differences in Picat.

  From https://twitter.com/SturnioloSimone/status/1214588457135329280
  """
  Suppose there's a certain set of numbers - x_1, x_2, ..., x_n. You don't know n. 
  Suppose you only know certain linear combinations of them, L_ij = x_i-x_j, but you DON'T know 
  i or j. What's an algorithm to find as many of the xs as possible just from looking 
  at the Ls? Ideas?

  ... 
  (a follow up tweet:)
  There are four numbers. I can tell you that 1, 3, and 6 are differences between pairs of 
  numbers chosen among them (you don't know which pairs). Can you tell me what the numbers are?

  Something like that, except you don't know it's four.
  """

  A better description of the problem and findings can be found at
  http://hakank.org/picat/Recover_origin_of_pair_differences.pdf


  DESCRIPTION OF THE PROBLEM AND SOME FINDINGS
  ============================================
  Note: This description is the one that I (hakank) has interpreted from
        the tweet above. Perhaps this exact model (with all its 
        assumptions and symmetry breakings) is not that useful.

  [This part mainly explains the variant where we assume that the list
  of differences don't contain any duplicates; these are tested via
  go/0 where AllowDuplicates = false. 

  Most of the same applies to where we allow duplicates, tested in 
  go2/9 with AllowDuplicates = true.
  ]

  The objective of this model is: Given a list of pairs differences
  we try to recover the list that generated this list of differences.
  Note that we don't know anything about this generating array:
   - not the length
   - not the maximum or minimum values
  The only thing we have is the list of pair differences.
  
  This model finds a solution that is of the shortest length by 
  looping though possible lengths and whenever it finds a solution
  it is the shortest (one of possible many solutions).


  Example:
  We are given the following list of pair of differences:
  
  diffs: [4, 7, 9, 11, 13, 15, 16, 18, 20, 22, 24, 29, 31, 33, 35, 38, 42, 44, 51]

  (It happens to be generated by the following array: [5,14,18,34,38,49,56] but the
  model don't know this.)

  The model outputs the following solutions:
   x: [1, 10, 14, 30, 34, 45, 52]
  
  See below for the full model, but some of the main constraints
  are:
   - it is sorted in increasing order and has unique numbers
   - it starts with 1
   - the last value is 1 + the maximum value of the difference list
  
  and with the more complex constraints we ensure:
   - the solution array must generate all the differences in the given
     difference list
   - the solution array must generate only the differences in the
     given difference list.

  Note: It is in general not possible to restore the original
  array that generated the differences, just a solution that
  satisfies the constraints above.

  The length of the solutions is - what we have seen - often
  the same length of the original array. There are cases where
  the solution is much shorter, e.g. for the extremal cases
  where the differences is a list of all number from 1..<some n>.

  Number of solutions
  - - - - - - - - - -
  The problem shown above has at least 2 solutions:
    x: [1, 8, 19, 23, 39, 43, 52]
    x: [1, 10, 14, 30, 34, 45, 52]

  These two solutions has the following list differences:
   x_diffs: [7, 11, 4, 16, 4, 9]
   x_diffs: [9, 4, 16, 4, 11, 7]
  i.e. it's the same list but they are reversed of each other.
  If we only want 1 solution, then the model adds an constraint
  that
     x[2]-x[1] < x[n]-x[n-1]

  Some more comments on the list differences:
  For the example above, i.e. the original array of [5, 14, 18, 34, 38, 49, 56],
  it has the following differences:
     [9, 4, 16, 4, 11, 7]

  which is exactly the same as the differences of the solution
    x: [1, 10, 14, 30, 34, 45, 52]
    x_diffs: [9, 4, 16, 4, 11, 7]

  The other solution has a difference list that is the reversal
    x: [1, 8, 19, 23, 39, 43, 52]
    x_diffs: [7, 11, 4, 16, 4, 9]

  See below for more about the number of solutions.


  Some further explorations: sum of the differences
  - - - - - - - - - - - - - - - - - - - - - - - - - 
  For other problem instances there might be a huge number of
  solutions. For example the problem that is generated by the array 
     [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15],
  has 130 solutions, e.g.
     x = [1,2,6,9,11,13,15]
     x = [1,2,4,7,11,14,15]
     x = [1,3,5,6,8,14,15]
     x = [1,3,5,7,10,14,15]
     x = [1,2,5,6,8,13,15]

  Interestingly, the sum of the list differences are all 14.

  The number of solutions for the generator array of 1..20 is 326,
  and the sum of list differences are all 19, i.e. one less than the 
  number of elements.

  A generator of [1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39]
  also give 326 solutions, with the sum of list differences of 38.
  However, for most random instances there seems to be mostly 2 solutions
  (with symmetric list differences of the solution array).

  Invariants: sum of list differences 
  - - - - - - - - - - - - - - - - - -
  It seems that the sum of list differences is an invariant:
   - the sum of the list differences in the solution is the same as
     the sum of the list differences of the generator array
  For example, the generator array of [16,41,43,54,56,62,66,67,84,95]
  has a sum of list differences of 79, which is the same as for the 2 
  solutions:
    l = [16,41,43,54,56,62,66,67,84,95]
    l_len = 10
    l_shifted = [1,26,28,39,41,47,51,52,69,80]
    l_diffs = [25,2,11,2,6,4,1,17,11]
    l_diff_sum = 79
    diffs = [1,2,4,5,6,8,10,11,12,13,15,17,18,19,21,22,23,24,25,26,27,28,29,30,33,38,39,40,41,43,46,50,51,52,54,68,79]
    diffLen = 37
    mode = all
    n = 9
    n = 10
    x = [1,12,29,30,34,40,42,53,55,80]
    x_shifted = [1,12,29,30,34,40,42,53,55,80]
    x_diffs = [11,17,1,4,6,2,11,2,25]
    sum_x_diffs = 79
    check = [1,2,4,5,6,8,10,11,12,13,15,17,18,19,21,22,23,24,25,26,27,28,29,30,33,38,39,40,41,43,46,50,51,52,54,68,79]

    x = [1,26,28,39,41,47,51,52,69,80]
    x_shifted = [1,26,28,39,41,47,51,52,69,80]
    x_diffs = [25,2,11,2,6,4,1,17,11]
    sum_x_diffs = 79
    check = [1,2,4,5,6,8,10,11,12,13,15,17,18,19,21,22,23,24,25,26,27,28,29,30,33,38,39,40,41,43,46,50,51,52,54,68,79]
 
  Relation between sum of differences of generating array, the differences and the solution
  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  Another invariant seems to be a relation between the sum of the differences of the generating array
  (which as we saw above is the same as the sum of differences of the solution) and the sum of
  the differences of the differences array (i.e. the input array of the problem).
    the sum of the differences of the generator == the sum of the differences of the difference array - diffs[1]
  For example (slightly edited):
     l = [27,29,39,44,50,52,61,64,67,86]
     l_len = 10
     l_diffs = [2,10,5,6,2,9,3,3,19]
     l_diff_sum = 59
     diffs = [2,3,5,6,8,9,10,11,12,13,14,15,17,19,20,21,22,23,25,28,32,34,35,36,37,38,40,42,47,57,59]
     diffs_diff_sum = 57
     difference_of_diff_sums = 2

     i.e. l_diff_sum = diffs_diff_sum + diffs[1] = 57 + 2 = 59
     which is the same as
          x_diff_sum + diffs[1] = 59


  Number of solutions for differences with 1..Diffs.len (perfect rulers)
  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  A perfect ruler is a sub-problem of this problem, i.e. when the
  differences are "extremal": a list of all numbers from 1..Diffs.len,
  and without duplicates.

  Here is the (empirical) number of solutions (for the shortest length of solution). 
  Note: The cp solver was used for these experiments.
    N  #solutions for diffs=1..N
    ----------------------------
    1   1
    2   1
    3   1
    4   2
    5   3
    6   4
    7   2
    8  12
    9   8
   10   4
   11  38
   12  30
   13  14
   14   6
   15 130
   16  80
   17  32
   18  12
   19 500
   20 326
   21 150
   22 66

  For OEIS: 1,1,2,3,4,2,12,8,4,38,30,14,6,130,80,32,12,500,326,150
  And not surprisingly it's: "Number of perfect rulers with length n (n>=0)"
  https://oeis.org/search?q=1%2C1%2C2%2C3%2C4%2C2%2C12%2C8%2C4%2C38%2C30%2C14%2C6%2C130%2C80%2C32%2C12%2C500%2C326%2C150&language=english&go=Search

  For more about perfect rulers, see:
    http://mathworld.wolfram.com/PerfectRuler.html

  A related problem is Golomb ruler: https://en.wikipedia.org/wiki/Golomb_ruler
  which requires that the solution is optimal; this is not handled in this model
  (the minimization mode is not really the same as Golomb ruler). Perhaps this will
  explored in the future (I don't want to side tracks too much now).

  And somewhere in thiw Wiki post might be some more good ideas: 
  https://oeis.org/wiki/User:Peter_Luschny/PerfectRulers


  Allowing duplicates
  - - - - - - - - - -
  In go2/0 we do some experiments where duplicates are allowed in the generating array,
  i.e. that the difference list contains a 0.

  Here in an example where the generating list has a couple of duplicates:
     - 7, 23, 37, and 38

  Note that there are only one 0 which "collects" all these four duplicates.
  As a result, the solution will also only contain one duplicate, here 37.
  When we don't allow duplicates, the length of the solutions is often near the
  length of the generating list (apart from some extreme cases where it's much shorter).

  With duplicates the solution list is always much shorter.
    
   l = [2,4,7,7,8,12,17,18,20,21,23,23,25,29,30,36,37,37,38,38]
   l_len = 20
   diffs = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]
   n = 9
   n = 10
   n = 11
   x = [1,2,4,7,14,21,28,32,36,37,37]
   x_shifted = [1,2,4,7,14,21,28,32,36,37,37]
   x_diffs = [1,2,3,7,7,7,4,4,1,0]
   sum_x_diffs = 36
   check = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]

  A little extreme example is the generating list of 10 pairs of duplicates:
     [1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10]
  with the difference list
     [0,1,2,3,4,5,6,7,8,9]
  The solution compresses down to the problem of a difference list of 0,1..9,
  and is just 6 elements with one duplicate:
     [1,2,5,8,10,10]


  AllowDuplicates = experimental
  -  -  -  -  -  -  -  -  -  -  

  The is a variant of AllowDuplicates, experimental, where we don't remove duplicates in the Diffs
  list. With AllowDuplicates=true, the only change is that we allow (one) 0 in the duplicate list.
  With experimental we are keeping all the differences from the gerating list. 
  The number of difference in this variant will thus be 
     NumDiffs = (Len*(Len-1)) div 2
  where Len is the length of the generating list. E.g. for Len=20 then NumDiffs = 20*19 div 2 = 190.

  This feature is tested in go3/0.

  Here is an example with AllowDuplicates=experimental:
  - Generating list: [7,10,11,12,12,12,13,15,15,17]
    Note that we have a couple of duplicates: 12,12,12 and 15,15
  - Differences: 
    [0,0,0,0,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,5,5,5,5,6,6,7,8,8,10]
  - Solution:
    [1,4,5,6,6,6,7,9,9,11]
    and it has the same differences as above:
    [0,0,0,0,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,5,5,5,5,6,6,7,8,8,10]

    The solution has as many duplicates as the generating list, though probably not exactly the 
    same numbers, here 6,6,6 and 9,9

  The additional constraint used here is that the difference list of the solution (X) 
  must have exactly the same number of occurrences as in the difference list, i.e.
     - 4 0s
     - 6 1s
     - 8 2s
     - etc.

  Note that this extra constraint might slow down the model.


  DESCRIPTION OF THE PICAT MODEL
  =============================
  Note: The PEQNP model (http://hakank.org/peqnp/linear_combinations.py ) use the same general
        approach, but the details (especially the syntax) varies little.

  What follows is a description of this program (model as we say regarding constraint modelling)
  as an introduction to constraint modelling. For about modelling in Picat, see our
  book "Constraint Solving and Planning with Picat" (available free from 
  http://picat-lang.org/picatbook2015.html )

  Please contact me (hakank@gmail.com) if you have more questions about this.

  Note: After some iterations of this model, I added support for allowing duplicates.
        The model is now handling this by checking the boolean variable AllowDuplicates.
        There are only a few cases where this addition of functionality have to been 
        taking care of.

  Picat (http://picat-lang.org/ ) is a multi-paradigm system supporting logic programming,
  constraint modelling (CP, SAT, MIP, and SAT solvers), tabling, functions as well as
  for loops, while loops, reassignments.

  This model is a constraint model, i.e. use a declarative paradigm of stating all the
  constraints that must be satisfied. 

  The predicate go/0 is the main program, it generates a random sample of a certain size, e.g
    [L, Diffs] = generate_problem(25)

  which first randomize the main list (L) and then takes the differences (+ sorting
  and removing duplicates). Note that L is not used in the solving part, it's only
  used for presentation.

  Then the (foreach) loop loops over different N's, starting with the minimum possible 
  length of the solution. Here we calculate this minimum length as

    floor(1+sqrt(1+8*Len)/2)

  found by solving the equation for how many pairs there are in a list of length Len
     Len == (N*(N-1)) div 2

  Note that for some difference lists the solution list can be quite small, e.g. 
  the solution for the differences of a list [1,2,3,...30] with the differences
  [1,2,..29] (i.e. a length of 29) the solution is much smaller, 9 values. Here is
  an example solution:

     x = [1,2,4,7,14,21,25,29,30]

  The maximum length of the solution is - completely arbitrary - set to Len*4+1.

  Note that this model do the following:
   - set the first value of x to 1
   - set the last value of x to max(diffs)+1

  For each loop it check if linear_comb/3 has a solution, it so it ends the program;
  if not it increases the length (N) by 1 and continue to run.


  Constraints
  -----------

  The constraints are collected in the linear_comb(N,Diffs, X) procedure where

  the inputs are
    - N, the length of the array of decision variables to create
    - Diffs, the list of differences that we try to find

  Terms used:
   - decision variable: A variable with a certain domain and contains a value
     that the model will assign to during the running of the model.
   - domain: Each decision variable must have a specifc domain, i.e. an interval
     of integers that are valid in the model. 
   - global constraint: a certain kind of constraint that is over a collection of
     decision variables. These are often highly optimized, and also they are
     a great tool when modeling a problem since they are "high level patterns".
    
  The output variable is

    - X, a list of decision variables of size N

   Here is a walktrough of the code.

   * X = new_list(N)
     Create the list of decision variables of size N

   * Max = max(Diffs)*2,
     X :: 1..Max, 

     Here the upper bound of the decision variables are calculated, Max. It 
     is - quite arbitrary - set to 2*the max value in the difference list. 
     It would be great with a better (=smaller) max value of this domain.
      
     The X :: 1..Max then ensures that all the decision variables in the list X
     must be contained in this interval.
      
   * all_distinct(X),
     increasing(X),

     These two global constraints ensure that
      - all_distinct(X): all the values in the list X must be distinct.
        There is a variant - all_different(X) - which is faster but don't propagate
        the constraints as strong as all_distinct(X).
      - increasing(X): ensure that the list X is ordered, in increasing order.
        This constraint is for symmetry breaking which often speeds up the solving.

   * X[1] #= 1
     This assigns the value of 1 to the first element of X. Noticed the special
     operator (#=) which means that it is a constraint (operates on the constraint store).
     This constraint is for symmetry breaking and to ensure that we have the number 1
     in the solution, which - hopefully - makes the model faster

   
  * % All the differences must be in the Diff list
    Ts = [],
    foreach(I in 1..N, J in I+1..N)
       T :: Diffs,  
       T #= abs(X[I]-X[J]),
       Ts := Ts ++ [T]
    end,

    This loop ensure that all the differences that can be construced from the list X
    must be in the list Diffs. I.e. we do not allow any difference that is not in
    the given difference list.

    In more detail: for each pair of values in X (I from 1..N, J from I+1..N) we
    calculate the absolute value of the difference (as a decision variable)
      T #= abs(X[I]-X[J]) 
    and ensure that this difference T is in the Diffs list
      T :: Diffs

    Also, we add the T's to the Ts list since the solver need to know all the 
    decision decision variables. The construct
       Ts := Ts ++ [T]
    is Picat's way of adding a value to a list. Note that we must use the 
    reassignment operator (':=') since otherwise (if we used '=') it will 
    check for equality or binds the variable (the Prolog way).

  * The next secion is the main part of the model, i.e. ensure that we conver 
    all the differences in the Diffs list. (The previous part ensured the other
    way around, i.e. that the pair differences in X must be in Diffs list.)

    % Ensure that we cover all the differences
    IJs = [],
    foreach(D in Diffs)
       % Create the indices I,J where I < J
       I :: 1..N,
       J :: 1..N,
       I #< J,
     
       element(I,X,Ix),
       element(J,X,Jx),
       % Ensure that this difference is covered
       D #= abs(Ix-Jx),
       IJs := IJs ++ [I,J]
    end,

    Here we loop over all the values (D) in the difference list (Diffs). The main point
    is that we want to find two indices (I and J where I < J)) such that the difference
      abs(X[I]-X[J]) 
    is the difference value D.

    Since I and J are decision variables, we have - as usual - define their domains
    
       I :: 1..N,
       J :: 1..N,

    i.e. that they have values as indices in the X list (which has size N).
    
    However, Picat don't handle when the indices of a list is a decision variable so we
    have to use the element/3 predicate:

       element(Index, List, Result) 

   is Picat's (and most other logic constraint programming systems') version of

      Result = List[Index]

   Thus
  
       element(I,X,Ix),
       element(J,X,Jx),
      
   means
       Ix = X[I]
   and
       Jx = X[J]
 
   The variables Ix and Jx are thus the two values and we then ensure that
      D #= abs(X[I]-X[J])
    
   As before, we collect the decision variables in a list to the solver.

 * The only thing that is now left is to call the solver:

     Vars = X ++ Ts ++ IJs,

   Collect all the decision variables into one big list.

   And call the solver:

    solve($[ffd,updown],Vars).  

   Note: the "ffd" and "updown" are search heuristics for the constraint programming solver
   (they are ignored by the other solvers). A big part of modelling is to find the best combinations 
   of search heuristics. But since we use the SAT solver (or perhaps the SMT solver) it's not
   interesting to go into more here.
    




  Some example runs of the program.

diffs = [2,4,6,8]
n = 5
[1,3,5,9]

l = [5,14,18,34,38,49]
diffs = [4,9,11,13,15,16,20,24,29,31,33,35,44]
n = 2
max = 88
n = 3
max = 88
n = 4
max = 88
n = 5
max = 88
n = 6
max = 88
[15,26,30,46,50,59]


l = [17,22,27,44,50,52,70,82,84,87]
diffs = [2,3,5,6,8,10,12,14,17,18,20,22,23,25,26,27,28,30,32,33,34,35,37,38,40,43,48,53,55,57,60,62,65,67,70]
n = 2
max = 140
...
n = 10
max = 140
[70,75,80,97,103,105,123,135,137,140]


l = [22,34,37,39,61,68,79,82,84]
diffs = [2,3,5,7,11,12,14,15,16,17,18,21,22,23,24,27,29,31,34,39,40,42,43,45,46,47,48,50,57,60,62]
...
n = 9
max = 124
[1,13,16,18,40,47,58,61,63]

l = [3,23,27,39,42,53,62,88,92,99]
diffs = [3,4,7,9,11,12,14,15,16,19,20,23,24,26,30,35,36,37,39,46,49,50,53,57,59,60,61,65,69,72,76,85,89,96]
...
n = 10
max = 192
[1,21,25,37,40,51,60,86,90,97]


l = [1,8,60,72,76,78,79,84,89]
diffs = [1,2,3,4,5,6,7,8,10,11,12,13,16,17,18,19,24,29,52,59,64,68,70,71,75,76,77,78,81,83,88]
n = 2
n = 3
n = 4
n = 5
n = 6
n = 7
n = 8
n = 9
x = [1,6,11,12,14,18,30,82,89]


l = [5,10,29,45,46,56,58,80,83,93]
diffs = [1,2,3,5,10,11,12,13,16,17,19,22,24,25,27,29,34,35,36,37,38,40,41,46,47,48,51,53,54,64,70,73,75,78,83,88]
n = 2
n = 3
n = 4
n = 5
n = 6
n = 7
n = 8
n = 9
n = 10
x = [1,6,25,41,42,52,54,76,79,89]


l = [10,11,15,17,27,29,50,75,90]
diffs = [1,2,4,5,6,7,10,12,14,15,16,17,18,19,21,23,25,33,35,39,40,46,48,58,60,61,63,64,65,73,75,79,80]
n = 2
n = 3
n = 4
n = 5
n = 6
n = 7
n = 8
n = 9
x = [1,2,6,8,18,20,41,66,81]
check = [1,2,4,5,6,7,10,12,14,15,16,17,18,19,21,23,25,33,35,39,40,46,48,58,60,61,63,64,65,73,75,79,80]


l = [16,22,28,36,50,60,63,64,78,84,87,91,93,94]
l_len = 14
diffs = [1,2,3,4,6,7,8,9,10,12,13,14,15,16,18,20,21,22,23,24,27,28,29,30,31,32,33,34,35,36,37,38,41,42,43,44,47,48,50,51,55,56,57,58,59,62,63,65,66,68,69,71,72,75,77,78]
n = 11
n = 12
n = 13
n = 14
x = [1,2,4,8,11,17,31,32,35,45,59,67,73,79]
check = [1,2,3,4,6,7,8,9,10,12,13,14,15,16,18,20,21,22,23,24,27,28,29,30,31,32,33,34,35,36,37,38,41,42,43,44,47,48,50,51,55,56,57,58,59,62,63,65,66,68,69,71,72,75,77,78]


SAT: 1min 47.32s
SMT: 56min 12.46s
CP: >4hours
l = [16,22,27,28,36,37,41,50,60,63,64,73,78,84,87,91,93,94]
l_len = 18
diffs = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,41,42,43,44,45,46,47,48,50,51,52,53,54,55,56,57,58,59,60,62,63,64,65,66,67,68,69,71,72,75,77,78]
diffLen = 69
n = 12
n = 13
n = 14
n = 15
x = [1,2,4,10,15,17,22,38,46,47,60,69,70,73,79]
check = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,41,42,43,44,45,46,47,48,50,51,52,53,54,55,56,57,58,59,60,62,63,64,65,66,67,68,69,71,72,75,77,78]


SAT: 4min25.59s
SMT: >5h
l = [12,16,22,27,28,30,36,37,41,50,60,63,64,68,69,73,78,83,84,87,91,93,94]
l_len = 23
diffs = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,71,72,75,77,78,79,81,82]
diffLen = 75
n = 13
n = 14
n = 15
n = 16
x = [1,2,4,5,6,8,16,28,30,43,52,61,62,68,73,83]
check = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,71,72,75,77,78,79,81,82]


l = [3,4,7,14,16,22,30,33,44,46,48,50,51,55,57,69,71,73,79,82,89,91,96]
l_len = 23
diffs = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,72,73,74,75,76,77,78,79,80,82,84,85,86,87,88,89,92,93]
diffLen = 88
n = 14
n = 15
n = 16
n = 17
x = [1,6,8,18,20,26,37,48,59,70,80,81,85,86,90,93,94]
check = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,72,73,74,75,76,77,78,79,80,82,84,85,86,87,88,89,92,93]


SAT: 51.5s
CP: 2s
l = [3,5,13,19,37,48,56,59,60,64,66,74,77,78,80,86,90,98]
l_len = 18
diffs = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,26,27,29,30,31,32,34,35,37,38,39,40,41,42,43,45,46,47,49,50,51,53,54,55,56,57,58,59,61,63,64,65,67,69,71,72,73,74,75,77,79,81,83,85,87,93,95]
diffLen = 71
n = 12
n = 13
n = 14
n = 15
n = 16
x = [1,3,11,17,35,46,57,58,62,64,72,75,76,84,88,96]
check = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,26,27,29,30,31,32,34,35,37,38,39,40,41,42,43,45,46,47,49,50,51,53,54,55,56,57,58,59,61,63,64,65,67,69,71,72,73,74,75,77,79,81,83,85,87,93,95]

SAT: 1min 41s
CP: 1.7s
L =  [40, 108, 122, 136, 141, 143, 146, 188, 192, 205, 259, 266, 287, 312, 330, 341, 349, 381, 386, 397, 429, 459, 475, 496]

SAT: 5.3s
CP: 12min 30s!
L = [4, 7, 9, 11, 13, 15, 16, 18, 20, 22, 24, 29, 31, 33, 35, 38, 42, 44, 51]

  Also see http://hakank.org/picat/linear_combinations2.pi for a 0/1 model which seems
  to be quite faster, at least when max(Diffs) are not too large. For larger
  max(Diffs) this current model might be faster.


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

*/


% import sat. % Seems to be faster than CP and SMT on larger problems.
import cp. % Very fast (often fastest) on small problem (up to say 30-40 difference values), then very slow.
% import smt. % Seems to be quite fast on small problems but often not as fast as cp. z3 is faster than cvc4
% import mip. % slow


main => go.


%
% In this experiment we assume that there are no duplicates in the generating list, i.e.
% that the duplicate list don't contains 0 (and has no duplicates).
%
% See go2/0 for experiments with duplicates.
% See go3/0 for experiments with duplicates in the difference lists.
% 
go ?=>
  nolog,
  _ = random2(), % comment this for a static random seed.

  % Generate a problem:
  % L is the (random) generated list, Diffs are the differences of L, sorted and with duplicates remove
  % L = generate_problem(30,200), % remove duplicates

  % L := [99,100,101],

  % L := [I : I in 1..10],
  % L := [I : I in 1..21],
  % L := [2,3,5,7,11,13,17],
  % L := [sum([J : J in 1..I]) : I in 1..10],

  % L := [4,6,7],
  % L := [16,22,27,28,36,37,41,50,60,63,64,73,78,84,87,91,93,94],

  % L := [5, 14, 18, 34, 38, 49],
  % L := [5, 13, 19, 37, 48, 60, 86, 90, 98],
  % L := [13,28,46,67,71,73,75,81,91,96],
  % L := [4,10,12,19,27,28,34,41,43,44,52,53,60,61,77,80,96],
  % L := [1, 2, 5, 10, 15, 26, 37, 48, 54, 60, 66, 67, 68],
  % L := [3, 5, 13, 19, 37, 48, 56, 59, 60, 64, 66, 74, 77, 78, 80, 86, 90, 98],
  % L := [4, 7, 9, 11, 13, 15, 16, 18, 20, 22, 24, 29, 31, 33, 35, 38, 42, 44, 51],
  % L := [40, 108, 122, 136, 141, 143, 146, 188, 192, 205, 259, 266, 287, 312, 330, 341, 349, 381, 386, 397, 429, 459, 475, 496],
  % L := [16,41,43,54,56,62,66,67,84,95],
  % L :=  [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,173,174,175,176,177,178,179,180,181,182,183,184,185,187,188,189,190],

  % L := [1,4,8,9],
  % L := [1,3,8],
  % L:= [16,41,43,54,56,62,66,67,84,95],
  L :=  [40,108,122,136,141,143,146,188,192],
  Diffs = diffs(L),
  Mode = all,
  AllowDuplicates = false,
  do_experiment(L, Diffs, Mode, AllowDuplicates),
  nl.

go => true.


%
% In this experiment we assume that there might be duplicates in the generating list, i.e.
% the difference list might include 0. However, there is no duplicates in the difference list. 
%
% See go/0 for experiments without duplicates.
% See go3/0 for experiments with duplicates in the difference list
% 
go2 ?=>
  nolog,
  _ = random2(), % comment this for a static random seed.

  % Generate a problem:
  % L is the (random) generated list, Diffs are the (sorted) differences of L. Duplicates are kept.
  % Note: To force duplicates: use a quite low max_val, e.g. just length*2  
  L = generate_problem2(20,30),

  % true: allow duplicates but there is no duplicates in Diffs
  % experiental: allow duplicates in Diffs
  % AllowDuplicates = true,
  AllowDuplicates = true,  
  Mode = first,
  if AllowDuplicates == experimental then
    Diffs = diffs2(L)
  else
    Diffs = diffs(L)
  end,
  
  do_experiment(L, Diffs, Mode, AllowDuplicates),
 
  nl.

go2 => true.


%
% In this experiment we assume that there might be duplicates in the generating list
% (the difference list might include 0) and thus that the solution might have duplicates.
% The AllowDuplicates = experimental yield _all_ the possible differences in the Diffs lists.
%
% See go/0 for experiments without duplicates.
% See go2/0 for experiments without duplicates in the difference list.
% 
go3 ?=>
  nolog,
  _ = random2(), % comment this for a static random seed.

  % Generate a problem:
  % L is the (random) generated list, Diffs are the (sorted) differences of L. Duplicates are kept.
  % Note: To force duplicates: use a quite low max_val, e.g. just length*2  
  L = generate_problem2(20,30),

  % true: allow duplicates but there is no duplicates in Diffs
  % experiental: allow duplicates in Diffs
  AllowDuplicates = experimental,  
  Mode = first, % or all
  if AllowDuplicates == experimental then
    Diffs = diffs2(L)
  else
    Diffs = diffs(L)
  end,
  do_experiment(L, Diffs, Mode, AllowDuplicates),
 
  nl.

go3 => true.


%
% Here we can do real experiment with differences lists without a known origin
%
go4 =>
  nolog,
  % Diffs = [4,9,11,13,15,16,20,24,29,31,33,35,44],
  % Diffs := [3,4,7,9,11,12,14,15,16,19,20,23,24,26,30,35,36,37,39,46,49,50,53,57,59,60,61,65,69,72,76,85,89,96],
   

  % Diffs :=  [1,3,6], % Not possible
  % Diffs := [1,2,4], % Not possible
  % Diffs := [1,2,3],
  %% Diffs := [1,3,6], % Note: this is the example given in the tweet. But it's not a possible difference list!

  % Diffs := [1,5,6,10,11,12], % This has only one solution
  % Diffs := [1,2,5,6,7,12], %
  % Diffs :=   [4, 7, 9, 11, 13, 15, 16, 18, 20, 22, 24, 29, 31, 33, 35, 38, 42, 44, 51],
  % Diffs := [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 61, 62, 63, 66], 


  % The example from whuber's example @ https://or.stackexchange.com/questions/3349/restoring-a-list-from-differences/3555
  %%%Diffs := [0,1,3,4,5,7,8,9],
  % Diffs := [1,3,4,5,7,8,9],

  Diffs:= [1,6,7,8,12,13,14,15,19,20,21,26,27,33,34],

  println(diffs=Diffs),
  Len = Diffs.len,
  println(diffsLen=Diffs.len),
  
  MinLen = floor(1+sqrt(1+8*Len)/2), % DiffLen = (N*(N-1)) div 2
  Mode = all,
  AllowDuplicates = all,
  Found = false,  
  foreach(N in MinLen..MinLen*2, Found == false)
    println(n=N),
    if linear_comb(Diffs,N,Mode,AllowDuplicates, Sol), Sol \= [] then
      if Mode == all then
        foreach(XX in Sol)
          print_solution(Diffs,XX,AllowDuplicates)
        end,
        println(num_solutions=Sol.len)
      else
        print_solution(Diffs,Sol,AllowDuplicates)
      end % ,
      % Found := true    
    end
  end,
  nl.


%
% Experiment: What is the lengths of solutions and total number of solutions
% for the "extremal" difference, i.e. the difference list 1..Len.
%
% Note: Here we use the Diffs = 1..Len (not Orig list = 1..Len).
%
% The total number of solutions for 1..10: 1,1,3,4,9,17,33,63,128,248
% And it's the expected https://oeis.org/A103295: "Number of complete rulers with length n."
%
% The number of solutions for a certain Len and N:
%
% [len =  1,counts = [1]
% [len =  2,counts = [1]
% [len =  3,counts = [2,1]
% [len =  4,counts = [3,1]
% [len =  5,counts = [4,4,1]
% [len =  6,counts = [2,9,5,1]
% [len =  7,counts = [12,14,6,1]
% [len =  8,counts = [8,27,20,7,1]
% [len =  9,counts = [4,40,48,27,8,1]
% [len = 10,counts = [38,90,75,35,9,1]
% [len = 11,counts = [30,134,166,110,44,10,1]]
% [len = 12,counts = [14,166,311,277,154,54,11,1]]
% [len = 13,counts = [6,166,490,590,431,208,65,12,1]]
% [len = 14,counts = [130,660,1096,1022,639,273,77,13,1]]
% [len = 15,counts = [80,800,1816,2132,1662,912,350,90,14,1]]
% [len = 16,counts = [32,790,2641,3980,3796,2574,1262,440,104,15,1]]
% [len = 17,counts = [12,710,3476,6721,7796,6371,3836,1702,544,119,16,1]]
% [len = 18,counts = [500,4140,10445,14655,14184,10208,5538,2246,663,135,17,1]]
% [len = 19,counts = [326,4380,14791,25322,28880,24394,15746,7784,2909,798,152,18,1]]
% [len = 20,counts = [150,4064,19107,40545,54419,53297,40141,23530,10693,3707,950,170,19,1]]
% [len = 21,counts = [66,3504,23236,60922,95823,107912,93458,63672,34223,14400,4657,1120,189,20,1]]
%
% They are all (lumped together) in sequence https://oeis.org/A103294
% "Triangle T, read by rows: T(n,k) = number of complete rulers with length n and k segments (n >= 0, k >= 0)."
% 
% Note that this is just about complete rulers, which is only a very special case of the
% problem we study, i.e. recovering a list based on its differences.
%
% Also: the cp solver is much faster on this.
go5 =>
   nolog,
   TotalSums = [],
   % Mode = all,
   AllowDuplicates = false,
   
   foreach(Len in 1..21)
     println(len=Len),
     % Diffs = diffs(1..Len),
     Diffs = 1..Len,     
     println(diffs=Diffs),
     MinLen = floor(1+sqrt(1+8*Len)/2), % DiffLen = (N*(N-1)) div 2
     MaxLen = Len+1, % 
     println(lens=[MinLen,MaxLen]),
     FoundAll = false,
     Counts = [],
     Total = 0,
     % There are no more solutions for Len+2, so we limit to Len+1.
     % The last solution is always 1 (for the complete 1..Len)
     foreach(N in MinLen..MaxLen,FoundAll==false)
        garbage_collect(200_000_000),

        % Later note: Using Model=all is now much faster in earlier version of
        %             this program. 
        % Sol = _,
        % if linear_comb(Diffs,N,all, AllowDuplicates, Sol) then
        %    % println(sols=Sol),
        %    Count = Sol.len,
        %    println(N=Count)
        %    if Count > 0 then
        %       println(N=Count),
        %       Counts := Counts ++ [Count],
        %       Total := Total + Count
        %    end
        % end,

        % Even though Mode=all now is faster, the following is quite 
        % faster than linear_comb(...,all,...).
        Sol2 = _,
        % Note that since we use findall/2, we should NOT use first or all in
        % Mode since they have special constraints/behaviour.
        All = findall(Sol2, linear_comb(Diffs,N,xxx, AllowDuplicates, Sol2)),
        % println(allLen=All.len),
        % H = new_map(),
        % foreach(A in All)
        %   H.put(A,1)
        % end,
        % AllH = H.keys.sort,
        % println(allH=AllH),
        % Count2 = AllH.len,
        Count2 = All.len,
        if Count2 > 0 then
          println(N=Count2),
          Counts := Counts ++ [Count2],
          Total := Total + Count2
        end,
        flush(stdout)
     end,
     println(counts=[len=Len, counts=Counts]),
     TotalSums := TotalSums ++ [Total],
     nl,nl
   end,
   println(totalSums=TotalSums),
   nl.

%
% Not all differences are possible.
% For example the list Diffs = [1,3,6] is not possible.
% Assumptions about the origin list:
%
%  - It should start with 1 and end with Diffs[-1]+1 (7)
%
% We know that the list has at least [1,7]

% We have the following differences
%   7-1 = 6
% Now we are after the numbers that yield the difference 3.
% Then we must have a number that
%   X-Y = 3
% 1) Try 7-Y = 3
%    give Y=4 which is not in the difference list.
% 2) Try X-1 = 3
%    give X=4 which is not in the difference list.
%
% And any other X and Y for which X-Y = 3
% must give at least one number X or Y not in the difference list.
%
% QED: The difference list [1,3,6] is not possible.
%
% 
% So, what difference lists are possible? Here we look at very small
% 
%
% These are the only possible difference lists of length 2 (domain 1..10)
% [diffs = [1,2],sol = [1,2,3]]
% [diffs = [2,4],sol = [1,3,5]]
% [diffs = [3,6],sol = [1,4,7]]
% [diffs = [4,8],sol = [1,5,9]]
% [diffs = [5,10],sol = [1,6,11]]
%
% Length 3 domain 1..10
% [diffs = [1,2,3],sol = [1,3,4]]
% [diffs = [1,3,4],sol = [1,4,5]]
% [diffs = [1,4,5],sol = [1,5,6]]
% [diffs = [1,5,6],sol = [1,6,7]]
% [diffs = [1,6,7],sol = [1,7,8]]
% [diffs = [1,7,8],sol = [1,8,9]]
% [diffs = [1,8,9],sol = [1,9,10]]
% [diffs = [1,9,10],sol = [1,10,11]]
% [diffs = [2,3,5],sol = [1,4,6]]
% [diffs = [2,4,6],sol = [1,5,7]]
% [diffs = [2,5,7],sol = [1,6,8]]
% [diffs = [2,6,8],sol = [1,7,9]]
% [diffs = [2,7,9],sol = [1,8,10]]
% [diffs = [2,8,10],sol = [1,9,11]]
% [diffs = [3,4,7],sol = [1,5,8]]
% [diffs = [3,5,8],sol = [1,6,9]]
% [diffs = [3,6,9],sol = [1,7,10]]
% [diffs = [3,7,10],sol = [1,8,11]]
% [diffs = [4,5,9],sol = [1,6,10]]
% [diffs = [4,6,10],sol = [1,7,11]]

%
% Length 4 (domain 1..10):
% [diffs = [1,2,3,4],sol = [1,3,4,5]]
% [diffs = [1,2,3,5],sol = [1,3,4,6]]
% [diffs = [1,3,4,5],sol = [1,2,5,6]]
% [diffs = [1,3,4,7],sol = [1,4,5,8]]
% [diffs = [1,4,5,6],sol = [1,2,6,7]]
% [diffs = [1,4,5,9],sol = [1,5,6,10]]
% [diffs = [1,5,6,7],sol = [1,2,7,8]]
% [diffs = [1,6,7,8],sol = [1,2,8,9]]
% [diffs = [1,7,8,9],sol = [1,2,9,10]]
% [diffs = [1,8,9,10],sol = [1,2,10,11]]
% [diffs = [2,3,5,7],sol = [1,3,6,8]]
% [diffs = [2,3,5,8],sol = [1,4,6,9]]
% [diffs = [2,4,6,8],sol = [1,5,7,9]]
% [diffs = [2,4,6,10],sol = [1,5,7,11]]
% [diffs = [2,5,7,9],sol = [1,3,8,10]]
% [diffs = [2,6,8,10],sol = [1,3,9,11]]
% [diffs = [3,4,7,10],sol = [1,4,8,11]]
%
% Length 5 (domain 1..10)
% [diffs = [1,2,3,4,5],sol = [1,3,5,6]]
% [diffs = [1,2,3,4,6],sol = [1,4,5,7]]
% [diffs = [1,2,3,5,6],sol = [1,4,6,7]]
% [diffs = [1,2,4,5,6],sol = [1,5,6,7]]
% [diffs = [1,2,5,6,7],sol = [1,6,7,8]]
% [diffs = [1,2,6,7,8],sol = [1,7,8,9]]
% [diffs = [1,2,7,8,9],sol = [1,8,9,10]]
% [diffs = [1,2,8,9,10],sol = [1,9,10,11]]
% [diffs = [1,3,4,5,8],sol = [1,5,6,9]]
% [diffs = [1,3,4,6,7],sol = [1,4,7,8]]
% [diffs = [1,3,4,7,8],sol = [1,5,8,9]]
% [diffs = [1,4,5,6,10],sol = [1,6,7,11]]
% [diffs = [1,4,5,8,9],sol = [1,5,9,10]]
% [diffs = [1,4,5,9,10],sol = [1,6,10,11]]
% [diffs = [2,3,4,5,7],sol = [1,4,6,8]]
% [diffs = [2,3,5,6,8],sol = [1,4,7,9]]
% [diffs = [2,3,5,7,10],sol = [1,6,8,11]]
% [diffs = [2,3,5,8,10],sol = [1,6,9,11]]
% [diffs = [2,4,5,7,9],sol = [1,6,8,10]]
% [diffs = [2,4,6,8,10],sol = [1,5,9,11]]
% [diffs = [3,4,6,7,10],sol = [1,5,8,11]]
%
% Number of lists and possible difference lists for certain list lengts and max domain (1..Max)
%
% Length 3:
% Len Max #Lists #Diffs #Pct
% ------------------------------
% 3   3      1     1   1.0000
% 3   4      4     2   0.5000
% 3   5     10     4   0.4000
% 3   6     20     6   0.3000
% 3   7     35     9   0.2571
% 3   8     56    12   0.2143
% 3   9     84    16   0.1905
% 3  10    120    20   0.1667
% 3  11    165    25   0.1515
% 3  12    220    30   0.1364
% 3  13    286    36   0.1259
% 3  14    364    42   0.1154
% 3  15    455    49   0.1077
% 3  16    560    56   0.1000
% 3  17    680    64   0.0941
% 3  18    816    72   0.0882
% 3  19    969    81   0.0836
% 3  20   1140    90   0.0789
% 3  30   4060   210   0.0517
% 3  40   9880   380   0.0385
% 3  50  19600   600   0.0306
% 3 100 161700  2450   0.0152
%
% Length 4:
% Len  Max  #Lists #Diffs  #Pct
% ------------------------------
% 4   4      1     1   1.0000
% 4   5      5     3   0.6000
% 4   6     15     4   0.2667
% 4   7     35     7   0.2000
% 4   8     70    10   0.1429
% 4   9    126    13   0.1032
% 4  10    210    17   0.0810
% 4  11    330    22   0.0667
% 4  12    495    26   0.0525
% 4  13    715    32   0.0448
% 4  14   1001    38   0.0380
% 4  15   1365    44   0.0322
% 4  16   1820    51   0.0280
% 4  17   2380    59   0.0248
% 4  18   3060    66   0.0216
% 4  19   3876    75   0.0193
% 4  20   4845    84   0.0173
% 4  30  27405   200   0.0073
% 4  40  91390   367   0.0040
% 4  50 230300   584   0.0025
% 4 100 3921225  2417   0.0006

% 

go6 ?=>
   Mode = first,
   AllowDuplicates = false,
   Len = 3,
   MaxList = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,30,40,50,100],
   % MaxList = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20],   
   AllSols = [],
   AllNumLists = [],   
   foreach(Max in MaxList,Max >= Len)
     go6_model(Len,Max,All),
     
     NumLists := All.len,
     NumSols = 0,     
     foreach(L in All)
       Found = false, % for first solution       
       foreach(N in 2..Len+1, Found==false)
         if linear_comb(L,N,Mode,AllowDuplicates,_Sol) then
            % println([diffs=L,sol=Sol]),
            NumSols := NumSols + 1
            % , Found := true
          end
       end
     end,
     printf("%2d %3d %6d %5d   %0.4f\n", Len,Max,NumLists,NumSols,NumSols/NumLists),
     AllNumLists := AllNumLists ++ [NumLists],
     AllSols := AllSols ++ [NumSols]
   end,
   println(allSols=AllSols),
   println(allNumLists=AllNumLists),   
   nl.

go6 => true.

go6_model(Len,Max, All) =>
     L = new_list(Len),
     L :: 1..Max,
     all_different(L),
     increasing_strict(L),
     solve_all(L)=All.

%
% Experiment: number of solutions for "near extremal"/"holes".
% Notation:
% I will call "extremal" lists of 1..Len with m holes (missing values)
% as holes(Len,M). Example
%   A list of 1..10 with 2 holes is called hole(10,2)
% A difference list of hole(Len,M) must have Len as its last element,
% otherwise it reduces to hole(Len-1,M) or hole(Len-2,M) etc.
%
% Some expertimental observations:
%
% * For random difference lists (that give a solution) it seems that there are
%   only 2 solutions (where the differences of the solution are reversed)
% 
% * It seems that the most solutions are for the "extremal" difference lists
%   of 1..Len. See experiment in go5. These has also the most number of
%   different lengths of the solutions.
% 
%   For example the difference list 1..9 give these solutions [4,40,48,27,8,1]
%   for N=5..10
%
% * An "almost extremal" difference list is a list which is an extremal
%   list but with a hole, e.g. the following length 8 difference lists
%   with the number of solutions for different N, i.e. these are hole(9,1).
%   Note that some are not solvable.
%    - [2,3,4,5,6,7,8,9] % missing 1. Not solvable
%    - [1,3,4,5,6,7,8,9] % missing 2  Not solvable
%    - [1,2,4,5,6,7,8,9] % missing 3  [2] for N=5
%    - [1,2,3,5,6,7,8,9] % missing 4  [4,3] for N=5..6
%    - [1,2,3,4,6,7,8,9] % missing 5  [2] for N=5
%    - [1,2,3,4,5,7,8,9] % missing 6  [4,2] for N=5..6
%    - [1,2,3,4,5,6,8,9] % missing 7  [6,8,2] for N=5..7
%    - [1,2,3,4,5,6,7,9] % missing 8  [8,13,6,1] for N=5..8
%
% For the full 1..10 difference list these are the solutions: [38,90,75,35,9,1].
%
% As we see, as more holes the fewer possible solutions. Also, that there are
% fewer and fewer solvable difference lists.
%
% * hole(11,1)

% *  hole(10,1).
%   The same one-hole runs with a single holes in 1..10.
%   - [2,3,4,5,6,7,8,9,10] % missing 1 Not solvable
%   - [1,3,4,5,6,7,8,9,10] % missing 2 Not solvable
%   - [1,2,4,5,6,7,8,9,10] % missing 3 [2,2] for N=5..6
%   - [1,2,3,5,6,7,8,9,10] % missing 4 [2,6,2] for N=5..7
%   - [1,2,3,4,6,7,8,9,10] % missing 5 [2,7] for N=5..6
%   - [1,2,3,4,5,7,8,9,10] % missing 6 [2,6] for N=5..6
%   - [1,2,3,4,5,6,8,9,10] % missing 7 [7,3] for N=6..7
%   - [1,2,3,4,5,6,7,9,10] % missing 8 [4,16,10,2] for N=5..8
%   - [1,2,3,4,5,6,7,8,10] % missing 9 [2,19,19,7,1] for N=5..9
%
% Here are all possible solutions for hole(9,M)
% * hole(9,2). For two-holes differences for 1..9
% [1,2,3,4,5,6,9] = [4,1]
% [1,2,3,4,5,7,9] = [2,1]
% [1,2,3,4,6,7,9] = [2]
% [1,2,3,5,6,8,9] = [2]
% [1,2,3,6,7,8,9] = [2]
% [1,2,5,6,7,8,9] = [2,1]
% [1,3,4,5,6,8,9] = [2]
% [1,3,4,5,7,8,9] = [2,1]
% [2,3,4,5,6,7,9] = [2]
% count_diffs = 28
% count_sols = 13
% max_count = 2


% * hole(9,3): 3-holes for 1..9
% [1,2,3,6,7,9] = [2]
% [1,2,3,6,8,9] = [2]
% [1,2,6,7,8,9] = [2,2]
% [1,3,4,5,6,9] = [2]
% [1,3,4,5,8,9] = [2,2]
% [1,3,5,6,8,9] = [2]
% [2,3,4,5,6,9] = [2]
% [2,3,4,5,7,9] = [2,2]
% [2,3,4,6,7,9] = [2]
% count_diffs = 56
% count_sols = 12
% max_count = 2
%
% * hole(9,4): 4 holes for 1..9
% [1,2,7,8,9] = [2]
% [1,4,5,8,9] = [2]
% [2,4,5,7,9] = [2]
% count_diffs = 70
% count_sols = 3
% max_count = 1
%
% * hole(9,5)
% [1,4,5,9] = [1]
% [1,7,8,9] = [1]
% [2,5,7,9] = [1]
% count_diffs = 56
% count_sols = 3
% max_count = 1

% * hole(9,6)
% [1,8,9] = [2]
% [2,7,9] = [2]
% [3,6,9] = [2,1]
% [4,5,9] = [2]
% count_diffs = 28
% count_sols = 5
% max_count = 2

% * hole(9,7)
% (No solution)
% count_diffs = 8
% count_sols = 0
% max_count = 0

% * hole(9,8)
% [9] = [1]
% count_diffs = 1
% count_sols = 1
% max_count = 1


% One would except be some kind of formula when a solution has
% only solutions of a single length.
%
% One problem with this is that difference lists that are "regular", e.g.
%   1..12-lists with 6-holes
% Notation: I call this hole(12,6), i.e. 1..12 with 6 holes.
%
%   e.g. [2,4,6,8,10,12] has four different "solution points": [2,9,5,1]
%
% For hole(12,6) are 462 possible difference lists of hole(12,6)
% but only these 19 are solvable.
% [1,2,3,9,10,12] = [2]
% [1,2,3,9,11,12] = [2]
% [1,2,5,6,7,12] = [1]
% [1,2,9,10,11,12] = [2,2]
% [1,3,4,8,9,12] = [2]
% [1,3,4,8,11,12] = [2]
% [1,3,8,9,11,12] = [2]
% [1,4,5,7,8,12] = [2]
% [1,4,5,7,11,12] = [2]
% [1,4,7,8,11,12] = [2]
% [1,5,6,7,11,12] = [2,2]
% [1,5,6,10,11,12] = [1]
% [2,3,5,7,9,12] = [2]
% [2,3,5,7,10,12] = [2,2]
% [2,3,7,9,10,12] = [2]
% [2,4,6,8,10,12] = [2,9,5,1]
% [3,4,5,7,8,12] = [2]
% [3,4,5,7,9,12] = [2]
% [3,4,5,8,9,12] = [2]
%
% Also, note the "anomalies of these two which has only one solution. Why?
%   [1,5,6,10,11,12]  (solution: [1,2,7,12,13])
%   [1,2,5,6,7,12]    (solution: [1,6,7,8,13])
% Ah, because the list difference of the solution are reverse of itself (palindromes), interestingly the list differences of both solution are [5,1,1,5]. 
%
% Another way is to study the differences of the difference list.
%
% 
% Here are the solutions for hole(11,1)
% [1,2,3,4,5,6,7,8,9,11] = (diffs = [1,1,1,1,1,1,1,1,2]) = [2,28,42,26,8,1]
% [1,2,3,4,5,6,7,8,10,11] = (diffs = [1,1,1,1,1,1,1,2,1]) = [20,28,12,2]
% [1,2,3,4,5,6,7,9,10,11] = (diffs = [1,1,1,1,1,1,2,1,1]) = [10,12,3]
% [1,2,3,4,5,6,8,9,10,11] = (diffs = [1,1,1,1,1,2,1,1,1]) = [10,6]
% [1,2,3,4,5,7,8,9,10,11] = (diffs = [1,1,1,1,2,1,1,1,1]) = [2,4]
% [1,2,3,4,6,7,8,9,10,11] = (diffs = [1,1,1,2,1,1,1,1,1]) = [12,6]
% [1,2,3,5,6,7,8,9,10,11] = (diffs = [1,1,2,1,1,1,1,1,1]) = [6,6,1]
% [1,2,4,5,6,7,8,9,10,11] = (diffs = [1,2,1,1,1,1,1,1,1]) = [4]
% count_diffs = 10
% count_sols = 23
% max_count = 6
%
% hole(10,1)
% [1,2,3,4,5,6,7,8,10] = (diffs = [1,1,1,1,1,1,1,2]) = [2,19,19,7,1]
% [1,2,3,4,5,6,7,9,10] = (diffs = [1,1,1,1,1,1,2,1]) = [4,16,10,2]
% [1,2,3,4,5,6,8,9,10] = (diffs = [1,1,1,1,1,2,1,1]) = [7,3]
% [1,2,3,4,5,7,8,9,10] = (diffs = [1,1,1,1,2,1,1,1]) = [2,6]
% [1,2,3,4,6,7,8,9,10] = (diffs = [1,1,1,2,1,1,1,1]) = [2,7]
% [1,2,3,5,6,7,8,9,10] = (diffs = [1,1,2,1,1,1,1,1]) = [2,6,2]
% [1,2,4,5,6,7,8,9,10] = (diffs = [1,2,1,1,1,1,1,1]) = [2,2]
% count_diffs = 9
% count_sols = 20
% max_count = 5
%
% * hole(9.1):
% [1,2,3,4,5,6,7,9] = (diffs = [1,1,1,1,1,1,2]) = [8,13,6,1]
% [1,2,3,4,5,6,8,9] = (diffs = [1,1,1,1,1,2,1]) = [6,8,2]
% [1,2,3,4,5,7,8,9] = (diffs = [1,1,1,1,2,1,1]) = [4,2]
% [1,2,3,4,6,7,8,9] = (diffs = [1,1,1,2,1,1,1]) = [2]
% [1,2,3,5,6,7,8,9] = (diffs = [1,1,2,1,1,1,1]) = [4,3]
% [1,2,4,5,6,7,8,9] = (diffs = [1,2,1,1,1,1,1]) = [2]
% count_diffs = 8
% count_sols = 13
% max_count = 4
%
% hole(12,2)
% [1,2,3,4,5,6,7,8,9,12] = (diffs = [1,1,1,1,1,1,1,1,3]) = [12,16,7,1]
% [1,2,3,4,5,6,7,8,10,12] = (diffs = [1,1,1,1,1,1,1,2,2]) = [16,19,7,1]
% [1,2,3,4,5,6,7,8,11,12] = (diffs = [1,1,1,1,1,1,1,3,1]) = [12,10,2]
% [1,2,3,4,5,6,7,9,10,12] = (diffs = [1,1,1,1,1,1,2,1,2]) = [12,10,2]
% [1,2,3,4,5,6,7,9,11,12] = (diffs = [1,1,1,1,1,1,2,2,1]) = [2,6,2]
% [1,2,3,4,5,6,7,10,11,12] = (diffs = [1,1,1,1,1,1,3,1,1]) = [4,2]
% [1,2,3,4,5,6,8,9,10,12] = (diffs = [1,1,1,1,1,2,1,1,2]) = [2,4,2]
% [1,2,3,4,5,6,8,9,11,12] = (diffs = [1,1,1,1,1,2,1,2,1]) = [6,2]
% [1,2,3,4,5,6,8,10,11,12] = (diffs = [1,1,1,1,1,2,2,1,1]) = [2]
% [1,2,3,4,5,6,9,10,11,12] = (diffs = [1,1,1,1,1,3,1,1,1]) = [2]
% [1,2,3,4,5,7,8,9,10,12] = (diffs = [1,1,1,1,2,1,1,1,2]) = [8,6]
% [1,2,3,4,5,7,8,9,11,12] = (diffs = [1,1,1,1,2,1,1,2,1]) = [2,6,2]
% [1,2,3,4,5,7,8,10,11,12] = (diffs = [1,1,1,1,2,1,2,1,1]) = [2]
% [1,2,3,4,5,7,9,10,11,12] = (diffs = [1,1,1,1,2,2,1,1,1]) = [4]
% [1,2,3,4,6,7,8,9,10,12] = (diffs = [1,1,1,2,1,1,1,1,2]) = [2,6,3]
% [1,2,3,4,6,7,8,10,11,12] = (diffs = [1,1,1,2,1,1,2,1,1]) = [2]
% [1,2,3,4,6,8,9,10,11,12] = (diffs = [1,1,1,2,2,1,1,1,1]) = [2,2]
% [1,2,3,4,7,8,9,10,11,12] = (diffs = [1,1,1,3,1,1,1,1,1]) = [6,2]
% [1,2,3,5,6,7,8,9,11,12] = (diffs = [1,1,2,1,1,1,1,2,1]) = [2,2]
% [1,2,3,5,6,7,9,10,11,12] = (diffs = [1,1,2,1,1,2,1,1,1]) = [2]
% [1,2,3,5,6,8,9,10,11,12] = (diffs = [1,1,2,1,2,1,1,1,1]) = [2,1]
% [1,2,3,6,7,8,9,10,11,12] = (diffs = [1,1,3,1,1,1,1,1,1]) = [4,4,1]
% [1,2,4,5,6,7,8,10,11,12] = (diffs = [1,2,1,1,1,1,2,1,1]) = [4,2]
% [1,2,4,5,6,7,9,10,11,12] = (diffs = [1,2,1,1,1,2,1,1,1]) = [4,2]
% [1,2,4,5,6,8,9,10,11,12] = (diffs = [1,2,1,1,2,1,1,1,1]) = [2,1]
% [1,2,4,6,7,8,9,10,11,12] = (diffs = [1,2,2,1,1,1,1,1,1]) = [2]
% [1,3,4,5,6,7,8,9,11,12] = (diffs = [2,1,1,1,1,1,1,2,1]) = [2]
% [1,3,4,5,6,7,8,10,11,12] = (diffs = [2,1,1,1,1,1,2,1,1]) = [4,2]
% [2,3,4,5,6,7,8,9,10,12] = (diffs = [1,1,1,1,1,1,1,1,2]) = [2]
% count_diffs = 55
% count_sols = 60
% max_count = 4

%
% They are quite regular: All 1s except one 2. The position of the 2 don't
% seems to yield any discernable pattern. Except that:
%  - the longest range of solutions are when the 2 is last
%
% 
% Finding the regularity in this is a bit like prime numbers and
% composite numbers in this irregularity: The relation between primes and
% the composite numbers is extremly regular but it's harder to find it
% with simple formulas...
% 
% And then there's the problem why a certain difference
% list don't have any solution...
%
% Perhaps the problem is that most of these experiments are done on
% difference lists without duplicates...
%
go7 ?=>
  nolog,
  _ = random2(), % comment this for a static random seed.
  Map := get_global_map(),
  Map.put(max,0),
  Map.put(count_diffs,0),
  Map.put(count_sols,0),    

  % L = [I : I in 0..3],
  % L := generate_problem2(10,10),
  % println(l=L),
  % println(l_len=L.len),
  
  Mode = all,  
  AllowDuplicates = false,
  
  % if AllowDuplicates == experimental then
  %   Diffs = diffs2(L)
  % else
  %   Diffs = diffs(L)
  % end,

  % Diffs := [1,2,3,4,5,6,7,8,9],
  Max = 15,
  M = 12, % length of diffs list
  Diffs = new_list(M),
  Diffs :: 1..Max,
  increasing_strict(Diffs),
  % increasing(Diffs),  % with duplicates
  Diffs[M] #= Max,
  solve(Diffs),
  %println(diffs=Diffs),
  Len = Diffs.len,
  % println(diffs_len=Len),
  MinLen = min_len(Len),
  MaxLen=MinLen*6,
  % println([minLen=MinLen,maxLen=MaxLen]),
  Found = false,
  % println("\nSolutions:"),
  Map.put(count_diffs,Map.get(count_diffs)+1),
  Counts = [],
  foreach(N in MinLen..MaxLen, Found == false)
    % println(n=N),
    Count = 0,
    if linear_comb(Diffs,N,Mode,AllowDuplicates, Sol), Sol \= [] then
      if Mode == all then
        % println(n=N),
        % foreach(XX in Sol)
        %    print_solution(Diffs,XX,AllowDuplicates)
        % end,
        % println(num_solutions=Sol.len),
        Count := Sol.len
      else
        % print_solution(Diffs,Sol,AllowDuplicates),
        Count := Count + 1
      end
       % , Found := true    
    end,
    % println(N=Count),
    if Count > 0 then
      Map.put(count_sols, Map.get(count_sols)+1),
      Counts := Counts ++ [Count]
    end,
    if length(Counts) > Map.get(max) then
      Map.put(max,length(Counts))
    end
  end,
  Counts2 = [C : C in Counts, C > 0],
  if Counts2.len > 0 then
    % DiffsDiffs = [Diffs[I]-Diffs[I-1] : I in 2..Diffs.len],
    % println(Diffs=(diffs=DiffsDiffs)=Counts2=(sum=Counts2.sum)),
    println(Diffs=Counts2=Counts2.sum)
    % println(pair_difference=pair_differences(Diffs)),
  end,
  fail,
  nl.

go7 =>
  println(count_diffs=get_global_map().get(count_diffs)),
  println(count_sols=get_global_map().get(count_sols)),  
  println(max_count=get_global_map().get(max)).


%
% A variant of go7/0:
% Investigate the distribution of pair differences
% for the difference lists of a certain solution lengthN
% given possible values of 1..Len.
%
% It seems that this the distribution is some kind of
% ordering of the occurences of 1s, 2s, etc there
% are in the pair differences. The number with the
% hightest number of solution (the extremal 1..8)
% has the most number of 1s, not fewer 2s than any
% other difference list etc. However, it's not a complete
% lexicographic ordered table, since for the
% difference list
%   [1,2,3,4,5,6,7,9,11] with 20 solutions
% there are 6 1s, whereas the difference list
%   [1,2,3,5,6,7,8,9,10] with 10 solutions
% has 7 1s.
% Another difference list
%   [2,3,4,5,6,7,8,9,11]  with only 2 solutions
% has also 7 1s.
%   
% 
% For Max=12 M=9 (i.e. length 9 with possible 1..12)
% #Sols  1   2   3   4   5   6   7   8   9  10  11  12
% -----------------------------------------------------
% 128    8   7   6   5   4   3   2   1   0   0   0   0 [1,2,3,4,5,6,7,8,9] = [4,40,48,27,8,1] = (sum = 128)
%
%  48    7   7   6   5   4   3   2   1   1   0   0   0 [1,2,3,4,5,6,7,8,10] = [2,19,19,7,1] = (sum = 48)
%
%  32    7   6   6   5   4   3   2   2   1   0   0   0 [1,2,3,4,5,6,7,9,10] = [4,16,10,2] = (sum = 32)
%
%  23    7   6   6   5   4   3   2   1   1   1   0   0 [1,2,3,4,5,6,7,8,11] = [4,12,6,1] = (sum = 23)
%
%  20    6   7   5   5   4   3   2   2   1   1   0   0 [1,2,3,4,5,6,7,9,11] = [2,11,6,1] = (sum = 20)
% 
%  10    7   6   5   5   4   3   3   2   1   0   0   0 [1,2,3,5,6,7,8,9,10] = [2,6,2] = (sum = 10)
%  10    7   6   5   5   4   3   3   2   1   0   0   0 [1,2,3,4,5,6,8,9,10] = [7,3] = (sum = 10)
%  10    7   5   5   5   4   3   2   2   2   1   0   0 [1,2,3,4,5,6,7,10,11] = [2,6,2] = (sum = 10)
%  10    7   5   4   4   4   3   3   3   2   1   0   0 [1,2,3,6,7,8,9,10,11] = [2,6,2] = (sum = 10)
%  10    7   5   3   3   3   3   3   3   3   2   1   0 [1,2,3,7,8,9,10,11,12] = [2,6,2] = (sum = 10)
%  10    6   6   6   4   4   3   3   2   1   1   0   0 [1,2,3,4,5,6,8,9,11] = [2,6,2] = (sum = 10)
%  10    6   5   5   5   4   3   3   2   1   1   1   0 [1,2,3,4,5,7,8,9,12] = [2,6,2] = (sum = 10)
%  10    5   5   5   6   3   3   3   3   1   1   1   0 [1,3,4,5,7,8,9,11,12] = [2,6,2] = (sum = 10)
%
%   9    7   6   5   4   4   4   3   2   1   0   0   0 [1,2,3,4,6,7,8,9,10] = [2,7] = (sum = 9)
%   8    7   6   5   4   4   4   3   2   1   0   0   0 [1,2,3,4,5,7,8,9,10] = [2,6] = (sum = 8)
%   7    7   6   5   5   4   3   2   1   1   1   1   0 [1,2,3,4,5,6,7,8,12] = [2,4,1] = (sum = 7)
%   5    6   6   6   4   4   3   2   2   1   1   1   0 [1,2,3,4,5,6,7,9,12] = [4,1] = (sum = 5)
%   5    6   5   6   5   3   3   3   2   1   1   1   0 [1,2,3,4,5,6,8,9,12] = [4,1] = (sum = 5)
%
%   4    7   6   6   5   4   3   2   2   1   0   0   0 [1,2,4,5,6,7,8,9,10] = [2,2] = (sum = 4)
%   4    7   5   4   4   4   3   3   3   2   1   0   0 [1,2,3,4,5,6,9,10,11] = [2,2] = (sum = 4)
%   4    6   7   5   5   4   3   2   2   1   1   0   0 [2,3,4,5,6,7,8,10,12] = [2,2] = (sum = 4)
%   4    6   6   6   5   4   3   3   1   1   1   0   0 [1,3,4,5,6,7,8,10,11] = [2,2] = (sum = 4)
%   4    6   6   5   5   3   4   3   2   1   1   0   0 [1,2,3,4,5,7,8,9,11] = [2,2] = (sum = 4)
%   4    6   6   5   4   5   3   3   2   1   1   0   0 [2,3,4,5,7,8,9,10,12] = [2,2] = (sum = 4)
%   4    6   6   5   4   4   3   3   2   2   1   0   0 [1,2,3,4,5,6,8,10,11] = [2,2] = (sum = 4)
%   4    6   6   5   4   4   3   2   2   2   1   1   0 [1,2,3,4,5,6,7,10,12] = [2,2] = (sum = 4)
%   4    6   5   5   5   5   3   2   2   2   1   0   0 [1,2,4,5,6,7,9,10,11] = [2,2] = (sum = 4)
%   4    6   5   5   5   5   3   2   2   2   1   0   0 [1,2,3,5,6,7,8,10,11] = [2,2] = (sum = 4)
%   4    6   5   5   5   4   3   3   2   2   1   0   0 [1,2,4,5,6,8,9,10,11] = [2,2] = (sum = 4)
%   4    6   5   5   5   4   3   3   2   2   1   0   0 [1,2,3,4,6,7,8,10,11] = [2,2] = (sum = 4)
%   4    6   5   5   5   4   3   3   2   1   1   1   0 [1,3,4,5,6,7,8,11,12] = [2,2] = (sum = 4)
%   4    6   5   5   4   4   3   3   3   2   1   0   0 [1,2,3,5,6,8,9,10,11] = [2,2] = (sum = 4)
%   4    6   4   5   5   4   4   3   1   1   2   1   0 [1,2,4,5,6,7,8,11,12] = [2,2] = (sum = 4)
%   4    6   4   4   4   5   4   2   2   2   2   1   0 [1,2,4,5,6,7,10,11,12] = [3,1] = (sum = 4)
%   4    6   4   4   4   4   3   3   3   2   2   1   0 [1,2,4,5,6,9,10,11,12] = [2,2] = (sum = 4)
%   4    5   6   5   5   4   4   2   2   1   1   1   0 [1,2,4,5,6,7,8,10,12] = [2,2] = (sum = 4)
%   4    5   6   5   4   4   3   3   2   2   1   1   0 [1,2,3,4,5,7,8,10,12] = [2,2] = (sum = 4)
%   4    5   6   5   3   4   3   3   3   2   1   1   0 [1,2,3,4,5,7,9,10,12] = [2,2] = (sum = 4)
%
%   2    7   7   6   5   4   3   2   1   1   0   0   0 [2,3,4,5,6,7,8,9,11] = [2] = (sum = 2)
%   2    7   5   4   3   3   4   4   3   2   1   0   0 [1,2,3,4,7,8,9,10,11] = [2] = (sum = 2)
%   2    7   5   3   2   2   3   4   4   3   2   1   0 [1,2,3,4,8,9,10,11,12] = [2] = (sum = 2)
%   2    6   6   5   5   4   4   2   2   1   1   0   0 [1,3,4,5,6,7,9,10,11] = [2] = (sum = 2)
%   2    6   6   5   4   4   3   3   2   2   1   0   0 [1,2,4,6,7,8,9,10,11] = [2] = (sum = 2)
%   2    6   5   6   4   3   4   3   2   2   1   0   0 [1,2,3,4,5,7,8,10,11] = [2] = (sum = 2)
%   2    6   5   5   4   3   3   3   2   2   2   1   0 [1,2,3,4,5,6,8,11,12] = [2] = (sum = 2)
%   2    6   5   5   3   3   3   3   3   2   2   1   0 [1,2,4,7,8,9,10,11,12] = [2] = (sum = 2)
%   2    6   5   4   3   3   3   3   3   3   2   1   0 [1,2,3,4,5,7,10,11,12] = [2] = (sum = 2)
%   2    5   6   4   4   4   3   4   2   2   1   1   0 [1,2,3,5,7,8,9,10,12] = [2] = (sum = 2)
%   2    5   6   4   4   3   4   3   3   2   1   1   0 [1,2,3,4,6,8,9,10,12] = [2] = (sum = 2)
%   2    5   6   4   4   3   3   3   3   2   2   1   0 [1,2,3,4,5,7,9,11,12] = [2] = (sum = 2)
%   2    5   5   6   3   4   4   2   3   2   1   1   0 [1,2,3,4,6,7,9,10,12] = [2] = (sum = 2)
%   2    5   5   5   3   4   3   3   3   2   2   1   0 [1,2,4,5,7,9,10,11,12] = [2] = (sum = 2)
%   2    5   5   4   5   4   4   2   2   2   2   1   0 [1,2,4,6,7,8,10,11,12] = [2] = (sum = 2)
%   2    5   4   6   4   3   4   3   2   2   2   1   0 [1,2,3,5,6,8,9,11,12] = [2] = (sum = 2)
%
%   1    6   6   5   5   4   4   2   2   1   1   0   0 [1,2,3,5,6,7,8,9,11] = [1] = (sum = 1)
%   1    6   6   4   4   3   4   3   3   2   1   0   0 [1,2,3,5,7,8,9,10,11] = [1] = (sum = 1)
%   1    6   5   4   4   4   4   3   2   2   1   1   0 [1,2,3,6,7,8,9,10,12] = [1] = (sum = 1)
%   1    6   5   4   3   3   3   3   3   3   2   1   0 [1,2,3,6,8,9,10,11,12] = [1] = (sum = 1)
%
%
% For Max=11 M=8 (order by total number of solution, decreasing)
%
% #Sols  1   2   3   4   5   6   7   8   9  10  11
% ------------------------------------------------
 % 63    7   6   5   4   3   2   1   0   0   0   0 [1,2,3,4,5,6,7,8] = [8,27,20,7,1] = (sum = 63)
 % 28    6   6   5   4   3   2   1   1   0   0   0 [1,2,3,4,5,6,7,9] = [8,13,6,1] = (sum = 28)
 % 16    6   5   5   4   3   2   2   1   0   0   0 [1,2,3,4,5,6,8,9] = [6,8,2] = (sum = 16)
 %  9    6   5   5   4   3   2   1   1   1   0   0 [1,2,3,4,5,6,7,10] = [4,4,1] = (sum = 9)
 %  9    5   6   4   4   3   2   2   1   1   0   0 [1,2,3,4,5,6,8,10] = [4,4,1] = (sum = 9)
 %  7    6   5   4   4   3   3   2   1   0   0   0 [1,2,3,5,6,7,8,9] = [4,3] = (sum = 7)
 %  6    6   5   4   4   3   3   2   1   0   0   0 [1,2,3,4,5,7,8,9] = [4,2] = (sum = 6)
 %  6    6   4   3   3   3   3   3   2   1   0   0 [1,2,3,6,7,8,9,10] = [4,2] = (sum = 6)
 %  6    6   4   2   2   2   3   3   3   2   1   0 [1,2,3,7,8,9,10,11] = [4,2] = (sum = 6)
 %  6    5   6   4   4   3   2   2   1   1   0   0 [2,3,4,5,6,7,9,11] = [4,2] = (sum = 6)
 %  6    5   5   4   4   3   3   2   1   1   0   0 [1,2,3,4,6,7,8,10] = [4,2] = (sum = 6)
 %  6    5   4   5   3   3   3   2   2   1   0   0 [1,2,3,4,6,7,9,10] = [4,2] = (sum = 6)
 %  6    4   4   5   3   3   3   2   2   1   1   0 [1,2,3,5,6,8,9,11] = [4,2] = (sum = 6)
 %  4    6   4   4   4   3   2   2   2   1   0   0 [1,2,3,4,5,6,9,10] = [2,2] = (sum = 4)
 %  4    5   5   5   3   3   3   2   1   1   0   0 [1,2,3,4,5,7,8,10] = [2,2] = (sum = 4)
 %  4    5   5   4   4   4   2   2   1   1   0   0 [1,3,4,5,6,8,9,10] = [2,2] = (sum = 4)
 %  4    5   5   4   4   4   2   2   1   1   0   0 [1,2,3,5,6,7,8,10] = [2,2] = (sum = 4)
 %  4    5   3   4   4   4   3   1   1   2   1   0 [1,2,4,5,6,7,10,11] = [2,2] = (sum = 4)
 %  4    4   4   5   4   3   2   3   1   1   1   0 [1,3,4,6,7,8,10,11] = [2,2] = (sum = 4)
 %  3    6   5   4   4   3   2   1   1   1   1   0 [1,2,3,4,5,6,7,11] = [2,1] = (sum = 3)
 %  3    4   4   5   3   3   3   2   2   1   1   0 [1,3,4,6,7,9,10,11] = [2,1] = (sum = 3)
 %  2    6   6   5   4   3   2   1   1   0   0   0 [2,3,4,5,6,7,8,10] = [2] = (sum = 2)
 %  2    6   5   5   4   3   2   2   1   0   0   0 [1,2,4,5,6,7,8,9] = [2] = (sum = 2)
 %  2    6   5   5   4   3   2   1   1   1   0   0 [2,3,4,5,6,7,8,11] = [2] = (sum = 2)
 %  2    6   5   4   3   4   3   2   1   0   0   0 [1,2,3,4,6,7,8,9] = [2] = (sum = 2)
 %  2    6   4   3   3   3   3   3   2   1   0   0 [1,2,3,4,5,8,9,10] = [2] = (sum = 2)
 %  2    5   5   5   4   3   3   1   1   1   0   0 [1,3,4,5,6,7,9,10] = [2] = (sum = 2)
 %  2    5   5   5   4   3   2   2   1   0   1   0 [1,3,4,5,6,7,8,11] = [2] = (sum = 2)
 %  2    5   5   5   3   3   3   2   1   1   0   0 [2,3,4,5,6,8,9,11] = [2] = (sum = 2)
 %  2    5   5   4   4   4   2   2   1   1   0   0 [2,3,4,6,7,8,9,11] = [2] = (sum = 2)
 %  2    5   5   4   3   3   3   2   2   1   0   0 [1,2,3,4,5,7,9,10] = [2] = (sum = 2)
 %  2    5   5   4   3   3   2   2   2   1   1   0 [1,2,3,4,5,6,9,11] = [2] = (sum = 2)
 %  2    5   5   3   3   3   3   3   2   1   0   0 [1,2,3,5,7,8,9,10] = [2] = (sum = 2)
 %  2    5   5   3   3   3   3   3   2   1   0   0 [1,2,3,4,6,8,9,10] = [2] = (sum = 2)
 %  2    5   4   4   5   3   2   2   2   1   0   0 [1,2,4,5,6,8,9,10] = [2] = (sum = 2)
 %  2    5   4   4   4   3   3   2   1   1   1   0 [1,3,4,5,6,7,10,11] = [2] = (sum = 2)
 %  2    5   4   3   3   4   3   2   2   1   1   0 [1,2,3,6,7,8,9,11] = [2] = (sum = 2)
 %  2    5   4   3   3   3   3   3   2   1   1   0 [1,2,3,4,7,8,9,11] = [2] = (sum = 2)
 %  2    5   4   3   2   3   2   3   3   2   1   0 [1,2,3,6,8,9,10,11] = [2] = (sum = 2)
 %  2    5   4   3   2   2   3   3   3   2   1   0 [1,2,3,4,7,9,10,11] = [2] = (sum = 2)
 %  2    5   3   4   3   2   3   3   2   2   1   0 [1,2,3,4,7,8,10,11] = [2] = (sum = 2)
 %  2    5   3   3   4   4   2   2   2   2   1   0 [1,2,4,5,6,9,10,11] = [2] = (sum = 2)
 %  2    4   5   4   4   4   2   2   1   1   1   0 [1,2,4,5,6,7,9,11] = [2] = (sum = 2)
 %  2    4   4   4   3   4   2   2   2   2   1   0 [1,2,3,5,6,8,10,11] = [2] = (sum = 2)
 %  2    4   4   4   3   3   3   2   2   2   1   0 [1,2,3,5,7,8,10,11] = [2] = (sum = 2)

%
go8 =>
  Mode=all,
  AllowDuplicates=false,

  nolog,
  _ = random2(), % comment this for a static random seed.
  Map := get_global_map(),
  Map.put(max,0),
  Map.put(count_diffs,0),
  Map.put(count_sols,0),    

  % L = [I : I in 0..3],
  % L := generate_problem2(10,10),
  % println(l=L),
  % println(l_len=L.len),
  
  Mode = all,  
  AllowDuplicates = false,
  
  % if AllowDuplicates == experimental then
  %   Diffs = diffs2(L)
  % else
  %   Diffs = diffs(L)
  % end,

  % Diffs := [1,2,3,4,5,6,7,8,9],
  Max = 15,
  M = 9, % Max-M holes
  Diffs = new_list(M),
  Diffs :: 1..Max,
  increasing_strict(Diffs),
  % increasing(Diffs),  % with duplicates
  % Diffs[M] #= Max,
  solve(Diffs),
  %println(diffs=Diffs),
  Len = Diffs.len,
  % println(diffs_len=Len),
  MinLen = min_len(Len),
  MaxLen=MinLen*6,
  % println([minLen=MinLen,maxLen=MaxLen]),
  Found = false,
  % println("\nSolutions:"),
  Map.put(count_diffs,Map.get(count_diffs)+1),
  Counts = [],
  foreach(N in MinLen..MaxLen, Found == false)
    % println(n=N),
    Count = 0,
    if linear_comb(Diffs,N,Mode,AllowDuplicates, Sol), Sol \= [] then
      if Mode == all then
        % println(n=N),
        % foreach(XX in Sol)
        %    print_solution(Diffs,XX,AllowDuplicates)
        % end,
        % println(num_solutions=Sol.len),
        Count := Sol.len
      else
        % print_solution(Diffs,Sol,AllowDuplicates),
        Count := Count + 1
      end
       % , Found := true    
    end,
    % println(N=Count),
    if Count > 0 then
      Map.put(count_sols, Map.get(count_sols)+1),
      Counts := Counts ++ [Count]
    end,
    if length(Counts) > Map.get(max) then
      Map.put(max,length(Counts))
    end
  end,
  Counts2 = [C : C in Counts, C > 0],
  if Counts2.len > 0 then
    % DiffsDiffs = [Diffs[I]-Diffs[I-1] : I in 2..Diffs.len],
    % println(Diffs=Counts2=(sum=Counts2.sum)),
    % println(Diffs=Counts2=Counts2.sum),
    PD = pair_differences(Diffs),
    % println(diffs=Diffs),
    % println(pair_difference=PD),
    PairDiffMap = new_map(),    
    foreach(P in PD)
      PairDiffMap.put(P, PairDiffMap.get(P,0)+1)
    end,
    printf("%3d  ", Counts2.sum),
    foreach(I in 1..Max)
       printf("%3d ", PairDiffMap.get(I,0))
    end,
    println(Diffs=Counts2=(sum=Counts2.sum))    
  end,
  
  fail,

  nl.


go8 =>
  println(count_diffs=get_global_map().get(count_diffs)),
  println(count_sols=get_global_map().get(count_sols)),  
  println(max_count=get_global_map().get(max)).


go9 =>
  % Orig = [0,1,3,4,5,7,8,9],
  Orig = [1,3,4,5,7,8,9],  
  println(orig=Orig),
  LL = [
      [0,1,4,5,8,9], [0,1,4,8,9], [0,1,5,8,9]
  ],
  foreach(L in LL)
    println(l=L),
    L2 = [L[I]+1 : I in 1..L.len],
    println(l2=L2),
    println(pair_differences=pair_differences(L).sort_remove_dups),
    println(pair_differences2=pair_differences(L2).sort_remove_dups),
    nl
  end,
  nl.



%
% The (sorted) pair differences of a list L.
%
pair_differences(L) = [abs(L[I]-L[J]) : I in 1..L.len, J in I+1..L.len].sort.

%
% (Theoretical) Minumum length of a origin list given a difference list of length Len .
%
min_len(Len) = floor(1+sqrt(1+8*Len)/2).


%
% Shifted (normalized) list of list L.
% A shifted/normalized is a list which starts with 1 and ends with last(L)-L[1]+1
%
shifted_list(L) = [L[I] - L[1] + 1 : I in 1..L.len].


%
% Run the experiment
%
do_experiment(L, Diffs, Mode, AllowDuplicates) =>
  println(l=L),
  println(l_len=L.len),
  println(l_shifted=shifted_list(L)),
  println(l_diffs=[L[I] - L[I-1] : I in 2..L.len]),
  LDiffSum=sum([L[I] - L[I-1] : I in 2..L.len]),
  println(l_diff_sum=LDiffSum),

  println(diffs=Diffs),
  Len = Diffs.len,
  println(diffLen=Len), 
  DiffsDiffSum = sum([Diffs[I] - Diffs[I-1] : I in 2..Diffs.len]),
  println(diffs_diff_sum=DiffsDiffSum),
  
  println(difference_of_diff_sums=(LDiffSum-DiffsDiffSum)),
  println(differences_of_diffs_sums2=(LDiffSum-DiffsDiffSum-Diffs[1])),

  % Numbers that are missing in Diffs from the "extremal" 1..Diffs.len
  MissingFromDiffs = [I : I in 1..max(Diffs), not member(I,Diffs)],
  println(missing_from_diffs=MissingFromDiffs),
  println(num_missing_from_diffs=MissingFromDiffs.len),  

  % println(diffs_shifted=[Diffs[I] - Diffs[1]+1 : I in 1..Len]),  

  % Assumption: Length of X can be from MinLen to MinLen*2
  MinLen = floor(1+sqrt(1+8*Len)/2), % DiffLen = (N*(N-1)) div 2
  
  println(mode=Mode),
  Found = false,
  nl,
  flush(stdout),
  foreach(N in MinLen..MinLen*2, Found==false)
    println(n=N),
    if linear_comb(Diffs,N,Mode,AllowDuplicates, Sol), Sol \= [] then
      if Mode == all then
        foreach(XX in Sol)
          print_solution(Diffs,XX,AllowDuplicates)
        end,
        println(num_solutions=Sol.len)
      else
        print_solution(Diffs,Sol,AllowDuplicates)
      end,
      if Mode != all then
        Found := true
      end
    end
  end,
  
  nl.

%
% Print the solution.
%
print_solution(Diffs,X,AllowDuplicates) =>
        println(x=X),
        println(x0=[X[I]-1 : I in 1..X.len]),        
        N = X.len,
        % println(x_shifted=shifted_list(X)),
        println(x_diffs=[X[I] - X[I-1] : I in 2..N]),
        println(sum_x_diffs=sum([X[I] - X[I-1] : I in 2..N])),
        
        % Check that we got a correct solution

        %% Experimental
        if AllowDuplicates == experimental then
          Check = [abs(X[I]-X[J]) : I in 1..N, J in I+1..N].sort % allow duplicates
        else 
          Check = [abs(X[I]-X[J]) : I in 1..N, J in I+1..N].sort_remove_dups
        end,

        println(check=Check),
        if Diffs != Check then
          println("Diffs != Check!")
        end,
        nl.


%
% Recover linear combinations (pairs of differences)
%
% Diffs   : The list of differences
% N       : The length of X, the solution
% Mode    : - minimize: minimize the sum of the list differences in X
%           - all: find all solutions. A tip: for smaller problems is might be better to use the CP solver for this.
%           - first: find the first solution. Use SAT for large lists.
% AllowDuplicates: Allow duplicates in Diffs if true.
%           Note: Even if AllowDuplicates is true, then we handle it as false if Diffs[1] > 0.
% X       : the solution, the "recovered" list
%
linear_comb(Diffs,N,Mode, AllowDuplicates, Sol)=>

  % Create the list X of decision variables of size N
  X = new_list(N),
  
  %% The min value is 1 and the max value is 1 + max(Diffs)
  Max = max(Diffs)+1,
  X :: 1..Max, % create the domain of the decision variables

  % Assumptions: The unknown list must have distinct values
  % Note: SAT solver seems to go faster without this.
  if member(cp,loaded_modules()), (AllowDuplicates == false ; Diffs[1] > 0) then  
      all_distinct(X)
      % all_different(X)
  end,
  
  % The values are sorted (symmetry breaking)
  if AllowDuplicates != false, Diffs[1] == 0 then
    increasing(X)
  else
    increasing_strict(X)
  end,


  % Symmetry breaking: 1 is the first number.
  X[1] #= 1,

  % Since we start with 1, the last (and maximum) number in X is
  % max(Diffs) + 1 (which is the max domain value as well)
  X[N] #= max(Diffs)+1,

  % Another symmetry breaking: It seems that there are always two
  % solutions for each problem instance. There is a symmetry here
  % where these two solution has the same list differences
  % (X[2]-X[1],...X[N]-X[N-1]) but one is the reversal of the other.
  % Hence we can ensure that that the first difference must be
  % less that the last.
  %
  if Mode == first then
     X[2]-X[1] #>= X[N]-X[N-1]
  end,

  
  % From the observation that the sum of the differences of the solutions
  % is sum of the list of the differences + the first element in (the sorted list) Diffs
  % sum([X[I]-X[I-1] : I in 2..N]) #= sum([Diffs[I]-Diffs[I-1] : I in 2..Diffs.len]) + Diffs[1],

  % Experimental
  if AllowDuplicates == experimental then
     foreach(T in 0..Max)
        sum([abs(X[I]-X[J])#=T : I in 1..N, J in I+1..N]) #= count(T,Diffs)
     end
  end,

  %
  % For large value of N Mode = all give a lot of symmetric IJs and Ts which
  % yield a lot of duplicate solutions. Hence it's better to skip these constraints
  % and use the sum approach, which is generally slower but does not give duplicates.
  %
  IJs = [],
  Ts = [],
  if Mode == first then
  
    % Ensure that we cover all the differences (in Diffs),
    % i.e. find some I and J (I<J) where abs(X[I]-X[J]) == D  
    foreach(D in Diffs)
       % Create the indices I,J where I < J
       I :: 1..N,
       J :: 1..N,
       I #< J,
     
       element(I,X,Ix), % Ix = X[I]
       element(J,X,Jx), % Jx = X[J]
       % % Ensure that this difference is covered
       D #= abs(Ix-Jx),
       IJs := IJs ++ [I,J]
    end,

    % All the differences from the pairs in the list X must be in the Diff list
    foreach(I in 1..N, J in I+1..N)
      T #= abs(X[I]-X[J]),
      T :: Diffs,       
      Ts := Ts ++ [T]
    end

  else
  
     % Ensure that we cover all the differences (in Diffs),
     % i.e. find some I and J (I<J) where abs(X[I]-X[J]) == D
     foreach(D in Diffs)
        sum([D #= abs(X[I]-X[J]) : I in 1..N, J in I+1..N]) #>= 1
     end,

     % All the differences from the pairs in the list X must be in the Diff list
     foreach(I in 1..N, J in I+1..N)
       sum([D #= abs(X[I]-X[J]) : D in Diffs]) #>= 1
       % sum([Diffs[D] #= abs(X[I]-X[J]) : D in 1..Diffs.len]) #>= 1     
     end  
      
  end,

  % Solve
  % Vars = X ++ Ts ++ IJs,
  Vars = X ++ IJs ++ Ts,
  if Mode = minimize then
    SumDiffs #= sum([abs(X[I]-X[J]) : I in 1..N, J in I+1..N]),
    solve($[min(SumDiffs),report(printf("X: %w\n", X))],Vars),
    Sol = X
  elseif Mode == all then
      % This might takes a long time.
      % Note: The variables must be in a list since we
      %       only care about the first element (X).
      %       We got a (huge) number of the same X
      %       since Ts and IJs are different.
      solve_all($[ffd,updown],[X,Ts,IJs])=AllSols,
      % println(allSols=AllSols),
      Sols = new_map(),
      foreach(S in AllSols)
        XX = S[1], % take the X
        if not Sols.has_key(XX) then
          Sols.put(XX,1)
        end
      end,
      Sol = Sols.keys()
  else
     % if Mode == all then
     %   solve($[ff,split],Vars)
     %else 
     solve($[ffc,split],Vars),
     %end, 

     Sol = X
  end.


%
% Return the list difference of a list L
%
diffs(L) = Diffs =>
  Len = L.len,
  Diffs1 = [],
  foreach(I in 1..Len, J in I+1..Len)
    Diffs1 := Diffs1 ++ [abs(L[I]-L[J])] 
  end,
  Diffs = Diffs1.sort_remove_dups.


% Don't remove duplicates
% Experimental
diffs2(L) = Diffs =>
  Len = L.len,
  Diffs1 = [],
  foreach(I in 1..Len, J in I+1..Len)
    Diffs1 := Diffs1 ++ [abs(L[I]-L[J])] 
  end,
  Diffs = Diffs1.sort.


%
% Generate a random problem instance of (at least) length N
%
% Sort and remove all duplicates
generate_problem(N,MaxVal) = [1+random() mod MaxVal : _ in 1..N].sort_remove_dups.

% Sort and keep duplicates.
generate_problem2(N,MaxVal) = [1+random() mod MaxVal : _ in 1..N].sort.