scripts/gen_occ_lookup.pl

   1 #!/usr/bin/perl -w
   2
   3 #
   4 # Generate lookup table that, given a packed DNA byte (four bases) and
   5 # a character (A, C, G or T), returns how many times that character
   6 # occurs in that packed byte.  Useful for quickly counting character
   7 # occurrences in long strings.  The LUT is indexed first by character
   8 # (0-3) then by byte (0-255).
   9 #
  10 # Larger lookup tables are also possible, though they seem
  11 # counterproductive.  E.g., looking up eight bases at a time yields a
  12 # 256K LUT, which doesn't fit in L1.  A four-base LUT is 1KB, easily
  13 # fitting in L1.
  14 #
  15 # See ebwt.h.
  16 #
  17
  18 my @as4 = (), @as3 = (), @as2 = (), @as1 = ();
  19 my @cs4 = (), @cs3 = (), @cs2 = (), @cs1 = ();
  20 my @gs4 = (), @gs3 = (), @gs2 = (), @gs1 = ();
  21 my @ts4 = (), @ts3 = (), @ts2 = (), @ts1 = ();
  22
  23 # Compile character arrays
  24 my $i;
  25 for($i = 0; $i < 256; $i++) {
  26         my $b01 = ($i >> 0) & 3;
  27         my $b23 = ($i >> 2) & 3;
  28         my $b45 = ($i >> 4) & 3;
  29         my $b67 = ($i >> 6) & 3;
  30
  31         my $a4 = ($b01 == 0) + ($b23 == 0) + ($b45 == 0) + ($b67 == 0);
  32         my $c4 = ($b01 == 1) + ($b23 == 1) + ($b45 == 1) + ($b67 == 1);
  33         my $g4 = ($b01 == 2) + ($b23 == 2) + ($b45 == 2) + ($b67 == 2);
  34         my $t4 = ($b01 == 3) + ($b23 == 3) + ($b45 == 3) + ($b67 == 3);
  35
  36         push @as4, $a4;
  37         push @cs4, $c4;
  38         push @gs4, $g4;
  39         push @ts4, $t4;
  40
  41         my $a3 = ($b01 == 0) + ($b23 == 0) + ($b45 == 0);
  42         my $c3 = ($b01 == 1) + ($b23 == 1) + ($b45 == 1);
  43         my $g3 = ($b01 == 2) + ($b23 == 2) + ($b45 == 2);
  44         my $t3 = ($b01 == 3) + ($b23 == 3) + ($b45 == 3);
  45
  46         push @as3, $a3;
  47         push @cs3, $c3;
  48         push @gs3, $g3;
  49         push @ts3, $t3;
  50
  51         my $a2 = ($b01 == 0) + ($b23 == 0);
  52         my $c2 = ($b01 == 1) + ($b23 == 1);
  53         my $g2 = ($b01 == 2) + ($b23 == 2);
  54         my $t2 = ($b01 == 3) + ($b23 == 3);
  55
  56         push @as2, $a2;
  57         push @cs2, $c2;
  58         push @gs2, $g2;
  59         push @ts2, $t2;
  60
  61         my $a1 = ($b01 == 0) + 0;
  62         my $c1 = ($b01 == 1) + 0;
  63         my $g1 = ($b01 == 2) + 0;
  64         my $t1 = ($b01 == 3) + 0;
  65
  66         push @as1, $a1;
  67         push @cs1, $c1;
  68         push @gs1, $g1;
  69         push @ts1, $t1;
  70 }
  71
  72 my $entsPerLine = 16;
  73
  74 print "#include <stdint.h>\n\n";
  75 print "/* Generated by gen_lookup_tables.pl */\n\n";
  76
  77 # Count occurrences in all 4 bit pairs
  78
  79 print "uint8_t cCntLUT_4[4][4][256] = {\n";
  80 print "\t/* All 4 bit pairs */ {\n";
  81
  82 # Print As array
  83 print "\t\t/* As */ {\n";
  84 for($i = 0; $i < 256; $i++) {
  85         print "\t\t\t" if(($i % $entsPerLine) == 0);
  86         print "$as4[$i], ";
  87         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
  88 }
  89 print "\t\t},\n";
  90
  91 # Print Cs array
  92 print "\t\t/* Cs */ {\n";
  93 for($i = 0; $i < 256; $i++) {
  94         print "\t\t\t" if(($i % $entsPerLine) == 0);
  95         print "$cs4[$i], ";
  96         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
  97 }
  98 print "\t\t},\n";
  99
 100 # Print Gs array
 101 print "\t\t/* Gs */ {\n";
 102 for($i = 0; $i < 256; $i++) {
 103         print "\t\t\t" if(($i % $entsPerLine) == 0);
 104         print "$gs4[$i], ";
 105         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 106 }
 107 print "\t\t},\n";
 108
 109 # Print Ts array
 110 print "\t\t/* Ts */ {\n";
 111 for($i = 0; $i < 256; $i++) {
 112         print "\t\t\t" if(($i % $entsPerLine) == 0);
 113         print "$ts4[$i], ";
 114         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 115 }
 116 print "\t\t}\n\t},\n";
 117
 118 # Count occurrences in low 1 bit pair
 119
 120 print "\t/* Least significant 1 bit pair */ {\n";
 121
 122 # Print As array
 123 print "\t\t/* As */ {\n";
 124 for($i = 0; $i < 256; $i++) {
 125         print "\t\t\t" if(($i % $entsPerLine) == 0);
 126         print "$as1[$i], ";
 127         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 128 }
 129 print "\t\t},\n";
 130
 131 # Print Cs array
 132 print "\t\t/* Cs */ {\n";
 133 for($i = 0; $i < 256; $i++) {
 134         print "\t\t\t" if(($i % $entsPerLine) == 0);
 135         print "$cs1[$i], ";
 136         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 137 }
 138 print "\t\t},\n";
 139
 140 # Print Gs array
 141 print "\t\t/* Gs */ {\n";
 142 for($i = 0; $i < 256; $i++) {
 143         print "\t\t\t" if(($i % $entsPerLine) == 0);
 144         print "$gs1[$i], ";
 145         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 146 }
 147 print "\t\t},\n";
 148
 149 # Print Ts array
 150 print "\t\t/* Ts */ {\n";
 151 for($i = 0; $i < 256; $i++) {
 152         print "\t\t\t" if(($i % $entsPerLine) == 0);
 153         print "$ts1[$i], ";
 154         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 155 }
 156 print "\t\t}\n\t},\n";
 157
 158 # Count occurrences in low 2 bit pairs
 159
 160 print "\t/* Least significant 2 bit pairs */ {\n";
 161
 162 # Print As array
 163 print "\t\t/* As */ {\n";
 164 for($i = 0; $i < 256; $i++) {
 165         print "\t\t\t" if(($i % $entsPerLine) == 0);
 166         print "$as2[$i], ";
 167         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 168 }
 169 print "\t\t},\n";
 170
 171 # Print Cs array
 172 print "\t\t/* Cs */ {\n";
 173 for($i = 0; $i < 256; $i++) {
 174         print "\t\t\t" if(($i % $entsPerLine) == 0);
 175         print "$cs2[$i], ";
 176         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 177 }
 178 print "\t\t},\n";
 179
 180 # Print Gs array
 181 print "\t\t/* Gs */ {\n";
 182 for($i = 0; $i < 256; $i++) {
 183         print "\t\t\t" if(($i % $entsPerLine) == 0);
 184         print "$gs2[$i], ";
 185         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 186 }
 187 print "\t\t},\n";
 188
 189 # Print Ts array
 190 print "\t\t/* Ts */ {\n";
 191 for($i = 0; $i < 256; $i++) {
 192         print "\t\t\t" if(($i % $entsPerLine) == 0);
 193         print "$ts2[$i], ";
 194         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 195 }
 196 print "\t\t}\n\t},\n";
 197
 198 # Count occurrences in low 3 bit pairs
 199
 200 print "\t/* Least significant 3 bit pairs */ {\n";
 201
 202 # Print As array
 203 print "\t\t/* As */ {\n";
 204 for($i = 0; $i < 256; $i++) {
 205         print "\t\t\t" if(($i % $entsPerLine) == 0);
 206         print "$as3[$i], ";
 207         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 208 }
 209 print "\t\t},\n";
 210
 211 # Print Cs array
 212 print "\t\t/* Cs */ {\n";
 213 for($i = 0; $i < 256; $i++) {
 214         print "\t\t\t" if(($i % $entsPerLine) == 0);
 215         print "$cs3[$i], ";
 216         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 217 }
 218 print "\t\t},\n";
 219
 220 # Print Gs array
 221 print "\t\t/* Gs */ {\n";
 222 for($i = 0; $i < 256; $i++) {
 223         print "\t\t\t" if(($i % $entsPerLine) == 0);
 224         print "$gs3[$i], ";
 225         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 226 }
 227 print "\t\t},\n";
 228
 229 # Print Ts array
 230 print "\t\t/* Ts */ {\n";
 231 for($i = 0; $i < 256; $i++) {
 232         print "\t\t\t" if(($i % $entsPerLine) == 0);
 233         print "$ts3[$i], ";
 234         print "\n" if(($i % $entsPerLine) == ($entsPerLine-1));
 235 }
 236 print "\t\t}\n\t}\n";
 237
 238 print "};\n";