```python
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.cpp_extension import load_inline

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/AccumulateType.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>

template <typename T>
__device__ __forceinline__ T activation(T x);

template <>
__device__ __forceinline__ float activation<float>(float x) {
    return tanhf(x);
}

template <>
__device__ __forceinline__ double activation<double>(double x) {
    return tanh(x);
}

template <typename scalar_t>
__global__ void fused_bias_tanh_kernel(
    scalar_t* __restrict__ output,
    const scalar_t* __restrict__ conv_bias,
    const scalar_t* __restrict__ subtract_bias,
    int64_t spatial,
    int64_t channels,
    int64_t batch_channels
) {
    using acc_t = at::acc_type<scalar_t, true>;

    for (int64_t bc = blockIdx.y; bc < batch_channels; bc += gridDim.y) {
        const int64_t channel = bc % channels;
        const acc_t add = static_cast<acc_t>(conv_bias[channel]);
        const acc_t sub = static_cast<acc_t>(subtract_bias[channel]);
        const int64_t base = bc * spatial;

        for (int64_t tile = static_cast<int64_t>(blockIdx.x) * 1024;
             tile < spatial;
             tile += static_cast<int64_t>(gridDim.x) * 1024) {
            #pragma unroll
            for (int k = 0; k < 4; ++k) {
                const int64_t position = tile + threadIdx.x + k * 256;
                if (position < spatial) {
                    // Preserve the rounding of the separate PyTorch operations.
                    const scalar_t biased = static_cast<scalar_t>(
                        static_cast<acc_t>(output[base + position]) + add);
                    const scalar_t shifted = static_cast<scalar_t>(
                        static_cast<acc_t>(biased) - sub);
                    output[base + position] = static_cast<scalar_t>(
                        activation<acc_t>(static_cast<acc_t>(shifted)));
                }
            }
        }
    }
}

torch::Tensor fused_bias_tanh_cuda(
    torch::Tensor output,
    torch::Tensor conv_bias,
    torch::Tensor subtract_bias
) {
    TORCH_CHECK(output.is_cuda() && output.is_contiguous(),
                "output must be a contiguous CUDA tensor");
    TORCH_CHECK(output.dim() == 4, "expected NCHW output");
    TORCH_CHECK(conv_bias.is_contiguous() && subtract_bias.is_contiguous(),
                "bias tensors must be contiguous");
    TORCH_CHECK(output.device() == conv_bias.device() &&
                output.device() == subtract_bias.device(),
                "all tensors must have the same device");
    TORCH_CHECK(output.scalar_type() == conv_bias.scalar_type() &&
                output.scalar_type() == subtract_bias.scalar_type(),
                "all tensors must have the same dtype");

    const c10::cuda::CUDAGuard guard(output.device());
    const int64_t channels = output.size(1);
    TORCH_CHECK(conv_bias.numel() == channels &&
                subtract_bias.numel() == channels,
                "expected one bias value per output channel");

    if (output.numel() == 0) return output;

    const int64_t spatial = output.size(2) * output.size(3);
    const int64_t batch_channels = output.size(0) * channels;
    const unsigned int gx = static_cast<unsigned int>(
        std::min<int64_t>((spatial + 1023) / 1024, 65535));
    const unsigned int gy = static_cast<unsigned int>(
        std::min<int64_t>(batch_channels, 65535));

    const auto stream = at::cuda::getCurrentCUDAStream();
    AT_DISPATCH_FLOATING_TYPES_AND2(
        at::ScalarType::Half, at::ScalarType::BFloat16,
        output.scalar_type(), "fused_bias_tanh_cuda", [&] {
            fused_bias_tanh_kernel<scalar_t><<<dim3(gx, gy), 256, 0, stream>>>(
                output.data_ptr<scalar_t>(),
                conv_bias.data_ptr<scalar_t>(),
                subtract_bias.data_ptr<scalar_t>(),
                spatial, channels, batch_channels);
        });
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = load_inline(
    name="conv_transpose_fused_bias_tanh_v1",
    cpp_sources=(
        "torch::Tensor fused_bias_tanh_cuda("
        "torch::Tensor output, torch::Tensor conv_bias, "
        "torch::Tensor subtract_bias);"
    ),
    cuda_sources=_cuda_source,
    functions=["fused_bias_tanh_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class _FusedBiasTanh(torch.autograd.Function):
    @staticmethod
    def forward(ctx, output, conv_bias, subtract_bias):
        ctx.mark_dirty(output)
        _extension.fused_bias_tanh_cuda(output, conv_bias, subtract_bias)
        ctx.save_for_backward(output)
        ctx.subtract_shape = subtract_bias.shape
        return output

    @staticmethod
    def backward(ctx, grad_output):
        (output,) = ctx.saved_tensors
        grad = grad_output * (1 - output.square())
        channel_grad = grad.sum(dim=(0, 2, 3))
        return grad, channel_grad, -channel_grad.reshape(ctx.subtract_shape)


class ModelNew(nn.Module):
    def __init__(
        self,
        in_channels,
        out_channels,
        kernel_size,
        bias_shape,
        stride=2,
        padding=1,
        output_padding=1,
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose2d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
            output_padding=output_padding,
        )
        self.bias = nn.Parameter(torch.randn(bias_shape))

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

        if (
            not x.is_cuda
            or self.bias.shape != (conv.out_channels, 1, 1)
            or self.bias.dtype != x.dtype
            or conv.weight.dtype != x.dtype
            or torch.is_autocast_enabled()
        ):
            return torch.tanh(conv(x) - self.bias)

        output = F.conv_transpose2d(
            x,
            conv.weight,
            None,
            conv.stride,
            conv.padding,
            conv.output_padding,
            conv.groups,
            conv.dilation,
        ).contiguous()

        return _FusedBiasTanh.apply(
            output,
            conv.bias.contiguous(),
            self.bias.contiguous(),
        )
```