import torch
import torch.nn as nn

class Model(nn.Module):
    """
    Simple model that performs RMS Normalization.
    """
    def __init__(self, num_features: int, eps: float = 1e-5):
        """
        Initializes the RMSNorm layer.

        Args:
            num_features (int): Number of features in the input tensor.
            eps (float, optional): A small value added to the denominator to avoid division by zero. Defaults to 1e-5.
        """
        super(Model, self).__init__()
        self.num_features = num_features
        self.eps = eps

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        Applies RMS Normalization to the input tensor.

        Args:
            x (torch.Tensor): Input tensor of shape (batch_size, num_features, *).

        Returns:
            torch.Tensor: Output tensor with RMS Normalization applied, same shape as input.
        """
        # Calculate the RMS along the feature dimension
        rms = torch.sqrt(torch.mean(x ** 2, dim=1, keepdim=True) + self.eps)

        # Normalize the input by dividing by the RMS
        return x / rms

batch_size = 112
features = 64
dim1 = 512
dim2 = 512

def get_inputs():
    x = torch.rand(batch_size, features, dim1, dim2)
    return [x]

def get_init_inputs():
    return [features]

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

_cpp_src = r'''
#include <torch/extension.h>
torch::Tensor rmsnorm_forward(torch::Tensor x, double eps);
'''

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

__global__ void rmsnorm64_512_kernel(const float* __restrict__ x, float* __restrict__ y, float eps) {
    int s = (int)blockIdx.x * 256 + (int)threadIdx.x;
    int base = ((int)blockIdx.y << 24) + s;  // b * 64 * 262144 + s

    float v[64];
    float sum = 0.0f;
    #pragma unroll
    for (int c = 0; c < 64; ++c) {
        float t = x[base + (c << 18)];
        v[c] = t;
        sum = fmaf(t, t, sum);
    }

    float inv = rsqrtf(sum * 0.015625f + eps);
    #pragma unroll
    for (int c = 0; c < 64; ++c) {
        y[base + (c << 18)] = v[c] * inv;
    }
}

__global__ void rmsnorm64_kernel(const float* __restrict__ x, float* __restrict__ y,
                                 int spatial, float eps) {
    int s = blockIdx.x * blockDim.x + threadIdx.x;
    if (s >= spatial) return;
    int base = ((int)blockIdx.y * 64) * spatial + s;

    float v[64];
    float sum = 0.0f;
    #pragma unroll
    for (int c = 0; c < 64; ++c) {
        float t = x[base + c * spatial];
        v[c] = t;
        sum = fmaf(t, t, sum);
    }

    float inv = rsqrtf(sum * 0.015625f + eps);
    #pragma unroll
    for (int c = 0; c < 64; ++c) {
        y[base + c * spatial] = v[c] * inv;
    }
}

__global__ void rmsnorm_generic_kernel(const float* __restrict__ x, float* __restrict__ y,
                                       int spatial, int C, float eps) {
    int s = blockIdx.x * blockDim.x + threadIdx.x;
    if (s >= spatial) return;
    int base = ((int)blockIdx.y * C) * spatial + s;

    float sum = 0.0f;
    for (int c = 0; c < C; ++c) {
        float t = x[base + c * spatial];
        sum = fmaf(t, t, sum);
    }
    float inv = rsqrtf(sum / (float)C + eps);

    for (int c = 0; c < C; ++c) {
        y[base + c * spatial] = x[base + c * spatial] * inv;
    }
}

torch::Tensor rmsnorm_forward(torch::Tensor x, double eps) {
    auto y = torch::empty_like(x);

    int C = (int)x.size(1);
    int B = (int)x.size(0);
    int64_t spatial64 = 1;
    for (int64_t i = 2; i < x.dim(); ++i) spatial64 *= x.size(i);
    int spatial = (int)spatial64;

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    if (C == 64 && spatial == 262144 && x.scalar_type() == at::kFloat) {
        rmsnorm64_512_kernel<<<dim3(1024, B), dim3(256), 0, stream>>>(
            x.data_ptr<float>(), y.data_ptr<float>(), (float)eps);
    } else if (C == 64 && x.scalar_type() == at::kFloat) {
        const int threads = 256;
        rmsnorm64_kernel<<<dim3((spatial + threads - 1) / threads, B), dim3(threads), 0, stream>>>(
            x.data_ptr<float>(), y.data_ptr<float>(), spatial, (float)eps);
    } else {
        const int threads = 256;
        rmsnorm_generic_kernel<<<dim3((spatial + threads - 1) / threads, B), dim3(threads), 0, stream>>>(
            x.data_ptr<float>(), y.data_ptr<float>(), spatial, C, (float)eps);
    }
    return y;
}
'''

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

class ModelNew(nn.Module):
    def __init__(self, num_features: int, eps: float = 1e-5):
        super(ModelNew, self).__init__()
        self.num_features = num_features
        self.eps = eps

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return _ext.rmsnorm_forward(x, self.eps)