```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_unfold_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_unfold_3x3_4x4(
    const float* __restrict__ x,
    float* __restrict__ out,
    int channels, int height, int width,
    int oh, int ow, int64_t count)
{
    const int spatial = oh * ow;
    const int64_t idx =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (idx >= count) return;

    const int s = idx % spatial;
    const int64_t nc = idx / spatial;
    const int n = nc / channels;
    const int c = nc % channels;
    const int y = (s / ow) * 4;
    const int z = (s % ow) * 4;

    const float* src =
        x + (nc * height + y) * width + z;

    float sums[9] = {};

    #pragma unroll
    for (int r = 0; r < 6; ++r) {
        const float* row = src + r * width;
        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 h0 = ((v0 + v1) + v2) + v3;
        const float h1 = ((v1 + v2) + v3) + v4;
        const float h2 = ((v2 + v3) + v4) + v5;

        #pragma unroll
        for (int ky = 0; ky < 3; ++ky) {
            if (r >= ky && r < ky + 4) {
                sums[ky * 3] += h0;
                sums[ky * 3 + 1] += h1;
                sums[ky * 3 + 2] += h2;
            }
        }
    }

    const int64_t base =
        (static_cast<int64_t>(n) * channels * 9 + c * 9) * spatial + s;

    #pragma unroll
    for (int k = 0; k < 9; ++k) {
        out[base + static_cast<int64_t>(k) * spatial] =
            sums[k] * 0.0625f;
    }
}

__global__ void pooled_unfold_generic(
    const float* __restrict__ x,
    float* __restrict__ out,
    int channels, int height, int width,
    int kh, int kw, int ph, int pw,
    int oh, int ow, int64_t count)
{
    const int64_t idx =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (idx >= count) return;

    const int spatial = oh * ow;
    const int s = idx % spatial;
    int64_t t = idx / spatial;
    const int kx = t % kw;
    t /= kw;
    const int ky = t % kh;
    t /= kh;
    const int c = t % channels;
    const int n = t / channels;

    const int y = (s / ow) * ph + ky;
    const int z = (s % ow) * pw + kx;
    const float* src =
        x + ((static_cast<int64_t>(n) * channels + c) * height + y)
                * width + z;

    float sum = 0.0f;
    for (int py = 0; py < ph; ++py) {
        for (int px = 0; px < pw; ++px) {
            sum += src[py * width + px];
        }
    }
    out[idx] = sum / static_cast<float>(ph * pw);
}

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

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

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

    value = threadIdx.x < 8 ? warp_sums[lane] : 0.0f;
    if (warp == 0) {
        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1) {
            value += __shfl_down_sync(0xffffffff, value, offset);
        }
    }
    return value;
}

__global__ void sigmoid_partial_sum(
    const float* __restrict__ x,
    const float* __restrict__ bias,
    float* __restrict__ partial,
    int spatial, int64_t per_batch, int chunks)
{
    const int n = blockIdx.y;
    const int chunk = blockIdx.x;
    const int64_t begin = static_cast<int64_t>(chunk) * 4096;
    const int64_t base = static_cast<int64_t>(n) * per_batch;

    float sum = 0.0f;
    #pragma unroll
    for (int j = 0; j < 16; ++j) {
        const int64_t i = begin + threadIdx.x + j * 256;
        if (i < per_batch) {
            const float value = x[base + i] + bias[i / spatial];
            sum += 1.0f / (1.0f + expf(-value));
        }
    }

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

__global__ void finish_sum(
    const float* __restrict__ partial,
    float* __restrict__ out,
    int chunks)
{
    float sum = 0.0f;
    for (int i = threadIdx.x; i < chunks; i += 256) {
        sum += partial[static_cast<int64_t>(blockIdx.x) * chunks + i];
    }
    sum = block_sum(sum);
    if (threadIdx.x == 0) out[blockIdx.x] = sum;
}

torch::Tensor pooled_unfold_cuda(
    torch::Tensor x, int64_t kh, int64_t kw,
    int64_t ph, int64_t pw)
{
    const c10::cuda::CUDAGuard guard(x.device());

    const int n = x.size(0);
    const int c = x.size(1);
    const int h = x.size(2);
    const int w = x.size(3);
    const int oh = (h - kh + 1) / ph;
    const int ow = (w - kw + 1) / pw;

    TORCH_CHECK(oh > 0 && ow > 0, "Invalid convolution/pooling dimensions");

    auto out = torch::empty({n, c * kh * kw, oh * ow}, x.options());
    auto stream = at::cuda::getCurrentCUDAStream();

    if (kh == 3 && kw == 3 && ph == 4 && pw == 4) {
        const int64_t count = static_cast<int64_t>(n) * c * oh * ow;
        pooled_unfold_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_unfold_generic<<<(count + 255) / 256, 256, 0, stream>>>(
            x.data_ptr<float>(), out.data_ptr<float>(),
            c, h, w, kh, kw, ph, pw, oh, ow, count);
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return out;
}

torch::Tensor sigmoid_sum_cuda(
    torch::Tensor x, torch::Tensor bias)
{
    const c10::cuda::CUDAGuard guard(x.device());

    const int n = x.size(0);
    const int spatial = x.size(2);
    const int64_t per_batch = x.size(1) * x.size(2);
    const int chunks = (per_batch + 4095) / 4096;

    auto partial = torch::empty({n, chunks}, x.options());
    auto out = torch::empty({n}, x.options());
    auto stream = at::cuda::getCurrentCUDAStream();

    sigmoid_partial_sum<<<dim3(chunks, n), 256, 0, stream>>>(
        x.data_ptr<float>(), bias.data_ptr<float>(),
        partial.data_ptr<float>(), spatial, per_batch, chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_sum<<<n, 256, 0, stream>>>(
        partial.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_unfold_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()
            or x.shape[0] == 0
        ):
            return torch.sigmoid(self.avg_pool(self.conv(x))).sum((1, 2, 3))

        # Average pooling commutes with the linear convolution:
        # pool each input patch first, then multiply by convolution weights.
        columns = _ops.pooled_unfold_cuda(
            x.contiguous(), *self._kernel, *self._pool
        )
        weights = self.conv.weight.reshape(self.conv.out_channels, -1)
        values = torch.matmul(weights.unsqueeze(0), columns)
        return _ops.sigmoid_sum_cuda(values, self.conv.bias)
```