doubledigest.jl

# Double digest problem
# Paul Lange

# length 10000 - 24 solutions
# a = [5976, 1543, 1319, 1120, 42];
# b = [4513, 2823, 2057, 607];
# c = [4513, 1543, 1319, 1120, 607, 514, 342, 42];

# length 20000 - 72 solutions
# a = [8479, 4868, 3696, 2646, 169, 142];
# b = [11968, 5026, 1081, 1050, 691, 184];
# c = [8479, 4167, 2646, 1081, 881, 859, 701, 691, 184, 169, 142];

# length 30000 - 4 solutions
# a = [13719, 6942, 3597, 3576, 2166];
# b = [10682, 5351, 4095, 3116, 2617, 2060, 1235, 844];
# c = [9299, 4420, 4095, 2962, 2617, 2060, 1383, 1235, 931, 614, 230, 154];

# length 40000 - 2?
a = [9979, 9348, 8022, 4020, 2693, 1892, 1714, 1371, 510, 451]
b = [9492, 8453, 7749, 7365, 2292, 2180, 1023, 959, 278, 124, 85]
c = [7042, 5608, 5464, 4371, 3884, 3121, 1901, 1768, 1590, 959, 899, 707, 702, 510, 451, 412, 278, 124, 124, 85]

# length 50000 -
# a = [8162, 7643, 6335, 6119, 4202, 4130, 4050, 2752, 2652, 2156, 1127, 672]
# b = [11374, 9801, 6341, 5671, 4949, 3189, 2689, 1687, 1274, 1195, 1121, 709]
# c = [7643, 6335, 4016, 3650, 3593, 2748, 2689, 2652, 2241, 2156, 1687, 1655, 1274, 1195, 1127, 1121, 975, 709, 672, 552, 537, 511, 262]


beta = 1e-7;
# steps = int(4.5e6);
steps = int(1e4);
iters = 10;

#using Debug;

function CalculateCurrentC(confa, confb, data_a, data_b)
  current_a = data_a[confa];
  current_b = data_b[confb];

  idxa = 1;
  idxb = 1;

  sa = 0;
  sb = 0;
  pos = 0;
  cc = 1;

  limit = length(current_a) + length(current_b) - 1;
  res = zeros(limit);

  while (cc <= limit)
    cura = current_a[idxa];
    curb = current_b[idxb];

    da = sa + cura;
    db = sb + curb;

    if (da <= db)
      sa = sa + cura;
      res[cc] = sa-pos;
      cc += 1;
      pos = sa;
      idxa += 1;
    else
      sb = sb + curb;
      res[cc] = sb-pos;
      cc += 1;
      pos = sb;
      idxb += 1;
    end
  end

  return sort(res, rev=true);
end


function Energy(original, current)
  return sum((original - current).^2 ./ original);
end


function Neighbor(confa, confb)
  split = rand(1:2);

  data = split == 1 ? confa : confb;
  idx = swap = randperm(length(data));

  data = copy(data);
  tmp = data[swap[1]];
  data[swap[1]] = data[swap[2]];
  data[swap[2]] = tmp;

  if split==1
    data, confb
  else
    confa, data
  end
end


# @debug function Step()
function Step()
  confa = randperm(length(a));
  confb = randperm(length(b));

  hes = zeros(steps)
  c_new = CalculateCurrentC(confa, confb, a, b);

  h_old = Energy(c, CalculateCurrentC(confa, confb, a, b));

  for t in 1:steps
    n1, n2 = Neighbor(confa, confb);
    c_new = CalculateCurrentC(n1, n2, a, b);
    h_new = Energy(c, c_new);

    if (h_new == 0)
      println("found something")
      return n1, n2
    end
    #@bp

    delta_h = h_new - h_old;
    hes[t] = h_old;

    if delta_h < 0
      confa = n1;
      confb = n2;
      h_old = h_new;
      hes[t] = h_new;
    else
      prob = exp(-beta * t * delta_h);
      r = rand();
      if r < prob
        confa = n1;
        confb = n2

        hes[t] = h_new;
        h_old = h_new;
      end
    end
    if mod(t,100000)==0
      println(h_old)
    end
  end
  return 0, 0;
end


function Iter()
  solutions = zeros(iters, length(a) + length(b));
  s=1;
  for i = 1:1
    #println(i)
    r1, r2 = Step();
    if r1 != 0 && r2 != 0
      solutions[s,:] = vcat(r1, r2)';
      s+=1;
    end
  end

  return unique(solutions, 1);
end

#results = Iter();
# file = matopen("matfile.mat", "w")
# write(file, "varname", results)
# close(file)