```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 finish(float value, float bias) {
    float sigmoid = __fdividef(1.0f, 1.0f + __expf(-(value + 3.0f)));
    float result = (value * sigmoid) / 6.0f - bias;
    result = result < -1.0f ? -1.0f : result;
    result = result > 1.0f ? 1.0f : result;
    return result;
}

__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 logsumexp16(const float4* values) {
    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);

    // Avoid the maximum reduction and sixteen subtractions when the
    // unshifted exponential sum is safely representable.
    if (sum >= 1.0e-20f && sum <= 1.0e30f)
        return __logf(sum);

    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 __logf((sx + sy) + (sz + sw)) + maximum;
}

__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, int PW>
__global__ __launch_bounds__(256)
void deconv_reduce_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 NX = PW ? 2 : 1;
    constexpr int TAPS = 3 * NZ * NY * NX;

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

    int tid = threadIdx.x;

    if (tid < TAPS * 4) {
        int group = tid & 3;
        int tap = tid >> 2;
        int xw = tap % NX;
        tap /= NX;
        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 kw = PW ? 2 * xw : 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();

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

    int qw = index % (W - PW);
    int remaining = index / (W - PW);
    int qh = remaining % (H - PH);
    remaining /= H - PH;
    int qd = remaining % (D - PD);
    int batch = remaining / (D - PD);
    int input_base = batch * 3 * SPATIAL;

    float4 values[4];
    #pragma unroll
    for (int group = 0; group < 4; ++group)
        values[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) {
            int id = qd + (PD ? 1 - zd : 0);

            #pragma unroll
            for (int yh = 0; yh < NY; ++yh) {
                int ih = qh + (PH ? 1 - yh : 0);

                #pragma unroll
                for (int xw = 0; xw < NX; ++xw) {
                    int iw = qw + (PW ? 1 - xw : 0);
                    float value = __ldg(
                        x + input_base + ic * SPATIAL
                          + (id * H + ih) * W + iw);

                    constexpr int GROUPS = 4;
                    int tap = ((ic * NZ + zd) * NY + yh) * NX + xw;

                    #pragma unroll
                    for (int group = 0; group < GROUPS; ++group)
                        accumulate4(
                            values[group], value, weights[tap * 4 + group]);
                }
            }
        }
    }

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

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

    output[out_index] = finish(logsumexp16(values), output_bias);
}

template<int PD, int PH, int PW>
void launch_parity(
    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 - PW);
    constexpr int threads = 256;
    deconv_reduce_kernel<PD, PH, PW>
        <<<(count + threads - 1) / threads, threads, 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_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_parity<0, 0, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<0, 0, 1>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<0, 1, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<0, 1, 1>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<1, 0, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<1, 0, 1>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<1, 1, 0>(xp, wp, cp, bp, op, batch, stream);
    launch_parity<1, 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;

    float value;
    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]);
        }

        value = logsumexp16(values);
    } 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);

        value = __logf(sum) + maximum;
    }

    output[index] = finish(value, 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_constant_geometry_v4",
    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)
```