import torch
import torch.nn as nn

class Model(nn.Module):
    """
    Simple model that performs a matrix multiplication, applies Mish, and applies Mish again.
    """
    def __init__(self, in_features, out_features):
        super(Model, self).__init__()
        self.linear = nn.Linear(in_features, out_features)

    def forward(self, x):
        x = self.linear(x)
        x = torch.nn.functional.mish(x)
        x = torch.nn.functional.mish(x)
        return x

batch_size = 1024
in_features = 8192
out_features = 8192

def get_inputs():
    return [torch.rand(batch_size, in_features)]

def get_init_inputs():
    return [in_features, out_features]

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

_cpp_src = r'''
#include <torch/extension.h>

// CUDA launcher (defined in cuda_sources)
void mish2_launcher(torch::Tensor y);

torch::Tensor linear_mish2(torch::Tensor x, torch::Tensor weight, torch::Tensor bias) {
    // Use PyTorch/ATen's optimized cuBLAS path for the large GEMM, then fuse the
    // two Mish applications into a single custom CUDA kernel.
    auto y = at::linear(x, weight, bias);
    mish2_launcher(y);
    return y;
}
'''

_cuda_src = r'''
#include <torch/extension.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <ATen/cuda/CUDAContext.h>

__device__ __forceinline__ float softplus_approx(float x) {
    // Matches PyTorch softplus beta=1 threshold=20 closely for fp32.
    if (x > 20.0f) return x;
    if (x < -20.0f) return expf(x);
    return log1pf(expf(x));
}

__device__ __forceinline__ float mish_once(float x) {
    return x * tanhf(softplus_approx(x));
}

__global__ void mish2_kernel(float* __restrict__ y, int64_t n) {
    int64_t idx = (int64_t)blockIdx.x * blockDim.x + threadIdx.x;
    int64_t stride = (int64_t)blockDim.x * gridDim.x;
    for (; idx < n; idx += stride) {
        float v = y[idx];
        v = mish_once(v);
        v = mish_once(v);
        y[idx] = v;
    }
}

void mish2_launcher(torch::Tensor y) {
    const int64_t n = y.numel();
    if (n == 0) return;
    const int threads = 256;
    int blocks = (int)((n + threads - 1) / threads);
    if (blocks > 65535) blocks = 65535;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    mish2_kernel<<<blocks, threads, 0, stream>>>(y.data_ptr<float>(), n);
}
'''

_ext = load_inline(
    name='linear_mish2_ext',
    cpp_sources=[_cpp_src],
    cuda_sources=[_cuda_src],
    functions=['linear_mish2'],
    extra_cflags=['-O3'],
    extra_cuda_cflags=['-O3', '--use_fast_math'],
    verbose=False,
)

class ModelNew(nn.Module):
    """
    Drop-in replacement: linear layer followed by Mish twice.  The GEMM uses
    ATen/cuBLAS and the two activations are fused in one CUDA kernel.
    """
    def __init__(self, in_features, out_features):
        super(ModelNew, self).__init__()
        self.linear = nn.Linear(in_features, out_features)

    def forward(self, x):
        return _ext.linear_mish2(x, self.linear.weight, self.linear.bias)