```python
# Copyright (c) Meta Platforms, Inc. and affiliates.
# All rights reserved.
#
# This source code is licensed under the license found in the
# LICENSE file in the root directory of this source tree.

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

# Fused CUDA kernel: (bias add) + (add_value add) + Swish + Tanh + GELU + Hardtanh
# The GEMM itself is delegated to cuBLAS via F.linear, which is already optimal
# for these sizes; all elementwise work is fused into a single kernel pass.
fused_act_cuda_source = """
#include <torch/extension.h>
#include <cuda_runtime.h>

__device__ __forceinline__ float act_chain(float v) {
    // Swish: sigmoid(v) * v
    float s = 1.0f / (1.0f + __expf(-v));
    v = s * v;
    // Tanh
    v = tanhf(v);
    // GELU (exact, erf-based): 0.5 * x * (1 + erf(x / sqrt(2)))
    v = 0.5f * v * (1.0f + erff(v * 0.70710678118654752440f));
    // Hardtanh in [-1, 1]
    v = fminf(fmaxf(v, -1.0f), 1.0f);
    return v;
}

__global__ void fused_act_vec4_kernel(const float4* __restrict__ x,
                                      const float4* __restrict__ bias,
                                      const float4* __restrict__ addv,
                                      float4* __restrict__ out,
                                      long n4, int width4) {
    long idx = (long)blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n4) return;
    int c = (int)(idx % width4);
    float4 xv = x[idx];
    float4 b = bias[c];
    float4 a = addv[c];
    float4 r;
    r.x = act_chain(xv.x + b.x + a.x);
    r.y = act_chain(xv.y + b.y + a.y);
    r.z = act_chain(xv.z + b.z + a.z);
    r.w = act_chain(xv.w + b.w + a.w);
    out[idx] = r;
}

__global__ void fused_act_scalar_kernel(const float* __restrict__ x,
                                        const float* __restrict__ bias,
                                        const float* __restrict__ addv,
                                        float* __restrict__ out,
                                        long n, int width) {
    long idx = (long)blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;
    int c = (int)(idx % width);
    out[idx] = act_chain(x[idx] + bias[c] + addv[c]);
}

torch::Tensor fused_act(torch::Tensor x, torch::Tensor bias, torch::Tensor addv) {
    TORCH_CHECK(x.is_cuda() && bias.is_cuda() && addv.is_cuda(),
                "All inputs must be CUDA tensors");
    TORCH_CHECK(x.scalar_type() == torch::kFloat32, "x must be float32");

    auto xc = x.contiguous();
    auto bc = bias.contiguous();
    auto ac = addv.contiguous();
    auto out = torch::empty_like(xc);

    long n = xc.numel();
    int width = xc.size(-1);
    const int threads = 256;

    if ((n % 4 == 0) && (width % 4 == 0)) {
        long n4 = n / 4;
        int width4 = width / 4;
        long blocks = (n4 + threads - 1) / threads;
        fused_act_vec4_kernel<<<blocks, threads>>>(
            reinterpret_cast<const float4*>(xc.data_ptr<float>()),
            reinterpret_cast<const float4*>(bc.data_ptr<float>()),
            reinterpret_cast<const float4*>(ac.data_ptr<float>()),
            reinterpret_cast<float4*>(out.data_ptr<float>()),
            n4, width4);
    } else {
        long blocks = (n + threads - 1) / threads;
        fused_act_scalar_kernel<<<blocks, threads>>>(
            xc.data_ptr<float>(),
            bc.data_ptr<float>(),
            ac.data_ptr<float>(),
            out.data_ptr<float>(),
            n, width);
    }
    return out;
}
"""

fused_act_cpp_source = (
    "torch::Tensor fused_act(torch::Tensor x, torch::Tensor bias, torch::Tensor addv);"
)

fused_act_module = load_inline(
    name="fused_swish_tanh_gelu_hardtanh",
    cpp_sources=fused_act_cpp_source,
    cuda_sources=fused_act_cuda_source,
    functions=["fused_act"],
    verbose=False,
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3", "--use_fast_math"],
)


class ModelNew(nn.Module):
    """
    Optimized model: cuBLAS GEMM (no epilogue) followed by a single fused CUDA
    kernel that applies linear-bias add, add_value add, Swish, Tanh, GELU and
    Hardtanh in one memory pass (float4 vectorized).
    """
    def __init__(self, in_features, out_features, add_value_shape):
        super(ModelNew, self).__init__()
        self.matmul = nn.Linear(in_features, out_features)
        self.add_value = nn.Parameter(torch.randn(add_value_shape))
        self.fused = fused_act_module

    def forward(self, x):
        if x.is_cuda and x.dtype == torch.float32:
            # GEMM without bias epilogue; bias and add_value are folded into
            # the fused activation kernel to avoid extra elementwise passes.
            y = F.linear(x, self.matmul.weight)
            return self.fused.fused_act(y, self.matmul.bias, self.add_value)

        # Fallback for non-CUDA / non-float32 inputs
        x = self.matmul(x)
        x = x + self.add_value
        x = torch.sigmoid(x) * x  # Swish
        x = torch.tanh(x)
        x = F.gelu(x)  # GELU
        x = F.hardtanh(x, min_val=-1, max_val=1)  # Hardtanh
        return x
```