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


_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor pooled_columns_cuda(
    torch::Tensor x, int64_t kh, int64_t kw, int64_t ph, int64_t pw);

torch::Tensor sigmoid_sum_cuda(
    torch::Tensor x, torch::Tensor bias);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <cstdint>

__global__ void pooled_columns_3x3_4x4(
    const float* __restrict__ x,
    float* __restrict__ out,
    int C, int H, int W, int OH, int OW,
    int64_t count)
{
    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (index >= count) return;

    const int L = OH * OW;
    const int pos = index % L;
    const int64_t nc = index / L;
    const int oy = pos / OW;
    const int ox = pos % OW;

    const float* src = x + nc * H * W + oy * 4 * W + ox * 4;

    float a00 = 0.f, a01 = 0.f, a02 = 0.f;
    float a10 = 0.f, a11 = 0.f, a12 = 0.f;
    float a20 = 0.f, a21 = 0.f, a22 = 0.f;

    #pragma unroll
    for (int r = 0; r < 6; ++r) {
        const float* row = src + r * W;
        const float v0 = row[0];
        const float v1 = row[1];
        const float v2 = row[2];
        const float v3 = row[3];
        const float v4 = row[4];
        const float v5 = row[5];

        const float q0 = ((v0 + v1) + v2) + v3;
        const float q1 = ((v1 + v2) + v3) + v4;
        const float q2 = ((v2 + v3) + v4) + v5;

        if (r < 4) {
            a00 += q0; a01 += q1; a02 += q2;
        }
        if (r >= 1 && r < 5) {
            a10 += q0; a11 += q1; a12 += q2;
        }
        if (r >= 2) {
            a20 += q0; a21 += q1; a22 += q2;
        }
    }

    float* dst = out + nc * 9 * L + pos;
    dst[0LL * L] = a00 * 0.0625f;
    dst[1LL * L] = a01 * 0.0625f;
    dst[2LL * L] = a02 * 0.0625f;
    dst[3LL * L] = a10 * 0.0625f;
    dst[4LL * L] = a11 * 0.0625f;
    dst[5LL * L] = a12 * 0.0625f;
    dst[6LL * L] = a20 * 0.0625f;
    dst[7LL * L] = a21 * 0.0625f;
    dst[8LL * L] = a22 * 0.0625f;
}

__global__ void pooled_columns_generic(
    const float* __restrict__ x,
    float* __restrict__ out,
    int C, int H, int W, int OH, int OW,
    int KH, int KW, int PH, int PW,
    int64_t count)
{
    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (index >= count) return;

    const int L = OH * OW;
    const int pos = index % L;
    const int64_t column = index / L;
    const int kernel_offset = column % (KH * KW);
    const int64_t nc = column / (KH * KW);

    const int iy = (pos / OW) * PH + kernel_offset / KW;
    const int ix = (pos % OW) * PW + kernel_offset % KW;
    const float* src = x + nc * H * W + iy * W + ix;

    float sum = 0.f;
    for (int py = 0; py < PH; ++py) {
        for (int px = 0; px < PW; ++px) {
            sum += src[py * W + px];
        }
    }
    out[index] = sum / static_cast<float>(PH * PW);
}

__device__ __forceinline__ float warp_sum(float value)
{
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        value += __shfl_down_sync(0xffffffff, value, offset);
    }
    return value;
}

__device__ __forceinline__ float block_sum(float value)
{
    __shared__ float warp_values[8];
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    value = warp_sum(value);
    if (lane == 0) warp_values[warp] = value;
    __syncthreads();

    value = (threadIdx.x < 8) ? warp_values[lane] : 0.f;
    if (warp == 0) value = warp_sum(value);
    return value;
}

__global__ void sigmoid_sum_partials(
    const float* __restrict__ x,
    const float* __restrict__ bias,
    float* __restrict__ partials,
    int spatial, int elements, int chunks)
{
    constexpr int CHUNK = 4096;
    const int chunk = blockIdx.x;
    const int batch = blockIdx.y;
    const int start = chunk * CHUNK;
    const int end = min(start + CHUNK, elements);
    const float* src = x + static_cast<int64_t>(batch) * elements;

    float sum = 0.f;
    for (int i = start + threadIdx.x; i < end; i += blockDim.x) {
        const float value = src[i] + bias[i / spatial];
        sum += 1.f / (1.f + __expf(-value));
    }

    sum = block_sum(sum);
    if (threadIdx.x == 0) {
        partials[static_cast<int64_t>(batch) * chunks + chunk] = sum;
    }
}

__global__ void finish_sum(
    const float* __restrict__ partials,
    float* __restrict__ out,
    int chunks)
{
    const int batch = blockIdx.x;
    const float* src = partials + static_cast<int64_t>(batch) * chunks;

    float sum = 0.f;
    for (int i = threadIdx.x; i < chunks; i += blockDim.x) {
        sum += src[i];
    }

    sum = block_sum(sum);
    if (threadIdx.x == 0) out[batch] = sum;
}

torch::Tensor pooled_columns_cuda(
    torch::Tensor x, int64_t kh, int64_t kw, int64_t ph, int64_t pw)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 4 && x.is_contiguous(),
                "Expected contiguous NCHW input");
    TORCH_CHECK(kh > 0 && kw > 0 && ph > 0 && pw > 0,
                "Kernel sizes must be positive");

    const c10::cuda::CUDAGuard guard(x.device());
    const int64_t N = x.size(0);
    const int C = x.size(1);
    const int H = x.size(2);
    const int W = x.size(3);

    TORCH_CHECK(H >= kh + ph - 1 && W >= kw + pw - 1,
                "Input is too small for convolution and pooling");

    const int OH = (H - kh + 1) / ph;
    const int OW = (W - kw + 1) / pw;
    const int64_t L = static_cast<int64_t>(OH) * OW;

    auto out = torch::empty({N, C * kh * kw, L}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    if (N == 0) return out;

    if (kh == 3 && kw == 3 && ph == 4 && pw == 4) {
        const int64_t count = N * C * L;
        pooled_columns_3x3_4x4<<<(count + 255) / 256, 256, 0, stream>>>(
            x.data_ptr<float>(), out.data_ptr<float>(),
            C, H, W, OH, OW, count);
    } else {
        const int64_t count = out.numel();
        pooled_columns_generic<<<(count + 255) / 256, 256, 0, stream>>>(
            x.data_ptr<float>(), out.data_ptr<float>(),
            C, H, W, OH, OW, kh, kw, ph, pw, count);
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return out;
}

torch::Tensor sigmoid_sum_cuda(torch::Tensor x, torch::Tensor bias)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat &&
                x.dim() == 3 && x.is_contiguous(),
                "Expected contiguous CUDA float32 input");
    TORCH_CHECK(bias.device() == x.device() &&
                bias.scalar_type() == at::kFloat &&
                bias.is_contiguous() && bias.numel() == x.size(1),
                "Invalid bias");

    const c10::cuda::CUDAGuard guard(x.device());
    const int N = x.size(0);
    const int spatial = x.size(2);
    const int elements = x.size(1) * spatial;
    const int chunks = (elements + 4095) / 4096;

    auto out = torch::empty({N}, x.options());
    if (N == 0) return out;

    auto partials = torch::empty({N, chunks}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    sigmoid_sum_partials<<<dim3(chunks, N), 256, 0, stream>>>(
        x.data_ptr<float>(), bias.data_ptr<float>(),
        partials.data_ptr<float>(), spatial, elements, chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_sum<<<N, 256, 0, stream>>>(
        partials.data_ptr<float>(), out.data_ptr<float>(), chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return out;
}
"""

_ops = load_inline(
    name="pooled_conv_sigmoid_sum_cuda_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["pooled_columns_cuda", "sigmoid_sum_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self, in_channels, out_channels, kernel_size, pool_kernel_size):
        super().__init__()
        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size)
        self.avg_pool = nn.AvgPool2d(pool_kernel_size)

        self._kernel = self.conv.kernel_size
        self._pool = (
            (pool_kernel_size, pool_kernel_size)
            if isinstance(pool_kernel_size, int)
            else tuple(pool_kernel_size)
        )

    def forward(self, x):
        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or self.conv.weight.dtype != torch.float32
            or torch.is_autocast_enabled()
        ):
            return torch.sigmoid(self.avg_pool(self.conv(x))).sum((1, 2, 3))

        kh, kw = self._kernel
        ph, pw = self._pool

        # Pooling is linear: average input patches before the weight product.
        columns = _ops.pooled_columns_cuda(
            x.contiguous(), kh, kw, ph, pw
        )
        weights = self.conv.weight.reshape(self.conv.out_channels, -1)
        products = torch.bmm(
            weights.unsqueeze(0).expand(x.shape[0], -1, -1),
            columns,
        )
        return _ops.sigmoid_sum_cuda(products, self.conv.bias)
```