module main;
import std.stdio, std.datetime, std.parallelism, std.algorithm, core.memory;
//Adjacency matrix version
//What can we cut from the matrix? First column - I think, first row will
//list all of the connections it makes - actually kept this to keep alignment
//with the sequence for bit hacking. Also the last row can be omitted, every
//connection to the last point will already have been listed. Finally every row can
//start listing points one later as the previous connections will already have been listed
//4.2 seconds SPEEDY!
//3409 matrixes with an unconstrained fold space, yet the constrained gives 3367 folds and
//the correct answer too
enum LEN = 15;
enum OFFSET = 2;
enum DIM = LEN / 2 + (LEN & 1? 1 : 0);
enum SUBOPTIMAL = 2;
enum TARGET = 263_916;
struct coord { int y, x; }
pure int bitRev(int n)
{ int rev = 0;
foreach(i;0..LEN)
if(n & 1 << i)
rev += 1 << LEN - 1 - i;
return rev;
}
pure int bitCount(int i) //Slower than lookup
{ i = i - ((i >> 1) & 0x55555555);
i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
return (((i + (i >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
}
//Recursive shape search
void fold( int position, // current term of H/P
ref int optimal,
ref coord[] path,
ref int[DIM][DIM] board,
ref bool[ushort[LEN - 1]] unique_matrix,
ref ushort[LEN - 1][] matrixes) // coordinates of groups
{ board[path[position].y][path[position].x] = position + 1; //Block the current square
if(position == 0) //Reached the end of the sequence
{ int hh_count = 0;
foreach(i, c;path)
{ if(c.x > 0 && board[c.y][c.x - 1]) //Left is full
++hh_count;
if(c.y > 0 && board[c.y - 1][c.x]) //Up is full
++hh_count;
}
//Check if this is the highest number of bonds arrangement
if(hh_count > optimal)
optimal = hh_count;
if(hh_count >= optimal - SUBOPTIMAL) //Add suboptimal folds to the path list
{ ushort[LEN - 1] matrix;
//i iteration number, c coordinate. We don't need to check the last coord
//as everything linking to it will already have been found
foreach(i, c;path[0..$ - 1])
{ if(c.x > 0 && board[c.y][c.x - 1] - 1 > i) //Left, test bound and that coord is later in path
matrix[i] ^= 1 << (board[c.y][c.x - 1] - 1);
if(c.x < DIM - 1 && board[c.y][c.x + 1] - 1 > i) //Right
matrix[i] ^= 1 << (board[c.y][c.x + 1] - 1);
if(c.y > 0 && board[c.y - 1][c.x] - 1 > i) //Up
matrix[i] ^= 1 << (board[c.y - 1][c.x] - 1);
if(c.y < DIM - 1 && board[c.y + 1][c.x] - 1 > i) //Down
matrix[i] ^= 1 << (board[c.y + 1][c.x] - 1);
}
if(unique_matrix.get(matrix, 0) == 0)
{ unique_matrix[matrix] = true;
matrixes ~= matrix;
}
}
}
else
{ //Otherwise- try moving in all 4 possible directions
void advanceProtein(int my, int mx)
{ if(board[my][mx] == 0)
{ path[position] = coord(my, mx);
fold(position, optimal, path, board, unique_matrix, matrixes);
board[my][mx] = 0;
}
}
const int x = path[position].x;
const int y = path[position].y;
--position;
if(x > 0) advanceProtein(y, x - 1);
if(x < board.length - 1) advanceProtein(y, x + 1);
if(y > 0) advanceProtein(y - 1, x);
if(y < board.length - 1) advanceProtein(y + 1, x);
}
}
pure int[] solve(const bool[] redundancy, const ushort[LEN - 1][] matrixes, const int[2^^LEN] bits, const int[] numbers)
{ int redun = true;
static immutable int[] opti_table = [0,0,1,2,4,5,7,8,10,12,13,15,17,18,20,22];
int[] max; //Max H-H bonds seen for a given sequence
max.length = 2^^LEN;
foreach(sequence;numbers)
{ int[] filled;
foreach(i;0..LEN - 1) //-1 because the last adjacency row is not used
if((sequence & 1 << i))
filled ~= i;
foreach(im, matrix;matrixes)
{ if(redundancy[im])
continue;
int hh_count = 0;
foreach(i;filled)
//Counts the hydrogens present in the sequence which are neighbours to the row
hh_count += bits[matrix[i] & sequence];
//Check if this is the highest number of bonds arrangement
if(hh_count > max[sequence])
{ max[sequence] = hh_count;
if(hh_count == opti_table[bits[sequence]]) //If this is the optimal number for n hydrogens break
break;
}
}
}
foreach(sequence;numbers)
max[bitRev(sequence)] = max[sequence]; //Copy bit reversed optimals
int temp = max.reduce!("a + b");
delete max;
//GC gc;
//gc.collect;
if(temp != TARGET) redun = false;
return [redun, temp];
}
void main()
{ StopWatch sw;
sw.start;
coord[] path;
path.length = LEN;
int[DIM][DIM] board;
bool[ushort[LEN - 1]] unique_matrix;
ushort[LEN - 1][] matrixes;
int optimal = 0;
path[LEN - 1] = coord(OFFSET, OFFSET); //y, then x
fold(LEN - 1, optimal, path, board, unique_matrix, matrixes);
sw.peek.msecs.writeln("msecs after tree search");
//Strict inferiority filter, removes folds that are strictly worse than others
ushort[LEN - 1][] matrixes_filtered;
bool[] superior;
superior.length = matrixes.length;
foreach(inum, i;matrixes.parallel)
{ int j;
for(j = 0;j < matrixes.length;++j)
{ if(j == inum)
continue;
int k;
for(k = 0;k < LEN - 1;++k)
if(i[k] &(~matrixes[j][k]))
break;
if(k == LEN - 1) //Got to the end with all matches, either a copy or inferior
break;
}
if(j == matrixes.length) //Matrix had no matches nor strictly superior
superior[inum] = true;
}
foreach(i, b;superior)
if(b) matrixes_filtered ~= matrixes[i];
matrixes = matrixes_filtered;
//Also could keep path information and use that at the end? May be faster
sw.peek.msecs.writeln("msecs after adjacency filtering");
int[] numbers; //Remove reverse bits
foreach(sequence;3..2^^LEN)
if(sequence <= bitRev(sequence))
numbers ~= sequence;
int[2^^LEN] bits; //Set bits for bit counter
bool[] redun; //Test for useless folds
redun.length = matrixes.length;
foreach(i;0..2^^LEN)
bits[i] = bitCount(i);
enum SKIP = 0;
foreach(i;0..SKIP)
redun[i] = true;
/*
foreach(i, rr;redun)
{ int[] temp = solve(redun, matrixes, bits, numbers);
temp[1].writeln;
//writefln("%.20f", (cast(double) temp[1]) / 2^^LEN);
//redun[i] = cast(bool) temp[0];
}
*/
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
solve(redun, matrixes, bits, numbers)[1].writeln;
int used = 0;
foreach(red;redun)
if(red == 0) ++used;
used.writeln;
sw.peek.msecs.writeln("msecs at the end");
}