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

__device__ __forceinline__ float clamp_result(float x) {
    x = x < -1.0f ? -1.0f : x;
    return x > 1.0f ? 1.0f : x;
}

__device__ __forceinline__ float finish(float x, float bias) {
    float sigmoid = __fdividef(1.0f, 1.0f + __expf(-(x + 3.0f)));
    return clamp_result((x * sigmoid) / 6.0f - bias);
}

__device__ __forceinline__ void madd(
    float4& a, float x, const float4& w)
{
    a.x = fmaf(x, w.x, a.x);
    a.y = fmaf(x, w.y, a.y);
    a.z = fmaf(x, w.z, a.z);
    a.w = fmaf(x, w.w, a.w);
}

__device__ __forceinline__ void add_bias(float4& a, const float4& b) {
    a.x += b.x;
    a.y += b.y;
    a.z += b.z;
    a.w += b.w;
}

// Two neighboring lanes collectively hold the sixteen output channels.
__device__ __forceinline__ float reduce_pair(
    const float4& a, const float4& b, float bias, int lane)
{
    float sx = __expf(a.x) + __expf(b.x);
    float sy = __expf(a.y) + __expf(b.y);
    float sz = __expf(a.z) + __expf(b.z);
    float sw = __expf(a.w) + __expf(b.w);
    float sum = (sx + sy) + (sz + sw);
    sum += __shfl_xor_sync(0xffffffffu, sum, 1);

    if (sum >= 1.0e-20f && sum <= 1.0e30f) {
        if (lane == 0) {
            float value = __logf(sum);
            float sigmoid =
                __fdividef(sum, sum + 0.049787068367863944f);
            return clamp_result((value * sigmoid) / 6.0f - bias);
        }
        return 0.0f;
    }

    float maximum = fmaxf(
        fmaxf(fmaxf(a.x, a.y), fmaxf(a.z, a.w)),
        fmaxf(fmaxf(b.x, b.y), fmaxf(b.z, b.w)));
    maximum = fmaxf(
        maximum, __shfl_xor_sync(0xffffffffu, maximum, 1));

    if (isinf(maximum))
        maximum = 0.0f;

    sx = __expf(a.x - maximum) + __expf(b.x - maximum);
    sy = __expf(a.y - maximum) + __expf(b.y - maximum);
    sz = __expf(a.z - maximum) + __expf(b.z - maximum);
    sw = __expf(a.w - maximum) + __expf(b.w - maximum);
    sum = (sx + sy) + (sz + sw);
    sum += __shfl_xor_sync(0xffffffffu, sum, 1);

    return lane == 0 ? finish(__logf(sum) + maximum, bias) : 0.0f;
}

__global__ void pack_weights(
    const float* __restrict__ weight,
    float* __restrict__ packed)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < 1296) {
        int channel = i & 15;
        int tap = i >> 4;
        int ic = tap % 3;
        int kernel_index = tap / 3;
        packed[i] = weight[(ic * 16 + channel) * 27 + kernel_index];
    }
}

template<int PD, int PH>
__global__ __launch_bounds__(256, 3)
void deconv_kernel(
    const float* __restrict__ input,
    const float4* __restrict__ packed,
    const float4* __restrict__ conv_bias,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int count)
{
    constexpr int NZ = PD ? 2 : 1;
    constexpr int NY = PH ? 2 : 1;
    constexpr int TAPS = 3 * NZ * NY;
    constexpr int SPATIAL = 16 * 32 * 32;

    __shared__ float4 weights[TAPS * 12];
    __shared__ float4 channel_bias[4];
    __shared__ float scalar_bias;

    int tid = threadIdx.x;

    if (tid < TAPS * 12) {
        int group = tid & 3;
        int t = tid >> 2;
        int kw = t % 3;
        t /= 3;
        int yh = t % NY;
        t /= NY;
        int zd = t % NZ;
        int ic = t / NZ;
        int kd = PD ? 2 * zd : 1;
        int kh = PH ? 2 * yh : 1;
        int k = kd * 9 + kh * 3 + kw;
        weights[tid] = __ldg(packed + (k * 3 + ic) * 4 + group);
    }

    if (tid < 4)
        channel_bias[tid] = __ldg(conv_bias + tid);
    if (tid == 0)
        scalar_bias = __ldg(bias);

    __syncthreads();

    int lane = tid & 1;
    int index = blockIdx.x * 128 + (tid >> 1);
    if (index >= count)
        return;

    int qw = index & 31;
    int remaining = index >> 5;
    int qh = remaining % (32 - PH);
    remaining /= 32 - PH;
    int qd = remaining % (16 - PD);
    int batch = remaining / (16 - PD);
    int base = batch * 3 * SPATIAL;
    bool have_odd = qw < 31;

    float4 even0 = make_float4(0.f, 0.f, 0.f, 0.f);
    float4 even1 = make_float4(0.f, 0.f, 0.f, 0.f);
    float4 odd0 = make_float4(0.f, 0.f, 0.f, 0.f);
    float4 odd1 = make_float4(0.f, 0.f, 0.f, 0.f);

    #pragma unroll
    for (int ic = 0; ic < 3; ++ic) {
        #pragma unroll
        for (int zd = 0; zd < NZ; ++zd) {
            int id = qd + (PD ? 1 - zd : 0);

            #pragma unroll
            for (int yh = 0; yh < NY; ++yh) {
                int ih = qh + (PH ? 1 - yh : 0);
                const float* ptr = input + base + ic * SPATIAL
                                 + (id * 32 + ih) * 32 + qw;
                float current = __ldg(ptr);
                float next = have_odd ? __ldg(ptr + 1) : 0.0f;
                int t = ((ic * NZ + zd) * NY + yh) * 12 + lane;

                madd(even0, current, weights[t + 4]);
                madd(odd0, next, weights[t]);
                madd(odd0, current, weights[t + 8]);

                madd(even1, current, weights[t + 6]);
                madd(odd1, next, weights[t + 2]);
                madd(odd1, current, weights[t + 10]);
            }
        }
    }

    float4 b0 = channel_bias[lane];
    float4 b1 = channel_bias[lane + 2];
    add_bias(even0, b0);
    add_bias(even1, b1);
    add_bias(odd0, b0);
    add_bias(odd1, b1);

    // All lanes participate in both reductions, including the row boundary.
    float even_result = reduce_pair(even0, even1, scalar_bias, lane);
    float odd_result = reduce_pair(odd0, odd1, scalar_bias, lane);

    if (lane == 0) {
        int out_index =
            ((batch * 31 + 2 * qd + PD) * 63 + 2 * qh + PH) * 63
            + 2 * qw;
        output[out_index] = even_result;
        if (have_odd)
            output[out_index + 1] = odd_result;
    }
}

template<int PD, int PH>
void launch_deconv(
    const float* x, const float4* weight, const float4* conv_bias,
    const float* bias, float* output, int batch, cudaStream_t stream)
{
    int count = batch * (16 - PD) * (32 - PH) * 32;
    deconv_kernel<PD, PH><<<(count + 127) / 128, 256, 0, stream>>>(
        x, weight, conv_bias, bias, output, count);
}

torch::Tensor fused_deconv3d(
    torch::Tensor x, torch::Tensor weight,
    torch::Tensor conv_bias, torch::Tensor bias)
{
    const c10::cuda::CUDAGuard guard(x.device());
    int batch = x.size(0);
    auto output = torch::empty({batch, 1, 31, 63, 63}, x.options());

    if (batch == 0)
        return output;

    auto packed = torch::empty({1296}, weight.options());
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    pack_weights<<<6, 256, 0, stream>>>(
        weight.data_ptr<float>(), packed.data_ptr<float>());

    const float* xp = x.data_ptr<float>();
    const float4* wp =
        reinterpret_cast<const float4*>(packed.data_ptr<float>());
    const float4* cp =
        reinterpret_cast<const float4*>(conv_bias.data_ptr<float>());
    const float* bp = bias.data_ptr<float>();
    float* op = output.data_ptr<float>();

    launch_deconv<0, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_deconv<0, 1>(xp, wp, cp, bp, op, batch, stream);
    launch_deconv<1, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_deconv<1, 1>(xp, wp, cp, bp, op, batch, stream);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}

__global__ void postprocess_kernel(
    const float* __restrict__ input,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int64_t spatial, int64_t total, int channels)
{
    int64_t index = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    if (index >= total)
        return;

    int64_t batch = index / spatial;
    int64_t position = index - batch * spatial;
    int64_t base = batch * channels * spatial + position;

    float maximum = -INFINITY;
    for (int c = 0; c < channels; ++c)
        maximum = fmaxf(maximum, input[base + c * spatial]);

    if (isinf(maximum))
        maximum = 0.0f;

    float sum = 0.0f;
    for (int c = 0; c < channels; ++c)
        sum += __expf(input[base + c * spatial] - maximum);

    output[index] = finish(__logf(sum) + maximum, bias[0]);
}

torch::Tensor fused_postprocess(torch::Tensor x, torch::Tensor bias)
{
    const c10::cuda::CUDAGuard guard(x.device());
    int channels = x.size(1);
    int64_t spatial = x.size(2) * x.size(3) * x.size(4);
    int64_t total = x.size(0) * spatial;
    auto output = torch::empty(
        {x.size(0), 1, x.size(2), x.size(3), x.size(4)}, x.options());

    if (total != 0) {
        cudaStream_t stream = at::cuda::getCurrentCUDAStream();
        postprocess_kernel<<<(total + 255) / 256, 256, 0, stream>>>(
            x.data_ptr<float>(), bias.data_ptr<float>(),
            output.data_ptr<float>(), spatial, total, channels);
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }
    return output;
}
"""

_ops = load_inline(
    name="fused_deconv3d_cooperative_pairs_v6",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_deconv3d", "fused_postprocess"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3", "--fmad=false"],
    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))

        conv = self.conv_transpose
        self._specialized = (
            in_channels == 3
            and out_channels == 16
            and conv.kernel_size == (3, 3, 3)
            and conv.stride == (2, 2, 2)
            and conv.padding == (1, 1, 1)
            and conv.dilation == (1, 1, 1)
            and conv.output_padding == (0, 0, 0)
            and conv.groups == 1
        )

    def forward(self, x):
        conv = self.conv_transpose
        if (
            self._specialized
            and x.is_cuda
            and x.ndim == 5
            and tuple(x.shape[1:]) == (3, 16, 32, 32)
            and x.dtype == torch.float32
            and conv.weight.dtype == torch.float32
            and conv.bias is not None
            and conv.bias.dtype == torch.float32
            and self.bias.dtype == torch.float32
            and conv.weight.device == x.device
            and conv.bias.device == x.device
            and self.bias.device == x.device
            and not torch.is_autocast_enabled()
        ):
            return _ops.fused_deconv3d(
                x.contiguous(),
                conv.weight.contiguous(),
                conv.bias.contiguous(),
                self.bias.contiguous(),
            )

        x = conv(x)
        if (
            x.is_cuda
            and x.ndim == 5
            and x.dtype == torch.float32
            and self.bias.dtype == torch.float32
            and self.bias.device == x.device
        ):
            return _ops.fused_postprocess(
                x.contiguous(), self.bias.contiguous()
            )

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