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


_cuda_source = r"""
#include <torch/extension.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAStream.h>
#include <c10/cuda/CUDAException.h>
#include <cuda.h>
#include <cuda_runtime.h>

__device__ __forceinline__ float max_propagate_nan(float a, float b) {
    return (isnan(a) || isnan(b)) ? CUDART_NAN_F : fmaxf(a, b);
}

// Specialization for ConvTranspose3d(3, 16, 3, stride=2,
// padding=1, output_padding=1), followed by MaxPool3d(2, 2).
// Each thread computes all eight convolution outputs in one pooling cell.
__global__ void conv_pool_kernel(
    const float* __restrict__ input,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ pooled,
    int D, int H, int W, int spatial, int total
) {
    __shared__ float weights[81];
    const int channel = blockIdx.y;

    if (threadIdx.x < 81) {
        const int c = threadIdx.x / 27;
        const int k = threadIdx.x % 27;
        weights[threadIdx.x] = weight[(c * 16 + channel) * 27 + k];
    }
    __syncthreads();

    const int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= total) return;

    const int batch = index / spatial;
    const int s = index - batch * spatial;
    const int z = s / (H * W);
    const int y = (s / W) % H;
    const int x = s % W;

    const bool has_x = x + 1 < W;
    const bool has_y = y + 1 < H;
    const bool has_z = z + 1 < D;
    const int plane = H * W;

    float a000 = bias[channel];
    float a001 = a000;
    float a010 = a000;
    float a011 = a000;
    float a100 = a000;
    float a101 = a000;
    float a110 = a000;
    float a111 = a000;

    #pragma unroll
    for (int c = 0; c < 3; ++c) {
        const float* p = input + (batch * 3 + c) * spatial + s;
        const float* w = weights + c * 27;

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

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

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

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

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

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

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

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

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

    float m0 = max_propagate_nan(a000, a001);
    float m1 = max_propagate_nan(a010, a011);
    float m2 = max_propagate_nan(a100, a101);
    float m3 = max_propagate_nan(a110, a111);
    m0 = max_propagate_nan(m0, m1);
    m2 = max_propagate_nan(m2, m3);

    pooled[(batch * 16 + channel) * spatial + s] =
        max_propagate_nan(m0, m2);
}

// Fuse channel softmax, subtraction, swish, and channel maximum.
__global__ void finish_kernel(
    const float* __restrict__ pooled,
    const float* __restrict__ subtract,
    float* __restrict__ output,
    int spatial,
    int total
) {
    const int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= total) return;

    const int batch = index / spatial;
    const int s = index - batch * spatial;
    const float* p = pooled + batch * 16 * spatial + s;

    float values[16];
    float maximum = -CUDART_INF_F;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        values[c] = p[c * spatial];
        maximum = max_propagate_nan(maximum, values[c]);
    }

    float denominator = 0.0f;
    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        values[c] = __expf(values[c] - maximum);
        denominator += values[c];
    }

    const float inverse = 1.0f / denominator;
    float result = -CUDART_INF_F;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        const float value = values[c] * inverse - subtract[c];
        const float sigmoid = 1.0f / (1.0f + __expf(-value));
        result = max_propagate_nan(result, value * sigmoid);
    }

    output[index] = result;
}

torch::Tensor fused_forward(
    torch::Tensor input,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor subtract
) {
    const c10::cuda::CUDAGuard device_guard(input.device());

    const int N = input.size(0);
    const int D = input.size(2);
    const int H = input.size(3);
    const int W = input.size(4);
    const int spatial = D * H * W;
    const int total = N * spatial;

    auto pooled = torch::empty({N, 16, D, H, W}, input.options());
    auto output = torch::empty({N, D, H, W}, input.options());
    if (total == 0) return output;

    constexpr int threads = 256;
    const int blocks = (total + threads - 1) / threads;
    cudaStream_t stream =
        c10::cuda::getCurrentCUDAStream(input.get_device()).stream();

    conv_pool_kernel<<<dim3(blocks, 16), threads, 0, stream>>>(
        input.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        pooled.data_ptr<float>(),
        D, H, W, spatial, total
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_kernel<<<blocks, threads, 0, stream>>>(
        pooled.data_ptr<float>(),
        subtract.data_ptr<float>(),
        output.data_ptr<float>(),
        spatial, total
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_extension = load_inline(
    name="fused_transpose3d_pool_swish_v1",
    cpp_sources="""
        torch::Tensor fused_forward(
            torch::Tensor input,
            torch::Tensor weight,
            torch::Tensor bias,
            torch::Tensor subtract);
    """,
    cuda_sources=_cuda_source,
    functions=["fused_forward"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


def _triple(value):
    if isinstance(value, (tuple, list)):
        return tuple(value)
    return (value, value, value)


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._specialized = (
            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):
        weight = self.conv_transpose.weight
        bias = self.conv_transpose.bias
        tensors = (x, weight, bias, self.subtract)

        if (
            self._specialized
            and x.ndim == 5
            and x.shape[1] == 3
            and all(size > 0 for size in x.shape[2:])
            and x.is_cuda
            and not torch.is_autocast_enabled()
            and all(
                t is not None
                and t.device == x.device
                and t.dtype == torch.float32
                and t.is_contiguous()
                for t in tensors
            )
            and x.numel() // 3 * 16 < 2**31
            and not (
                torch.is_grad_enabled()
                and any(t.requires_grad for t in tensors)
            )
        ):
            return _extension.fused_forward(x, weight, bias, self.subtract)

        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]
```