```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 channel_max(float x) {
    #pragma unroll
    for (int delta = 8; delta > 0; delta >>= 1) {
        x = maximum(x, __shfl_xor_sync(0xffffffffu, x, delta, 16));
    }
    return x;
}

__device__ __forceinline__ float channel_sum(float x) {
    #pragma unroll
    for (int delta = 8; delta > 0; delta >>= 1) {
        x += __shfl_xor_sync(0xffffffffu, x, delta, 16);
    }
    return x;
}

__global__ __launch_bounds__(256)
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)
{
    // Adjacent lanes compute adjacent channels of the same output.
    __shared__ float weights[3 * 27 * 16];
    __shared__ float biases[16];
    __shared__ float offsets[16];

    for (int i = threadIdx.x; i < 3 * 27 * 16; i += 256) {
        const int c = i & 15;
        const int k = (i >> 4) % 27;
        const int ic = i / (27 * 16);
        weights[i] = weight[(ic * 16 + c) * 27 + k];
    }
    if (threadIdx.x < 16) {
        biases[threadIdx.x] = bias[threadIdx.x];
        offsets[threadIdx.x] = subtract[threadIdx.x];
    }
    __syncthreads();

    const int c = threadIdx.x & 15;
    const int group = threadIdx.x >> 4;
    const float b = biases[c];
    const float offset = offsets[c];
    const int64_t plane = static_cast<int64_t>(H) * W;

    // Reuse the shared weights across eight positions per channel group.
    #pragma unroll 1
    for (int iteration = 0; iteration < 8; ++iteration) {
        const int64_t index =
            static_cast<int64_t>(blockIdx.x) * 128 +
            iteration * 16 + group;
        const bool valid = index < total;

        float pooled = 0.0f;
        if (valid) {
            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 / plane;

            const bool next_w = w + 1 < W;
            const bool next_h = h + 1 < H;
            const bool next_d = d + 1 < D;

            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* p =
                    x + (batch * 3 + ic) * spatial + position;
                const float* q = weights + ic * 27 * 16 + c;

                const float v000 = p[0];
                const float v001 = next_w ? p[1] : 0.0f;
                const float v010 = next_h ? p[W] : 0.0f;
                const float v011 =
                    next_h && next_w ? p[W + 1] : 0.0f;
                const float v100 = next_d ? p[plane] : 0.0f;
                const float v101 =
                    next_d && next_w ? p[plane + 1] : 0.0f;
                const float v110 =
                    next_d && next_h ? p[plane + W] : 0.0f;
                const float v111 =
                    next_d && next_h && next_w
                    ? p[plane + W + 1] : 0.0f;

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

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

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

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

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

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

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

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

            pooled = maximum(
                maximum(maximum(a000, a001), maximum(a010, a011)),
                maximum(maximum(a100, a101), maximum(a110, a111))
            ) + b;
        }

        // Every lane participates, including lanes for a partial final tile.
        const float m = channel_max(pooled);
        const float e = __expf(pooled - m);
        const float denominator = channel_sum(e);
        const float probability = e / denominator;
        const float z = __fsub_rn(probability, offset);
        const float activated = z * (1.0f / (1.0f + __expf(-z)));
        const float result = channel_max(activated);

        if (valid && c == 0) {
            output[index] = result;
        }
    }
}

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) {
        const int64_t blocks = (total + 127) / 128;
        TORCH_CHECK(blocks <= 2147483647LL, "Input is too large");

        fused_transpose_pool_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_max_v4",
    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]
```