#include "cnn.h"

u8 isrun;

void _cnn_run(){
	isrun = 1;
}

/*调试用打印函数*/
/// @brief 打印图像矩阵
/// @param array 输入的图像
/// @param array_num 输入图像总的个数
/// @param elements_per_line 每一行的个数，个数到了自动换行
void PrintfArray(float *array, int array_num, int elements_per_line)
{
    for (int i = 0; i < array_num; i++)
    {
        printf("%f ", array[i]); // 打印当前元素
        // 每打印完指定数量的元素后换行
        if ((i + 1) % elements_per_line == 0)
        {
            printf("\n");
        }
    }
    // 如果最后一行元素不足，手动换行
    if (array_num % elements_per_line != 0)
    {
        printf("\n");
    }
}



void PrintfResArray(float *array, int array_num, int elements_per_line)
{
    for (int i = 0; i < array_num; i++)
    {
        printf("%f ", array[i]); // 打印当前元素
        // 每打印完指定数量的元素后换行
        if ((i + 1) % elements_per_line == 0)
        {
            printf("\t");
        }
    }
}



/*二维数组指针*/
/// @brief 动态申请二维数组，类型为elementSize
/// @param depth 二维数组深度，即a[depth][]
/// @param num 长度，即a[][num]
/// @param elementSize 二维数组类型
void **allocate2DArray(int depth, int num, size_t elementSize)
{
    void **array = (void **)mymalloc(SRAMEX, depth * sizeof(void *));
    for (int d = 0; d < depth; d++)
    {
        array[d] = mymalloc(SRAMEX, num * elementSize);
    }
    return array;
}

/// @brief 释放通过allocate2DArray申请的二维数组
/// @param array 二维数组
/// @param depth 深度
void free2DArray(float **array, int depth)
{
    for (int d = 0; d < depth; d++)
    {
        myfree(SRAMEX, array[d]);
    }
    myfree(SRAMEX, array);
}

/*卷积相关函数*/
/// @brief 对输入的图像进行四周0填充
/// @param inputArray 输入的图像
/// @param input_size 输入的单行的元素个数，例如：100即一行有100个元素
/// @param outputArray 保存输出的图像
void Full(float *inputArray, int input_size, float *outputArray)
{
    int i, j;
    for (i = 0; i < ((input_size + 2) * (input_size + 2)); i++)
    { 
        outputArray[i] = 0;
    }
    for (i = 0; i < input_size; i++)
    {
        for (j = 0; j < input_size; j++)
        {
            outputArray[(i + 1) * (input_size + 2) + (j + 1)] = inputArray[i * input_size + j];
        }
    }
}


/// @brief 对输入的图像进行池化，选取kernel_size*kernel_size的范围的最大值，步长为step
/// @param inputArray 输入的图像
/// @param input_size 输入图像的单行的元素个数
/// @param kernel_size 池化核的单行的元素个数
/// @param step 步进值
/// @param outputArray 保存输出的图像
void Pooling(float *inputArray, int input_size,
             int kernel_size, unsigned int step,
             float *outputArray)
{
    int output_size = (input_size - kernel_size) / step + 1;
    
    for (int i = 0; i < output_size; i++)
    {
        for (int j = 0; j < output_size; j++)
        {
            float max_value = 0;
            
            for (int m = 0; m < kernel_size; m++)
            {
                for (int n = 0; n < kernel_size; n++)
                {
                    int input_row = i * step + m; 
                    int input_col = j * step + n;
                    int input_idx = input_row * input_size + input_col;
                    
                    if (inputArray[input_idx] > max_value)
                    {
                        max_value = inputArray[input_idx];
                    }
                }
            }
            
            int output_idx = i * output_size + j;
            outputArray[output_idx] = max_value;
        }
    }
}

/// @brief 对输入图像进行卷积，选取kernel_size*kernel_size的范围，这里步长为1
/// @param inputArray 输入的图像
/// @param input_size 输入图像的单行的元素个数
/// @param kernel 输入的卷积核，即kernel_size * kernel_size的卷积核参数
/// @param kernel_size 池化核的单行的元素个数
/// @param outputArray 保存输出的图像
void Convolution(float *inputArray, int input_size,
                 float *kernel, int kernel_size,
                 float *outputArray)
{
    int i, j, m, n;
    int half_k = kernel_size / 2;
    int output_size = input_size - 2 * half_k;

    for (i = half_k; i < input_size - half_k; i++)
    {
        for (j = half_k; j < input_size - half_k; j++)
        {
            float sum = 0;

            for (m = 0; m < kernel_size; m++)
            {
                for (n = 0; n < kernel_size; n++)
                {
                    int input_row = i + m - half_k;
                    int input_col = j + n - half_k;
                    int input_idx = input_row * input_size + input_col;
                    int kernel_idx = m * kernel_size + n;
                    sum += inputArray[input_idx] * kernel[kernel_idx];
                }
            }
            int output_idx = (i - half_k) * output_size + (j - half_k);
            outputArray[output_idx] = sum;
        }
    }
}

/// @brief 对二维数组做合并处理，对应位相加
/// @param inputArray 输入的图像
/// @param input_depth 输入图像的深度，即a[input_depth][]
/// @param input_size 输入图像的单行的元素个数
/// @param outputArray 保存输出的图像
void Combine(float **inputArray, int input_depth, int input_size, float *outputArray)
{
    int i, j, k; 
    int input_idx; 

    for (i = 0; i < input_size; i++)
    {
        
        for (j = 0; j < input_size; j++)
        {
            float sum = 0; 
            
            input_idx = i * input_size + j;
            for (k = 0; k < input_depth; k++)
            {
                sum += inputArray[k][input_idx];
            }
            outputArray[i * input_size + j] = sum;
        }
    }
}

/// @brief 展开二维数组成一维数组
/// @param inputArray 输入的图像
/// @param input_depth 输入图像的深度，即a[input_depth][]
/// @param input_size 输入图像的单行的元素个数
/// @param outputArray 保存输出的图像
void Flatten2D(float **inputArray, int input_depth, int input_size, float *outputArray)
{
    int i, j, k;
    for (k = 0; k < input_depth; k++)
    {
        for (i = 0; i < input_size; i++)
        {
            for (j = 0; j < input_size; j++)
            {
                int input_idx = i * input_size + j;
                outputArray[k * input_size * input_size + input_idx] = inputArray[k][input_idx];
            }
        }
    }
}

/// @brief 对输入图像的每一位添加偏置
/// @param inputArray 输入的图像
/// @param input_num 输入图像总的元素个数
/// @param bias 偏置核
/// @param outputArray 保存输出的图像
void AddBias(float *inputArray, int input_num, float bias, float *outputArray)
{
    for (int i = 0; i < input_num; i++)
    {
        outputArray[i] = inputArray[i] + bias;
    }
}

/// @brief 全连接层操作：对一个神经元做全连接
/// @param inputArray 输入的图像
/// @param input_num 输入图像总的元素个数
/// @param input_w 输入的权重核
/// @param input_b 输入的偏置核
/// @return 一个神经元的结果，类型为float
float ConnectedLayer(float *inputArray, int input_num,
                      float *input_w, float input_b)
{
    int i;
    float sum = 0;
    for (i = 0; i < input_num; i++)
    {
        sum += inputArray[i] * input_w[i];
    }
    sum = sum + input_b;
    return sum;
}

/// @brief ReLU函数激活，此函数是对一整个数组做激活：若值大于等于0，则保留值，否则置0
/// @param inputArray 输入的图像
/// @param num 输入图像总的元素个数
/// @param outputArray 保存输出的矩阵
void ReLU1(float *inputArray, int num, float *outputArray) 
{
    for (int i = 0; i < num; i++) {
        outputArray[i] = (inputArray[i] > 0) ? inputArray[i] : 0;
    }
}

/// @brief ReLU函数激活，此函数是对一个数做激活：若值大于等于0，则保留值，否则置0
/// @param data 输入的数
/// @return 激活完的数，类型为float
float ReLU2(float data)
{
    if (data > 0) {
        return data;
    }
    else {
        return 0;
    }
}


//求数值

/// @brief 对一个输入数组求平均值、最大值和标准差
/// @param arr 输入的数组
/// @param size 数组的长度
/// @param max_val 最大值保存
/// @param mean 平均值保存
/// @param std_dev 标准差保存
void calculate_statistics(float arr[], int size, float *max_val, float *mean, float *std_dev) 
{
    *max_val = fabs(arr[0]);
    float sum = 0.0;
    float sum_sq = 0.0;

    for (int i = 0; i < size; i++) {
        float abs_val = fabs(arr[i]); 

        if (abs_val > *max_val) {
            *max_val = abs_val; 
        }
        sum += abs_val;
        sum_sq += abs_val * abs_val; 
    }

    *mean = sum / size;

    float variance = (sum_sq / size) - (*mean * *mean);
    *std_dev = sqrt(variance);
}
//判断放电

/// @brief 判断是否放电
/// @param arr 输入数组
/// @param size 数组的长度
/// @param mean 平均值
/// @param std_dev 标准差
/// @return 如果超过阈值则返回1，否则返回0
int check_threshold(float arr[], int size, float *mean, float *std_dev) 
{
    const float threshold = 20.0;

    for (int i = 0; i < size; i++) {
        float K = (arr[i] - *mean) / *std_dev;
        if (K > threshold) {
            return 1;
        }
    }
    return 0;
}

/// @brief 生成100*100矩阵
/// @param get_data 输入的数组
/// @param Max_value 输入最大值
/// @param totalPoints 总的数组元素个数
/// @param CNN_data 输出的数组（二维）
void generateMatrix(float *get_data, float Max_value, int totalPoints, float CNN_data[100*100]) 
{
    for (int i = 0; i < 100; i++) {
        for (int j = 0; j < 100; j++) {
            CNN_data[i*100+j] = 0;
        }
    }

    int pointsPerInterval = totalPoints / 100; 
    float amplitudeStep = Max_value / 100;

    for (int i = 0; i < totalPoints; i++) {
        float amplitudeValue = fabsf(get_data[i]);

        if (amplitudeValue == 0.0f) {
            continue;
        }

        int intervalIndex = i / pointsPerInterval;
        if (intervalIndex >= 100) intervalIndex = 99;

        int amplitudeIndex = 99 - (int)(amplitudeValue / amplitudeStep);
        if (amplitudeIndex < 0) amplitudeIndex = 0; 

        CNN_data[amplitudeIndex*100+intervalIndex]++;
    }
}

/// @brief 计算概率
/// @param input_array 输入数组 
/// @param output_array 保存输出的数组
/// @param input_num 输入元素个数
/// @return 返回概率中最大值的索引（默认+1，即如果是数组的a[0]最大，则返回1）
int calculate_probabilities(float *input_array, float *output_array, int input_num)
{
    float sum = 0.0;
    float temp[input_num];
    for (int i = 0; i < input_num; i++)
    {
        temp[i] = exp(input_array[i]);
        sum = sum + temp[i];
    }
    for (int j = 0; j < input_num; j++)
    {
        output_array[j] = temp[j] / sum;
    }
    int max_index = 0; 
    float max_value = output_array[0];
    for (int k = 1; k < input_num; k++)
    {
        if (output_array[k] > max_value)
        {
            max_value = output_array[k];
            max_index = k; 
        }
    }
    return max_index + 1;
}


/*这个程序仅是基本验证程序能跑，还无优化*/

/*原始数据个数
*目前不保证能否在数目不对的情况下计算，有可能程序出错。
*目前没法再没有提前给出元素个数的情况下计算。
*原始数据个数仅印象预处理的计算，后面的卷积程序用不到。
*/

/// @brief 完整的程序
/// @param data 输入的一维数组原始数据 --> data.array代替
/// @param size 原始数据个数 --> data.maxlength代替
void cnn_run()
{
    float maxvalue = 0.0;          //最大值
    float meanvalue = 0.0;         //平均值
    float std_devvalue = 0.0;      //标准差
    int index = 0;                  //索引（忘了干啥的）
	float *Matrix_data = (float*)mymalloc(SRAMEX, sizeof(float) * 100 * 100);  //生成的100*100的矩阵保存在这

	printf("CNN模型正在运行！！！\r\nData数据集为：%s\r\n",data.dname);
	DEBUG_PRINTF("data[%d]: %f\r\n",0,data.array[0]);
	DEBUG_PRINTF("data[%d]: %f\r\n",1,data.array[1]);
	DEBUG_PRINTF("data[%d]: %f\r\n",2,data.array[2]);
	DEBUG_PRINTF("data[%d]: %f\r\n",299,data.array[299]);
	DEBUG_PRINTF("data[%d]: %f\r\n",300,data.array[300]);
	DEBUG_PRINTF("data[%d]: %f\r\n",301,data.array[301]);
	DEBUG_PRINTF("data[%d]: %f\r\n",1249997,data.array[1249997]);
	DEBUG_PRINTF("data[%d]: %f\r\n",1249998,data.array[1249998]);
	DEBUG_PRINTF("data[%d]: %f\r\n",1249999,data.array[1249999]);

    /*计算数据*/
    calculate_statistics(data.array, data.maxlength, &maxvalue, &meanvalue,&std_devvalue);

    /*判断数据是否超过阈值*/
    int x = check_threshold(data.array, data.maxlength, &meanvalue, &std_devvalue);
    if (x == 1) {
        printf("Start\r\n");
        //pre 开始预处理
        /*生成矩阵*/
        generateMatrix(data.array, maxvalue, data.maxlength, Matrix_data);
        float *CNN_data = (float*)mymalloc(SRAMEX, sizeof(float) * 100 * 100);
        /*二维数组转一维数组，以便后续计算*/
        for (int i = 0; i < 100; i++) {
            for (int j = 0; j < 100; j++) {
                CNN_data[9999 - index++] = Matrix_data[i*100+j];
            }
        }
		DEBUG_PRINTF("开始第一层！\r\n");
        //1 第一层卷积
        /*数组资源申请*/
        float* Full_output1 = (float*)mymalloc(SRAMEX, sizeof(float) * 102 * 102);
        /*填充矩阵*/
        Full(CNN_data, 100, Full_output1);
        float** Convolution_result1_before = (float**)allocate2DArray(32, 100 * 100, sizeof(float));
        float** Convolution_result1_relu = (float**)allocate2DArray(32, 100 * 100, sizeof(float));
        float** Convolution_result1 = (float**)allocate2DArray(32, 100 * 100, sizeof(float));
        /*卷积操作*/
        for (int i = 0; i < 32; i++) {
            /*提取卷积核（还没优化）*/
            float conv1_weight_new[9] = {
                conv1_weight.array[i*9+0],conv1_weight.array[i*9+1],
                conv1_weight.array[i*9+2],conv1_weight.array[i*9+3],
                conv1_weight.array[i*9+4],conv1_weight.array[i*9+5],
                conv1_weight.array[i*9+6],conv1_weight.array[i*9+7],
                conv1_weight.array[i*9+8] };
            Convolution(Full_output1, 102, conv1_weight_new, 3, Convolution_result1_before[i]);
        }
        /*添加偏置*/
        for (int i = 0; i < 32; i++) {
            AddBias(Convolution_result1_before[i], 100 * 100, conv1_bias.array[i], Convolution_result1_relu[i]);
        }
        /*做ReLU激活*/
        for (int i = 0; i < 32; i++) {
            ReLU1(Convolution_result1_relu[i], 100 * 100, Convolution_result1[i]);
        }
        /*池化操作*/
        float ** Pooling_result1 = (float**)allocate2DArray(32, 50 * 50, sizeof(float));
        for (int i = 0; i < 32; i++) {
            Pooling(Convolution_result1[i], 100, 2, 2, Pooling_result1[i]);
        }
		DEBUG_PRINTF("第一层完毕！\r\n");
        //2 第二层卷积
        /*释放上一层的资源*/
        myfree(SRAMEX, Full_output1);
        free2DArray(Convolution_result1_relu,32);
        free2DArray(Convolution_result1_before,32);
        free2DArray(Convolution_result1,32);
        /*填充操作*/
        float** Full_output2 = (float**)allocate2DArray(32, 52 * 52, sizeof(float));
        for (int i = 0; i < 32; i++) {
            Full(Pooling_result1[i], 50, Full_output2[i]);
        }
        /*申请本层资源*/
        float** Convolution_result_temp_2 = (float**)allocate2DArray(32, 50 * 50, sizeof(float));
        float** Convolution_result2_before = (float**)allocate2DArray(64, 50 * 50, sizeof(float));
        float** Convolution_result2_relu = (float**)allocate2DArray(64, 50 * 50, sizeof(float));
        float** Convolution_result2 = (float**)allocate2DArray(64, 50 * 50, sizeof(float));
        /*卷积操作*/
        for (int i = 0; i < 64; i++) {
            for (int j = 0; j < 32; j++) {
                /*提取卷积核*/
                float conv2_weight_new[9] = {
                    conv2_weight.array[i*32*9+9*j+0],conv2_weight.array[i*32*9+9*j+1],
                    conv2_weight.array[i*32*9+9*j+2],conv2_weight.array[i*32*9+9*j+3],
                    conv2_weight.array[i*32*9+9*j+4],conv2_weight.array[i*32*9+9*j+5],
                    conv2_weight.array[i*32*9+9*j+6],conv2_weight.array[i*32*9+9*j+7],
                    conv2_weight.array[i*32*9+9*j+8]
                };
                Convolution(Full_output2[j], 52, conv2_weight_new, 3, Convolution_result_temp_2[j]);
            }
            /*每一个子图像合并*/
            Combine(Convolution_result_temp_2, 32, 50, Convolution_result2_before[i]);
        }
        /*添加偏置*/
        for (int i = 0; i < 64; i++) {
            AddBias(Convolution_result2_before[i], 50 * 50, conv2_bias.array[i], Convolution_result2_relu[i]);
        }
        /*数组ReLU激活*/
        for (int i = 0; i < 64; i++) {
            ReLU1(Convolution_result2_relu[i], 50 * 50, Convolution_result2[i]);
        }
        /*池化操作*/
        float** Pooling_result2 = (float**)allocate2DArray(64, 25 * 25, sizeof(float));
        for (int i = 0; i < 64; i++) {
            Pooling(Convolution_result2[i], 50, 2, 2, Pooling_result2[i]);
        }
		DEBUG_PRINTF("第二层完毕！\r\n");
        //3 第三层卷积
        /*释放上一层的资源*/
        free2DArray(Full_output2,32);
        free2DArray(Pooling_result1,32);
        free2DArray(Convolution_result_temp_2,32);
        free2DArray(Convolution_result2_relu,64);
        free2DArray(Convolution_result2_before,64);
        free2DArray(Convolution_result2,64);
        /*填充操作*/
        float** Full_output3 = (float**)allocate2DArray(64, 27 * 27, sizeof(float));
        for (int i = 0; i < 64; i++) {
            Full(Pooling_result2[i], 25, Full_output3[i]);
        }
        /*申请本层资源*/
        float** Convolution_result_temp_3 = (float**)allocate2DArray(64, 25 * 25, sizeof(float));
        float** Convolution_result3_before = (float**)allocate2DArray(128, 25 * 25, sizeof(float));
        float** Convolution_result3_relu = (float**)allocate2DArray(128, 25 * 25, sizeof(float));
        float** Convolution_result3 = (float**)allocate2DArray(128, 25 * 25, sizeof(float));
        /*卷积操作*/
        for (int i = 0; i < 128; i++) {
            for (int j = 0; j < 64; j++) {
                float conv3_weight_new[9] = {
                    conv3_weight.array[i*64*9+9*j+0],conv3_weight.array[i*64*9+9*j+1],
                    conv3_weight.array[i*64*9+9*j+2],conv3_weight.array[i*64*9+9*j+3],
                    conv3_weight.array[i*64*9+9*j+4],conv3_weight.array[i*64*9+9*j+5],
                    conv3_weight.array[i*64*9+9*j+6],conv3_weight.array[i*64*9+9*j+7],
                    conv3_weight.array[i*64*9+9*j+8]
                };
                Convolution(Full_output3[j], 27, conv3_weight_new, 3, Convolution_result_temp_3[j]);
            }
            /*子图像合并*/
            Combine(Convolution_result_temp_3, 64, 25, Convolution_result3_before[i]);
        }
        /*添加偏置*/
        for (int i = 0; i < 128; i++) {
            AddBias(Convolution_result3_before[i], 25 * 25, conv3_bias.array[i], Convolution_result3_relu[i]);
        }
        /*ReLU激活*/
        for (int i = 0; i < 128; i++) {
            ReLU1(Convolution_result3_relu[i], 25 * 25, Convolution_result3[i]);
        }
        /*池化操作*/
        float** Pooling_result3 = (float**)allocate2DArray(128, 12 * 12, sizeof(float));
        for (int i = 0; i < 128; i++) {
            Pooling(Convolution_result3_before[i], 25, 2, 2, Pooling_result3[i]);
        }
        /*池化的输出展开成一维数组*/
        float* xi = (float *)mymalloc(SRAMEX, sizeof(float) * 128 * 12 * 12);
        Flatten2D(Pooling_result3, 128, 12, xi);
		DEBUG_PRINTF("第三层完毕！\r\n");
        //4 第一连接层
        /*释放上一层的资源*/
        free2DArray(Full_output3,64);
        free2DArray(Pooling_result2,64);
        free2DArray(Convolution_result_temp_3,64);
        free2DArray(Convolution_result3_relu,128);
        free2DArray(Convolution_result3_before,128);
        free2DArray(Convolution_result3,128);
        /*连接层操作*/
        //float yi[128] = {0};
		DEBUG_PRINTF("开始第四层！\r\n");
		float* yi = (float *)mymalloc(SRAMEX, 128 * sizeof(float));
		memset(yi, 0, 128 * sizeof(float));
        for (int i = 0; i < 128; i++) {
            float sum = 0;
            //float fc1_weight_new[128*12*12];
			float* fc1_weight_new = (float*)mymalloc(SRAMEX, sizeof(float) * 128 * 12 * 12);
            /*提取卷积核*/
            memcpy(fc1_weight_new,&fc1_weight.array[i*128*12*12],128*12*12 * sizeof(float));
            sum = ConnectedLayer(xi, 128 * 12 * 12, fc1_weight_new, fc1_bias.array[i]);
            /*单个数做ReLU激活*/
            yi[i] = ReLU2(sum);
			myfree(SRAMEX, fc1_weight_new);
        }
		DEBUG_PRINTF("第四层完毕！\r\n");
        //5 最后连接层
        /*释放上一层的资源*/
        free2DArray(Pooling_result3,128);
        /*连接层操作*/
        float zi[7] = {0};
//		float *zi = (float *)mymalloc(SRAMEX, 7 * sizeof(float));
//		memset(zi, 0, 7 * sizeof(float));
        for (int i = 0; i < 7; i++) {
            //float fc2_weight_new[128];
			float* fc2_weight_new = (float*)mymalloc(SRAMEX, sizeof(float) * 128);
            /*提取卷积核*/
            memcpy(fc2_weight_new,&fc2_weight.array[i*128],128 * sizeof(float));
            /*单个数做ReLU激活*/
            zi[i] = ConnectedLayer(yi, 128, fc2_weight_new, fc2_bias.array[i]);
			myfree(SRAMEX, fc2_weight_new);
        }
		DEBUG_PRINTF("第五层完毕！\r\n");
        //end 输出处理层
        float result[7];
        int max_probability_idx = calculate_probabilities(zi,result,7);
        PrintfResArray(zi,7,7);
        PrintfResArray(result,7,7);
        printf("%f, Label %d\t", result[max_probability_idx - 1] * 100, max_probability_idx);
		myfree(SRAMEX, xi);
		myfree(SRAMEX, yi);
		myfree(SRAMEX, CNN_data);
		myfree(SRAMEX, Full_output1);
		myfree(SRAMEX, Matrix_data);
    }
    else printf("No\r\n");/*如果没放电打印“No”*/
    printf("maxvalue is:%f\t",maxvalue);
    printf("meanvalue is:%f\t",meanvalue);
    printf("std_devvalue is:%f\r\n",std_devvalue);
	isrun = 0;
}

