```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 that performs, in a single memory pass:
#   out = hardtanh(gelu(tanh(swish(x + add_value))))
# The matmul (nn.Linear) is left to cuBLAS via F.linear, which already fuses
# the linear bias with an addmm call.
fused_act_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda_runtime.h>
#include <cstdint>

__device__ __forceinline__ float act_chain(float v) {
    // Swish: v * sigmoid(v)
    float sig = 1.0f / (1.0f + expf(-v));
    v = v * sig;
    // Tanh
    v = tanhf(v);
    // GELU (exact / erf formulation, matches PyTorch default)
    v = 0.5f * v * (1.0f + erff(v * 0.70710678118654752440f));
    // Hardtanh clamp to [-1, 1]
    v = fminf(fmaxf(v, -1.0f), 1.0f);
    return v;
}

__global__ void fused_act_vec4_kernel(const float4* __restrict__ in,
                                      const float4* __restrict__ addv,
                                      float4* __restrict__ out,
                                      int n4, int row4) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n4) {
        float4 v = in[i];
        float4 a = addv[i % row4];
        v.x = act_chain(v.x + a.x);
        v.y = act_chain(v.y + a.y);
        v.z = act_chain(v.z + a.z);
        v.w = act_chain(v.w + a.w);
        out[i] = v;
    }
}

__global__ void fused_act_scalar_kernel(const float* __restrict__ in,
                                        const float* __restrict__ addv,
                                        float* __restrict__ out,
                                        int n, int row) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        out[i] = act_chain(in[i] + addv[i % row]);
    }
}

torch::Tensor fused_act_cuda(torch::Tensor x, torch::Tensor addv) {
    TORCH_CHECK(x.is_cuda() && addv.is_cuda(), "Inputs must be CUDA tensors.");
    TORCH_CHECK(x.scalar_type() == at::kFloat && addv.scalar_type() == at::kFloat,
                "Only float32 is supported.");

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

    const int n = (int)xc.numel();
    const int row = (int)ac.numel();
    TORCH_CHECK(row > 0 && (n % row) == 0, "add_value must divide the row size.");

    const int threads = 256;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    bool aligned =
        ((reinterpret_cast<std::uintptr_t>(xc.data_ptr<float>()) & 0xF) == 0) &&
        ((reinterpret_cast<std::uintptr_t>(ac.data_ptr<float>()) & 0xF) == 0) &&
        ((reinterpret_cast<std::uintptr_t>(out.data_ptr<float>()) & 0xF) == 0);

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

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

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


class ModelNew(nn.Module):
    """
    Optimized version:
      - nn.Linear (matmul + bias) executed with a single cuBLAS addmm call.
      - The bias add, Swish, Tanh, GELU and Hardtanh are fused into one
        custom CUDA kernel (single memory pass, vectorized float4 accesses).
    """
    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))

    def forward(self, x):
        if x.is_cuda and x.dtype == torch.float32:
            h = F.linear(x, self.matmul.weight, self.matmul.bias)
            return fused_act_module.fused_act_cuda(h, self.add_value)
        # Fallback (CPU / unsupported dtype): original reference computation.
        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.0, max_val=1.0)  # Hardtanh
        return x
```