#include "cnn.h"



void print_rslt(float* rslt, u8 input_matrix_length, u32 length){
    int _tmp = 0;
    printf("[0:0]");
    for (int i = 0; i < length; i++) {
        printf("%f ",rslt[i]);
        if ((i + 1) % input_matrix_length == 0) {
            printf("\n[%d:%d]",++_tmp,i+1);
        }
    }
    printf("\r\n");
}

// 将原始矩阵复制到填充后的矩阵中央
float* expand(const float* old_matrix, int old_matrix_length, int layer){
    float* new_matrix = (float *)malloc(sizeof(float)*layer*(old_matrix_length+2)*(old_matrix_length+2));
    memset(new_matrix, 0, sizeof(float)*layer*(old_matrix_length+2)*(old_matrix_length+2));
    for(int l=0; l < layer; l++){
        for (int i = 0; i < old_matrix_length; i++) {
            for (int j = 0; j < old_matrix_length; j++) {
                new_matrix[(i + 1) * (old_matrix_length+2) + (j + 1) +
                           l * (old_matrix_length+2) * (old_matrix_length+2)]
                = old_matrix[i * old_matrix_length + j +
                             l * (old_matrix_length) * (old_matrix_length)];
            }
        }
    }
    return new_matrix;
}

//model 模型名字
//input_matrix 输入图像
//input_matrix_length 输入图像的边长：102
//c_rl 输出图像的边长：100
//返回卷积的结果
float* convolution(Model model_w, Model model_b, const float* input_matrix, int input_matrix_length){
    // 初始化卷积层参数
    int im_l = input_matrix_length;
    int cr_l = input_matrix_length - 2;
    float conv_temp; // 临时变量，用于存储卷积计算的中间结果

    //用于合并前的数组，具有32*64*50*50(第二层)的大小
    float* _conv_rlst = (float *) malloc(sizeof (float) * model_w.num_kernels * (cr_l * cr_l));
    memset(_conv_rlst, 0, sizeof (float) * model_w.num_kernels * (cr_l * cr_l));
    //子图合并后的数组
    float* conv_rlst = (float *) malloc(sizeof (float) * model_w.num_kernels * (cr_l * cr_l));
    memset(conv_rlst, 0, sizeof (float) * model_w.num_kernels * (cr_l * cr_l));

    // 遍历30个卷积核（假设有30个通道）
    for(int c=0; c < model_w.channel; c++){
        for(int k=0; k < model_w.num_kernels; k++){
            for(int row = 0; row < cr_l; row++) {
                for (int col = 0; col < cr_l; col++) {
                    conv_temp = 0; // 每个输出像素初始化为0
                    // 进行3x3的卷积操作
                    for (int x = 0; x < 3; x++) {
                        for (int y = 0; y < 3; y++) {
                            // 将输入图像的对应像素与卷积核权重相乘，并累加到conv_temp
                            conv_temp += input_matrix[(c*im_l*im_l) + (row*(im_l)+col) + (x*(im_l)+y)]
                                    * model_w.array[((c+k*model_w.channel)*3*3) + (x*3+y)];
                        }
                    }
                    _conv_rlst[(k*cr_l*cr_l) + (row*cr_l+col)] = conv_temp;
                }
           }
        }
        //合并子图
        for(int k=0; k < model_w.num_kernels; k++) {
            for (int row = 0; row < cr_l; row++) {
                for (int col = 0; col < cr_l; col++) {
                    conv_rlst[(k*cr_l*cr_l) + (row*cr_l+col)] += _conv_rlst[(k*cr_l*cr_l) + (row*cr_l+col)];
                }
            }
        }
    }
    for(int k=0; k < model_w.num_kernels; k++) {
        for (int row = 0; row < cr_l; row++) {
            for (int col = 0; col < cr_l; col++) {
                conv_temp = conv_rlst[(k * cr_l * cr_l) + (row * cr_l + col)];
                // 加上对应卷积核的偏置
                conv_temp += model_b.array[k];
                // 激活函数：ReLU（将小于0的值设为0）
                if (conv_temp > 0)
                    conv_rlst[(k * (cr_l * cr_l)) + (row * cr_l) + (col)] = conv_temp; // 如果卷积结果大于0，存入结果数组
                else
                    conv_rlst[(k * (cr_l * cr_l)) + (row * cr_l) + (col)] = 0; // 否则存入0
            }
        }
    }
    free(_conv_rlst);
    _conv_rlst = NULL;
    return conv_rlst;
}

//num_kernels 卷积核的个数：32
//area 池化的面积：2*2
//input_matrix 输入图像
//input_matrix_length 输入图像的边长：100
//输出图像的边长：50
//返回池化的结果
float* pooling(Model model_w, const float* input_matrix, u8 input_matrix_length){
    u8 im_l = input_matrix_length;
    float pool_temp = 0; // 临时变量，用于存储池化操作的最大值
    float* pool_rslt = (float *) malloc(sizeof (float)*model_w.num_kernels*im_l*im_l);
    memset(pool_rslt, 0, sizeof (float)*model_w.num_kernels*im_l*im_l);
    // 遍历30个通道（与卷积核数量相同）
    for(u8 n=0; n<model_w.num_kernels; n++)
    {
        // 遍历输入图像的每一行，步长为2（2x2的池化窗口）
        for(u8 row=0; row<im_l; row=row+2)
        {
            // 遍历输入图像的每一列，步长为2
            for(u8 col=0; col<im_l; col=col+2)
            {
                pool_temp = 0; // 每个池化区域的最大值初始化为0
                // 进行2x2的最大池化操作
                for(u8 x=0; x<2; x++)
                {
                    for(u8 y=0; y<2; y++)
                    {
                        // 更新当前池化区域的最大值
                        if(pool_temp <= input_matrix[row*im_l+col+x*im_l+y+n*(im_l*im_l)])
                            pool_temp = input_matrix[row*im_l+col+x*im_l+y+n*(im_l*im_l)];
                    }
                }
                // 将最大值存入池化结果数组
                pool_rslt[(row/2)*(im_l/2)+col/2+n*((im_l/2)*(im_l/2))] = pool_temp;
            }
        }
    }
    return pool_rslt;
}

float* hidden(const float* input_matrix){
    float affine1_temp; // 临时变量，用于存储全连接层的中间结果
    float *affine1_rslt = (float *) malloc(sizeof(float)*128);
    memset(affine1_rslt, 0, sizeof(float)*128);

    // 遍历128个神经元（假设隐藏层有128个神经元）
    for(u8 n=0; n<128; n++)
    {
        affine1_temp = 0; // 每个神经元的输出初始化为0

        // 进行矩阵乘法，将池化层输出展平为一维向量后，与全连接层权重进行点积
        for(int i=0; i<(128*12*12); i++)
        {
            affine1_temp = affine1_temp + input_matrix[i] * fc1_weight.array[i+(128*12*12)*n];
        }

        // 加上对应神经元的偏置
        affine1_temp = affine1_temp + fc1_bias.array[n];

        // 激活函数：ReLU（将小于0的值设为0）
        if(affine1_temp > 0)
            affine1_rslt[n] = affine1_temp; // 如果结果大于0，存入结果数组
        else
            affine1_rslt[n] = 0; // 否则存入0
    }

    //    print_rslt(affine1_rslt,1,128);

    return affine1_rslt;
}

float* output(Model model_w, const float* input_matrix){
    u8 num = model_w.num_kernels;
    float affine2_temp; // 临时变量，用于存储输出层的中间结果
    float *affine2_rslt = (float *) malloc(sizeof(float)*num);
    memset(affine2_rslt, 0, sizeof(float)*num);

// 遍历10个输出神经元（假设有10个类别）
    for(int n=0; n<num; n++)
    {
        affine2_temp = 0; // 当前神经元的输出初始化为0

        // 进行矩阵乘法，将隐藏层的输出与输出层权重进行点积
        for(int i=0; i<128; i++)
        {
            affine2_temp = affine2_temp + fc2_weight.array[i+128*n] * input_matrix[i];
        }

        // 加上对应神经元的偏置
        affine2_temp = affine2_temp + fc2_weight.array[n];
        affine2_rslt[n] = affine2_temp; // 存储输出层的结果
    }

    return affine2_rslt;
}



void calculate_statistics(Model model, float* value)
{
    value[0] = fabsf(model.array[0]);
    float sum = 0;
    float sum_sq = 0;

    for (int i = 0; i < model.maxlength; i++) {
        float abs_val = fabsf(model.array[i]);

        if (abs_val > value[0]) {
            value[0] = abs_val;
        }
        sum += abs_val;
        sum_sq += abs_val * abs_val;
    }

    value[1] = sum / (float)model.maxlength;

    float variance = (sum_sq / (float)model.maxlength) - (value[1] * value[1]);
    value[2] = sqrtf(variance);
}

u8 check_threshold(Model model, const float* value)
{
    const float threshold = 20;

    for (int i = 0; i < model.maxlength; i++) {
        float K = (fabsf(model.array[i]) - value[1]) / value[2];
        if (K > threshold) {
            return 1;
        }
    }
    return 0;
}

float* generateMatrix(Model model, const float* value)
{
    float* CNN_data = (float*) malloc(sizeof(float)*100*100);
    memset(CNN_data, 0, sizeof(float)*100*100);

    u16 x = model.maxlength / 100;
    float y = value[0] / 100;

    for (int i = 0; i < model.maxlength; i++) {
        float absolutevalue = fabsf(model.array[i]);

        if (!absolutevalue) {
            continue;
        }

        int xIndex = i / x;
        if (xIndex >= 100) xIndex = 99;

        int yIndex = (int)(absolutevalue / y);
        if (yIndex < 0) yIndex = 0;

        CNN_data[yIndex * 100 + xIndex]++;
    }
    return CNN_data;
}

float calculate_probabilities(Model model_w, float *input_array)
{
    float sum = 0;
    u8 input_num = model_w.num_kernels;
    float *result = (float *) malloc(sizeof(float)*input_num);
    memset(result, 0, sizeof(float)*input_num);

    float *temp = (float *) malloc(sizeof(float)*input_num);
    memset(temp, 0, sizeof(float)*input_num);

    for (int i = 0; i < input_num; i++)
    {
        temp[i] = expf(input_array[i]);
        sum = sum + temp[i];
    }
    for (int j = 0; j < input_num; j++)
    {
        result[j] = temp[j] / sum;
        if(isnan(result[j]))result[j] = 1;
    }

    int max_index = 0;
    float max_value = result[0];
    for (int k = 1; k < input_num; k++)
    {
        if (result[k] > max_value)
        {
            max_value = result[k];
            max_index = k;
        }
    }

    float _tmp = result[max_index] * 100;
    free(temp);
    temp = NULL;
    free(result);
    result = NULL;
    return _tmp;
}


u8 calculate_layer(Model model_w, float *input_array){
    u8 input_num = model_w.num_kernels;
    u8 predict_num = 0;
    float max_temp = -100;
    for(int n=0; n<input_num; n++)
    {
        if(max_temp <= input_array[n])
        {
            max_temp = input_array[n]; // 更新最大值
            predict_num = n; // 记录最大值对应的类别索引
        }
    }
    print_rslt(input_array,input_num,input_num);
    return predict_num+0;
}


void cnn_run(){
#if is1300000
    float value[3] = {0};
    calculate_statistics(data,&value[0]);
    if (check_threshold(data,&value[0])){
        float* _data = generateMatrix(data,&value[0]);

        //第一层：填充102 * 102
        float* expand_matrix_1 = expand(_data, 100, 1);
        float* conv_rlst_1 = convolution(conv1_weight,conv1_bias,expand_matrix_1, 102);
        float* pool_rslt_1 = pooling(conv1_weight, conv_rlst_1, 100);

        free(_data);
        _data = NULL;
        free(expand_matrix_1);
        expand_matrix_1 = NULL;
        free(conv_rlst_1);
        conv_rlst_1 = NULL;

        //第二层：填充32 * 52 * 52
        float* expand_matrix_2 = expand(pool_rslt_1, 50, 32);
        float* conv_rlst_2 = convolution(conv2_weight,conv2_bias,expand_matrix_2, 52);
        float* pool_rslt_2 = pooling(conv2_weight, conv_rlst_2, 50);

        free(pool_rslt_1);
        pool_rslt_1 = NULL;
        free(expand_matrix_2);
        expand_matrix_2 = NULL;
        free(conv_rlst_2);
        conv_rlst_2 = NULL;

        //第三层：填充 64 * 27 * 27
        float* expand_matrix_3 = expand(pool_rslt_2, 25, 64);
        float* conv_rlst_3 = convolution(conv3_weight,conv3_bias,expand_matrix_3, 27);
        float* pool_rslt_3 = pooling(conv3_weight, conv_rlst_3, 25);

        free(pool_rslt_2);
        pool_rslt_2 = NULL;
        free(expand_matrix_3);
        expand_matrix_3 = NULL;
        free(conv_rlst_3);
        conv_rlst_3 = NULL;

        float* affine1_rslt = hidden(pool_rslt_3);
        float* affine2_rslt = output(fc2_weight, affine1_rslt);

        printf("概率：%f\r\n",calculate_probabilities(fc2_weight, affine2_rslt));
        printf("Label is:%d\r\n",calculate_layer(fc2_weight, affine2_rslt));

        free(pool_rslt_3);
        pool_rslt_3 = NULL;
        free(affine1_rslt);
        affine1_rslt = NULL;
        free(affine2_rslt);
        affine2_rslt = NULL;

    } else{
        printf("未放电！最大值:%f 平均值:%f 标准差:%f\r\n",value[0],value[1],value[2]);
    }
#else
        float* expand_matrix_1 = expand(data.array, 100, 1);
        float* conv_rlst_1 = convolution(conv1_weight,conv1_bias,expand_matrix_1, 102);
        float* pool_rslt_1 = pooling(conv1_weight, conv_rlst_1, 100);

        free(expand_matrix_1);
        expand_matrix_1 = NULL;
        free(conv_rlst_1);
        conv_rlst_1 = NULL;

        //第二层：填充32 * 52 * 52
        float* expand_matrix_2 = expand(pool_rslt_1, 50, 32);
        float* conv_rlst_2 = convolution(conv2_weight,conv2_bias,expand_matrix_2, 52);
        float* pool_rslt_2 = pooling(conv2_weight, conv_rlst_2, 50);

        free(pool_rslt_1);
        pool_rslt_1 = NULL;
        free(expand_matrix_2);
        expand_matrix_2 = NULL;
        free(conv_rlst_2);
        conv_rlst_2 = NULL;

        //第三层：填充 64 * 27 * 27
        float* expand_matrix_3 = expand(pool_rslt_2, 25, 64);
        float* conv_rlst_3 = convolution(conv3_weight,conv3_bias,expand_matrix_3, 27);
        float* pool_rslt_3 = pooling(conv3_weight, conv_rlst_3, 25);

        free(pool_rslt_2);
        pool_rslt_2 = NULL;
        free(expand_matrix_3);
        expand_matrix_3 = NULL;
        free(conv_rlst_3);
        conv_rlst_3 = NULL;

        float* affine1_rslt = hidden(pool_rslt_3);
        float* affine2_rslt = output(fc2_weight, affine1_rslt);

        printf("概率：%f\r\n",calculate_probabilities(fc2_weight, affine2_rslt));
        printf("Label is:%d\r\n",calculate_layer(fc2_weight, affine2_rslt));

        free(pool_rslt_3);
        pool_rslt_3 = NULL;
        free(affine1_rslt);
        affine1_rslt = NULL;
        free(affine2_rslt);
        affine2_rslt = NULL;
#endif
}