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

# C++ Source
cpp_source = """
#include <torch/extension.h>

torch::Tensor smelu_cuda_forward(const torch::Tensor& input, float beta);
"""

# CUDA Source
cuda_source = """
#include <torch/extension.h>
#include <cuda_runtime.h>
#include <math.h>

// Vectorized type
struct __align__(16) Float4 {
    float x, y, z, w;
};

// Core computation logic (Device Inline)
__device__ __forceinline__ float smelu_op(float x, float beta, float inv_4beta) {
    // Branch 1: Right linear region
    if (x >= beta) {
        return x;
    } 
    // Branch 2: Left zero region
    else if (x <= -beta) {
        return 0.0f;
    } 
    // Branch 3: Middle quadratic region
    else {
        float tmp = x + beta;
        return tmp * tmp * inv_4beta;
    }
}

__global__ void smelu_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    const float beta,
    const float inv_4beta, // Precomputed 1 / (4 * beta)
    const int n_elements) 
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;

    // 1. Vectorized Loop
    int vec_loops = n_elements / 4;
    const Float4* vec_input = reinterpret_cast<const Float4*>(input);
    Float4* vec_output = reinterpret_cast<Float4*>(output);

    for (int i = idx; i < vec_loops; i += stride) {
        Float4 in_val = vec_input[i];
        Float4 out_val;
        
        out_val.x = smelu_op(in_val.x, beta, inv_4beta);
        out_val.y = smelu_op(in_val.y, beta, inv_4beta);
        out_val.z = smelu_op(in_val.z, beta, inv_4beta);
        out_val.w = smelu_op(in_val.w, beta, inv_4beta);
        
        vec_output[i] = out_val;
    }

    // 2. Scalar Loop for tail
    int tail_start = vec_loops * 4;
    for (int i = tail_start + idx; i < n_elements; i += stride) {
        output[i] = smelu_op(input[i], beta, inv_4beta);
    }
}

// Wrapper
torch::Tensor smelu_cuda_forward(const torch::Tensor& input, float beta) {
    TORCH_CHECK(input.is_cuda(), "Input tensor must be a CUDA tensor");
    TORCH_CHECK(input.is_contiguous(), "Input tensor must be contiguous");
    TORCH_CHECK(beta > 0, "Beta must be positive");

    auto output = torch::empty_like(input);
    const int n_elements = input.numel();
    
    // Precompute constant inverse to replace division with multiplication in kernel
    float inv_4beta = 1.0f / (4.0f * beta);

    const int block_size = 256;
    int grid_size = (n_elements + block_size * 4 - 1) / (block_size * 4);
    if (grid_size > 65535) grid_size = 65535;

    smelu_kernel<<<grid_size, block_size>>>(
        input.data_ptr<float>(),
        output.data_ptr<float>(),
        beta,
        inv_4beta,
        n_elements
    );

    return output;
}
"""

smelu_op = load_inline(
    name='smelu_op',
    cpp_sources=cpp_source,
    cuda_sources=cuda_source,
    functions=['smelu_cuda_forward'],
    verbose=False,
    extra_cuda_cflags=['-O3']
)

class SmeLUNew(nn.Module):
    def __init__(self, beta=2.0):
        super(SmeLUNew, self).__init__()
        self.beta = beta

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return smelu_op.smelu_cuda_forward(x, self.beta)

class ModelNew(nn.Module):
    def __init__(self, beta=2.0):
        super(ModelNew, self).__init__()
        self.act = SmeLUNew(beta=beta)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return self.act(x)