The turn-4 attempt was truncated mid-source (missing the Python `ModelNew` entirely). This version restores the working fusion and makes the dominant eval kernel substantially faster: one block per (batch, qz) with a contiguous 3×3-plane shared tile, register-blocked accumulation over all 16 output channels, and the effective kernel pre-folded into `__constant__` memory as `float4`s so inner-loop coefficient loads are uniform vector broadcasts instead of shared-memory traffic.

```python
# Copyright (c) Meta Platforms, Inc. and affiliates.
# All rights reserved.
#
# This source code is licensed under the license found in the
# LICENSE file in the root directory of this source tree.

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.cpp_extension import load_inline

# ---------------------------------------------------------------------------
# Math behind the fusion (kernel_size=3, stride=2, padding=1 transposed conv,
# followed by two stride-2 2x2x2 average pools):
#
#   convT output p receives x[i] with weight W[k] iff p = 2*i + k - 1.
#   Two avg pools with odd sizes (63 -> 31 -> 15) only use convT outputs
#   p in [0,59] per dim, i.e. final[q] = (1/64) * sum_{d in [0,3]^3} y[4q+d].
#   Gathering the 64 sub-terms gives a *valid* stride-2 kernel-3 convolution
#   on x[0..30] with an "effective" kernel E:
#       final[oc][q] = sum_ic sum_{j in [0,2]^3} E[oc][ic][j] * x[ic][2q+j]
#   where E folds W over the 64 pooled positions (and the 1/64 factor).
#   BatchNorm is per-channel affine, so it commutes with the average pool and
#   collapses into a per-channel scale/shift applied at the end.
#
# Performance notes:
#   * E is pre-folded into __constant__ memory as float4 groups
#     (layout [ic][j][oc]) so the inner loop issues 1 shared load + 4
#     uniform LDC.128 + 16 FFMA per tap.
#   * Each block handles one output z-plane (qz) for one batch entry: a
#     3(ic) x 3(gz) x 31 x 31 input tile is staged in shared memory with
#     near-fully-coalesced loads, then each of the 256 threads produces one
#     (qy, qx) output for all 16 output channels with 16 register
#     accumulators (16-way ILP).
#   * Train mode: batch statistics of the *full* 63^3 convT output are
#     computed on the fly (y is never materialized) with per-thread float
#     accumulation reduced in double.
# ---------------------------------------------------------------------------

cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda_runtime.h>

#define SLICE_ELEMS 250047   // 63*63*63

// Effective kernel, transposed layout E2[(ic*27 + j)*16 + oc], as float4.
__constant__ float4 Ec4[3 * 27 * 4];   // 3*27*16 floats = 5184 B

// Build effective kernel from conv-transpose W[ic][oc][27].
// Valid (d,j)->k mapping per dim: j=0 -> k={1,2}, j=1 -> k={0,1,2}, j=2 -> k={0}.
// grid = (IC*27) blocks, block = 16 threads (one per oc).
__global__ void build_eff_kernel(const float* __restrict__ W,
                                 float* __restrict__ E2,
                                 int IC, int OC) {
    const int ic = blockIdx.x / 27;
    const int j  = blockIdx.x % 27;
    const int oc = threadIdx.x;
    if (oc >= OC) return;
    const int jz = j / 9, jy = (j / 3) % 3, jx = j % 3;

    int kzl[3], kyl[3], kxl[3], nz, ny, nx;
    if (jz == 0)      { kzl[0]=1; kzl[1]=2; nz=2; }
    else if (jz == 1) { kzl[0]=0; kzl[1]=1; kzl[2]=2; nz=3; }
    else              { kzl[0]=0; nz=1; }
    if (jy == 0)      { kyl[0]=1; kyl[1]=2; ny=2; }
    else if (jy == 1) { kyl[0]=0; kyl[1]=1; kyl[2]=2; ny=3; }
    else              { kyl[0]=0; ny=1; }
    if (jx == 0)      { kxl[0]=1; kxl[1]=2; nx=2; }
    else if (jx == 1) { kxl[0]=0; kxl[1]=1; kxl[2]=2; nx=3; }
    else              { kxl[0]=0; nx=1; }

    const float* wb = W + ((long)ic * OC + oc) * 27;
    float acc = 0.f;
    for (int a = 0; a < nz; ++a)
        for (int b = 0; b < ny; ++b)
            for (int c = 0; c < nx; ++c)
                acc += wb[kzl[a] * 9 + kyl[b] * 3 + kxl[c]];
    E2[(ic * 27 + j) * OC + oc] = acc * (1.0f / 64.0f);
}

// Fully fused: effective stride-2 valid conv (k=3) + BN affine + both pools.
// x: [N][3][32][32][32], out: [N][16][15][15][15].
// grid = (15, N), block = 256 threads; thread -> (qy, qx).
__global__ void __launch_bounds__(256) fused_conv_bn_pool_kernel(
        const float* __restrict__ x,
        const float* __restrict__ bn_w,
        const float* __restrict__ bn_b,
        const float* __restrict__ bn_mean,   // running_mean (eval) or batch mean (train)
        const float* __restrict__ bn_var,    // running_var  (eval) or batch biased var (train)
        const float* __restrict__ conv_bias,
        float* __restrict__ out,
        float eps) {
    __shared__ float xs[3 * 3 * 31 * 31];    // [ic][gz:3][gy:31][xx:31]

    const int n   = blockIdx.y;
    const int qz  = blockIdx.x;
    const int tid = threadIdx.x;
    const int qy  = tid >> 4;                // 0..15
    const int qx  = tid & 15;                // 0..15

    // stage tile: x[n][ic][2qz+dz][0..30][0..30]
    const long xb = (long)n * 98304 + (long)(2 * qz) * 1024;
    #pragma unroll
    for (int c = 0; c < 9; ++c) {
        const float* src = x + xb + (long)(c / 3) * 32768 + (long)(c % 3) * 1024;
        float* dst = xs + c * 961;
        #pragma unroll
        for (int i = tid; i < 961; i += 256) {
            const int gy = i / 31;
            const int xx = i - gy * 31;
            dst[i] = src[gy * 32 + xx];
        }
    }
    __syncthreads();

    if (qy >= 15 || qx >= 15) return;

    float acc[16];
    #pragma unroll
    for (int oc = 0; oc < 16; ++oc) acc[oc] = 0.f;

    #pragma unroll
    for (int ic = 0; ic < 3; ++ic) {
        #pragma unroll
        for (int j = 0; j < 27; ++j) {
            const int jz = j / 9, jy = (j / 3) % 3, jx = j % 3;
            const float v = xs[(((ic * 3 + jz) * 31) + (2 * qy + jy)) * 31 + 2 * qx + jx];
            const float4* e4 = Ec4 + (ic * 27 + j) * 4;
            #pragma unroll
            for (int g = 0; g < 4; ++g) {
                const float4 w4 = e4[g];
                acc[4*g+0] = fmaf(w4.x, v, acc[4*g+0]);
                acc[4*g+1] = fmaf(w4.y, v, acc[4*g+1]);
                acc[4*g+2] = fmaf(w4.z, v, acc[4*g+2]);
                acc[4*g+3] = fmaf(w4.w, v, acc[4*g+3]);
            }
        }
    }

    const long ob = ((long)n * 240 + qz) * 225 + qy * 15 + qx;
    #pragma unroll
    for (int oc = 0; oc < 16; ++oc) {
        const float scale = bn_w[oc] * rsqrtf(bn_var[oc] + eps);
        const float shift = bn_b[oc] + scale * (conv_bias[oc] - bn_mean[oc]);
        out[ob + (long)oc * 3375] = fmaf(acc[oc], scale, shift);
    }
}

// Train-mode batch statistics of the full convT output y (63^3 per slice),
// computed on the fly: y is never written to global memory.
//   y[n][oc][pz][py][px] = bias[oc]
//       + sum_ic sum_{valid k taps} x[n][ic][iz][iy][ix] * W[ic][oc][k]
// For each dim: p even -> single tap (k=1, i=p/2);
//               p odd  -> two taps (k=0, i=(p+1)/2) and (k=2, i=(p-1)/2).
// block = (64, 4) threads -> px = threadIdx.x, py = blockIdx.x*4 + threadIdx.y,
// loop over pz.  Per-thread partials in float, block reduction in double,
// atomically accumulated per output channel.
__global__ void convt_stats_kernel(const float* __restrict__ x,
                                   const float* __restrict__ W,
                                   const float* __restrict__ bias,
                                   double* __restrict__ sum,
                                   double* __restrict__ sumsq,
                                   int IC, int OC) {
    __shared__ float Ws[3 * 27];
    __shared__ double bs[256];
    __shared__ double bs2[256];

    const int oc  = blockIdx.z;
    const int n   = blockIdx.y;
    const int tid = threadIdx.x + threadIdx.y * blockDim.x;

    for (int i = tid; i < IC * 27; i += 256) {
        int ic = i / 27, k = i - ic * 27;
        Ws[i] = W[((long)ic * OC + oc) * 27 + k];
    }
    __syncthreads();

    const int px = threadIdx.x;
    const int py = blockIdx.x * 4 + threadIdx.y;

    float s = 0.f, s2 = 0.f;
    if (px < 63 && py < 63) {
        int ix[2], kx[2], iy[2], ky[2], nx, ny;
        if (px & 1) { ix[0]=(px+1)>>1; kx[0]=0; ix[1]=(px-1)>>1; kx[1]=2; nx=2; }
        else        { ix[0]=px>>1;     kx[0]=1;                       nx=1; }
        if (py & 1) { iy[0]=(py+1)>>1; ky[0]=0; iy[1]=(py-1)>>1; ky[1]=2; ny=2; }
        else        { iy[0]=py>>1;     ky[0]=1;                       ny=1; }

        const float b = bias[oc];
        const long xn = (long)n * IC * 32768;

        for (int pz = 0; pz < 63; ++pz) {
            int iz[2], kz[2], nz;
            if (pz & 1) { iz[0]=(pz+1)>>1; kz[0]=0; iz[1]=(pz-1)>>1; kz[1]=2; nz=2; }
            else        { iz[0]=pz>>1;     kz[0]=1;                       nz=1; }

            float acc = b;
            for (int ic = 0; ic < IC; ++ic) {
                const float* xp = x + xn + (long)ic * 32768;
                const float* wp = Ws + ic * 27;
                for (int a = 0; a < nz; ++a) {
                    const float* xz = xp + (long)iz[a] * 1024;
                    const int kza = kz[a] * 9;
                    for (int bb = 0; bb < ny; ++bb) {
                        const float* xy = xz + iy[bb] * 32;
                        const int kb = kza + ky[bb] * 3;
                        #pragma unroll
                        for (int c = 0; c < 2; ++c) {
                            if (c < nx) acc += xy[ix[c]] * wp[kb + kx[c]];
                        }
                    }
                }
            }
            s  += acc;
            s2 += acc * acc;
        }
    }

    bs[tid] = (double)s; bs2[tid] = (double)s2;
    __syncthreads();
    for (int off = 128; off > 0; off >>= 1) {
        if (tid < off) {
            bs[tid]  += bs[tid + off];
            bs2[tid] += bs2[tid + off];
        }
        __syncthreads();
    }
    if (tid == 0) {
        atomicAdd(&sum[oc], bs[0]);
        atomicAdd(&sumsq[oc], bs2[0]);
    }
}

// compute batch mean/var and update running stats like aten batch_norm
__global__ void bn_finalize_kernel(const double* __restrict__ sum,
                                   const double* __restrict__ sumsq,
                                   float* __restrict__ running_mean,
                                   float* __restrict__ running_var,
                                   float* __restrict__ out_mean,
                                   float* __restrict__ out_var,
                                   long count, double momentum, float eps, int C) {
    int c = blockIdx.x * blockDim.x + threadIdx.x;
    if (c >= C) return;
    double m  = sum[c] / (double)count;
    double sq = sumsq[c] / (double)count;
    double v  = sq - m * m;
    if (v < 0.0) v = 0.0;
    out_mean[c] = (float)m;
    out_var[c]  = (float)v;
    double unb = (count > 1) ? v * ((double)count / (double)(count - 1)) : 0.0;
    running_mean[c] = (float)((1.0 - momentum) * (double)running_mean[c] + momentum * m);
    running_var[c]  = (float)((1.0 - momentum) * (double)running_var[c]  + momentum * unb);
}

static inline void upload_eff(const torch::Tensor& W, cudaStream_t stream) {
    auto E2 = torch::empty({3 * 27 * 16}, W.options());
    build_eff_kernel<<<dim3(81), 16, 0, stream>>>(
        W.data_ptr<float>(), E2.data_ptr<float>(), 3, 16);
    cudaMemcpyToSymbolAsync(Ec4, E2.data_ptr<float>(),
                            3 * 27 * 16 * sizeof(float), 0,
                            cudaMemcpyDeviceToDevice, stream);
}

torch::Tensor fused_eval(torch::Tensor x, torch::Tensor W,
                         torch::Tensor bn_w, torch::Tensor bn_b,
                         torch::Tensor bn_mean, torch::Tensor bn_var,
                         torch::Tensor conv_bias,
                         int64_t N, double eps) {
    auto stream = at::cuda::getCurrentCUDAStream();
    upload_eff(W, stream.stream());
    auto out = torch::empty({N, 16, 15, 15, 15}, x.options());
    dim3 grid(15, (unsigned)N);
    fused_conv_bn_pool_kernel<<<grid, 256, 0, stream>>>(
        x.data_ptr<float>(),
        bn_w.data_ptr<float>(), bn_b.data_ptr<float>(),
        bn_mean.data_ptr<float>(), bn_var.data_ptr<float>(),
        conv_bias.data_ptr<float>(), out.data_ptr<float>(), (float)eps);
    return out;
}

// train-mode: on-the-fly convT batch stats (y never materialized),
// then fused effective-conv + BN + pools with the batch stats.
torch::Tensor fused_train(torch::Tensor x, torch::Tensor W,
                          torch::Tensor bn_w, torch::Tensor bn_b,
                          torch::Tensor running_mean, torch::Tensor running_var,
                          torch::Tensor conv_bias,
                          int64_t N, double eps, double momentum) {
    auto stream = at::cuda::getCurrentCUDAStream();
    upload_eff(W, stream.stream());

    auto opts_d = torch::TensorOptions().dtype(torch::kFloat64).device(x.device());
    auto sum    = torch::zeros({16}, opts_d);
    auto sumsq  = torch::zeros({16}, opts_d);

    dim3 sgrid(16, (unsigned)N, 16);
    dim3 sblock(64, 4);
    convt_stats_kernel<<<sgrid, sblock, 0, stream>>>(
        x.data_ptr<float>(), W.data_ptr<float>(), conv_bias.data_ptr<float>(),
        sum.data_ptr<double>(), sumsq.data_ptr<double>(), 3, 16);

    auto bmean = torch::empty({16}, x.options());
    auto bvar  = torch::empty({16}, x.options());
    long count = (long)N * SLICE_ELEMS;
    bn_finalize_kernel<<<1, 64, 0, stream>>>(
        sum.data_ptr<double>(), sumsq.data_ptr<double>(),
        running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
        bmean.data_ptr<float>(), bvar.data_ptr<float>(),
        count, momentum, (float)eps, 16);

    auto out = torch::empty({N, 16, 15, 15, 15}, x.options());
    dim3 grid(15, (unsigned)N);
    fused_conv_bn_pool_kernel<<<grid, 256, 0, stream>>>(
        x.data_ptr<float>(),
        bn_w.data_ptr<float>(), bn_b.data_ptr<float>(),
        bmean.data_ptr<float>(), bvar.data_ptr<float>(),
        conv_bias.data_ptr<float>(), out.data_ptr<float>(), (float)eps);
    return out;
}
"""

cpp_source = r"""
torch::Tensor fused_eval(torch::Tensor x, torch::Tensor W,
                         torch::Tensor bn_w, torch::Tensor bn_b,
                         torch::Tensor bn_mean, torch::Tensor bn_var,
                         torch::Tensor conv_bias,
                         int64_t N, double eps);
torch::Tensor fused_train(torch::Tensor x, torch::Tensor W,
                          torch::Tensor bn_w, torch::Tensor bn_b,
                          torch::Tensor running_mean, torch::Tensor running_var,
                          torch::Tensor conv_bias,
                          int64_t N, double eps, double momentum);
"""

_fused = load_inline(
    name="fused_convT_bn_avgpool_v3",
    cpp_sources=cpp_source,
    cuda_sources=cuda_source,
    functions=["fused_eval", "fused_train"],
    verbose=False,
    extra_cuda_cflags=["-O3", "--use_fast_math"],
)


class ModelNew(nn.Module):
    """
    Optimized version of Model.  The conv-transpose (k=3, s=2, p=1), batch-norm
    and both 2x2x2 average pools are fused into custom CUDA kernels by
    algebraically collapsing the whole chain into one valid stride-2 kernel-3
    convolution with an effective kernel (held in constant memory), plus a
    per-channel affine transform.  In training mode the batch statistics of
    the transposed-conv output are computed on the fly (the 63^3-per-channel
    intermediate is never materialized), which removes the dominant cuDNN
    dgrad + memset + stats passes entirely.
    """

    def __init__(self, in_channels, out_channels, kernel_size, stride, padding, bias_shape):
        super(ModelNew, self).__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size, stride=stride, padding=padding
        )
        self.batch_norm = nn.BatchNorm3d(out_channels)
        self._in_channels = in_channels
        self._out_channels = out_channels

    def _fallback(self, x):
        x = self.conv_transpose(x)
        x = self.batch_norm(x)
        x = F.avg_pool3d(x, kernel_size=2)
        x = F.avg_pool3d(x, kernel_size=2)
        return x

    def forward(self, x):
        ic, oc = self._in_channels, self._out_channels
        bn = self.batch_norm
        ct = self.conv_transpose

        fast_ok = (
            x.is_cuda and x.dtype == torch.float32 and x.dim() == 5
            and x.size(1) == ic and x.size(2) == 32 and x.size(3) == 32 and x.size(4) == 32
            and ic == 3 and oc == 16
            and ct.bias is not None and bn.affine and bn.track_running_stats
        )
        if not fast_ok:
            return self._fallback(x)

        x = x.contiguous()
        w = ct.weight.contiguous()
        conv_bias = ct.bias.contiguous()
        N = x.shape[0]

        if not self.training:
            return _fused.fused_eval(
                x, w,
                bn.weight.contiguous(), bn.bias.contiguous(),
                bn.running_mean.contiguous(), bn.running_var.contiguous(),
                conv_bias, N, float(bn.eps),
            )
        else:
            return _fused.fused_train(
                x, w,
                bn.weight.contiguous(), bn.bias.contiguous(),
                bn.running_mean, bn.running_var,
                conv_bias, N, float(bn.eps), float(bn.momentum),
            )


batch_size = 64
in_channels = 3
out_channels = 16
depth, height, width = 32, 32, 32
kernel_size = 3
stride = 2
padding = 1
bias_shape = (out_channels, 1, 1, 1)


def get_inputs():
    return [torch.rand(batch_size, in_channels, depth, height, width)]


def get_init_inputs():
    return [in_channels, out_channels, kernel_size, stride, padding, bias_shape]
```