Write a custom CUDA kernel to replace PyTorch's InstanceNorm + ReLU implementation for CNN layers.

You are given the following PyTorch architecture:

python
import torch
import torch.nn as nn

class Model(nn.Module):
"""
Simple model that performs InstanceNorm + ReLU.
"""
def init(self, num_features=64, eps=1e-5, affine=True, track_running_stats=False):
super(Model, self).init()
self.instance_norm = nn.InstanceNorm2d(
    num_features=num_features,
    eps=eps,
    affine=affine,
    track_running_stats=track_running_stats
)
self.relu = nn.ReLU(inplace=True)

def forward(self, x: torch.Tensor) -> torch.Tensor:  
    """  
    Applies InstanceNorm and then ReLU to the input tensor.  

    Args:  
        x (torch.Tensor): Input tensor of shape [B, C, H, W].  

    Returns:  
        torch.Tensor: ReLU(InstanceNorm(x)), same shape as input.  
    """  
    x_normalized = self.instance_norm(x)
    return self.relu(x_normalized)

batch_size = 256
num_features = 64
height = 128
width = 128

def get_inputs():
x = torch.randn(batch_size, num_features, height, width)
return [x]

def get_init_inputs():
return [num_features]


Your task is to optimize this InstanceNorm + ReLU implementation by:

1. **Operator Fusion**: Combine InstanceNorm computation (mean, variance, normalization) and ReLU activation into a single CUDA kernel to eliminate intermediate tensor storage and reduce memory bandwidth overhead.

2. **Memory Access Optimization**: Minimize global memory access by keeping intermediate computations in registers, and ensure coalesced memory access patterns for the [B, C, H, W] tensor layout.

3. **Shared Memory Optimization**: Use shared memory for efficient parallel reduction when computing mean and variance across spatial dimensions (H*W) within each instance and channel.

4. **Activation Fusion**: Apply ReLU activation (max(0, x)) directly within the kernel after normalization to avoid separate activation computation.

5. **Thread Configuration**: Use optimal block size (e.g., 256 threads) and compute grid dimensions based on batch_size and num_features to maximize GPU utilization.

6. **Numerical Stability**: Ensure proper epsilon handling in InstanceNorm computation to avoid division by zero and maintain numerical precision.

The optimized CUDA kernel should:
- Take input tensor x, weight, and bias as input (all float32)
- Compute InstanceNorm (mean, variance, normalization) and apply ReLU in a single kernel
- Apply ReLU activation (max(0, x)) directly within the kernel
- Use shared memory for efficient mean and variance computation
- Maintain numerical stability with proper epsilon handling
- Achieve significant speedup over PyTorch's separate InstanceNorm + ReLU implementation
- Support both affine and non-affine modes

Follow the inline CUDA extension syntax example provided in reference. The kernel should be optimized for GPU architectures and demonstrate performance improvements through reduced memory access, fused computation, and efficient parallel reduction.
