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


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

torch::Tensor fused_channel_reduce_cuda(
    torch::Tensor x, torch::Tensor subtract, bool pool);
"""

_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>

template <bool Pool>
__global__ void fused_channel_reduce_kernel(
    const float* __restrict__ input,
    const float* __restrict__ subtract,
    float* __restrict__ output,
    int64_t count,
    int depth,
    int height,
    int width,
    int out_depth,
    int out_height,
    int out_width)
{
    const int64_t idx =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (idx >= count) return;

    const int64_t out_spatial =
        static_cast<int64_t>(out_depth) * out_height * out_width;
    const int64_t in_spatial =
        static_cast<int64_t>(depth) * height * width;

    const int64_t batch = idx / out_spatial;
    const int64_t spatial = idx - batch * out_spatial;

    int64_t input_offset;
    if (Pool) {
        const int w = spatial % out_width;
        const int h = (spatial / out_width) % out_height;
        const int d = spatial / (static_cast<int64_t>(out_width) * out_height);
        input_offset =
            (static_cast<int64_t>(2 * d) * height + 2 * h) * width + 2 * w;
    } else {
        input_offset = spatial;
    }

    const float* base = input + batch * 16 * in_spatial + input_offset;

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

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        const float* p = base + c * in_spatial;
        float value;

        if (Pool) {
            value = -CUDART_INF_F;
            #pragma unroll
            for (int dz = 0; dz < 2; ++dz) {
                #pragma unroll
                for (int dy = 0; dy < 2; ++dy) {
                    const int64_t offset =
                        static_cast<int64_t>(dz) * height * width + dy * width;
                    const float a = p[offset];
                    const float b = p[offset + 1];
                    value = (isnan(a) || a > value) ? a : value;
                    value = (isnan(b) || b > value) ? b : value;
                }
            }
        } else {
            value = p[0];
        }

        values[c] = value;
        maximum = fmaxf(maximum, value);
    }

    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_denominator = 1.0f / denominator;
    float result = -CUDART_INF_F;

    #pragma unroll
    for (int c = 0; c < 16; ++c) {
        const float z = values[c] * inverse_denominator - subtract[c];
        const float activated = z * (1.0f / (1.0f + expf(-z)));
        result = (isnan(activated) || activated > result) ? activated : result;
    }

    output[idx] = result;
}

torch::Tensor fused_channel_reduce_cuda(
    torch::Tensor x, torch::Tensor subtract, bool pool)
{
    TORCH_CHECK(x.is_cuda() && subtract.is_cuda(),
                "Expected CUDA tensors");
    TORCH_CHECK(x.device() == subtract.device(),
                "Tensors must be on the same device");
    TORCH_CHECK(x.scalar_type() == at::kFloat &&
                subtract.scalar_type() == at::kFloat,
                "Expected float32 tensors");
    TORCH_CHECK(x.dim() == 5 && x.size(1) == 16,
                "Expected NCDHW input with 16 channels");
    TORCH_CHECK(subtract.numel() == 16,
                "Expected 16 subtraction parameters");

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

    const int depth = x.size(2);
    const int height = x.size(3);
    const int width = x.size(4);
    const int od = pool ? depth / 2 : depth;
    const int oh = pool ? height / 2 : height;
    const int ow = pool ? width / 2 : width;

    TORCH_CHECK(od > 0 && oh > 0 && ow > 0,
                "Spatial dimensions are too small");

    auto output = torch::empty({x.size(0), od, oh, ow}, x.options());
    const int64_t count = output.numel();
    if (count == 0) return output;

    constexpr int threads = 256;
    const int blocks = static_cast<int>((count + threads - 1) / threads);
    const cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (pool) {
        fused_channel_reduce_kernel<true><<<blocks, threads, 0, stream>>>(
            x.data_ptr<float>(), subtract.data_ptr<float>(),
            output.data_ptr<float>(), count,
            depth, height, width, od, oh, ow);
    } else {
        fused_channel_reduce_kernel<false><<<blocks, threads, 0, stream>>>(
            x.data_ptr<float>(), subtract.data_ptr<float>(),
            output.data_ptr<float>(), count,
            depth, height, width, od, oh, ow);
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="fused_pool_softmax_subtract_swish_max16",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_channel_reduce_cuda"],
    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))

        def triple(value):
            return (value, value, value) if isinstance(value, int) else tuple(value)

        actual_stride = pool_kernel_size if pool_stride is None else pool_stride
        self._fuse_pool = (
            triple(pool_kernel_size) == (2, 2, 2)
            and triple(actual_stride) == (2, 2, 2)
            and triple(pool_padding) == (0, 0, 0)
        )

    def forward(self, x):
        x = self.conv_transpose(x)

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and self.subtract.dtype == torch.float32
            and x.shape[1] == 16
        ):
            if not self._fuse_pool:
                x = self.max_pool(x)
            return _fused_ops.fused_channel_reduce_cuda(
                x, self.subtract, self._fuse_pool
            )

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