import math

import torch
import torch.nn as nn


class Model(nn.Module):
    """
    Causal grouped-query attention with KV-head expansion.

    Query heads outnumber KV heads. Keys and values are expanded by repeating
    each KV head across its assigned query-head group before causal attention.
    """

    def __init__(self):
        super().__init__()

    def forward(
        self,
        q: torch.Tensor,
        k: torch.Tensor,
        v: torch.Tensor,
    ) -> torch.Tensor:
        batch_size, seq_len, num_q_heads, head_dim = q.shape
        num_kv_heads = k.shape[2]
        group_size = num_q_heads // num_kv_heads

        expanded_k = k.repeat_interleave(group_size, dim=2)
        expanded_v = v.repeat_interleave(group_size, dim=2)

        out = torch.zeros_like(q)
        scale = 1.0 / math.sqrt(head_dim)

        for b in range(batch_size):
            for h in range(num_q_heads):
                for t in range(seq_len):
                    q_t = q[b, t, h]
                    k_hist = expanded_k[b, : t + 1, h]
                    v_hist = expanded_v[b, : t + 1, h]
                    scores = (k_hist * q_t.unsqueeze(0)).sum(dim=-1) * scale
                    attn = torch.softmax(scores, dim=0)
                    out[b, t, h] = (attn.unsqueeze(-1) * v_hist).sum(dim=0)

        return out


batch_size = 2
seq_len = 80
num_q_heads = 8
num_kv_heads = 2
head_dim = 32


def get_inputs():
    p = popcorn_pri
    q = torch.randn(p.jitter_int(batch_size), p.jitter_int(seq_len), p.jitter_int(num_q_heads), p.jitter_int(head_dim, align=8), dtype=torch.float32)
    k = torch.randn(p.jitter_int(batch_size), p.jitter_int(seq_len), p.jitter_int(num_kv_heads), p.jitter_int(head_dim, align=8), dtype=torch.float32)
    v = torch.randn(p.jitter_int(batch_size), p.jitter_int(seq_len), p.jitter_int(num_kv_heads), p.jitter_int(head_dim, align=8), dtype=torch.float32)
    return [q, k, v]


def get_init_inputs():
    return []

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

cuda_source = '''
#include <torch/extension.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>

__global__ void gqa_causal_attn_kernel(
    const float* __restrict__ q,
    const float* __restrict__ k,
    const float* __restrict__ v,
    float* __restrict__ out,
    int batch_size,
    int seq_len,
    int num_q_heads,
    int num_kv_heads,
    int head_dim,
    int group_size
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int total = batch_size * num_q_heads * seq_len;
    if (idx >= total) return;

    int t = idx % seq_len;
    int h = (idx / seq_len) % num_q_heads;
    int b = idx / (seq_len * num_q_heads);

    int kv_h = h / group_size;
    float scale = 1.0f / sqrtf((float)head_dim);

    float scores[256]; // assume seq_len <= 256
    float max_score = -1e30f;
    for (int s = 0; s <= t; s++) {
        float dot = 0.0f;
        for (int d = 0; d < head_dim; d++) {
            float q_val = q[((b * seq_len + t) * num_q_heads + h) * head_dim + d];
            float k_val = k[((b * seq_len + s) * num_kv_heads + kv_h) * head_dim + d];
            dot += q_val * k_val;
        }
        scores[s] = dot * scale;
        if (scores[s] > max_score) max_score = scores[s];
    }

    float sum_exp = 0.0f;
    for (int s = 0; s <= t; s++) {
        scores[s] = expf(scores[s] - max_score);
        sum_exp += scores[s];
    }

    for (int d = 0; d < head_dim; d++) {
        float val = 0.0f;
        for (int s = 0; s <= t; s++) {
            float v_val = v[((b * seq_len + s) * num_kv_heads + kv_h) * head_dim + d];
            val += (scores[s] / sum_exp) * v_val;
        }
        out[((b * seq_len + t) * num_q_heads + h) * head_dim + d] = val;
    }
}

torch::Tensor gqa_causal_attn_cuda(
    torch::Tensor q,
    torch::Tensor k,
    torch::Tensor v
) {
    int batch_size = q.size(0);
    int seq_len = q.size(1);
    int num_q_heads = q.size(2);
    int head_dim = q.size(3);
    int num_kv_heads = k.size(2);
    int group_size = num_q_heads / num_kv_heads;

    auto out = torch::zeros_like(q);

    int total = batch_size * num_q_heads * seq_len;
    int threads = 256;
    int blocks = (total + threads - 1) / threads;

    gqa_causal_attn_kernel<<<blocks, threads>>>(
        q.data_ptr<float>(),
        k.data_ptr<float>(),
        v.data_ptr<float>(),
        out.data_ptr<float>(),
        batch_size,
        seq_len,
        num_q_heads,
        num_kv_heads,
        head_dim,
        group_size
    );

    return out;
}
'''

cpp_source = '''
torch::Tensor gqa_causal_attn_cuda(torch::Tensor q, torch::Tensor k, torch::Tensor v);
'''

module = load_inline(
    name='gqa_attn',
    cpp_sources=cpp_source,
    cuda_sources=cuda_source,
    functions=['gqa_causal_attn_cuda'],
    verbose=False
)

class ModelNew(nn.Module):
    def __init__(self):
        super().__init__()
        self.module = module

    def forward(self, q: torch.Tensor, k: torch.Tensor, v: torch.Tensor) -> torch.Tensor:
        return self.module.gqa_causal_attn_cuda(q, k, v)
```