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


_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 <cmath>
#include <cstdint>

// Each thread computes all 16 channels at one output position, then
// immediately reduces them and applies the remaining operations.
template<int PD, int PH, int PW>
__device__ __forceinline__ void compute_voxel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ conv_bias,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int D, int H, int W)
{
    const int nd = D - PD;
    const int nh = H - PH;
    const int nw = W - PW;

    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    const int64_t count = static_cast<int64_t>(nd) * nh * nw;
    if (index >= count)
        return;

    const int iw = index % nw;
    const int ih = (index / nw) % nh;
    const int id = index / (static_cast<int64_t>(nw) * nh);

    float values[16];

    #pragma unroll
    for (int oc = 0; oc < 16; ++oc)
        values[oc] = conv_bias[oc];

    const int64_t input_volume = static_cast<int64_t>(D) * H * W;
    const int64_t batch_base =
        static_cast<int64_t>(blockIdx.y) * 3 * input_volume;

    #pragma unroll
    for (int ic = 0; ic < 3; ++ic) {
        #pragma unroll
        for (int zd = 0; zd < 1 + PD; ++zd) {
            const int kz = PD ? 2 * zd : 1;
            const int xd = id + PD - zd;

            #pragma unroll
            for (int yh = 0; yh < 1 + PH; ++yh) {
                const int ky = PH ? 2 * yh : 1;
                const int xh = ih + PH - yh;

                #pragma unroll
                for (int xw = 0; xw < 1 + PW; ++xw) {
                    const int kx = PW ? 2 * xw : 1;
                    const int xx = iw + PW - xw;

                    const int64_t input_index =
                        batch_base + ic * input_volume +
                        (static_cast<int64_t>(xd) * H + xh) * W + xx;

                    const float input_value = x[input_index];

                    #pragma unroll
                    for (int oc = 0; oc < 16; ++oc) {
                        const int weight_index =
                            (ic * 16 + oc) * 27 + kz * 9 + ky * 3 + kx;
                        values[oc] = fmaf(
                            input_value, weight[weight_index], values[oc]);
                    }
                }
            }
        }
    }

    float maximum = values[0];
    #pragma unroll
    for (int oc = 1; oc < 16; ++oc)
        maximum = fmaxf(maximum, values[oc]);

    float sum = 0.0f;
    #pragma unroll
    for (int oc = 0; oc < 16; ++oc)
        sum += expf(values[oc] - maximum);

    const float lse = isinf(maximum) ? maximum : maximum + logf(sum);
    const float sigmoid = 1.0f / (1.0f + expf(-(lse + 3.0f)));
    float result = (lse * sigmoid) / 6.0f - bias[0];

    if (result < -1.0f) result = -1.0f;
    if (result >  1.0f) result =  1.0f;

    const int od = 2 * id + PD;
    const int oh = 2 * ih + PH;
    const int ow = 2 * iw + PW;
    const int Do = 2 * D - 1;
    const int Ho = 2 * H - 1;
    const int Wo = 2 * W - 1;

    const int64_t output_index =
        ((static_cast<int64_t>(blockIdx.y) * Do + od) * Ho + oh) * Wo + ow;
    output[output_index] = result;
}

__global__ void fused_transpose3d_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ conv_bias,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int D, int H, int W)
{
    // Parity is uniform across each block. It determines exactly which
    // convolution taps contribute, without inserting or reading zeros.
    #define RUN(PD, PH, PW) \
        compute_voxel<PD, PH, PW>( \
            x, weight, conv_bias, bias, output, D, H, W)

    switch (blockIdx.z) {
        case 0: RUN(0, 0, 0); break;
        case 1: RUN(0, 0, 1); break;
        case 2: RUN(0, 1, 0); break;
        case 3: RUN(0, 1, 1); break;
        case 4: RUN(1, 0, 0); break;
        case 5: RUN(1, 0, 1); break;
        case 6: RUN(1, 1, 0); break;
        case 7: RUN(1, 1, 1); break;
    }

    #undef RUN
}

torch::Tensor fused_transpose3d_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor conv_bias,
    torch::Tensor bias)
{
    const c10::cuda::CUDAGuard device_guard(x.device());

    const int64_t N = x.size(0);
    const int D = static_cast<int>(x.size(2));
    const int H = static_cast<int>(x.size(3));
    const int W = static_cast<int>(x.size(4));

    auto output = torch::empty(
        {N, 1, 2LL * D - 1, 2LL * H - 1, 2LL * W - 1},
        x.options());

    if (N == 0)
        return output;

    constexpr int threads = 128;
    const int64_t largest_partition = static_cast<int64_t>(D) * H * W;
    const dim3 grid(
        static_cast<unsigned>((largest_partition + threads - 1) / threads),
        static_cast<unsigned>(N),
        8);

    fused_transpose3d_kernel<<<
        grid, threads, 0, at::cuda::getCurrentCUDAStream()
    >>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        conv_bias.data_ptr<float>(),
        bias.data_ptr<float>(),
        output.data_ptr<float>(),
        D, H, W);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""


_fused_ops = load_inline(
    name="fused_transpose3d_lse_sigmoid_clamp_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_transpose3d_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, bias_shape
    ):
        super().__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))

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

        use_fused = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 5
            and x.shape[1] == 3
            and x.shape[0] <= 65535
            and all(s > 0 for s in x.shape[2:])
            and x.is_contiguous()
            and conv.in_channels == 3
            and conv.out_channels == 16
            and conv.kernel_size == (3, 3, 3)
            and conv.stride == (2, 2, 2)
            and conv.padding == (1, 1, 1)
            and conv.output_padding == (0, 0, 0)
            and conv.dilation == (1, 1, 1)
            and conv.groups == 1
            and conv.weight.dtype == torch.float32
            and conv.weight.device == x.device
            and conv.weight.is_contiguous()
            and conv.bias is not None
            and conv.bias.dtype == torch.float32
            and conv.bias.device == x.device
            and conv.bias.is_contiguous()
            and self.bias.dtype == torch.float32
            and self.bias.device == x.device
            and self.bias.numel() == 1
            and not torch.is_autocast_enabled()
        )

        # Preserve autograd when gradients are requested.
        needs_grad = torch.is_grad_enabled() and (
            x.requires_grad
            or conv.weight.requires_grad
            or (conv.bias is not None and conv.bias.requires_grad)
            or self.bias.requires_grad
        )

        if use_fused and not needs_grad:
            return _fused_ops.fused_transpose3d_cuda(
                x, conv.weight, conv.bias, self.bias
            )

        x = conv(x)
        x = torch.logsumexp(x, dim=1, keepdim=True)
        x = x * torch.sigmoid(x + 3) / 6
        return torch.clamp(x - self.bias, min=-1, max=1)
```