       tx = min(31, max(0, t28 - imin)); vtmp = __shfl(vv, tx); if (t28 >= imin && t28 < imax) totals[threadIdx.x][0] += vtmp;
       tx = min(31, max(0, t29 - imin)); vtmp = __shfl(vv, tx); if (t29 >= imin && t29 < imax) totals[threadIdx.x][0] += vtmp;
       tx = min(31, max(0, t30 - imin)); vtmp = __shfl(vv, tx); if (t30 >= imin && t30 < imax) totals[threadIdx.x][0] += vtmp;
-      tx = min(31, max(0, t31 - imin)); vtmp = __shfl(vv, tx); if (t31 >= imin && t31 < imax) totals[threadIdx.x][0] += vtmp;
-    }
-
+      tx = min(31, max(0, t31 - imin)); vtmp = __shfl(vv, tx); if (t31 >= imin && t31 < imax) totals[threadIdx.x][0] += vtmp; */
+
+      // Check each feature index in turn, and see if it matches an input feature. Skip t0 which is actually the thresholds.
+      int nrem = min(32, nrows - tid);
+      vtmp = __shfl(vv, t1); if (t1 < nrem) totals[threadIdx.x][0] += vtmp; t1 -= 32;
+      vtmp = __shfl(vv, t2); if (t2 < nrem) totals[threadIdx.x][0] += vtmp; t2 -= 32;
+      vtmp = __shfl(vv, t3); if (t3 < nrem) totals[threadIdx.x][0] += vtmp; t3 -= 32;
+      vtmp = __shfl(vv, t4); if (t4 < nrem) totals[threadIdx.x][0] += vtmp; t4 -= 32;
+      vtmp = __shfl(vv, t5); if (t5 < nrem) totals[threadIdx.x][0] += vtmp; t5 -= 32;
+      vtmp = __shfl(vv, t6); if (t6 < nrem) totals[threadIdx.x][0] += vtmp; t6 -= 32;
+      vtmp = __shfl(vv, t7); if (t7 < nrem) totals[threadIdx.x][0] += vtmp; t7 -= 32;
+      vtmp = __shfl(vv, t8); if (t8 < nrem) totals[threadIdx.x][0] += vtmp; t8 -= 32;
+      vtmp = __shfl(vv, t9); if (t9 < nrem) totals[threadIdx.x][0] += vtmp; t9 -= 32;
+      vtmp = __shfl(vv, t10); if (t10 < nrem) totals[threadIdx.x][0] += vtmp; t10 -= 32;
+      vtmp = __shfl(vv, t11); if (t11 < nrem) totals[threadIdx.x][0] += vtmp; t11 -= 32;
+      vtmp = __shfl(vv, t12); if (t12 < nrem) totals[threadIdx.x][0] += vtmp; t12 -= 32;
+      vtmp = __shfl(vv, t13); if (t13 < nrem) totals[threadIdx.x][0] += vtmp; t13 -= 32;
+      vtmp = __shfl(vv, t14); if (t14 < nrem) totals[threadIdx.x][0] += vtmp; t14 -= 32;
+      vtmp = __shfl(vv, t15); if (t15 < nrem) totals[threadIdx.x][0] += vtmp; t15 -= 32;
+      vtmp = __shfl(vv, t16); if (t16 < nrem) totals[threadIdx.x][0] += vtmp; t16 -= 32;
+      vtmp = __shfl(vv, t17); if (t17 < nrem) totals[threadIdx.x][0] += vtmp; t17 -= 32;
+      vtmp = __shfl(vv, t18); if (t18 < nrem) totals[threadIdx.x][0] += vtmp; t18 -= 32;
+      vtmp = __shfl(vv, t19); if (t19 < nrem) totals[threadIdx.x][0] += vtmp; t19 -= 32;
+      vtmp = __shfl(vv, t20); if (t20 < nrem) totals[threadIdx.x][0] += vtmp; t20 -= 32;
+      vtmp = __shfl(vv, t21); if (t21 < nrem) totals[threadIdx.x][0] += vtmp; t21 -= 32;
+      vtmp = __shfl(vv, t22); if (t22 < nrem) totals[threadIdx.x][0] += vtmp; t22 -= 32;
+      vtmp = __shfl(vv, t23); if (t23 < nrem) totals[threadIdx.x][0] += vtmp; t23 -= 32;
+      vtmp = __shfl(vv, t24); if (t24 < nrem) totals[threadIdx.x][0] += vtmp; t24 -= 32;
+      vtmp = __shfl(vv, t25); if (t25 < nrem) totals[threadIdx.x][0] += vtmp; t25 -= 32;
+      vtmp = __shfl(vv, t26); if (t26 < nrem) totals[threadIdx.x][0] += vtmp; t26 -= 32;
+      vtmp = __shfl(vv, t27); if (t27 < nrem) totals[threadIdx.x][0] += vtmp; t27 -= 32;
+      vtmp = __shfl(vv, t28); if (t28 < nrem) totals[threadIdx.x][0] += vtmp; t28 -= 32;
+      vtmp = __shfl(vv, t29); if (t29 < nrem) totals[threadIdx.x][0] += vtmp; t29 -= 32;
+      vtmp = __shfl(vv, t30); if (t30 < nrem) totals[threadIdx.x][0] += vtmp; t30 -= 32;
+      vtmp = __shfl(vv, t31); if (t31 < nrem) totals[threadIdx.x][0] += vtmp; t31 -= 32;
+
+
+    }
+
+    // Compare the total for tree = threadIdx.x with its threshold. Save right child index if its bigger, else left child. 
     newt = 2 * tp;
     if (totals[threadIdx.x][0] > fthresh) {
       newt++;

    t30 = totals[threadIdx.x][30];
    t31 = totals[threadIdx.x][31];

    fthresh = *((float *)&t0);

    // Clear totals for each tree
    totals[threadIdx.x][0] = 0;
  
    // Now read the column and update totals
    for (tid = 0; tid + threadIdx.x < nrows; tid += blockDim.x) {
      vv = feats[tid + threadIdx.x + bd * nrows];
      imin = tid;
      imax = tid + blockDim.x;
      tx = min(31, max(0, t1 - imin)); vtmp = __shfl(vv, tx); if (t1 >= imin && t1 < imax) totals[threadIdx.x][0] += vtmp;
      tx = min(31, max(0, t2 - imin)); vtmp = __shfl(vv, tx); if (t2 >= imin && t2 < imax) totals[threadIdx.x][0] += vtmp;
      tx = min(31, max(0, t3 - imin)); vtmp = __shfl(vv, tx); if (t3 >= imin && t3 < imax) totals[threadIdx.x][0] += vtmp;
      tx = min(31, max(0, t4 - imin)); vtmp = __shfl(vv, tx); if (t4 >= imin && t4 < imax) totals[threadIdx.x][0] += vtmp;

      tx = min(31, max(0, t13 - imin)); vtmp = __shfl(vv, tx); if (t13 >= imin && t13 < imax) totals[threadIdx.x][0] += vtmp;
      tx = min(31, max(0, t14 - imin)); vtmp = __shfl(vv, tx); if (t14 >= imin && t14 < imax) totals[threadIdx.x][0] += vtmp;
      tx = min(31, max(0, t15 - imin)); vtmp = __shfl(vv, tx); if (t15 >= imin && t15 < imax) totals[threadIdx.x][0] += vtmp;
    }

    newt = 2 * tp;
