239. Sliding Window Maximum

本文介绍了一种高效求解滑动窗口最大值问题的方法,使用优先队列和双向队列实现,时间复杂度分别为O(NlogK)和O(N),适用于处理大规模数据集。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

转载作者:ethammmli

转载地址:https://segmentfault.com/a/1190000003903509

注意时间复杂度为:线性

Sliding Window Maximum

Given an array nums, there is a sliding window of size k which is moving from the very left of the array to the very right. You can only see the k numbers in the window. Each time the sliding window moves right by one position.

For example, Given nums = [1,3,-1,-3,5,3,6,7], and k = 3.

Window position                Max
---------------               -----
[1  3  -1] -3  5  3  6  7       3
 1 [3  -1  -3] 5  3  6  7       3
 1  3 [-1  -3  5] 3  6  7       5
 1  3  -1 [-3  5  3] 6  7       5
 1  3  -1  -3 [5  3  6] 7       6
 1  3  -1  -3  5 [3  6  7]      7 

Therefore, return the max sliding window as [3,3,5,5,6,7].

Note: You may assume k is always valid, ie: 1 ≤ k ≤ input array's size for non-empty array.

Follow up: Could you solve it in linear time?

Hint:

How about using a data structure such as deque (double-ended queue)? The queue size need not be the same as the window’s size. Remove redundant elements and the queue should store only elements that need to be considered.

优先队列

复杂度

时间 O(NlogK) 空间 O(K)

思路

维护一个大小为K的最大堆,依此维护一个大小为K的窗口,每次读入一个新数,都把堆中窗口最左边的数扔掉,再把新数加入堆中,这样堆顶就是这个窗口内最大的值。

注意

-结果数组的大小是nums.length + 1 - k, 赋值时下标也是i + 1 - k

代码

public class Solution {
    public int[] maxSlidingWindow(int[] nums, int k) {
        if(nums == null || nums.length == 0) return new int[0];
        PriorityQueue<Integer> pq = new PriorityQueue<Integer>(Collections.reverseOrder());
        int[] res = new int[nums.length + 1 - k];
        for(int i = 0; i < nums.length; i++){
            // 把窗口最左边的数去掉
            if(i >= k) pq.remove(nums[i - k]);
            // 把新的数加入窗口的堆中
            pq.offer(nums[i]);
            // 堆顶就是窗口的最大值
            if(i + 1 >= k) res[i + 1 - k] = pq.peek();
        }
        return res;
    }
}

双向队列

复杂度

时间 O(N) 空间 O(K)

思路

我们用双向队列可以在O(N)时间内解决这题。当我们遇到新的数时,将新的数和双向队列的末尾比较,如果末尾比新数小,则把末尾扔掉,直到该队列的末尾比新数大或者队列为空的时候才住手。这样,我们可以保证队列里的元素是从头到尾降序的,由于队列里只有窗口内的数,所以他们其实就是窗口内第一大,第二大,第三大...的数。保持队列里只有窗口内数的方法和上个解法一样,也是每来一个新的把窗口最左边的扔掉,然后把新的加进去。然而由于我们在加新数的时候,已经把很多没用的数给扔了,这样队列头部的数并不一定是窗口最左边的数。这里的技巧是,我们队列中存的是那个数在原数组中的下标,这样我们既可以直到这个数的值,也可以知道该数是不是窗口最左边的数。这里为什么时间复杂度是O(N)呢?因为每个数只可能被操作最多两次,一次是加入队列的时候,一次是因为有别的更大数在后面,所以被扔掉,或者因为出了窗口而被扔掉。

代码

public class Solution {
    public int[] maxSlidingWindow(int[] nums, int k) {
        if(nums == null || nums.length == 0) return new int[0];
        LinkedList<Integer> deque = new LinkedList<Integer>();
        int[] res = new int[nums.length + 1 - k];
        for(int i = 0; i < nums.length; i++){
            // 每当新数进来时,如果发现队列头部的数的下标,是窗口最左边数的下标,则扔掉
            if(!deque.isEmpty() && deque.peekFirst() == i - k) deque.poll();
            // 把队列尾部所有比新数小的都扔掉,保证队列是降序的
            while(!deque.isEmpty() && nums[deque.peekLast()] < nums[i]) deque.removeLast();
            // 加入新数
            deque.offerLast(i);
            // 队列头部就是该窗口内第一大的
            if((i + 1) >= k) res[i + 1 - k] = nums[deque.peek()];
        }
        return res;
    }
}

内容概要:《中文大模型基准测评2025年上半年报告》由SuperCLUE团队发布,详细评估了2025年上半年中文大模型的发展状况。报告涵盖了大模型的关键进展、国内外大模型全景图及差距、专项测评基准介绍等。通过SuperCLUE基准,对45个国内外代表性大模型进行了六大任务(数学推理、科学推理、代码生成、智能体Agent、精确指令遵循、幻觉控制)的综合测评。结果显示,海外模型如o3、o4-mini(high)在推理任务上表现突出,而国内模型如Doubao-Seed-1.6-thinking-250715在智能体Agent和幻觉控制任务上表现出色。此外,报告还分析了模型性价比、效能区间分布,并对代表性模型如Doubao-Seed-1.6-thinking-250715、DeepSeek-R1-0528、GLM-4.5等进行了详细介绍。整体来看,国内大模型在特定任务上已接近国际顶尖水平,但在综合推理能力上仍有提升空间。 适用人群:对大模型技术感兴趣的科研人员、工程师、产品经理及投资者。 使用场景及目标:①了解2025年上半年中文大模型的发展现状与趋势;②评估国内外大模型在不同任务上的表现差异;③为技术选型和性能优化提供参考依据。 其他说明:报告提供了详细的测评方法、评分标准及结果分析,确保评估的科学性和公正性。此外,SuperCLUE团队还发布了多个专项测评基准,涵盖多模态、文本、推理等多个领域,为业界提供全面的测评服务。
<think>好的,我现在需要检查并修正用户提供的高斯滤波函数。首先,我得仔细阅读代码,理解每个部分的功能,然后找出可能的问题。 首先看函数定义,函数名为cupy_gaussian_filter,参数是data和window_size,返回np.ndarray。看起来是用CuPy在GPU上加速的高斯滤波实现。接下来分步骤分析: 1. 数据标准化部分:将数据归一化到0-1范围,防止除零加了1e-8,这里应该没问题。 2. 创建滑动窗口视图:使用sliding_window_view生成窗口视图,这一步可能是为了后续向量化操作。但要注意,滑动窗口的形状是否正确,比如window_size是否为1D,或者是否需要处理二维情况。但用户的问题可能是在一维数据上的滤波。 3. 计算方差:用cp.var计算每个窗口的方差,axis=1可能有问题。假设输入数据是一维的,滑动窗口后,win_view的形状应该是(n_windows, window_size),所以方差应该在axis=1上计算,这没问题。但如果是多维数据,可能需要调整轴。不过根据函数名和上下文,可能假设是一维的。 接下来,处理零方差的部分被注释掉了,用户可能暂时不需要,或者可能后面sigma为零导致问题。但注释掉的话,如果方差为零,计算sigma时会得到零,可能导致高斯核除零错误。 4. 计算标准差:得到sigmas,这里没问题,但需要考虑sigma为零的情况。 然后进入并行高斯滤波部分,使用CuPy的ElementwiseKernel。这里的问题可能较大。用户提供的核函数内部逻辑是空的,需要实现高斯计算。当前代码中,gauss_kernel的参数是x和sigma,但实际传入的是win_view和sigmas[:, None]。这里可能存在维度不匹配的问题,因为win_view的形状可能是(n_windows, window_size),而sigmas[:, None]的形状是(n_windows, 1)。ElementwiseKernel是按元素处理的,但这里的每个窗口需要处理整个窗口的数据,并应用高斯权重。 这里可能误解了ElementwiseKernel的用法。ElementwiseKernel适用于每个元素独立计算,而高斯滤波需要对每个窗口内的元素进行加权求和,这可能需要更复杂的核函数,比如每个线程处理一个窗口,或者使用二维索引。当前的设置可能无法正确计算每个窗口的高斯加权。 另外,高斯核函数的实现逻辑缺失,用户注释里需要在此处实现,但实际代码中没有。需要补全高斯计算的逻辑,比如根据每个窗口的sigma,计算每个位置x的权重,然后加权平均。 接下来,计算中值并反标准化。这里使用cp.median取高斯加权后的结果的中值,这似乎不符合高斯滤波的常规做法。通常高斯滤波是加权平均,而这里用了中值,可能用户想实现的是自适应高斯加权中值滤波?或者可能混淆了高斯滤波和中值滤波的概念? 另外,处理边界填充的部分,代码中使用concatenate来拼接原始数据的前pad部分和final以及后面的部分。但final的长度可能为len(data) - window_size + 1,所以填充方式是否正确?比如,假设原数据长度为N,窗口大小为w,pad为w//2。则滑动窗口视图生成的数据长度是N - w + 1。所以最终medians的长度也是N - w + 1。为了保持输出长度与原数据相同,需要在前后填充pad个样本。但当前的拼接方式可能有误,比如data[:pad]可能正确,但后半部分data[len(final)+pad - len(data):]可能需要更准确的索引计算。 总结问题点: 1. ElementwiseKernel的使用不正确,无法实现窗口内的高斯加权计算。 2. 高斯核函数内部逻辑未实现,导致计算结果错误。 3. 使用中值而非加权平均,可能不符合高斯滤波的需求。 4. 边界填充部分可能有错误,导致输出尺寸不正确。 5. 未处理sigma为零的情况,导致可能的除零错误。 6. 数据标准化后的还原是否正确? 接下来需要逐步修正这些问题: 首先,高斯滤波的核心是对每个窗口内的数据应用高斯权重,然后进行加权平均。而当前代码中使用ElementwiseKernel可能无法正确实现这一点,因为每个元素的操作需要访问整个窗口的数据。可能需要重新设计核函数,例如使用每个线程处理一个窗口,计算该窗口的高斯权重,然后求加权平均。 或者,可以考虑向量化的方式计算高斯权重矩阵,然后进行矩阵运算。例如,为每个窗口生成对应的高斯核,然后与窗口数据点乘求和。但这可能需要更多的内存,尤其是窗口较大时。 另一个思路是,对每个窗口计算其中心点的高斯权重,但可能用户是想根据窗口的方差自适应调整高斯核的sigma,实现自适应滤波。例如,每个窗口的sigma由该窗口的方差决定,然后根据该sigma生成对应的高斯权重,再对窗口内的数据进行加权平均。 当前的代码中,gauss_kernel的参数是x和sigma,其中x可能是窗口内的每个元素,sigma是该窗口对应的标准差。然后对每个元素应用高斯计算,但如何将这些计算结果聚合起来得到每个窗口的加权平均或中值? 可能原来的意图是对窗口内的每个元素计算高斯权重,然后将权重与数据相乘,最后求和或取中值。但当前的实现可能存在问题,因为每个元素单独处理,没有聚合步骤。例如,gauss_results可能是每个窗口内元素的高斯权重值,然后取中值,这不符合高斯滤波的要求。正确的高斯滤波应该是加权平均,而中值是另一种滤波方法,两者不应混淆。 因此,需要将median改为加权平均。或者,如果用户确实需要高斯加权中值,那可能需要不同的处理,但通常高斯滤波是加权平均。 修正步骤: 1. 高斯核函数的实现:每个窗口的每个元素的高斯权重计算。高斯核的公式是:$w_i = \exp(- (x_i - \mu)^2 / (2 \sigma^2))$,其中$\mu$是窗口的中心值,或者这里可能假设窗口的均值?或者用户可能将窗口的位置作为x?需要明确。 假设每个窗口的数据已经标准化到0-1,对于每个窗口,计算其高斯权重,其中中心点作为均值。例如,窗口的中心位置是窗口的中间点,但如果是数据的值,可能需要不同的处理。或者,这里的x是窗口内的数据点,而高斯核的均值是窗口的均值,sigma是该窗口的标准差。这样,每个数据点的权重是基于该点与窗口均值的差异,由sigma决定。 例如,对于一个窗口内的数据点x_i,权重为exp(- (x_i - mu)^2 / (2 sigma^2)),其中mu是窗口的均值,sigma是窗口的标准差。然后加权平均就是sum(w_i * x_i) / sum(w_i)。 但当前代码中,计算的是方差,所以方差是sigma平方。每个窗口的mu可以用窗口的均值来计算,即mu = cp.mean(win_view, axis=1)。但目前代码中没有计算mu,所以需要添加这个步骤。 所以,修正步骤可能包括: - 计算每个窗口的均值mu。 - 在核函数中,根据mu和sigma计算每个元素的高斯权重。 - 对每个窗口的权重和数据进行加权平均。 但当前的代码结构中使用ElementwiseKernel可能不够,因为需要每个窗口处理多个元素,并计算加权和。或许应该改用每个线程处理一个窗口,使用RawKernel而不是ElementwiseKernel。 或者,可以考虑将窗口数据、mu、sigma作为输入,然后在核函数中为每个窗口计算权重并求和。 这可能需要使用二维网格,其中每个线程块处理一个窗口,或者每个线程处理一个窗口。例如,定义一个核函数,每个线程计算一个窗口的加权平均。 另外,当前代码中,win_view的形状是(n_windows, window_size),假设输入数据是一维的。对于每个窗口i,有window_size个元素,需要计算对应的mu和sigma,然后对每个元素计算权重,求和得到加权平均。 因此,可能需要一个核函数,每个线程处理一个窗口,计算该窗口的高斯加权平均。 修改后的步骤可能包括: - 计算每个窗口的mu(均值)和sigma(标准差)。 - 对于每个窗口,生成高斯权重,然后计算加权平均。 或者,可以预先生成一个高斯权重矩阵,然后进行矩阵乘法。但由于每个窗口的sigma不同,权重矩阵会是三维的,这可能难以向量化。 另一个方法是使用CuPy的广播功能,为每个窗口计算对应的权重,然后求和。 例如: mu = cp.mean(win_view, axis=1, keepdims=True) sigma = cp.sqrt(variances)[:, None] weights = cp.exp(-0.5 * ((win_view - mu) / sigma)**2) / (sigma * cp.sqrt(2 * cp.pi)) weighted_sum = cp.sum(weights * win_view, axis=1) sum_weights = cp.sum(weights, axis=1) result = weighted_sum / sum_weights 这可以在不使用核函数的情况下,完全向量化实现。但可能内存消耗较大,因为需要存储weights数组,形状与win_view相同。 如果window_size很大,或者数据量很大,这可能导致内存不足。此时,使用自定义核函数可能更高效。 回到用户代码,他们可能想通过核函数来并行处理每个窗口,以减少内存占用。 因此,需要重新设计核函数部分。例如,使用RawKernel,每个线程处理一个窗口: 代码的大致思路: 计算每个窗口的mu和sigma。 然后,在核函数中,每个线程读取一个窗口的所有元素,计算高斯权重,然后加权平均。 但CuPy的RawKernel需要明确指定线程和块的布局。例如,总共有n_windows个窗口,每个线程处理一个窗口。每个线程需要遍历窗口内的所有元素,计算加权和。 示例代码: 计算mu和sigma: mu = cp.mean(win_view, axis=1) sigma = cp.sqrt(variances) 然后定义核函数: code = ''' extern "C" __global__ void gaussian_filter(const float* windows, const float* mu, const float* sigma, float* output, int window_size) { int tid = blockIdx.x * blockDim.x + threadIdx.x; if (tid >= {n_windows}) return; float sum = 0.0; float weight_sum = 0.0; float current_mu = mu[tid]; float current_sigma = sigma[tid]; for (int i = 0; i < window_size; i++) { float x = windows[tid * window_size + i]; float exponent = -0.5 * pow((x - current_mu) / current_sigma, 2); float weight = exp(exponent) / (current_sigma * sqrt(2 * M_PI)); sum += x * weight; weight_sum += weight; } output[tid] = sum / weight_sum; } ''' 但需要动态替换{n_windows}和window_size的值。不过CuPy的RawKernel可能不支持直接传入变量到代码中,可能需要使用格式化字符串。 然后,编译核函数并执行: n_windows = win_view.shape[0] window_size = win_view.shape[1] block_size = 256 grid_size = (n_windows + block_size - 1) // block_size kernel = cp.RawKernel(code, 'gaussian_filter') output_gpu = cp.empty(n_windows, dtype=cp.float32) kernel((grid_size,), (block_size,), (win_view, mu, sigma, output_gpu, window_size)) 然后,output_gpu即为每个窗口的加权平均结果。 但用户原来的代码中使用的是ElementwiseKernel,这可能不适合,因为每个元素操作无法访问整个窗口的数据。因此,必须改用更合适的核函数类型。 回到用户代码,原来的gauss_kernel的参数是x(可能是窗口中的每个元素)和sigma(每个窗口的sigma),但这样每个元素只能处理自己的数据,无法聚合整个窗口的信息。因此,原来的方法不正确,需要改为每个窗口一个线程,处理整个窗口的数据。 所以,原来的代码中,gauss_results的计算是错误的,因为无法通过ElementwiseKernel实现窗口内的聚合操作。 此外,用户代码中使用medians = cp.median(gauss_results, axis=1),这里假设gauss_results是每个窗口的高斯处理后的结果,但原来的意图可能不是这样。如果高斯滤波是对每个窗口进行加权平均,那么medians这一步应该是多余的,或者用户可能混淆了步骤。 现在总结问题: 1. 核函数未正确实现高斯加权逻辑,且ElementwiseKernel无法处理窗口级别的聚合。 2. 使用median代替加权平均,不符合高斯滤波定义。 3. 边界填充可能有误。 4. 未处理sigma为零的情况。 修正步骤: 1. 计算每个窗口的mu和sigma。 2. 使用自定义核函数(如RawKernel)对每个窗口计算加权平均。 3. 替换median为加权平均的结果。 4. 处理sigma为零的情况,例如设置最小sigma阈值。 5. 修正边界填充逻辑,确保输出尺寸与输入一致。 具体代码修改: 首先,计算mu和sigma: variances = cp.var(win_view, axis=1) # 处理零方差 variances = cp.maximum(variances, 1e-10) # 避免sigma为零 sigmas = cp.sqrt(variances) mu = cp.mean(win_view, axis=1) 然后,定义核函数: kernel_code = ''' extern "C" __global__ void gaussian_weighted_average(const float* windows, const float* mu, const float* sigma, float* output, int window_size) { int tid = blockIdx.x * blockDim.x + threadIdx.x; int total_windows = gridDim.x * blockDim.x; if (tid >= total_windows) return; float current_mu = mu[tid]; float current_sigma = sigma[tid]; float sum = 0.0; float weight_sum = 0.0; for (int i = 0; i < window_size; ++i) { float x = windows[tid * window_size + i]; float diff = x - current_mu; float exponent = - (diff * diff) / (2.0f * current_sigma * current_sigma); float weight = exp(exponent); sum += x * weight; weight_sum += weight; } if (weight_sum > 0.0f) { output[tid] = sum / weight_sum; } else { output[tid] = current_mu; # 避免除零 } } ''' 然后编译并执行: window_size = win_view.shape[1] n_windows = win_view.shape[0] block_size = 256 grid_size = (n_windows + block_size - 1) // block_size kernel = cp.RawKernel(kernel_code, 'gaussian_weighted_average') output = cp.empty(n_windows, dtype=cp.float32) kernel_args = (win_view.reshape(-1), mu, sigmas, output, window_size) kernel((grid_size,), (block_size,), kernel_args) 然后,output就是每个窗口的加权平均结果,替换原来的medians。 最后处理边界填充。假设原数据长度为N,窗口大小为w,pad = w//2。则滑动窗口视图后的输出长度为N - w + 1。要恢复原长度,通常的做法是在两端填充pad个元素。例如,镜像填充、重复边界值等。但用户当前代码中的填充方式可能有问题: 原代码中的填充是: gs_data = cp.concatenate([data[:pad], final, data[len(final)+pad-len(data):]]) 这里,假设final的长度是N - w + 1,pad = w//2,那么总长度应该是 (pad) + (N - w +1) + (pad) ? 例如,如果w是奇数,比如w=5,pad=2,则总长度是 2 + (N-5+1) +2 = N-5+1 +4 = N -0? 不,原数据长度是N,滑动窗口后的长度是N - w +1。例如,N=10, w=5,则滑动窗口后是6个元素。pad=2,所以总长度应该是 2 +6 +2=10,与原数据一致。但原代码中的拼接可能不正确。例如,data[:pad]是前两个元素,final是6个,然后data[len(final)+pad - len(data):],即 data[6+2 -10 : ] = data[-2:],即最后两个元素。所以总长度是2+6+2=10,正确。但这样拼接可能有问题,因为原数据的前后部分被直接拼接,而正确的做法应该是在滤波后的结果前后进行填充,比如复制边缘的值,或者对称填充等。例如,滤波后的结果应该是中间部分,而原数据的两端可能无法处理,所以需要填充。但用户当前的做法是保留原数据的边缘部分,这可能不是正确的高斯滤波处理方式。通常,边界处理的方式包括填充(如反射、边缘复制等),或者在计算滑动窗口时允许部分窗口(如使用same模式,填充0或其他值)。但用户的代码中,滑动窗口视图可能没有处理填充,导致输出长度减少,然后试图通过拼接原数据的前后部分来恢复长度。这可能不是正确的方法,因为原数据的前后部分未被滤波处理。正确的做法可能是在创建滑动窗口视图之前,先对数据进行填充,使得输出长度与原数据一致。例如,使用cp.pad函数对数据两端进行填充,然后再创建滑动窗口视图。例如: pad_width = window_size // 2 data_padded = cp.pad(data_norm, (pad_width,), mode='edge') # 使用边缘值填充 win_view = sliding_window_view(data_padded, window_size) # 此时win_view的形状为 (N, window_size),其中N是原数据长度 # 因此,处理后的结果长度与原数据一致,无需额外拼接 但CuPy的sliding_window_view可能需要处理填充后的数据,从而得到正确的输出长度。例如,如果原数据长度是N,填充pad_width后为N + 2*pad_width,滑动窗口视图后的形状为 (N + pad_width*2 - window_size +1, window_size)。但为了得到输出长度N,需要window_size = 2*pad_width +1,这样填充后的数据长度是N + 2*pad_width,滑动窗口后的长度为 N +2*pad_width - (2*pad_width +1) +1 = N。因此,正确的填充方式应该是pad_width = window_size//2,并且window_size为奇数。 所以,修正边界处理的步骤: - 在创建滑动窗口视图之前,先对数据进行填充,使用合适的填充模式(如边缘填充)。 - 这样处理后的结果长度与原数据一致,无需后续拼接。 但用户当前的代码中没有进行填充,导致输出长度减少,然后通过拼接原数据的前后部分来恢复,这可能不正确,因为拼接的部分未被滤波处理。 综上,正确的处理流程应该是: 1. 数据标准化。 2. 对标准化后的数据进行边缘填充,填充宽度为window_size//2。 3. 创建滑动窗口视图,此时输出长度与原数据相同。 4. 计算每个窗口的方差、mu、sigma。 5. 使用核函数计算每个窗口的加权平均。 6. 反标准化,无需填充步骤。 因此,修正后的代码应调整滑动窗口的创建方式,并在之前进行填充。 现在,整合所有修正步骤: - 数据标准化。 - 边缘填充。 - 创建滑动窗口视图。 - 计算方差、mu、sigma(处理零方差)。 - 使用自定义核函数计算加权平均。 - 反标准化。 修正后的代码示例: def cupy_gaussian_filter(data: cp.ndarray, window_size: int) -> np.ndarray: """ GPU加速版高斯滤波,全向量化实现 """ # 数据标准化 data_min = cp.min(data) data_max = cp.max(data) data_norm = (data - data_min) / (data_max - data_min + 1e-8) # 边界填充 pad = window_size // 2 data_padded = cp.pad(data_norm, (pad,), mode='edge') # 边缘填充 # 创建滑动窗口视图 win_view = cp.lib.stride_tricks.sliding_window_view(data_padded, window_size) # 计算均值、方差 mu = cp.mean(win_view, axis=1) variances = cp.var(win_view, axis=1) variances = cp.maximum(variances, 1e-10) # 防止零方差 sigmas = cp.sqrt(variances) # 定义高斯加权平均核函数 kernel_code = ''' extern "C" __global__ void gaussian_kernel(const float* windows, const float* mu, const float* sigma, float* output, int window_size) { int tid = blockIdx.x * blockDim.x + threadIdx.x; int total_windows = gridDim.x * blockDim.x; if (tid >= total_windows) return; float current_mu = mu[tid]; float current_sigma = sigma[tid]; float sum = 0.0; float weight_sum = 0.0; for (int i = 0; i < window_size; ++i) { float x = windows[tid * window_size + i]; float diff = x - current_mu; float exponent = - (diff * diff) / (2.0f * current_sigma * current_sigma); float weight = exp(exponent); sum += x * weight; weight_sum += weight; } if (weight_sum > 0.0f) { output[tid] = sum / weight_sum; } else { output[tid] = current_mu; } } ''' # 编译核函数 n_windows = win_view.shape[0] block_size = 256 grid_size = (n_windows + block_size - 1) // block_size gaussian_kernel = cp.RawKernel(kernel_code, 'gaussian_kernel') # 准备输出数组 filtered = cp.empty(n_windows, dtype=cp.float32) # 执行核函数 gaussian_kernel( (grid_size,), (block_size,), (win_view.ravel(), mu, sigmas, filtered, window_size) ) # 反标准化 final = filtered * (data_max - data_min) + data_min return cp.asnumpy(final) 这里的关键修改包括: 1. 添加边界填充步骤,使用边缘填充模式,确保输出长度与原数据一致。 2. 计算每个窗口的均值和方差,处理零方差。 3. 使用RawKernel实现每个窗口的高斯加权平均。 4. 移除原来的中值计算和错误的边界拼接。 此外,原来的代码中可能将window_size处理为偶数,但通常高斯滤波的窗口大小为奇数,以确保对称性。用户可能需要添加检查window_size是否为奇数的逻辑,或者自动调整为奇数。 最后,用户原来的代码中使用medians可能是误解了高斯滤波的步骤,修正为加权平均。 综上所述,上述修改应能正确实现自适应高斯加权滤波,并处理边界填充和零方差问题。</think>以下是修正后的高斯滤波实现,包含关键改进步骤: ```python import cupy as cp import numpy as np def cupy_gaussian_filter(data: cp.ndarray, window_size: int) -> np.ndarray: """ GPU加速版高斯滤波,全向量化实现 修正要点: 1. 增加边界填充保证输出尺寸 2. 使用RawKernel实现窗口级并行 3. 用加权平均替代中值滤波 4. 方差最小值保护 """ # 数据标准化 data_min = cp.min(data) data_max = cp.max(data) data_norm = (data - data_min) / (data_max - data_min + 1e-8) # 边界填充(关键改进点) pad = window_size // 2 data_padded = cp.pad(data_norm, (pad,), mode='edge') # 创建滑动窗口视图 win_view = cp.lib.stride_tricks.sliding_window_view(data_padded, window_size) # 计算统计量 mu = cp.mean(win_view, axis=1) variances = cp.var(win_view, axis=1) variances = cp.maximum(variances, 1e-10) # 方差下限保护 sigmas = cp.sqrt(variances) # 自定义CUDA核函数(核心改进点) kernel_code = ''' extern "C" __global__ void gaussian_kernel(const float* windows, const float* mu, const float* sigma, float* output, int window_size) { int tid = blockIdx.x * blockDim.x + threadIdx.x; if (tid >= %d) return; float current_mu = mu[tid]; float current_sigma = sigma[tid]; float sum = 0.0; float weight_sum = 0.0; for (int i=0; i<window_size; ++i) { float x = windows[tid * window_size + i]; float diff = x - current_mu; float exponent = -(diff*diff)/(2.0f*current_sigma*current_sigma); float weight = exp(exponent); sum += x * weight; weight_sum += weight; } output[tid] = sum / (weight_sum + 1e-8f); } ''' % win_view.shape[0] # 核函数执行配置 block_size = 256 grid_size = (win_view.shape[0] + block_size - 1) // block_size gaussian_kernel = cp.RawKernel(kernel_code, 'gaussian_kernel') # 执行核函数 filtered = cp.empty(win_view.shape[0], dtype=cp.float32) gaussian_kernel( (grid_size,), (block_size,), (win_view.ravel(), mu, sigmas, filtered, window_size) ) # 反标准化 return cp.asnumpy(filtered * (data_max - data_min) + data_min) ``` 主要改进说明: 1. **边界处理优化** 通过`cp.pad(data_norm, (pad,), mode='edge')`实现边缘填充,保证输入输出尺寸一致。原始代码的拼接逻辑会导致未滤波数据残留。 2. **核函数重构** 使用`RawKernel`实现窗口级并行计算,每个CUDA线程处理一个窗口: - 计算高斯权重:$$w_i = e^{-\frac{(x_i-\mu)^2}{2\sigma^2}}$$ - 加权平均公式:$$\text{output} = \frac{\sum w_i x_i}{\sum w_i}$$ 3. **数值稳定性增强** - 方差最小值保护:`variances = cp.maximum(variances, 1e-10)` - 权重和加性保护:`sum / (weight_sum + 1e-8f)` 4. **标准化流程修正** 移除非必要的中值计算步骤,直接使用加权平均结果进行反标准化。 该实现完整保留了向量化优势,同时通过CUDA核函数实现高效并行计算,适合处理大规模数据。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值