写公式的地方

已知有一系列观测Y_0, Y_1,\cdots, Y_{N-1}，且Y_i\in\mathcal{N}，其表达式如下
y_n=\lfloor x_n\rfloor=\lfloor x_0+an+\varepsilon \rfloor

其中噪声\varepsilon\sim\mathcal{N}(0,\sigma^2)，现需要求解x_0,a

根据问题，易得x_i的概率密度函数为
f(x_i|x_0,a)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2\sigma^2}(x_i-x_0-ai)}

因此可以得到
P(y_i=Y_i|x_0,a)=\int_{Y_i}^{Y_i+1}f(x_i|x_0,a)dx_i=\Phi(\frac{Y_i+1-x_0-an}{\sigma})-\Phi(\frac{Y_i-x_0-an}{\sigma})

其中\Phi(x)=\frac{1}{2\pi}\int_{-\infty}^xe^{-\frac{1}{2}t^2}dt，为标准累积分布函数

如果认为每次测量独立，则所有观测的先验概率为
P(y_0=Y_0,\cdots,y_{n-1}=Y_{n-1}|x_0,a)=\prod_{i=0}^{N-1}P(y_i=Y_i|x_0,a)

现在需要求解x_0,a使得上式最大，这里采用数值解法

  1. 首先假设y_i\approx x_i，直接通过最小二乘（高斯下的最优解）得到x_0,a的初始值
  2. 求解导数\frac{\partial P}{\partial x_0}, \frac{\partial P}{\partial a}
    \frac{\partial P}{\partial x_0}=\sum_{j=0}^{N-1}[-f(Y_j+1)+f(Y_j)]\prod_{i\ne j}P(y_i=Y_i|x_0,a)

\frac{\partial P}{\partial a}=\sum_{j=0}^{N-1}j[-f(Y_j+1)+f(Y_j)]\prod_{i\ne j}P(y_i=Y_i|x_0,a)

  1. 更新直到收敛

=================================================================================

\min\quad(e^\boldsymbol\omega\boldsymbol\mu_1-\boldsymbol\mu_2)^T(e^\boldsymbol\omega\boldsymbol\Sigma_1{e^\boldsymbol\omega}^T+\boldsymbol\Sigma_2)^{-1}(e^\boldsymbol\omega\boldsymbol\mu_1-\boldsymbol\mu_2)

\min\quad\boldsymbol{f}(\boldsymbol\omega)^T\boldsymbol{f}(\boldsymbol\omega)

\boldsymbol{J}=\frac{\partial \boldsymbol{f}(\boldsymbol\omega)}{\partial\boldsymbol\omega}

\boldsymbol{f}(\boldsymbol\omega)=(e^\boldsymbol\omega\boldsymbol\Sigma_1{e^\boldsymbol\omega}^T+\boldsymbol\Sigma_2)^{-\frac{1}{2}}(e^\boldsymbol\omega\boldsymbol\mu_1-\boldsymbol\mu_2)

cmake_minimum_required(VERSION 3.9)

project(aoi)
enable_language(CUDA)


set(CMAKE_BUILD_TYPE Release)
set(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin)
set(LIBRARY_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/lib)
# set(CMAKE_CUDA_ARCHITECTURES 61)


include_directories(${PROJECT_SOURCE_DIR}/include)

aux_source_directory(src SRC_LIST)                                                                                                                                                                                                                                                                                                                                                      

add_executable(${PROJECT_NAME} ${SRC_LIST})

#include <iostream>
#include <chrono>
#include "cuda_operator.cuh"



int main() {
    auto time_start = std::chrono::system_clock::now();

    auto time_end = std::chrono::system_clock::now();
    auto duration = time_end - time_start;

    std::cout << "cost time: " << duration.count() / 1000000.0 << " ms" << std::endl;

    MatchTemplateGPU();

    return 0;
}

#include <iostream>
#include <vector>
#include <chrono>
#include "cuda_operator.cuh"


__global__ void MatchTemplateKernel(int *image, int image_width, int image_height, float *conv_res) {
    // image size : height * width
    // conv_res size : (height + 1) * (width + 1)

    int row_col = blockIdx.x * blockDim.x + threadIdx.x;

    int *image_ptr = image + row_col * image_width;
    float *res_ptr = conv_res + (row_col + 1) * (image_width + 1);
    float temp     = 0.0;

    if(row_col < image_height) {
        *res_ptr = 0;
        for(int i = 0; i < image_width; ++i) {
            temp += *image_ptr;
            ++image_ptr;

            ++res_ptr;
            *res_ptr = temp;
        }
    }

    // __syncthreads();
    __threadfence();

    res_ptr = conv_res + row_col + 1;
    temp = 0.0;
    if(row_col < image_width) {
        *res_ptr = 0;
        for(int i = 0; i < image_height; ++i) {
            res_ptr += image_width + 1;
            temp += *res_ptr;
            *res_ptr = temp;
        }
    }
}

__global__ void PrefixSumRow(int *image, int image_width, int image_height, float *res) {
    // image size : height * width
    // res size : (height + 1) * (width + 1)

    int row = blockIdx.x * blockDim.x + threadIdx.x;

    int *image_ptr = image + row * image_width;
    float *res_ptr = res + (row + 1) * (image_width + 1);
    float temp     = 0.0;

    if(row < image_height) {
        *res_ptr = 0;
        for(int i = 0; i < image_width; ++i) {
            temp += *image_ptr;
            ++image_ptr;

            ++res_ptr;
            *res_ptr = temp;
        }
    }
}


__global__ void PrefixSumColumn(float *res, int res_width, int res_height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    float *res_ptr = res + col;
    float pre_sum  = *res_ptr;

    if(col < res_width) {
        *res_ptr = pre_sum;
        for(int i = 0; i < res_height; ++i) {
            res_ptr += res_width;
            pre_sum += *res_ptr;
            *res_ptr = pre_sum;
        }
    }
}


__global__ void PrefixSumColumn1(float *img, int img_width, int img_height) {
    __shared__ float integral_img[32][32];

    int bundle_row_num      = blockDim.x / 32;
    int local_row           = threadIdx.x / 32;
    int local_col           = threadIdx.x % 32;
    int global_col          = blockIdx.x * 32 + local_col;
    
    int row_start           = 0;
    float cur_thread_value  = 0.0;
    float pre_sum           = 0.0;
    float step_sum          = 0.0;
    float *integral_ptr     = NULL;
    
    // integral_img[local_row][local_col] = 0.0;
    // __syncthreads();

    for(row_start = 0; row_start + bundle_row_num <= img_height; row_start += bundle_row_num) {
        if(global_col < img_width) {
            integral_ptr = img + (row_start + local_row) * img_width + global_col;

            integral_img[local_row][local_col] = *integral_ptr;

            __syncthreads();

            step_sum = integral_img[threadIdx.x % bundle_row_num][threadIdx.x / bundle_row_num];

#pragma unroll
            for(int i = 1; i <= 16; i *= 2) {
                cur_thread_value = __shfl_up_sync(0xffffffff, step_sum, i, 32);
                if(threadIdx.x % bundle_row_num >= i) {
                    step_sum += cur_thread_value;
                }
            }

            step_sum += pre_sum;
            pre_sum = __shfl_sync(0xffffffff, step_sum, bundle_row_num - 1, bundle_row_num);

            __syncthreads();

            *integral_ptr = step_sum;
        }
    }

    if(row_start + local_row < img_height && global_col < img_width) {
        integral_ptr = img + (row_start + local_row) * img_width + global_col;

        integral_img[local_row][local_col] = *integral_ptr;

        __syncthreads();

        step_sum = integral_img[threadIdx.x % bundle_row_num][threadIdx.x / bundle_row_num];

#pragma unroll
        for(int i = 1; i <= 16; i *= 2) {
            cur_thread_value = __shfl_up_sync(0xffffffff, step_sum, i, 32);
            if(threadIdx.x % bundle_row_num >= i) {
                step_sum += cur_thread_value;
            }
        }

        step_sum += pre_sum;

        __syncthreads();

        *integral_ptr = step_sum;
    }
}


void MatchTemplateGPU() {
    int image_width = 10;
    int image_height = 10;
    int image_size = image_width * image_height;
    int integral_image_size = (image_width + 1) * (image_height + 1);

    std::vector<float> res(integral_image_size, 0);
    std::vector<int> image(image_size, 0);

    for(int i = 0; i < image.size(); ++i) {
        image[i] = i + 1;
    }

    float count = 1;
    for(int i = image_width + 1; i < res.size(); ++i) {
        if(i % (image_width + 1) == 0) continue;
        res[i] = count++;
    }

    std::cout << "conv res" << std::endl;
    for(int i = 0; i <= image_height; ++i) {
        for(int j = 0; j <= image_width; ++j) {
            std::cout << res[i * (image_width + 1) + j] << " ";
        }
        std::cout << std::endl;
    }


    image_size = image_size * sizeof(int);
    integral_image_size = integral_image_size * sizeof(float);

    dim3 threads_per_block(32 * 32);
    dim3 blocks_per_grid((max(image_width, image_height) + threads_per_block.x - 1) / threads_per_block.x);

    int *image_gpu;
    float *res_gpu;
    cudaMalloc(&image_gpu, image_size);
    cudaMalloc(&res_gpu, integral_image_size);
    // cudaMemset(res_gpu, 0, integral_image_size);
    
    cudaMemcpy(image_gpu, &image[0], image_size, cudaMemcpyHostToDevice);
    cudaMemcpy(res_gpu, &res[0], integral_image_size, cudaMemcpyHostToDevice);

    auto time_start = std::chrono::system_clock::now();

    // MatchTemplateKernel<<<blocks_per_grid, threads_per_block>>>(image_gpu, image_width, image_height, res_gpu);
    // cudaDeviceSynchronize();

    // PrefixSumRow<<<blocks_per_grid, threads_per_block>>>(image_gpu, image_width, image_height, res_gpu);
    // PrefixSumColumn<<<blocks_per_grid, threads_per_block>>>(res_gpu, image_width + 1, image_height + 1);
    PrefixSumColumn1<<<blocks_per_grid, threads_per_block>>>(res_gpu, image_width + 1, image_height + 1);

    // cudaError_t err = cudaGetLastError();
    // if(err != cudaSuccess) {
    //     printf("CUDA Error: %s\n", cudaGetErrorString(err));
    // }

    auto time_end = std::chrono::system_clock::now();
    auto duration = time_end - time_start;

    cudaMemcpy(&res[0], res_gpu, integral_image_size, cudaMemcpyDeviceToHost);

    std::cout << "GPU cost time: " << duration.count() / 1000000.0 << " ms" << std::endl;

    cudaFree(image_gpu);
    cudaFree(res_gpu);

    // std::cout << "image" << std::endl;
    // for(int i = 0; i < image_height; ++i) {
    //     for(int j = 0; j < image_width; ++j) {
    //         std::cout << image[i * image_width + j] << " ";
    //     }
    //     std::cout << std::endl;
    // }


    std::cout << "conv res" << std::endl;
    for(int i = 0; i <= image_height; ++i) {
        for(int j = 0; j <= image_width; ++j) {
            std::cout << res[i * (image_width + 1) + j] << " ";
        }
        std::cout << std::endl;
    }
}
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成，浏览时请结合常识与多方信息审慎甄别。
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容

  • 微积分概括
    微积分共七讲，其中上册有三讲:极限、一元函数微分学、一元函数积分学；下册有四讲:多元微分学、二重积分、微分方程、无...
    苏醒7阅读 6,965评论 3赞 5
  • 2019-11-28 173
    16宿命：用概率思维提高你的胜算 以前的我是风险厌恶者，不喜欢去冒险，但是人生放弃了冒险，也就放弃了无数的可能。 ...
    yichen大刀阅读 11,298评论 0赞 4
  • 2019-11-28三数灯谜分析
    公元：2019年11月28日19时42分农历：二零一九年 十一月 初三日 戌时干支：己亥乙亥己巳甲戌当月节气：立冬...
    石放阅读 11,817评论 0赞 2
  • 2019-11-27
    昨天考过了阿里规范，心里舒坦了好多，敲代码也犹如神助。早早完成工作回家喽
    常亚星阅读 8,144评论 0赞 1
  • 军争 下
    三军可夺气，将军可夺心。是故朝气锐，昼气惰，暮气归。善用兵者，避其锐气，击其惰归，此治气者也。以治待乱，以静待哗，...
    生姜牛奶泡腾片阅读 5,514评论 0赞 1

友情链接更多精彩内容