Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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");
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement