```python
import torch
import torch.nn as nn
from torch.nn.modules.utils import _triple
from torch.utils.cpp_extension import load_inline


_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor fused_forward(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor subtract);
"""

# Each thread computes two neighboring spatial outputs, sharing both
# convolution weights and overlapping input neighborhoods.
_terms = [
    [(13, 0)],
    [(14, 0), (12, 1)],
    [(16, 0), (10, 3)],
    [(17, 0), (15, 1), (11, 3), (9, 4)],
    [(22, 0), (4, 6)],
    [(23, 0), (21, 1), (5, 6), (3, 7)],
    [(25, 0), (19, 3), (7, 6), (1, 9)],
    [(26, 0), (24, 1), (20, 3), (18, 4),
     (8, 6), (6, 7), (2, 9), (0, 10)],
]

_accumulate = "\n".join(
    f"""
    {{
        const float coefficient = q[{k}];
        a[{j}].x = fmaf(inputs[ic][{v}], coefficient, a[{j}].x);
        a[{j}].y = fmaf(inputs[ic][{v + 1}], coefficient, a[{j}].y);
    }}
    """
    for j, terms in enumerate(_terms)
    for k, v in terms
)

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <stdint.h>
#include <math.h>

__device__ __forceinline__ float maximum(float a, float b) {
    return (isnan(b) || b > a) ? b : a;
}

__device__ __forceinline__ float minimum(float a, float b) {
    return (isnan(b) || b < a) ? b : a;
}

__device__ __forceinline__ float swish(float x) {
    return x * __fdividef(1.0f, 1.0f + __expf(-x));
}

__global__ __launch_bounds__(128, 3)
void fused_pair_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ subtract,
    float* __restrict__ output,
    int D,
    int H,
    int W,
    int64_t spatial,
    int64_t total_pairs)
{
    __shared__ float weights[3 * 16 * 27];
    __shared__ float biases[16];
    __shared__ float offsets[16];

    for (int i = threadIdx.x; i < 3 * 16 * 27; i += 128) {
        weights[i] = weight[i];
    }
    if (threadIdx.x < 16) {
        biases[threadIdx.x] = bias[threadIdx.x];
        offsets[threadIdx.x] = subtract[threadIdx.x];
    }
    __syncthreads();

    const int64_t pair_index =
        static_cast<int64_t>(blockIdx.x) * 128 + threadIdx.x;
    if (pair_index >= total_pairs) {
        return;
    }

    const int64_t index = pair_index * 2;
    const int64_t batch = index / spatial;
    const int64_t position = index - batch * spatial;
    const int64_t plane = static_cast<int64_t>(H) * W;

    const int w = position % W;
    const int h = (position / W) % H;
    const int d = position / plane;
    const bool next_w = w + 2 < W;

    float inputs[3][12];

    #pragma unroll
    for (int ic = 0; ic < 3; ++ic) {
        const float* base =
            x + (batch * 3 + ic) * spatial + position;

        #pragma unroll
        for (int dz = 0; dz < 2; ++dz) {
            #pragma unroll
            for (int dy = 0; dy < 2; ++dy) {
                constexpr int row_width = 3;
                const int row = dz * 2 + dy;
                const bool valid = d + dz < D && h + dy < H;

                float2 values = make_float2(0.0f, 0.0f);
                float third = 0.0f;
                if (valid) {
                    const float* p = base + dz * plane + dy * W;
                    values = *reinterpret_cast<const float2*>(p);
                    if (next_w) {
                        third = p[2];
                    }
                }

                inputs[ic][row * row_width] = values.x;
                inputs[ic][row * row_width + 1] = values.y;
                inputs[ic][row * row_width + 2] = third;
            }
        }
    }

    float pooled0[16];
    float pooled1[16];
    float max0 = -INFINITY;
    float max1 = -INFINITY;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        float2 a[8];

        #pragma unroll
        for (int j = 0; j < 8; ++j) {
            a[j] = make_float2(0.0f, 0.0f);
        }

        #pragma unroll
        for (int ic = 0; ic < 3; ++ic) {
            const float* q = weights + (ic * 16 + c) * 27;
""" + _accumulate + r"""
        }

        const float m0 = maximum(
            maximum(maximum(a[0].x, a[1].x),
                    maximum(a[2].x, a[3].x)),
            maximum(maximum(a[4].x, a[5].x),
                    maximum(a[6].x, a[7].x))
        ) + biases[c];

        const float m1 = maximum(
            maximum(maximum(a[0].y, a[1].y),
                    maximum(a[2].y, a[3].y)),
            maximum(maximum(a[4].y, a[5].y),
                    maximum(a[6].y, a[7].y))
        ) + biases[c];

        pooled0[c] = m0;
        pooled1[c] = m1;
        max0 = maximum(max0, m0);
        max1 = maximum(max1, m1);
    }

    float sum0 = 0.0f;
    float sum1 = 0.0f;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        pooled0[c] = __expf(pooled0[c] - max0);
        pooled1[c] = __expf(pooled1[c] - max1);
        sum0 += pooled0[c];
        sum1 += pooled1[c];
    }

    const float inverse0 = __fdividef(1.0f, sum0);
    const float inverse1 = __fdividef(1.0f, sum1);

    float lo0 = INFINITY;
    float lo1 = INFINITY;
    float hi0 = -INFINITY;
    float hi1 = -INFINITY;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        const float z0 = __fsub_rn(pooled0[c] * inverse0, offsets[c]);
        const float z1 = __fsub_rn(pooled1[c] * inverse1, offsets[c]);
        lo0 = minimum(lo0, z0);
        lo1 = minimum(lo1, z1);
        hi0 = maximum(hi0, z0);
        hi1 = maximum(hi1, z1);
    }

    // Swish has one minimum; its maximum occurs at an endpoint.
    const float2 result = make_float2(
        maximum(swish(lo0), swish(hi0)),
        maximum(swish(lo1), swish(hi1))
    );
    reinterpret_cast<float2*>(output)[pair_index] = result;
}

torch::Tensor fused_forward(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor subtract)
{
    TORCH_CHECK(
        x.is_cuda() && weight.is_cuda() &&
        bias.is_cuda() && subtract.is_cuda(),
        "Expected CUDA tensors");
    TORCH_CHECK(
        x.device() == weight.device() &&
        x.device() == bias.device() &&
        x.device() == subtract.device(),
        "All tensors must be on the same device");
    TORCH_CHECK(
        x.scalar_type() == at::kFloat &&
        weight.scalar_type() == at::kFloat &&
        bias.scalar_type() == at::kFloat &&
        subtract.scalar_type() == at::kFloat,
        "Expected float32 tensors");
    TORCH_CHECK(
        x.is_contiguous() && weight.is_contiguous() &&
        bias.is_contiguous() && subtract.is_contiguous(),
        "Expected contiguous tensors");
    TORCH_CHECK(
        x.dim() == 5 && x.size(1) == 3 &&
        x.size(2) > 0 && x.size(3) > 0 &&
        x.size(4) > 0 && x.size(4) % 2 == 0,
        "Expected [N, 3, D, H, W] with positive dimensions and even W");
    TORCH_CHECK(
        weight.dim() == 5 &&
        weight.size(0) == 3 && weight.size(1) == 16 &&
        weight.size(2) == 3 && weight.size(3) == 3 &&
        weight.size(4) == 3 &&
        bias.numel() == 16 && subtract.numel() == 16,
        "Unsupported parameter shapes");
    TORCH_CHECK(
        reinterpret_cast<uintptr_t>(x.data_ptr<float>()) % 8 == 0,
        "Input must be aligned to eight bytes");

    const c10::cuda::CUDAGuard guard(x.device());

    const int D = x.size(2);
    const int H = x.size(3);
    const int W = x.size(4);
    const int64_t spatial = static_cast<int64_t>(D) * H * W;
    const int64_t total_pairs = x.size(0) * spatial / 2;

    auto output = torch::empty({x.size(0), D, H, W}, x.options());

    if (total_pairs > 0) {
        const int64_t blocks = (total_pairs + 127) / 128;
        TORCH_CHECK(blocks <= 2147483647LL, "Input is too large");

        fused_pair_kernel<<<
            static_cast<unsigned int>(blocks),
            128,
            0,
            at::cuda::getCurrentCUDAStream()
        >>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            subtract.data_ptr<float>(),
            output.data_ptr<float>(),
            D, H, W, spatial, total_pairs);

        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    return output;
}
"""

_fused_ops = load_inline(
    name="fused_transpose3d_pool_softmax_swish_pairs_v5",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_forward"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self,
        in_channels,
        out_channels,
        kernel_size,
        stride,
        padding,
        output_padding,
        pool_kernel_size,
        pool_stride,
        pool_padding,
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
            output_padding=output_padding,
        )
        self.max_pool = nn.MaxPool3d(
            kernel_size=pool_kernel_size,
            stride=pool_stride,
            padding=pool_padding,
        )
        self.subtract = nn.Parameter(torch.randn(out_channels))

        effective_pool_stride = (
            pool_kernel_size if pool_stride is None else pool_stride
        )
        self._use_fused = (
            in_channels == 3
            and out_channels == 16
            and _triple(kernel_size) == (3, 3, 3)
            and _triple(stride) == (2, 2, 2)
            and _triple(padding) == (1, 1, 1)
            and _triple(output_padding) == (1, 1, 1)
            and _triple(pool_kernel_size) == (2, 2, 2)
            and _triple(effective_pool_stride) == (2, 2, 2)
            and _triple(pool_padding) == (0, 0, 0)
        )

    def forward(self, x):
        if (
            self._use_fused
            and x.is_cuda
            and x.ndim == 5
            and x.shape[-1] % 2 == 0
            and x.dtype == torch.float32
            and self.conv_transpose.weight.dtype == torch.float32
            and self.conv_transpose.bias.dtype == torch.float32
            and self.subtract.dtype == torch.float32
            and not torch.is_autocast_enabled()
        ):
            x = x.contiguous()
            if x.storage_offset() % 2:
                x = x.clone()
            return _fused_ops.fused_forward(
                x,
                self.conv_transpose.weight.contiguous(),
                self.conv_transpose.bias.contiguous(),
                self.subtract.contiguous(),
            )

        x = self.conv_transpose(x)
        x = self.max_pool(x)
        x = torch.softmax(x, dim=1)
        x = x - self.subtract.view(1, -1, 1, 1, 1)
        x = torch.sigmoid(x) * x
        return torch.max(x, dim=1)[0]
```