```python
# Copyright (c) Meta Platforms, Inc. and affiliates.
# All rights reserved.
#
# This source code is licensed under the license found in the
# LICENSE file in the root directory of this source tree.

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

# Fully fused CUDA kernel:
#   ConvTranspose3d (direct gather formulation, skipping stride-invalid taps)
#   + LogSumExp over channels
#   + x * sigmoid(x + 3) / 6
#   + bias subtraction
#   + clamp(-1, 1)
# Each thread computes all output channels for one spatial location in registers,
# so the huge (N, C_out, D, H, W) conv output never touches global memory.
_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>

template <int C_IN, int C_OUT, int K>
__global__ __launch_bounds__(256) void fused_ct3d_lse_hswish_clamp_kernel(
    const float* __restrict__ x,          // (N, C_IN, D_IN, H_IN, W_IN)
    const float* __restrict__ w,          // (C_IN, C_OUT, K, K, K)  [conv-transpose layout]
    const float* __restrict__ conv_bias,  // (C_OUT)
    const float* __restrict__ sub_bias,   // scalar
    float* __restrict__ out,              // (N, 1, D_OUT, H_OUT, W_OUT)
    const int N,
    const int D_IN, const int H_IN, const int W_IN,
    const int D_OUT, const int H_OUT, const int W_OUT,
    const int stride, const int pad)
{
    constexpr int CHUNK = C_OUT * C_IN;      // weights per kernel offset
    constexpr int OFFSETS = K * K * K;
    constexpr int W_SIZE = CHUNK * OFFSETS;
    __shared__ float sh[W_SIZE + C_OUT];

    // Load + transpose weights: (C_IN, C_OUT, kd, kh, kw) -> [offset][co][ci]
    const int total_w = C_IN * C_OUT * OFFSETS;
    for (int i = threadIdx.x; i < total_w; i += blockDim.x) {
        int tmp = i;
        const int kw = tmp % K; tmp /= K;
        const int kh = tmp % K; tmp /= K;
        const int kd = tmp % K; tmp /= K;
        const int co = tmp % C_OUT;
        const int ci = tmp / C_OUT;
        const int off = (kd * K + kh) * K + kw;
        sh[off * CHUNK + co * C_IN + ci] = w[i];
    }
    if (threadIdx.x < C_OUT) sh[W_SIZE + threadIdx.x] = conv_bias[threadIdx.x];
    __syncthreads();

    const long total = (long)N * D_OUT * H_OUT * W_OUT;
    const long in_ch_stride = (long)D_IN * H_IN * W_IN;
    const long s = (long)blockIdx.x * blockDim.x + threadIdx.x;
    if (s >= total) return;

    int ow = (int)(s % W_OUT);
    long t = s / W_OUT;
    int oh = (int)(t % H_OUT); t /= H_OUT;
    int od = (int)(t % D_OUT);
    const long n = t / D_OUT;

    float acc[C_OUT];
    #pragma unroll
    for (int co = 0; co < C_OUT; ++co) acc[co] = sh[W_SIZE + co];

    const int bd = od + pad;
    const int bh = oh + pad;
    const int bw = ow + pad;
    const float* xin = x + n * ((long)C_IN * in_ch_stride);

    // Gather form of transposed convolution: only taps with
    // (out_pos + pad - k) divisible by stride contribute.
    #pragma unroll
    for (int kd = 0; kd < K; ++kd) {
        const int td = bd - kd;
        if ((td % stride) != 0) continue;
        const int id = td / stride;
        if ((unsigned)id >= (unsigned)D_IN) continue;
        #pragma unroll
        for (int kh = 0; kh < K; ++kh) {
            const int th = bh - kh;
            if ((th % stride) != 0) continue;
            const int ih = th / stride;
            if ((unsigned)ih >= (unsigned)H_IN) continue;
            #pragma unroll
            for (int kw = 0; kw < K; ++kw) {
                const int tw = bw - kw;
                if ((tw % stride) != 0) continue;
                const int iw = tw / stride;
                if ((unsigned)iw >= (unsigned)W_IN) continue;

                const long in_off = ((long)id * H_IN + ih) * W_IN + iw;
                float xv[C_IN];
                #pragma unroll
                for (int ci = 0; ci < C_IN; ++ci)
                    xv[ci] = xin[in_off + (long)ci * in_ch_stride];

                const float* wsh = sh + ((kd * K + kh) * K + kw) * CHUNK;
                #pragma unroll
                for (int co = 0; co < C_OUT; ++co) {
                    float a = acc[co];
                    #pragma unroll
                    for (int ci = 0; ci < C_IN; ++ci)
                        a = fmaf(wsh[co * C_IN + ci], xv[ci], a);
                    acc[co] = a;
                }
            }
        }
    }

    // LogSumExp over channels (numerically stable).
    float m = acc[0];
    #pragma unroll
    for (int co = 1; co < C_OUT; ++co) m = fmaxf(m, acc[co]);
    float sum = 0.f;
    #pragma unroll
    for (int co = 0; co < C_OUT; ++co) sum += expf(acc[co] - m);
    const float lse = m + logf(sum);

    // x * sigmoid(x + 3) / 6
    const float sig = 1.f / (1.f + expf(-(lse + 3.f)));
    float v = (lse * sig) / 6.f;
    // subtract scalar bias, then clamp to [-1, 1]
    v -= *sub_bias;
    v = fminf(fmaxf(v, -1.f), 1.f);
    out[s] = v;
}

torch::Tensor fused_forward(torch::Tensor x, torch::Tensor w, torch::Tensor conv_bias,
                            torch::Tensor sub_bias, int64_t stride, int64_t pad)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat, "x must be CUDA float32");
    TORCH_CHECK(x.dim() == 5, "x must be 5D (N, C, D, H, W)");
    const c10::cuda::CUDAGuard guard(x.device());

    auto xc = x.contiguous();
    auto wc = w.contiguous();
    auto bc = conv_bias.contiguous();
    auto sb = sub_bias.contiguous();

    const int N    = (int)xc.size(0);
    const int D_IN = (int)xc.size(2), H_IN = (int)xc.size(3), W_IN = (int)xc.size(4);
    const int D_OUT = (int)((D_IN - 1) * stride - 2 * pad + wc.size(2));
    const int H_OUT = (int)((H_IN - 1) * stride - 2 * pad + wc.size(3));
    const int W_OUT = (int)((W_IN - 1) * stride - 2 * pad + wc.size(4));

    auto out = at::empty({N, 1, D_OUT, H_OUT, W_OUT}, xc.options());
    const long total = (long)N * D_OUT * H_OUT * W_OUT;
    if (total == 0) return out;

    const int threads = 256;
    const unsigned blocks = (unsigned)((total + threads - 1) / threads);
    auto stream = at::cuda::getCurrentCUDAStream();
    fused_ct3d_lse_hswish_clamp_kernel<3, 16, 3><<<blocks, threads, 0, stream>>>(
        xc.data_ptr<float>(), wc.data_ptr<float>(), bc.data_ptr<float>(),
        sb.data_ptr<float>(), out.data_ptr<float>(),
        N, D_IN, H_IN, W_IN, D_OUT, H_OUT, W_OUT, (int)stride, (int)pad);
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return out;
}
"""

_cpp_source = (
    "torch::Tensor fused_forward(torch::Tensor x, torch::Tensor w, torch::Tensor conv_bias, "
    "torch::Tensor sub_bias, int64_t stride, int64_t pad);"
)

_fused = load_inline(
    name="fused_ct3d_logsumexp_hswish_clamp_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_forward"],
    verbose=False,
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
)


class ModelNew(nn.Module):
    """
    Optimized model: the whole forward pass (3D transposed convolution,
    LogSumExp, HardSwish-like activation, bias subtraction and clamp) is
    executed by a single fused CUDA kernel for the target configuration.
    Falls back to the original PyTorch ops for any other configuration.
    """
    def __init__(self, in_channels, out_channels, kernel_size, stride, padding, bias_shape):
        super(ModelNew, self).__init__()
        self.conv_transpose = nn.ConvTranspose3d(in_channels, out_channels, kernel_size,
                                                 stride=stride, padding=padding)
        self.bias = nn.Parameter(torch.randn(1, 1, 1, 1))

        ct = self.conv_transpose
        self._use_fast = (
            ct.groups == 1
            and ct.bias is not None
            and ct.in_channels == 3
            and ct.out_channels == 16
            and tuple(ct.kernel_size) == (3, 3, 3)
            and tuple(ct.stride) == (2, 2, 2)
            and tuple(ct.padding) == (1, 1, 1)
            and tuple(ct.output_padding) == (0, 0, 0)
            and tuple(ct.dilation) == (1, 1, 1)
        )
        self._stride = ct.stride[0]
        self._pad = ct.padding[0]

    def forward(self, x):
        if (self._use_fast and x.is_cuda and x.dtype == torch.float32
                and self.conv_transpose.weight.is_cuda and self.bias.is_cuda):
            return _fused.fused_forward(
                x, self.conv_transpose.weight, self.conv_transpose.bias,
                self.bias, self._stride, self._pad)
        x = self.conv_transpose(x)
        x = torch.logsumexp(x, dim=1, keepdim=True)
        x = x * torch.sigmoid(x + 3) / 6
        x = x - self.bias
        return torch.clamp(x, min=-1, max=1)
```