         setup(hmm_index);
 
 //        const uint32 read_index    = thrust::get<0>(hmm_index);
-//        fprintf(stderr, "read %d: hmm_glocal(l_ref=%d qstart=%d, l_query=%d)\n", read_index, referenceLength, queryStart, queryLen);
-//        fprintf(stderr, "read %d: ref = { ", read_index);
+//        printf("read %d: hmm_glocal(l_ref=%d qstart=%d, l_query=%d)\n", read_index, referenceLength, queryStart, queryLen);
+//        printf("read %d: ref = { ", read_index);
 //        for(int c = 0; c < referenceLength; c++)
 //        {
-//            fprintf(stderr, "%c ", from_nvbio::iupac16_to_char(referenceBases[c]));
+//            printf("%c ", from_nvbio::iupac16_to_char(referenceBases[c]));
 //        }
-//        fprintf(stderr, "\n");
+//        printf("\n");
 //
-//        fprintf(stderr, "read %d: que = { ", read_index);
+//        printf("read %d: que = { ", read_index);
 //        for(int c = 0; c < queryLen; c++)
 //        {
-//            fprintf(stderr, "%c ", from_nvbio::iupac16_to_char(queryBases[c]));
+//            printf("%c ", from_nvbio::iupac16_to_char(queryBases[c]));
 //        }
-//        fprintf(stderr, "\n");
-
-//        fprintf(stderr, "read %d: _iqual = { % 3d % 3d % 3d % 3d % 3d ... % 3d % 3d % 3d % 3d % 3d }\n", read_index,
+//        printf("\n");
+//
+//        printf("read %d: _iqual = { % 3d % 3d % 3d % 3d % 3d ... % 3d % 3d % 3d % 3d % 3d }\n", read_index,
 //                inputQualities[0], inputQualities[1], inputQualities[2], inputQualities[3], inputQualities[4],
 //                inputQualities[queryLen - 5], inputQualities[queryLen - 4], inputQualities[queryLen - 3], inputQualities[queryLen - 2], inputQualities[queryLen - 1]);
-//        fprintf(stderr, "read %d: c->bw = %d, bw = %d, l_ref = %d, l_query = %d\n", read_index, MIN_BAND_WIDTH, bandWidth, referenceLength, queryLen);
+//        printf("read %d: c->bw = %d, bw = %d, l_ref = %d, l_query = %d\n", read_index, MIN_BAND_WIDTH, bandWidth, referenceLength, queryLen);
 
         /*** forward ***/
         // f[0]


                sum += fi[u] + fi[u+1] + fi[u+2];

    //            fprintf(stderr, "(%d,%d;%d): %.4f,%.4f,%.4f\n", i, k, u, fi[u], fi[u+1], fi[u+2]);
    //            fprintf(stderr, " .. u = %d v11 = %d v01 = %d v10 = %d e = %f\n", u, v11, v01, v10, e);
            }

            // rescale

            {
                int u;
                double e = calcEpsilon(referenceBases[k-1], queryBases[queryStart], inputQualities[queryStart]);
//                printf("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,
//                       inputQualities[queryStart],
//                       queryStart,
//                       from_nvbio::iupac16_to_char(queryBases[queryStart]), e);

                u = set_u(bandWidth, 1, k);

