```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);
"""

_terms = [
    [(13, 0)],
    [(14, 0), (12, 1)],
    [(16, 0), (10, 2)],
    [(17, 0), (15, 1), (11, 2), (9, 3)],
    [(22, 0), (4, 4)],
    [(23, 0), (21, 1), (5, 4), (3, 5)],
    [(25, 0), (19, 2), (7, 4), (1, 6)],
    [(26, 0), (24, 1), (20, 2), (18, 3),
     (8, 4), (6, 5), (2, 6), (0, 7)],
]

_weight_use = {
    coefficient: (accumulator, value)
    for accumulator, terms in enumerate(_terms)
    for coefficient, value in terms
}

# Pad each filter to seven aligned float4 loads. Descending coefficient
# order preserves the accumulation order within every convolution output.
_accumulate = ""
for _group in range(6, -1, -1):
    _accumulate += f"""
        {{
            const float4 coefficients =
                reinterpret_cast<const float4*>(q)[{_group}];
    """
    for _component in range(3, -1, -1):
        _coefficient = 4 * _group + _component
        if _coefficient in _weight_use:
            _accumulator, _value = _weight_use[_coefficient]
            _field = "xyzw"[_component]
            _accumulate += (
                f"\n            a[{_accumulator}] = "
                f"fmaf(inputs[ic][{_value}], coefficients.{_field}, "
                f"a[{_accumulator}]);"
            )
    _accumulate += "\n        }\n"


_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__(256, 2)
void fused_forward_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)
{
    __shared__ __align__(16) float weights[3 * 16 * 28];
    __shared__ float biases[16];
    __shared__ float offsets[16];

    for (int i = threadIdx.x; i < 3 * 16 * 28; i += 256) {
        const int filter = i / 28;
        const int coefficient = i % 28;
        weights[i] = coefficient < 27
            ? weight[filter * 27 + coefficient]
            : 0.0f;
    }

    if (threadIdx.x < 16) {
        biases[threadIdx.x] = bias[threadIdx.x];
        offsets[threadIdx.x] = subtract[threadIdx.x];
    }
    __syncthreads();

    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * 256 + threadIdx.x;
    if (index >= total) {
        return;
    }

    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 + 1 < W;
    const bool next_h = h + 1 < H;
    const bool next_d = d + 1 < D;

    float inputs[3][8];

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

        inputs[ic][0] = p[0];
        inputs[ic][1] = next_w ? p[1] : 0.0f;
        inputs[ic][2] = next_h ? p[W] : 0.0f;
        inputs[ic][3] = next_h && next_w ? p[W + 1] : 0.0f;
        inputs[ic][4] = next_d ? p[plane] : 0.0f;
        inputs[ic][5] = next_d && next_w ? p[plane + 1] : 0.0f;
        inputs[ic][6] = next_d && next_h ? p[plane + W] : 0.0f;
        inputs[ic][7] =
            next_d && next_h && next_w ? p[plane + W + 1] : 0.0f;
    }

    float pooled[16];
    float channel_max = -INFINITY;

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

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

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

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

        pooled[c] = m;
        channel_max = maximum(channel_max, m);
    }

    float denominator = 0.0f;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        pooled[c] = __expf(pooled[c] - channel_max);
        denominator += pooled[c];
    }

    const float inverse = __fdividef(1.0f, denominator);
    float z_min = INFINITY;
    float z_max = -INFINITY;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        const float z = __fsub_rn(pooled[c] * inverse, offsets[c]);
        z_min = minimum(z_min, z);
        z_max = maximum(z_max, z);
    }

    // Swish has a single minimum, so a maximum over channel values
    // is attained at either the smallest or the largest channel value.
    output[index] = maximum(swish(z_min), swish(z_max));
}

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,
        "Expected [N, 3, D, H, W] with positive spatial dimensions");
    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(
        x.size(2) <= 2147483647LL &&
        x.size(3) <= 2147483647LL &&
        x.size(4) <= 2147483647LL,
        "Spatial dimensions are too large");

    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 = x.size(0) * spatial;

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

    if (total > 0) {
        const int64_t blocks = (total + 255) / 256;
        TORCH_CHECK(blocks <= 2147483647LL, "Input is too large");

        fused_forward_kernel<<<
            static_cast<unsigned int>(blocks),
            256,
            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);

        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    return output;
}
"""

_fused_ops = load_inline(
    name="fused_transpose3d_pool_softmax_swish_vecweights_v6",
    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.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()
        ):
            return _fused_ops.fused_forward(
                x.contiguous(),
                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]
```