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_forward_cuda(torch::Tensor x, torch::Tensor weight, torch::Tensor bias, torch::Tensor mean, torch::Tensor var, double eps);
torch::Tensor bn_softmax_forward(torch::Tensor x, torch::Tensor weight, torch::Tensor bias, torch::Tensor mean, torch::Tensor var, double eps) { return bn_softmax_forward_cuda(x, weight, bias, mean, var, eps); }
'''
_cuda_src = r'''
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>
__inline__ __device__ float warp_reduce_max(float val){for(int o=16;o>0;o>>=1) val=fmaxf(val,__shfl_down_sync(0xffffffff,val,o)); return val;}
__inline__ __device__ float warp_reduce_sum(float val){for(int o=16;o>0;o>>=1) val+=__shfl_down_sync(0xffffffff,val,o); return val;}
__inline__ __device__ float block_reduce_max(float val){__shared__ float sh[32]; int lane=threadIdx.x&31,wid=threadIdx.x>>5; val=warp_reduce_max(val); if(lane==0) sh[wid]=val; __syncthreads(); float b=(threadIdx.x<(blockDim.x+31)/32)?sh[lane]:-INFINITY; if(wid==0) b=warp_reduce_max(b); if(threadIdx.x==0) sh[0]=b; __syncthreads(); return sh[0];}
__inline__ __device__ float block_reduce_sum(float val){__shared__ float sh[32]; int lane=threadIdx.x&31,wid=threadIdx.x>>5; val=warp_reduce_sum(val); if(lane==0) sh[wid]=val; __syncthreads(); float b=(threadIdx.x<(blockDim.x+31)/32)?sh[lane]:0.0f; if(wid==0) b=warp_reduce_sum(b); if(threadIdx.x==0) sh[0]=b; __syncthreads(); return sh[0];}
__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 c=(row/H)%C; int base=row*W; float scale=gamma[c]*rsqrtf(var[c]+eps); float shift=beta[c]-mean[c]*scale; float m=-INFINITY;
 for(int w=threadIdx.x; w<W; w+=blockDim.x){float v=x[base+w]*scale+shift; m=fmaxf(m,v);} m=block_reduce_max(m); float s=0.0f;
 for(int w=threadIdx.x; w<W; w+=blockDim.x){float v=expf((x[base+w]*scale+shift)-m); s+=v;} s=block_reduce_sum(s); float inv=1.0f/s;
 for(int w=threadIdx.x; w<W; w+=blockDim.x){out[base+w]=expf((x[base+w]*scale+shift)-m)*inv;}
}
torch::Tensor bn_softmax_forward_cuda(torch::Tensor x, torch::Tensor weight, torch::Tensor bias, torch::Tensor mean, torch::Tensor var, double eps){auto y=torch::empty_like(x); int N=x.size(0),C=x.size(1),H=x.size(2),W=x.size(3); int rows=N*C*H; bn_softmax_kernel<<<rows,256,0,at::cuda::getCurrentCUDAStream()>>>(x.data_ptr<float>(),weight.data_ptr<float>(),bias.data_ptr<float>(),mean.data_ptr<float>(),var.data_ptr<float>(),y.data_ptr<float>(),rows,C,H,W,(float)eps); return y;}
'''
_ext = load_inline(name='unet_bnsm_ext_final', cpp_sources=_cpp_src, cuda_sources=_cuda_src, functions=['bn_softmax_forward'], with_cuda=True, extra_cuda_cflags=['-O3'], 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)
        )

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):
        if self.training:
            return torch.softmax(bn(x), dim=-1)
        return _ext.bn_softmax_forward(x.contiguous(), bn.weight, bn.bias, bn.running_mean, bn.running_var, bn.eps)

    def _dc(self, x, m):
        s = m.double_conv
        x = F.conv2d(x, s[0].weight, s[0].bias, stride=1, padding=1)
        x = self._bn_sm(x, s[1])
        x = F.conv2d(x, s[3].weight, s[3].bias, stride=1, padding=1)
        x = self._bn_sm(x, s[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)
        dec4 = F.conv_transpose2d(bottleneck, self.upconv4.weight, self.upconv4.bias, stride=2)
        dec4 = self._dc(torch.cat((dec4, enc4), dim=1), self.decoder4)
        dec3 = F.conv_transpose2d(dec4, self.upconv3.weight, self.upconv3.bias, stride=2)
        dec3 = self._dc(torch.cat((dec3, enc3), dim=1), self.decoder3)
        dec2 = F.conv_transpose2d(dec3, self.upconv2.weight, self.upconv2.bias, stride=2)
        dec2 = self._dc(torch.cat((dec2, enc2), dim=1), self.decoder2)
        dec1 = F.conv_transpose2d(dec2, self.upconv1.weight, self.upconv1.bias, stride=2)
        dec1 = self._dc(torch.cat((dec1, enc1), dim=1), self.decoder1)
        return F.conv2d(dec1, self.final_conv.weight, self.final_conv.bias)