```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_pairs(
    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_output(float value) {
    value = value < -1.0f ? -1.0f : value;
    return value > 1.0f ? 1.0f : value;
}

__device__ __forceinline__ float finish_lse(float value, float bias) {
    float sigmoid = __fdividef(
        1.0f, 1.0f + __expf(-(value + 3.0f)));
    return clamp_output((value * sigmoid) / 6.0f - bias);
}

__device__ __forceinline__ void accumulate4(
    float4& accumulator, float value, const float4& weight)
{
    accumulator.x = fmaf(value, weight.x, accumulator.x);
    accumulator.y = fmaf(value, weight.y, accumulator.y);
    accumulator.z = fmaf(value, weight.z, accumulator.z);
    accumulator.w = fmaf(value, weight.w, accumulator.w);
}

__device__ __forceinline__ float reduce_finish16(
    const float4* values, float bias)
{
    float sx = 0.0f;
    float sy = 0.0f;
    float sz = 0.0f;
    float sw = 0.0f;

    #pragma unroll
    for (int g = 0; g < 4; ++g) {
        sx += __expf(values[g].x);
        sy += __expf(values[g].y);
        sz += __expf(values[g].z);
        sw += __expf(values[g].w);
    }

    float sum = (sx + sy) + (sz + sw);

    if (sum >= 1.0e-20f && sum <= 1.0e30f) {
        float value = __logf(sum);
        // sigmoid(log(sum) + 3) = sum / (sum + exp(-3)).
        float sigmoid = __fdividef(sum, sum + 0.049787068367863944f);
        return clamp_output((value * sigmoid) / 6.0f - bias);
    }

    float maximum = -INFINITY;
    #pragma unroll
    for (int g = 0; g < 4; ++g) {
        maximum = fmaxf(maximum, values[g].x);
        maximum = fmaxf(maximum, values[g].y);
        maximum = fmaxf(maximum, values[g].z);
        maximum = fmaxf(maximum, values[g].w);
    }

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

    sx = sy = sz = sw = 0.0f;
    #pragma unroll
    for (int g = 0; g < 4; ++g) {
        sx += __expf(values[g].x - maximum);
        sy += __expf(values[g].y - maximum);
        sz += __expf(values[g].z - maximum);
        sw += __expf(values[g].w - maximum);
    }

    return finish_lse(
        __logf((sx + sy) + (sz + sw)) + maximum, bias);
}

__global__ void pack_weights_kernel(
    const float* __restrict__ weight,
    float* __restrict__ packed)
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= 1296)
        return;

    int channel = index & 15;
    int tap = index >> 4;
    int ic = tap % 3;
    int kernel_index = tap / 3;
    packed[index] = weight[(ic * 16 + channel) * 27 + kernel_index];
}

template<int PD, int PH>
__global__ __launch_bounds__(256)
void deconv_pairs_kernel(
    const float* __restrict__ x,
    const float4* __restrict__ packed_weight,
    const float4* __restrict__ conv_bias,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int count)
{
    constexpr int D = 16;
    constexpr int H = 32;
    constexpr int W = 32;
    constexpr int OD = 31;
    constexpr int OH = 63;
    constexpr int OW = 63;
    constexpr int SPATIAL = D * H * W;
    constexpr int NZ = PD ? 2 : 1;
    constexpr int NY = PH ? 2 : 1;
    constexpr int TAPS = 3 * NZ * NY;

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

    const int tid = threadIdx.x;

    if (tid < TAPS * 12) {
        int group = tid & 3;
        int tap = tid >> 2;
        int kw = tap % 3;
        tap /= 3;
        int yh = tap % NY;
        tap /= NY;
        int zd = tap % NZ;
        int ic = tap / NZ;

        int kd = PD ? 2 * zd : 1;
        int kh = PH ? 2 * yh : 1;
        int kernel_index = kd * 9 + kh * 3 + kw;

        weights[tid] =
            __ldg(packed_weight + (kernel_index * 3 + ic) * 4 + group);
    }

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

    __syncthreads();

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

    const int qw = index & 31;
    int remaining = index >> 5;
    const int qh = remaining % (H - PH);
    remaining /= H - PH;
    const int qd = remaining % (D - PD);
    const int batch = remaining / (D - PD);
    const int input_base = batch * 3 * SPATIAL;
    const bool have_odd = qw < 31;

    float4 even[4];
    float4 odd[4];

    #pragma unroll
    for (int group = 0; group < 4; ++group) {
        even[group] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        odd[group] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
    }

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

            #pragma unroll
            for (int yh = 0; yh < NY; ++yh) {
                const int ih = qh + (PH ? 1 - yh : 0);
                const float* ptr =
                    x + input_base + ic * SPATIAL
                      + (id * H + ih) * W + qw;

                const float current = __ldg(ptr);
                const float next = have_odd ? __ldg(ptr + 1) : 0.0f;
                const int tap = (ic * NZ + zd) * NY + yh;

                #pragma unroll
                for (int group = 0; group < 4; ++group) {
                    accumulate4(
                        even[group], current,
                        weights[tap * 12 + 4 + group]);

                    accumulate4(
                        odd[group], next,
                        weights[tap * 12 + group]);

                    accumulate4(
                        odd[group], current,
                        weights[tap * 12 + 8 + group]);
                }
            }
        }
    }

    #pragma unroll
    for (int group = 0; group < 4; ++group) {
        const float4 b = channel_bias[group];
        even[group].x += b.x;
        even[group].y += b.y;
        even[group].z += b.z;
        even[group].w += b.w;
        odd[group].x += b.x;
        odd[group].y += b.y;
        odd[group].z += b.z;
        odd[group].w += b.w;
    }

    const int out_index =
        ((batch * OD + 2 * qd + PD) * OH + 2 * qh + PH) * OW
        + 2 * qw;

    output[out_index] = reduce_finish16(even, output_bias);
    if (have_odd)
        output[out_index + 1] = reduce_finish16(odd, output_bias);
}

template<int PD, int PH>
void launch_pairs(
    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;
    constexpr int threads = 256;

    deconv_pairs_kernel<PD, PH>
        <<<(count + threads - 1) / threads, threads, 0, stream>>>(
            x, weight, conv_bias, bias, output, count);
}

torch::Tensor fused_deconv3d_pairs(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor conv_bias,
    torch::Tensor bias)
{
    const c10::cuda::CUDAGuard guard(x.device());
    const 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_kernel<<<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_pairs<0, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_pairs<0, 1>(xp, wp, cp, bp, op, batch, stream);
    launch_pairs<1, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_pairs<1, 1>(xp, wp, cp, bp, op, batch, stream);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}

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

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

    if constexpr (CHANNELS == 16) {
        float4 values[4];

        #pragma unroll
        for (int g = 0; g < 4; ++g) {
            values[g] = make_float4(
                input[base + (4 * g + 0) * spatial],
                input[base + (4 * g + 1) * spatial],
                input[base + (4 * g + 2) * spatial],
                input[base + (4 * g + 3) * spatial]);
        }

        output[index] = reduce_finish16(values, bias[0]);
    } else {
        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_lse(__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);
    int spatial = x.size(2) * x.size(3) * x.size(4);
    int 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)
        return output;

    constexpr int threads = 256;
    int blocks = (total + threads - 1) / threads;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (channels == 16) {
        postprocess_kernel<16><<<blocks, threads, 0, stream>>>(
            x.data_ptr<float>(), bias.data_ptr<float>(),
            output.data_ptr<float>(), spatial, total, channels);
    } else {
        postprocess_kernel<0><<<blocks, threads, 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_width_pairs_v5",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_deconv3d_pairs", "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_pairs(
                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)
```