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

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

__global__ void softmax_width_kernel(const float* __restrict__ x, float* __restrict__ y, int rows, int W) {
    int row = blockIdx.x;
    if (row >= rows) return;
    int tid = threadIdx.x;
    const float* xr = x + ((long long)row) * W;
    float* yr = y + ((long long)row) * W;

    float m = -FLT_MAX;
    for (int j = tid; j < W; j += blockDim.x) {
        float v = xr[j];
        m = fmaxf(m, v);
    }
    __shared__ float smem[256];
    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();
    }
    m = smem[0];

    float sum = 0.0f;
    for (int j = tid; j < W; j += blockDim.x) {
        float e = expf(xr[j] - m);
        sum += e;
        yr[j] = e;
    }
    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 j = tid; j < W; j += blockDim.x) {
        yr[j] *= inv;
    }
}

__global__ void bn_softmax_width_kernel(const float* __restrict__ x,
                                        const float* __restrict__ mean,
                                        const float* __restrict__ var,
                                        const float* __restrict__ gamma,
                                        const float* __restrict__ beta,
                                        float* __restrict__ y,
                                        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;
    float g = gamma ? gamma[c] : 1.0f;
    float b = beta ? beta[c] : 0.0f;
    float scale = g * rsqrtf(var[c] + eps);
    float shift = b - mean[c] * scale;
    const float* xr = x + ((long long)row) * W;
    float* yr = y + ((long long)row) * W;

    float m = -FLT_MAX;
    for (int j = tid; j < W; j += blockDim.x) {
        float v = xr[j] * scale + shift;
        m = fmaxf(m, v);
    }
    __shared__ float smem[256];
    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();
    }
    m = smem[0];

    float sum = 0.0f;
    for (int j = tid; j < W; j += blockDim.x) {
        float e = expf((xr[j] * scale + shift) - m);
        sum += e;
        yr[j] = e;
    }
    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 j = tid; j < W; j += blockDim.x) {
        yr[j] *= inv;
    }
}

torch::Tensor softmax_width(torch::Tensor x) {
    x = x.contiguous();
    TORCH_CHECK(x.is_cuda(), "x must be CUDA");
    TORCH_CHECK(x.scalar_type() == at::kFloat, "only float32 supported");
    int64_t ndim = x.dim();
    TORCH_CHECK(ndim >= 1, "bad dim");
    int W = (int)x.size(ndim - 1);
    int rows = (int)(x.numel() / W);
    auto y = torch::empty_like(x);
    if (rows > 0) {
        softmax_width_kernel<<<rows, 256, 0, at::cuda::getCurrentCUDAStream()>>>(x.data_ptr<float>(), y.data_ptr<float>(), rows, W);
    }
    return y;
}

torch::Tensor bn_softmax_width(torch::Tensor x, torch::Tensor mean, torch::Tensor var,
                               torch::Tensor gamma, torch::Tensor beta, double eps) {
    x = x.contiguous();
    TORCH_CHECK(x.is_cuda(), "x must be CUDA");
    TORCH_CHECK(x.scalar_type() == at::kFloat, "only float32 supported");
    int C = (int)x.size(1);
    int H = (int)x.size(2);
    int W = (int)x.size(3);
    int rows = (int)(x.numel() / W);
    auto y = torch::empty_like(x);
    if (rows > 0) {
        bn_softmax_width_kernel<<<rows, 256, 0, at::cuda::getCurrentCUDAStream()>>>(
            x.data_ptr<float>(), mean.data_ptr<float>(), var.data_ptr<float>(),
            gamma.defined() ? gamma.data_ptr<float>() : nullptr,
            beta.defined() ? beta.data_ptr<float>() : nullptr,
            y.data_ptr<float>(), rows, C, H, W, (float)eps);
    }
    return y;
}
'''

_ext = load_inline(
    name='unet_bn_softmax_ext_v1',
    cpp_sources='',
    cuda_sources=_cuda_src,
    functions=['softmax_width', 'bn_softmax_width'],
    extra_cuda_cflags=['-O3', '--use_fast_math'],
    extra_cflags=['-O3'],
    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)

    @staticmethod
    def _bn_sm_eval(y, bn):
        return _ext.bn_softmax_width(y, bn.running_mean, bn.running_var, bn.weight, bn.bias, bn.eps)

    @staticmethod
    def _bn_sm_train(y, bn):
        y = F.batch_norm(y, bn.running_mean, bn.running_var, bn.weight, bn.bias,
                         True, bn.momentum, bn.eps)
        return _ext.softmax_width(y)

    def _conv_bn_sm(self, x, conv, bn):
        y = F.conv2d(x, conv.weight, conv.bias, conv.stride, conv.padding, conv.dilation, conv.groups)
        if bn.training:
            return self._bn_sm_train(y, bn)
        else:
            return self._bn_sm_eval(y, bn)

    def _double(self, x, block):
        s = block.double_conv
        x = self._conv_bn_sm(x, s[0], s[1])
        x = self._conv_bn_sm(x, s[3], s[4])
        return x

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

        dec4 = F.conv_transpose2d(bottleneck, self.upconv4.weight, self.upconv4.bias, stride=2)
        dec4 = self._double(torch.cat((dec4, enc4), dim=1), self.decoder4)
        dec3 = F.conv_transpose2d(dec4, self.upconv3.weight, self.upconv3.bias, stride=2)
        dec3 = self._double(torch.cat((dec3, enc3), dim=1), self.decoder3)
        dec2 = F.conv_transpose2d(dec3, self.upconv2.weight, self.upconv2.bias, stride=2)
        dec2 = self._double(torch.cat((dec2, enc2), dim=1), self.decoder2)
        dec1 = F.conv_transpose2d(dec2, self.upconv1.weight, self.upconv1.bias, stride=2)
        dec1 = self._double(torch.cat((dec1, enc1), dim=1), self.decoder1)
        return F.conv2d(dec1, self.final_conv.weight, self.final_conv.bias, stride=1, padding=0)