             {
                 int u, v11, v01, v10;
                 double e = calcEpsilon(referenceBases[k-1], qyi, inputQualities[queryStart+i-1]);
-//                printf("referenceBases[%d-1] = %c inputQualities[%d+%d-1] = %d qyi = %c -> e = %.4f\n",
+//                printf("read %d: referenceBases[%d-1] = %c inputQualities[%d+%d-1] = %d qyi = %c -> e = %.4f\n",
+//                        read_index,
 //                        k,
 //                        from_nvbio::iupac16_to_char(referenceBases[k-1]),
 //                        queryStart,

                k = (int)(double(-4.343) * log(double(1.0) - double(max)) + double(.499)); // = 10*log10(1-max)
                outputQualities[queryStart+i-1] = (char)(k > 100? 99 : (k < MIN_BASE_QUAL ? MIN_BASE_QUAL : k));

    //            printf("outputQualities[%d]: max = %.16f l = %.4f dk = %.4f k = %d -> %d\n", i, max, l, dk, k, outputQualities[queryStart+i-1]);
            }

    //        printf("(%.4f,%.4f) (%d,%d,%d,%.4f)\n", pb, sum, (i-1), (max_k>>2), (max_k&3), max);

            {
                int u;
                double e = calcEpsilon(referenceBases[k-1], queryBases[queryStart], inputQualities[queryStart]);
//                printf("read %d: referenceBases[%d-1] = %c inputQualities[%d] = %d queryBases[%d] = %c -> e = %.4f\n",
//                        read_index,
//                        k,
//                        from_nvbio::iupac16_to_char(referenceBases[k-1]),
//                        queryStart,
