Advertisement
Guest User

Strange Exception

a guest
May 29th, 2012
683
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 7.68 KB | None | 0 0
  1. module main;
  2. import std.stdio, std.datetime, std.parallelism, std.algorithm, core.memory;
  3.  
  4. //Adjacency matrix version
  5. //What can we cut from the matrix? First column - I think, first row will
  6. //list all of the connections it makes - actually kept this to keep alignment
  7. //with the sequence for bit hacking. Also the last row can be omitted, every
  8. //connection to the last point will already have been listed. Finally every row can
  9. //start listing points one later as the previous connections will already have been listed
  10. //4.2 seconds SPEEDY!
  11. //3409 matrixes with an unconstrained fold space, yet the constrained gives 3367 folds and
  12. //the correct answer too  
  13. enum LEN = 15;
  14. enum OFFSET = 2;
  15. enum DIM = LEN / 2 + (LEN & 1? 1 : 0);
  16. enum SUBOPTIMAL = 2;
  17. enum TARGET = 263_916;
  18.  
  19. struct coord { int y, x; }
  20.  
  21.  
  22. pure int bitRev(int n)
  23. {   int rev = 0;
  24.     foreach(i;0..LEN)
  25.         if(n & 1 << i)
  26.             rev += 1 << LEN - 1 - i;
  27.     return rev;
  28. }
  29.  
  30.  
  31. pure int bitCount(int i) //Slower than lookup
  32. {   i = i - ((i >> 1) & 0x55555555);
  33.     i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
  34.     return (((i + (i >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
  35. }
  36.  
  37.  
  38. //Recursive shape search
  39. void fold(  int position,     // current term of H/P
  40.             ref int optimal,
  41.             ref coord[] path,
  42.             ref int[DIM][DIM] board,
  43.             ref bool[ushort[LEN - 1]] unique_matrix,
  44.             ref ushort[LEN - 1][] matrixes)   // coordinates of groups
  45. {   board[path[position].y][path[position].x] = position + 1; //Block the current square
  46.     if(position == 0) //Reached the end of the sequence
  47.     {   int hh_count = 0;
  48.         foreach(i, c;path)
  49.         {   if(c.x > 0 && board[c.y][c.x - 1]) //Left is full
  50.                 ++hh_count;
  51.             if(c.y > 0 && board[c.y - 1][c.x]) //Up is full
  52.                 ++hh_count;
  53.         }
  54.  
  55.         //Check if this is the highest number of bonds arrangement
  56.         if(hh_count > optimal)
  57.             optimal = hh_count;
  58.         if(hh_count >= optimal - SUBOPTIMAL) //Add suboptimal folds to the path list
  59.         {   ushort[LEN - 1] matrix;
  60.             //i iteration number, c coordinate. We don't need to check the last coord
  61.             //as everything linking to it will already have been found
  62.             foreach(i, c;path[0..$ - 1])
  63.             {   if(c.x > 0 && board[c.y][c.x - 1] - 1 > i) //Left, test bound and that coord is later in path
  64.                     matrix[i] ^= 1 << (board[c.y][c.x - 1] - 1);
  65.                 if(c.x < DIM - 1 && board[c.y][c.x + 1] - 1 > i) //Right
  66.                     matrix[i] ^= 1 << (board[c.y][c.x + 1] - 1);
  67.                 if(c.y > 0 && board[c.y - 1][c.x] - 1 > i) //Up
  68.                     matrix[i] ^= 1 << (board[c.y - 1][c.x] - 1);
  69.                 if(c.y < DIM - 1 && board[c.y + 1][c.x] - 1 > i) //Down
  70.                     matrix[i] ^= 1 << (board[c.y + 1][c.x] - 1);
  71.             }
  72.             if(unique_matrix.get(matrix, 0) == 0)
  73.             {   unique_matrix[matrix] = true;
  74.                 matrixes ~= matrix;
  75.             }
  76.         }
  77.     }
  78.     else
  79.     {   //Otherwise- try moving in all 4 possible directions
  80.         void advanceProtein(int my, int mx)
  81.         {   if(board[my][mx] == 0)
  82.             {   path[position] = coord(my, mx);
  83.                 fold(position, optimal, path, board, unique_matrix, matrixes);
  84.                 board[my][mx] = 0;
  85.             }
  86.         }
  87.  
  88.         const int x = path[position].x;
  89.         const int y = path[position].y;
  90.         --position;
  91.         if(x > 0) advanceProtein(y, x -  1);
  92.         if(x < board.length - 1) advanceProtein(y, x + 1);
  93.         if(y > 0) advanceProtein(y - 1, x);
  94.         if(y < board.length - 1) advanceProtein(y + 1, x);
  95.     }
  96. }
  97.  
  98. pure int[] solve(const bool[] redundancy, const ushort[LEN - 1][] matrixes, const int[2^^LEN] bits, const int[] numbers)
  99. {   int redun = true;
  100.     static immutable int[] opti_table = [0,0,1,2,4,5,7,8,10,12,13,15,17,18,20,22];  
  101.     int[] max; //Max H-H bonds seen for a given sequence
  102.     max.length = 2^^LEN;
  103.     foreach(sequence;numbers)
  104.     {   int[] filled;
  105.         foreach(i;0..LEN - 1) //-1 because the last adjacency row is not used
  106.             if((sequence & 1 << i))
  107.                 filled ~= i;
  108.  
  109.         foreach(im, matrix;matrixes)
  110.         {   if(redundancy[im])
  111.                 continue;
  112.             int hh_count = 0;
  113.             foreach(i;filled)
  114.                 //Counts the hydrogens present in the sequence which are neighbours to the row
  115.                 hh_count += bits[matrix[i] & sequence];
  116.  
  117.             //Check if this is the highest number of bonds arrangement
  118.             if(hh_count > max[sequence])
  119.             {   max[sequence] = hh_count;
  120.                 if(hh_count == opti_table[bits[sequence]]) //If this is the optimal number for n hydrogens break
  121.                     break;
  122.             }
  123.         }
  124.     }
  125.     foreach(sequence;numbers)
  126.         max[bitRev(sequence)] = max[sequence]; //Copy bit reversed optimals
  127.     int temp = max.reduce!("a + b");
  128.     delete max;
  129.     //GC gc;
  130.     //gc.collect;
  131.  
  132.     if(temp != TARGET) redun = false;
  133.     return [redun, temp];
  134. }
  135.  
  136.  
  137. void main()
  138. {   StopWatch sw;
  139.     sw.start;
  140.    
  141.     coord[] path;
  142.     path.length = LEN;
  143.     int[DIM][DIM] board;
  144.     bool[ushort[LEN - 1]] unique_matrix;
  145.     ushort[LEN - 1][] matrixes;
  146.     int optimal = 0;
  147.  
  148.     path[LEN - 1] = coord(OFFSET, OFFSET); //y, then x
  149.     fold(LEN - 1, optimal, path, board, unique_matrix, matrixes);
  150.  
  151.     sw.peek.msecs.writeln("msecs after tree search");
  152.  
  153.     //Strict inferiority filter, removes folds that are strictly worse than others
  154.     ushort[LEN - 1][] matrixes_filtered;
  155.     bool[] superior;
  156.     superior.length = matrixes.length;
  157.     foreach(inum, i;matrixes.parallel)
  158.     {   int j;
  159.         for(j = 0;j < matrixes.length;++j)
  160.         {   if(j == inum)
  161.                 continue;
  162.             int k;
  163.             for(k = 0;k < LEN - 1;++k)
  164.                 if(i[k] &(~matrixes[j][k]))
  165.                     break;
  166.             if(k == LEN - 1) //Got to the end with all matches, either a copy or inferior
  167.                 break;
  168.         }
  169.         if(j == matrixes.length) //Matrix had no matches nor strictly superior
  170.             superior[inum] = true;
  171.     }
  172.     foreach(i, b;superior)
  173.         if(b) matrixes_filtered ~= matrixes[i];
  174.  
  175.     matrixes = matrixes_filtered;
  176.    
  177.     //Also could keep path information and use that at the end? May be faster
  178.     sw.peek.msecs.writeln("msecs after adjacency filtering");
  179.  
  180.     int[] numbers; //Remove reverse bits
  181.     foreach(sequence;3..2^^LEN)
  182.         if(sequence <= bitRev(sequence))
  183.             numbers ~= sequence;
  184.  
  185.     int[2^^LEN] bits; //Set bits for bit counter
  186.     bool[] redun; //Test for useless folds
  187.     redun.length = matrixes.length;
  188.     foreach(i;0..2^^LEN)
  189.         bits[i] = bitCount(i);
  190.  
  191.     enum SKIP = 0;
  192.     foreach(i;0..SKIP)
  193.         redun[i] = true;
  194.  
  195.     /*
  196.     foreach(i, rr;redun)
  197.     {   int[] temp = solve(redun, matrixes, bits, numbers);
  198.         temp[1].writeln;
  199.         //writefln("%.20f", (cast(double) temp[1]) / 2^^LEN);
  200.         //redun[i] = cast(bool) temp[0];
  201.     }
  202.     */
  203.     solve(redun, matrixes, bits, numbers)[1].writeln;
  204.     solve(redun, matrixes, bits, numbers)[1].writeln;
  205.     solve(redun, matrixes, bits, numbers)[1].writeln;
  206.     solve(redun, matrixes, bits, numbers)[1].writeln;
  207.     solve(redun, matrixes, bits, numbers)[1].writeln;
  208.     solve(redun, matrixes, bits, numbers)[1].writeln;
  209.     solve(redun, matrixes, bits, numbers)[1].writeln;
  210.     solve(redun, matrixes, bits, numbers)[1].writeln;
  211.  
  212.     int used = 0;
  213.     foreach(red;redun)
  214.         if(red == 0) ++used;
  215.  
  216.     used.writeln;
  217.     sw.peek.msecs.writeln("msecs at the end");
  218. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement