```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_transpose_pool(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor subtract);
"""

_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 * (1.0f / (1.0f + __expf(-x)));
}

__global__ __launch_bounds__(256, 2)
void fused_transpose_pool_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__ float weights[3 * 16 * 27];
    __shared__ float biases[16];
    __shared__ float offsets[16];

    for (int i = threadIdx.x; i < 3 * 16 * 27; i += blockDim.x) {
        weights[i] = weight[i];
    }
    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) * blockDim.x + threadIdx.x;
    if (index >= total) return;

    const int64_t batch = index / spatial;
    const int64_t position = index - batch * spatial;
    const int w = position % W;
    const int h = (position / W) % H;
    const int d = position / (static_cast<int64_t>(H) * W);
    const int64_t plane = static_cast<int64_t>(H) * W;

    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 a000 = 0.0f;
        float a001 = 0.0f;
        float a010 = 0.0f;
        float a011 = 0.0f;
        float a100 = 0.0f;
        float a101 = 0.0f;
        float a110 = 0.0f;
        float a111 = 0.0f;

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

            const float v000 = inputs[ic][0];
            const float v001 = inputs[ic][1];
            const float v010 = inputs[ic][2];
            const float v011 = inputs[ic][3];
            const float v100 = inputs[ic][4];
            const float v101 = inputs[ic][5];
            const float v110 = inputs[ic][6];
            const float v111 = inputs[ic][7];

            a000 = fmaf(v000, q[13], a000);

            a001 = fmaf(v000, q[14], a001);
            a001 = fmaf(v001, q[12], a001);

            a010 = fmaf(v000, q[16], a010);
            a010 = fmaf(v010, q[10], a010);

            a011 = fmaf(v000, q[17], a011);
            a011 = fmaf(v001, q[15], a011);
            a011 = fmaf(v010, q[11], a011);
            a011 = fmaf(v011, q[9], a011);

            a100 = fmaf(v000, q[22], a100);
            a100 = fmaf(v100, q[4], a100);

            a101 = fmaf(v000, q[23], a101);
            a101 = fmaf(v001, q[21], a101);
            a101 = fmaf(v100, q[5], a101);
            a101 = fmaf(v101, q[3], a101);

            a110 = fmaf(v000, q[25], a110);
            a110 = fmaf(v010, q[19], a110);
            a110 = fmaf(v100, q[7], a110);
            a110 = fmaf(v110, q[1], a110);

            a111 = fmaf(v000, q[26], a111);
            a111 = fmaf(v001, q[24], a111);
            a111 = fmaf(v010, q[20], a111);
            a111 = fmaf(v011, q[18], a111);
            a111 = fmaf(v100, q[8], a111);
            a111 = fmaf(v101, q[6], a111);
            a111 = fmaf(v110, q[2], a111);
            a111 = fmaf(v111, q[0], a111);
        }

        const float m0 = maximum(a000, a001);
        const float m1 = maximum(a010, a011);
        const float m2 = maximum(a100, a101);
        const float m3 = maximum(a110, a111);
        const float m = maximum(maximum(m0, m1), maximum(m2, m3))
                      + 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 = 1.0f / denominator;
    float z_min = INFINITY;
    float z_max = -INFINITY;

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

    // Swish decreases to its unique minimum, then increases.
    // Its maximum over these channel values occurs at an endpoint.
    output[index] = maximum(swish(z_min), swish(z_max));
}

torch::Tensor fused_transpose_pool(
    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,
                "Expected input shape [N, 3, D, H, 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(x.size(2) > 0 && x.size(3) > 0 && x.size(4) > 0,
                "Spatial dimensions must be positive");

    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) {
        constexpr int threads = 256;
        const int blocks = static_cast<int>((total + threads - 1) / threads);

        fused_transpose_pool_kernel<<<
            blocks, threads, 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_max_v3",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_transpose_pool"],
    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.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_transpose_pool(
                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]
```