import torch
import torch.nn as nn

# U-Net Implementation
class DoubleConv(nn.Module):
    def __init__(self, in_channels, out_channels):
        super().__init__()
        self.double_conv = nn.Sequential(
            nn.Conv2d(in_channels, out_channels, kernel_size=3, padding=1),
            nn.BatchNorm2d(out_channels),
            nn.Softmax(dim=-1),
            nn.Conv2d(out_channels, out_channels, kernel_size=3, padding=1),
            nn.BatchNorm2d(out_channels),
            nn.Softmax(dim=-1)
        )

    def forward(self, x):
        return self.double_conv(x)

class Model(nn.Module):
    def __init__(self, in_channels, out_channels, features):
        """
        :param in_channels: Number of input channels
        :param out_channels: Number of output channels
        :param features: Number of base features (will be doubled in each layer)
        """
        super(Model, self).__init__()
        self.encoder1 = DoubleConv(in_channels, features)
        self.pool1 = nn.MaxPool2d(kernel_size=2, stride=2)
        self.encoder2 = DoubleConv(features, features * 2)
        self.pool2 = nn.MaxPool2d(kernel_size=2, stride=2)
        self.encoder3 = DoubleConv(features * 2, features * 4)
        self.pool3 = nn.MaxPool2d(kernel_size=2, stride=2)
        self.encoder4 = DoubleConv(features * 4, features * 8)
        self.pool4 = nn.MaxPool2d(kernel_size=2, stride=2)

        self.bottleneck = DoubleConv(features * 8, features * 16)

        self.upconv4 = nn.ConvTranspose2d(features * 16, features * 8, kernel_size=2, stride=2)
        self.decoder4 = DoubleConv(features * 16, features * 8)
        self.upconv3 = nn.ConvTranspose2d(features * 8, features * 4, kernel_size=2, stride=2)
        self.decoder3 = DoubleConv(features * 8, features * 4)
        self.upconv2 = nn.ConvTranspose2d(features * 4, features * 2, kernel_size=2, stride=2)
        self.decoder2 = DoubleConv(features * 4, features * 2)
        self.upconv1 = nn.ConvTranspose2d(features * 2, features, kernel_size=2, stride=2)
        self.decoder1 = DoubleConv(features * 2, features)

        self.final_conv = nn.Conv2d(features, out_channels, kernel_size=1)

    def forward(self, x):
        """
        :param x: Input tensor, shape (batch_size, in_channels, height, width)
        :return: Output tensor, shape (batch_size, out_channels, height, width)
        """
        enc1 = self.encoder1(x)
        enc2 = self.encoder2(self.pool1(enc1))
        enc3 = self.encoder3(self.pool2(enc2))
        enc4 = self.encoder4(self.pool3(enc3))

        bottleneck = self.bottleneck(self.pool4(enc4))

        dec4 = self.upconv4(bottleneck)
        dec4 = torch.cat((dec4, enc4), dim=1)
        dec4 = self.decoder4(dec4)
        dec3 = self.upconv3(dec4)
        dec3 = torch.cat((dec3, enc3), dim=1)
        dec3 = self.decoder3(dec3)
        dec2 = self.upconv2(dec3)
        dec2 = torch.cat((dec2, enc2), dim=1)
        dec2 = self.decoder2(dec2)
        dec1 = self.upconv1(dec2)
        dec1 = torch.cat((dec1, enc1), dim=1)
        dec1 = self.decoder1(dec1)

        return self.final_conv(dec1)
    
batch_size = 8
in_channels = 8
out_channels = 4
height = 64
width = 512
features = 64
# Test code for UNet
def get_inputs():
    return [torch.rand(batch_size, in_channels, height, width)]

def get_init_inputs():
    return [in_channels, out_channels, features]

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

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

torch::Tensor bn_softmax_cuda(torch::Tensor x, torch::Tensor weight, torch::Tensor bias,
                              torch::Tensor running_mean, torch::Tensor running_var, double eps);
torch::Tensor softmax_cuda(torch::Tensor x);

torch::Tensor bn_softmax_forward(torch::Tensor x, torch::Tensor weight, torch::Tensor bias,
                                 torch::Tensor running_mean, torch::Tensor running_var, double eps) {
    return bn_softmax_cuda(x, weight, bias, running_mean, running_var, eps);
}

torch::Tensor softmax_forward(torch::Tensor x) {
    return softmax_cuda(x);
}
'''

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

__global__ void bn_softmax_kernel(const float* __restrict__ x,
                                  const float* __restrict__ gamma,
                                  const float* __restrict__ beta,
                                  const float* __restrict__ mean,
                                  const float* __restrict__ var,
                                  float* __restrict__ out,
                                  int rows, int C, int H, int W, float eps) {
    int row = blockIdx.x;
    if (row >= rows) return;
    int tid = threadIdx.x;
    int c = (row / H) % C;
    int base = row * W;
    float scale = gamma[c] * rsqrtf(var[c] + eps);
    float shift = beta[c] - mean[c] * scale;

    __shared__ float smem[256];
    float m = -FLT_MAX;
    for (int i = tid; i < W; i += blockDim.x) {
        float v = x[base + i] * scale + shift;
        m = fmaxf(m, v);
    }
    smem[tid] = m;
    __syncthreads();
    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
        if (tid < s) smem[tid] = fmaxf(smem[tid], smem[tid + s]);
        __syncthreads();
    }
    float mx = smem[0];
    float sum = 0.0f;
    for (int i = tid; i < W; i += blockDim.x) {
        float v = x[base + i] * scale + shift;
        sum += __expf(v - mx);
    }
    smem[tid] = sum;
    __syncthreads();
    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
        if (tid < s) smem[tid] += smem[tid + s];
        __syncthreads();
    }
    float inv = 1.0f / smem[0];
    for (int i = tid; i < W; i += blockDim.x) {
        float v = x[base + i] * scale + shift;
        out[base + i] = __expf(v - mx) * inv;
    }
}

__global__ void softmax_kernel(const float* __restrict__ x, float* __restrict__ out,
                               int rows, int W) {
    int row = blockIdx.x;
    if (row >= rows) return;
    int tid = threadIdx.x;
    int base = row * W;
    __shared__ float smem[256];
    float m = -FLT_MAX;
    for (int i = tid; i < W; i += blockDim.x) m = fmaxf(m, x[base + i]);
    smem[tid] = m;
    __syncthreads();
    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
        if (tid < s) smem[tid] = fmaxf(smem[tid], smem[tid + s]);
        __syncthreads();
    }
    float mx = smem[0];
    float sum = 0.0f;
    for (int i = tid; i < W; i += blockDim.x) sum += __expf(x[base + i] - mx);
    smem[tid] = sum;
    __syncthreads();
    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
        if (tid < s) smem[tid] += smem[tid + s];
        __syncthreads();
    }
    float inv = 1.0f / smem[0];
    for (int i = tid; i < W; i += blockDim.x) out[base + i] = __expf(x[base + i] - mx) * inv;
}

torch::Tensor bn_softmax_cuda(torch::Tensor x, torch::Tensor weight, torch::Tensor bias,
                              torch::Tensor running_mean, torch::Tensor running_var, double eps) {
    auto out = torch::empty_like(x);
    int64_t N = x.size(0), C = x.size(1), H = x.size(2), W = x.size(3);
    int rows = (int)(N * C * H);
    bn_softmax_kernel<<<rows, 256>>>(x.data_ptr<float>(), weight.data_ptr<float>(), bias.data_ptr<float>(),
                                     running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
                                     out.data_ptr<float>(), rows, (int)C, (int)H, (int)W, (float)eps);
    return out;
}

torch::Tensor softmax_cuda(torch::Tensor x) {
    auto out = torch::empty_like(x);
    int64_t rows64 = x.numel() / x.size(3);
    int W = (int)x.size(3);
    softmax_kernel<<<(int)rows64, 256>>>(x.data_ptr<float>(), out.data_ptr<float>(), (int)rows64, W);
    return out;
}
'''

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

class DoubleConv(nn.Module):
    def __init__(self, in_channels, out_channels):
        super().__init__()
        self.double_conv = nn.Sequential(
            nn.Conv2d(in_channels, out_channels, kernel_size=3, padding=1),
            nn.BatchNorm2d(out_channels),
            nn.Softmax(dim=-1),
            nn.Conv2d(out_channels, out_channels, kernel_size=3, padding=1),
            nn.BatchNorm2d(out_channels),
            nn.Softmax(dim=-1)
        )

    def forward(self, x):
        return self.double_conv(x)

class ModelNew(nn.Module):
    def __init__(self, in_channels, out_channels, features):
        super(ModelNew, self).__init__()
        self.encoder1 = DoubleConv(in_channels, features)
        self.pool1 = nn.MaxPool2d(kernel_size=2, stride=2)
        self.encoder2 = DoubleConv(features, features * 2)
        self.pool2 = nn.MaxPool2d(kernel_size=2, stride=2)
        self.encoder3 = DoubleConv(features * 2, features * 4)
        self.pool3 = nn.MaxPool2d(kernel_size=2, stride=2)
        self.encoder4 = DoubleConv(features * 4, features * 8)
        self.pool4 = nn.MaxPool2d(kernel_size=2, stride=2)

        self.bottleneck = DoubleConv(features * 8, features * 16)

        self.upconv4 = nn.ConvTranspose2d(features * 16, features * 8, kernel_size=2, stride=2)
        self.decoder4 = DoubleConv(features * 16, features * 8)
        self.upconv3 = nn.ConvTranspose2d(features * 8, features * 4, kernel_size=2, stride=2)
        self.decoder3 = DoubleConv(features * 8, features * 4)
        self.upconv2 = nn.ConvTranspose2d(features * 4, features * 2, kernel_size=2, stride=2)
        self.decoder2 = DoubleConv(features * 4, features * 2)
        self.upconv1 = nn.ConvTranspose2d(features * 2, features, kernel_size=2, stride=2)
        self.decoder1 = DoubleConv(features * 2, features)

        self.final_conv = nn.Conv2d(features, out_channels, kernel_size=1)

    def _bn_sm(self, x, bn):
        # eval path fuses affine batchnorm and softmax into one CUDA kernel.
        # training path preserves BatchNorm's exact running-stat update semantics.
        if bn.training:
            y = F.batch_norm(x, bn.running_mean, bn.running_var, bn.weight, bn.bias,
                             True, bn.momentum, bn.eps)
            return _ext.softmax_forward(y.contiguous())
        else:
            return _ext.bn_softmax_forward(x.contiguous(), bn.weight, bn.bias,
                                           bn.running_mean, bn.running_var, bn.eps)

    def _dc(self, x, dc):
        m = dc.double_conv
        c0 = m[0]
        x = F.conv2d(x, c0.weight, c0.bias, c0.stride, c0.padding, c0.dilation, c0.groups)
        x = self._bn_sm(x, m[1])
        c1 = m[3]
        x = F.conv2d(x, c1.weight, c1.bias, c1.stride, c1.padding, c1.dilation, c1.groups)
        x = self._bn_sm(x, m[4])
        return x

    def forward(self, x):
        enc1 = self._dc(x, self.encoder1)
        enc2 = self._dc(F.max_pool2d(enc1, 2, 2), self.encoder2)
        enc3 = self._dc(F.max_pool2d(enc2, 2, 2), self.encoder3)
        enc4 = self._dc(F.max_pool2d(enc3, 2, 2), self.encoder4)

        bottleneck = self._dc(F.max_pool2d(enc4, 2, 2), self.bottleneck)

        u = self.upconv4
        dec4 = F.conv_transpose2d(bottleneck, u.weight, u.bias, u.stride, u.padding, u.output_padding, u.groups, u.dilation)
        dec4 = self._dc(torch.cat((dec4, enc4), dim=1), self.decoder4)

        u = self.upconv3
        dec3 = F.conv_transpose2d(dec4, u.weight, u.bias, u.stride, u.padding, u.output_padding, u.groups, u.dilation)
        dec3 = self._dc(torch.cat((dec3, enc3), dim=1), self.decoder3)

        u = self.upconv2
        dec2 = F.conv_transpose2d(dec3, u.weight, u.bias, u.stride, u.padding, u.output_padding, u.groups, u.dilation)
        dec2 = self._dc(torch.cat((dec2, enc2), dim=1), self.decoder2)

        u = self.upconv1
        dec1 = F.conv_transpose2d(dec2, u.weight, u.bias, u.stride, u.padding, u.output_padding, u.groups, u.dilation)
        dec1 = self._dc(torch.cat((dec1, enc1), dim=1), self.decoder1)

        c = self.final_conv
        return F.conv2d(dec1, c.weight, c.bias, c.stride, c.padding, c.dilation, c.groups)