/*
** 计算激活函数对加权输入的导数，并乘以delta，得到当前层最终的delta（敏感度图）
** 输入：x    当前层的所有输出（维度为l.batch * l.out_c * l.out_w * l.out_h）
**       n    l.output的维度，即为l.batch * l.out_c * l.out_w * l.out_h（包含整个batch的）
**       ACTIVATION    激活函数类型
**       delta     当前层敏感度图（与当前成输出x维度一样）
** 说明1：该函数不但计算了激活函数对于加权输入的导数，还将该导数乘以了之前完成大部分计算的敏感度图delta（对应元素相乘），因此调用改函数之后，将得到该层最终的敏感度图
** 说明2：这里直接利用输出值求激活函数关于输入的导数值是因为神经网络中所使用的绝大部分激活函数，其关于输入的导数值都可以描述为输出值的函数表达式，
          比如对于Sigmoid激活函数（记作f(x)），其导数值为f(x)'=f(x)*(1-f(x)),因此如果给出y=f(x)，那么f(x)'=y*(1-y)，只需要输出值y就可以了，不需要输入x的值，
          （暂时不确定darknet中有没有使用特殊的激活函数，以致于必须要输入值才能够求出导数值，在activiation.c文件中，有几个激活函数暂时没看懂，也没在网上查到）。
** 说明3：关于l.delta的初值，可能你有注意到在看某一类型网络层的时候，比如卷积层中的backward_convolutional_layer()函数，没有发现在此之前对l.delta赋初值的语句，
**        只是用calloc为其动态分配了内存，这样的l.delta其所有元素的值都为0,那么这里使用*=运算符得到的值将恒为0。是的，如果只看某一层，或者说某一类型的层，的确有这个疑惑，
**        但是整个网络是有很多层的，且有多种类型，一般来说，不会以卷积层为最后一层，而回以COST或者REGION为最后一层，这些层中，会对l.delta赋初值，又由于l.delta是由后
**        网前逐层传播的，因此，当反向运行到某一层时，l.delta的值将都不会为0.
*/
void gradient_array(const float *x, const int n, const ACTIVATION a, float *delta)
{
    int i;
    #pragma omp parallel for
    for(i = 0; i < n; ++i){
        delta[i] *= gradient(x[i], a);
    }
}

/*
** 根据不同的激活函数求取对输入的梯度（导数）
** 输入：x    激活函数接收的输入值
**       a    激活函数类型，包括的激活函数类型见activations.h中枚举类型ACTIVATION的定义
** 输出：激活函数关于输入x的导数值
*/
float gradient(float x, ACTIVATION a)
{
	// 以下分别求取各种激活函数对输入的导数值，详见各个导数求取函数的内部注释
    switch(a){
        case LINEAR:
            return linear_gradient(x);
        case LOGISTIC:
            return logistic_gradient(x);
        case LOGGY:
            return loggy_gradient(x);
        case RELU:
            return relu_gradient(x);
        case NORM_CHAN:
            //return relu_gradient(x);
        case NORM_CHAN_SOFTMAX_MAXVAL:
            //...
        case NORM_CHAN_SOFTMAX:
            printf(" Error: should be used custom NORM_CHAN or NORM_CHAN_SOFTMAX-function for gradient \n");
            exit(0);
            return 0;
        case ELU:
            return elu_gradient(x);
        case SELU:
            return selu_gradient(x);
        case RELIE:
            return relie_gradient(x);
        case RAMP:
            return ramp_gradient(x);
        case LEAKY:
            return leaky_gradient(x);
        case TANH:
            return tanh_gradient(x);
        case PLSE:
            return plse_gradient(x);
        case STAIR:
            return stair_gradient(x);
        case HARDTAN:
            return hardtan_gradient(x);
        case LHTAN:
            return lhtan_gradient(x);
    }
    return 0;
}

/*
** 计算每个卷积核的偏置更新值，所谓偏置更新值，就是bias = bias - alpha * bias_update中的bias_update
** 输入：bias_updates     当前层所有偏置的更新值，维度为l.n（即当前层卷积核的个数）
**       delta            当前层的敏感度图（即l.delta）
**       batch            一个batch含有的图片张数（即l.batch）
**       n                当前层卷积核个数（即l.n）
**       k                当前层输入特征图尺寸（即l.out_w*l.out_h）
** 原理：当前层的敏感度图l.delta是误差函数对加权输入的导数，也就是偏置更新值，只是其中每l.out_w*l.out_h个元素都对应同一个
**      偏置，因此需要将其加起来，得到的和就是误差函数对当前层各偏置的导数（l.delta的维度为l.batch*l.n*l.out_h*l.out_w,
**      可理解成共有l.batch行，每行有l.n*l.out_h*l.out_w列，而这一大行又可以理解成有l.n，l.out_h*l.out_w列，这每一小行就
**      对应同一个卷积核也即同一个偏置）
*/
void backward_bias(float *bias_updates, float *delta, int batch, int n, int size)
{
    int i,b;
	// 遍历batch中每张输入图片
    // 注意，最后的偏置更新值是所有输入图片的总和（多张图片无非就是重复一张图片的操作，求和即可）。
    // 总之：一个卷积核对应一个偏置更新值，该偏置更新值等于batch中所有输入图片累积的偏置更新值，
    // 而每张图片也需要进行偏置更新值求和（因为每个卷积核在每张图片多个位置做了卷积运算，这都对偏置更新值有贡献）以得到每张图片的总偏置更新值。
    for(b = 0; b < batch; ++b){
        for(i = 0; i < n; ++i){
            bias_updates[i] += sum_array(delta+size*(i+b*n), size);
        }
    }
}

/*
** 卷积神经网络反向传播核心函数
** 主要流程：1） 调用gradient_array()计算当前层l所有输出元素关于加权输入的导数值（也即激活函数关于输入的导数值），
**             并乘上上一次调用backward_convolutional_layer()还没计算完的l.delta，得到当前层最终的敏感度图；
**          2） 如果网络进行了BN，则backward_batchnorm_layer。
**          3） 如果网络没有进行BN，则直接调用 backward_bias()计算当前层所有卷积核的偏置更新值；
**          4） 依次调用im2col_cpu()，gemm_nt()函数计算当前层权重系数更新值；
**          5） 如果上一层的delta已经动态分配了内存，则依次调用gemm_tn(), col2im_cpu()计算上一层的敏感度图（并未完成所有计算，还差一个步骤）；
** 强调：每次调用本函数会计算完成当前层的敏感度计算，同时计算当前层的偏置、权重更新值，除此之外，还会计算上一层的敏感度图，但是要注意的是，
**      并没有完全计算完，还差一步：乘上激活函数对加权输入的导数值。这一步在下一次调用本函数时完成。
*/
void backward_convolutional_layer(convolutional_layer l, network_state state)
{
    int i, j;
	// 卷积核个数，考虑到分组卷积
    int m = l.n / l.groups;
	// 每一个卷积核元素个数（包括l.c（l.c为该层网络接受的输入图片的通道数）个通道上的卷积核元素个数总数，比如卷积核尺寸为3*3,
    // 输入图片有3个通道，因为要同时作用于输入的3个通道上，所以实际上这个卷积核是一个立体的，共有3*3*3=27个元素，这些元素都是要训练的参数），同样需要考虑分组数
    int n = l.size*l.size*l.c / l.groups;
	// 每张输出特征图的元素个数：out_w，out_h是输出特征图的宽高
    int k = l.out_w*l.out_h;

	// 计算当前层激活函数对加权输入的导数值并乘以l.delta相应元素，从而彻底完成当前层敏感度图的计算，得到当前层的敏感度图l.delta。
    // l.output存储了该层网络的所有输出：该层网络接受一个batch的输入图片，其中每张图片经卷积处理后得到的特征图尺寸为：l.out_w,l.out_h，
    // 该层卷积网络共有l.n个卷积核，因此一张输入图片共输出l.n张宽高为l.out_w,l.out_h的特征图（l.output为一张图所有输出特征图的总元素个数），
    // 所以所有输入图片也即l.output中的总元素个数为：l.n*l.out_w*l.out_h*l.batch；
    // l.activation为该卷积层的激活函数类型，l.delta就是gradient_array()函数计算得到的l.output中每一个元素关于激活函数函数输入的导数值，
    // 注意，这里直接利用输出值求得激活函数关于输入的导数值是因为神经网络中所使用的绝大部分激活函数关于输入的导数值都可以描述为输出值的函数表达式，
    // 比如对于Sigmoid激活函数（记作f(x)），其导数值为f(x)'=f(x)*(1-f(x)),因此如果给出y=f(x)，那么f(x)'=y*(1-y)，只需要输出值y就可以了，不需要输入x的值，
    // （暂时不确定darknet中有没有使用特殊的激活函数，以致于必须要输入值才能够求出导数值，在activiation.c文件中，有几个激活函数暂时没看懂，也没在网上查到）。
    // l.delta是一个一维数组，长度为l.batch * l.outputs（其中l.outputs = l.out_h * l.out_w * l.out_c），在make_convolutional_layer()动态分配内存；
    // 再强调一次：gradient_array()不单单是完成激活函数对输入的求导运算，还完成计算当前层敏感度图的最后一步：l.delta中每个元素乘以激活函数对输入的导数（注意gradient_arry中使用的是*=运算符）。
    // 每次调用backward_convolutional_laye时，都会完成当前层敏感度图的计算，同时会计算上一层的敏感度图，但对于上一层，其敏感度图并没有完全计算完成，还差一步，
    // 需要等到下一次调用backward_convolutional_layer()时来完成，诚如col2im_cpu()中注释一样。
    if (l.activation == SWISH) gradient_array_swish(l.output, l.outputs*l.batch, l.activation_input, l.delta);
    else if (l.activation == MISH) gradient_array_mish(l.outputs*l.batch, l.activation_input, l.delta);
    else if (l.activation == NORM_CHAN_SOFTMAX || l.activation == NORM_CHAN_SOFTMAX_MAXVAL) gradient_array_normalize_channels_softmax(l.output, l.outputs*l.batch, l.batch, l.out_c, l.out_w*l.out_h, l.delta);
    else if (l.activation == NORM_CHAN) gradient_array_normalize_channels(l.output, l.outputs*l.batch, l.batch, l.out_c, l.out_w*l.out_h, l.delta);
    else gradient_array(l.output, l.outputs*l.batch, l.activation, l.delta);

    if (l.batch_normalize) {
		// 之后单独讲BN层的前向和反向传播
        backward_batchnorm_layer(l, state);
    }
    else {
		// 计算偏置的更新值：每个卷积核都有一个偏置，偏置的更新值也即误差函数对偏置的导数，这个导数的计算很简单，实际所有的导数已经求完了，都存储在l.delta中，
        // 接下来只需把l.delta中对应同一个卷积核的项加起来就可以（卷积核在图像上逐行逐列跨步移动做卷积，每个位置处都有一个输出，共有l.out_w*l.out_h个，
        // 这些输出都与同一个偏置关联，因此将l.delta中对应同一个卷积核的项加起来即得误差函数对这个偏置的导数）
        backward_bias(l.bias_updates, l.delta, l.batch, l.n, k);
    }
	// 遍历batch中的每张照片，对于l.delta来说，每张照片是分开存的，因此其维度会达到：l.batch*l.n*l.out_w*l.out_h，
    // 对于l.weights,l.weight_updates以及上面提到的l.bias,l.bias_updates，是将所有照片对应元素叠加起来
    // （循环的过程就是叠加的过程，注意gemm()这系列函数含有叠加效果，不是覆盖输入C的值，而是叠加到之前的C上），
    // 因此l.weights与l.weight_updates维度为l.n*l.size*l.size，l.bias与l.bias_updates的维度为l.h，都与l.batch无关
    for (i = 0; i < l.batch; ++i) {
        for (j = 0; j < l.groups; ++j) {
			float *a = l.delta + (i*l.groups + j)*m*k;
			// net.workspace的元素个数为所有层中最大的l.workspace_size（在make_convolutional_layer()计算得到workspace_size的大小，在parse_network_cfg()中动态分配内存，此值对应未使用gpu时的情况）,
			// net.workspace充当一个临时工作空间的作用，存储临时所需要的计算参数，比如每层单张图片重排后的结果（这些参数马上就会参与卷积运算），一旦用完，就会被马上更新（因此该变量的值的更新频率比较大）
            float *b = state.workspace;
			
            float *c = l.weight_updates + j*l.nweights / l.groups;
			// 进入本函数之前，在backward_network()函数中，已经将net.input赋值为prev.output，也即若当前层为第l层，net.input此时已经是第l-1层的输出
			// 注意反向传播从后往前来看
            float *im = state.input + (i*l.groups + j)* (l.c / l.groups)*l.h*l.w;

            //im2col_cpu(im, l.c / l.groups, l.h, l.w, l.size, l.stride, l.pad, b);
			// 下面两步：im2col_cpu()与gemm()是为了计算当前层的权重更新值（其实也就是误差函数对当前层权重的导数）
			// 将多通道二维图像net.input变成按一定存储规则排列的数组b，以方便、高效地进行矩阵（卷积）计算，详细查看该函数注释（比较复杂），
			// im2col_cpu_ext每次仅处理net.input（包含整个batch）中的一张输入图片（对于第一层，则就是读入的图片，对于之后的层，这些图片都是上一层的输出，通道数等于上一层卷积核个数）。
			// 最终重排的b为l.c * l.size * l.size行，l.out_h * l.out_w列。
			// 你会发现在前向forward_convolutional_layer()函数中，也为每层的输入进行了重排，但是很遗憾的是，并没有一个l.workspace把每一层的重排结果保存下来，而是统一存储到net.workspace中，
			// 并被不断擦除更新，那为什么不保存呢？保存下来不是省掉一大笔额外重复计算开销？原因有两个：1）net.workspace中只存储了一张输入图片的重排结果，所以重排下张图片时，马上就会被擦除，
			// 当然你可能会想，那为什么不弄一个l.worspaces将每层所有输入图片的结果保存呢？这引出第二个原因；2）计算成本是降低了，但存储空间需求急剧增加，想想每一层都有l.batch张图，且每张都是多通道的，
			// 重排后其元素个数还会增多，这个存储量搁谁都受不了，如果一个batch有128张图，输入图片尺寸为400*400，3通道，网络有16层（假设每层输入输出尺寸及通道数都一样），那么单单为了存储这些重排结果，
			// 就需要128*400*400*3*16*4/1024/1024/1024 = 3.66G，所以为了权衡，只能重复计算！
            im2col_cpu_ext(
                im,                 // input
                l.c / l.groups,     // input channels
                l.h, l.w,           // input size (h, w)
                l.size, l.size,     // kernel size (h, w)
                l.pad, l.pad,       // padding (h, w)
                l.stride_y, l.stride_x, // stride (h, w)
                l.dilation, l.dilation, // dilation (h, w)
                b);                 // output
			// 下面计算当前层的权重更新值，所谓权重更新值就是weight = weight - alpha * weight_update中的weight_update，
			// 权重更新值等于当前层敏感度图中每个元素乘以相应的像素值，因为一个权重跟当前层多个输出有关联（权值共享，即卷积核在图像中跨步移动做卷积，每个位置卷积得到的值
			// 都与该权值相关），所以对每一个权重更新值来说，需要在l.delta中找出所有与之相关的敏感度，乘以相应像素值，再求和，具体实现的方式依靠im2col_cpu()与gemm_nt()完成。
			// （backward_convolutional_layer整个函数的代码非常重要，仅靠文字没有公式与图表辅助说明可能很难说清，所以这部分更为清晰详细的说明，请参考个人博客！）
			// GEneral Matrix to Matrix Multiplication
			// 此处在im2col_cpu操作基础上，利用矩阵乘法c=alpha*a*b+beta*c完成对图像卷积的操作；
			// 0表示不对输入a进行转置，1表示对输入b进行转置；
			// m是输入a,c的行数，具体含义为卷积核的个数(l.n)；
			// n是输入b,c的列数，具体含义为每个卷积核元素个数乘以输入图像的通道数(l.size*l.size*l.c)；
			// k是输入a的列数也是b的行数，具体含义为每个输出特征图的元素个数（l.out_w*l.out_h）；
			// a,b,c即为三个参与运算的矩阵（用一维数组存储）,alpha=beta=1为常系数；
			// a为l.delta的一大行。l.delta为本层所有输出元素（包含整个batch中每张图片的所有输出特征图）关于加权输入的导数（即激活函数的导数值）集合,
			// 元素个数为l.batch * l.out_h * l.out_w * l.out_c（l.out_c = l.n），按行存储，共有l.batch行，l.out_c * l.out_h * l.out_w列，
			// 即l.delta中每行包含一张图的所有输出图，故这么一大行，又可以视作有l.out_c（l.out_c=l.n）小行，l.out_h*l*out_w小列，而一次循环就是处理l.delta的一大行，
			// 故可以将a视作l.out_c行，l.out_h*l*out_w列的矩阵；
			// b为单张输入图像经过im2col_cpu重排后的图像数据；
			// c为输出，按行存储，可视作有l.n行，l.c*l.size*l.size列（l.c是输入图像的通道数，l.n是卷积核个数），
			// 即c就是所谓的误差项（输出关于加权输入的导数），或者敏感度（强烈推荐：https://www.zybuluo.com/hanbingtao/note/485480）（一个核有l.c*l.size*l.size个权重，共有l.n个核）。
			// 由上可知：
			// a: (l.out_c) * (l.out_h*l*out_w)
			// b: (l.c * l.size * l.size) * (l.out_h * l.out_w)
			// c: (l.n) * (l.c*l.size*l.size)（注意：l.n = l.out_c）
			// 故要进行a * b + c计算，必须对b进行转置（否则行列不匹配），因故调用gemm_nt()函数
            gemm(0, 1, m, n, k, 1, a, k, b, k, 1, c, n);

			// 接下来，用当前层的敏感度图l.delta以及权重l.weights（还未更新）来获取上一层网络的敏感度图，BP算法的主要流程就是依靠这种层与层之间敏感度反向递推传播关系来实现。
			// 在network.c的backward_network()中，会从最后一层网络往前遍循环历至第一层，而每次开始遍历某一层网络之前，都会更新net.input为这一层网络前一层的输出，即prev.output,
			// 同时更新net.delta为prev.delta，因此，这里的net.delta是当前层前一层的敏感度图。
			// 已经强调很多次了，再说一次：下面得到的上一层的敏感度并不完整，完整的敏感度图是损失函数对上一层的加权输入的导数，
			// 而这里得到的敏感度图是损失函数对上一层输出值的导数，还差乘以一个输出值也即激活函数对加权输入的导数。
            if (state.delta) {
				// 当前层还未更新的权重
                a = l.weights + j*l.nweights / l.groups;
				
				// 每次循环仅处理一张输入图，注意移位（l.delta的维度为l.batch * l.out_c * l.out_w * l.out_h）（注意l.n = l.out_c，另外提一下，对整个网络来说，每一层的l.batch其实都是一样的）
                b = l.delta + (i*l.groups + j)*m*k;
				
				// net.workspace和上面一样，还是一张输入图片的重排，不同的是，此处我们只需要这个容器，而里面存储的值我们并不需要，在后面的处理过程中，
				// 会将其中存储的值一一覆盖掉（尺寸维持不变，还是(l.c * l.size * l.size) * (l.out_h * l.out_w）
                c = state.workspace;
				
				 // 相比上一个gemm，此处的a对应上一个的c,b对应上一个的a，c对应上一个的b，即此处a,b,c的行列分别为：
				// a: (l.n) * (l.c*l.size*l.size)，表示当前层所有权重系数
				// b: (l.out_c) * (l.out_h*l*out_w)（注意：l.n = l.out_c），表示当前层的敏感度图
				// c: (l.c * l.size * l.size) * (l.out_h * l.out_w)，表示上一层的敏感度图（其元素个数等于上一层网络单张输入图片的所有输出元素个数），
				// 此时要完成a * b + c计算，必须对a进行转置（否则行列不匹配），因故调用gemm_tn()函数。
				// 此操作含义是用：用当前层还未更新的权重值对敏感度图做卷积，得到包含上一层所有敏感度信息的矩阵，但这不是上一层最终的敏感度图，
				// 因为此时的c，也即net.workspace的尺寸为(l.c * l.size * l.size) * (l.out_h * l.out_w)，明显不是上一层的输出尺寸l.c*l.w*l.h，
				// 接下来还需要调用col2im_cpu()函数将其恢复至l.c*l.w*l.h（可视为l.c行，l.w*l.h列），这才是上一层的敏感度图（实际还差一个环节，
				// 这个环节需要等到下一次调用backward_convolutional_layer()才完成：将net.delta中每个元素乘以激活函数对加权输入的导数值）。
				// 完成gemm这一步，如col2im_cpu()中注释，是考虑了多个卷积核导致的一对多关系（上一层的一个输出元素会流入到下一层多个输出元素中），
				// 接下来调用col2im_cpu()则是考虑卷积核重叠（步长较小）导致的一对多关系。
                gemm(1, 0, n, k, m, 1, a, n, b, k, 0, c, k);

                //col2im_cpu(state.workspace, l.c / l.groups, l.h, l.w, l.size, l.stride,
                //     l.pad, state.delta + (i*l.groups + j)*l.c / l.groups*l.h*l.w);
				
				// 对c也即state.workspace进行重排，得到的结果存储在state.delta中，每次循环只会处理一张输入图片，因此，此处只会得到一张输入图产生的敏感图（注意net.delta的移位）,
				// 整个循环结束后，net.delta的总尺寸为l.batch * l.h * l.w * l.c，这就是上一层网络整个batch的敏感度图，可视为有l.batch行，l.h*l.w*l.c列，
				// 每行存储了一张输入图片所有输出特征图的敏感度
				// col2im_cpu()函数中会调用col2im_add_pixel()函数，该函数中使用了+=运算符，也即该函数要求输入的net.delta的初始值为0,而在gradient_array()中注释到l.delta的元素是不为0（也不能为0）的，
				// 看上去是矛盾的，实则不然，gradient_array()使用的l.delta是当前层的敏感度图，而在col2im_cpu()使用的net.delta是上一层的敏感度图，正如gradient_array()中所注释的，
				// 当前层l.delta之所以不为0,是因为从后面层反向传播过来的，对于上一层，显然还没有反向传播到那，因此net.delta的初始值都是为0的（注意，每一层在构建时，就为其delta动态分配了内存，
				// 且在前向传播时，为每一层的delta都赋值为0,可以参考network.c中forward_network()函数）
                col2im_cpu_ext(
                    state.workspace,        // input
                    l.c / l.groups,         // input channels (h, w)
                    l.h, l.w,               // input size (h, w)
                    l.size, l.size,         // kernel size (h, w)
                    l.pad, l.pad,           // padding (h, w)
                    l.stride_y, l.stride_x,     // stride (h, w)
                    l.dilation, l.dilation, // dilation (h, w)
                    state.delta + (i*l.groups + j)* (l.c / l.groups)*l.h*l.w); // output (delta)
            }
        }
    }
}

// 注释来自https://github.com/hgpvision/darknet/blob/master/src/col2im.c
/*
*  将输入图像im的channel通道上的第row行，col列像素灰度值加上val（直接修改im的值，因此im相当于是返回值）
** 输入：im         输入图像
**       channels   输入图像的im通道数（这个参数没用。。。）
**       height     输入图像im的高度（行）
**       width      输入图像im的宽度（列）
**       row        需要加上val的像素所在的行数（补零之后的行数，因此需要先减去pad才能得到真正在im中的行数）
**       col        需要加上val的像素所在的列数（补零之后的列数，因此需要先减去pad才能得到真正在im中的列数）
**       channel    需要加上val的像素所在的通道数
**       pad        四周补0长度
**       val        像素灰度添加值
*/
void col2im_add_pixel(float *im, int height, int width, int channels,
                        int row, int col, int channel, int pad, float val)
{
    row -= pad;
    col -= pad;

    if (row < 0 || col < 0 ||
        row >= height || col >= width) return;
    im[col + width*(row + height*channel)] += val;
}
//This one might be too, can't remember.
// 注释来自：https://github.com/hgpvision/darknet/blob/master/src/col2im.c

/*
** 此函数与im2col_cpu()函数的流程相反，目地是将im2col_cpu()函数重排得到的图片data_col恢复至正常的图像矩阵排列，并与data_im相加，最终data_im相当于是输出值，
** 要注意的是，data_im的尺寸是在函数外确定的，且并没有显示的将data_col转为一个与data_im尺寸相同的矩阵，而是将其中元素直接加在data_im对应元素上（data_im初始所有元素值都为0）。
** 得到的data_im尺寸为l.c*l.h*l.w，即为当前层的输入图像尺寸，上一层的输出图像尺寸，按行存储，可视为l.c行，l.h*l.w列，即其中每行对应一张输出特征图的敏感度图（实际上这还不是最终的敏感度，
** 还差一个环节：乘以激活函数对加权输入的导数，这将在下一次调用backward_convolutional_laye时完成）。
**
** 举个例子：第L-1层每张输入图片（本例子只分析单张输入图片）的输出为5*5*3（3为输出通道数），第L层共有2个卷积核，每个卷积核的尺寸为3*3，stride = 2,
**         第L-1层的输出是第L层的输入，第L层的2个卷积核同时对上一层3个通道的输出做卷积运算，为了做到这一点，需要调用im2col_cpu()函数将
**         上一层的输出，也就是本层的输入重排为27行4列的图，也就是由5*5*3变换至27*4，你会发现总的元素个数变多了（75增多到了98），
**         这是因为卷积核stride=2,小于卷积核的尺寸3,因此卷积在两个连续位置做卷积，会有重叠部分，而im2col_cpu()函数为了便于卷积运算，完全将其
**         铺排开来，并没有在空间上避免重复元素，因此像素元素会增多。此外，之所以是27行，是因为卷积核尺寸为3*3，而上一层的输出即本层输入有3个通道，
**         为了同时给3个通道做卷积运算，需要将3个通道上的输入一起考虑，即得到3*3*3行，4列是因为对于对于5*5的图像，使用3*3的卷积核，stride=2的卷积跨度，
**         最终会得到2*2的特征图，也就是4个元素。除了调用im2col_cpu()对输入图像做重排，相应的，也要将所有卷积核重排成一个2*27的矩阵，为什么是2呢？
**         因为有两个卷积核，为了做到同时将两个卷积核作用到输入图像上，需要将两个核合到一个矩阵中，每个核对应一行，因此有2行，那为什么是27呢？每个核
**         元素个数不是3*3=9吗？是的，但是考虑到要同时作用到3个通道上，所以实际一个卷积核有9*3=27个元素。综述，得到2*27的卷积核矩阵与27*4的输入图像矩阵，
**         两个矩阵相乘，即可完成将2个卷积核同时作用于3通道的输入图像上（非常方便，不枉前面非这么大劲的重排！），最终得到2*4的矩阵，这2*4矩阵又代表这什么呢？
**         2代表这有两个输出图（对应2个卷积核，即l.out_c=2），每个输出图占一行，4代表这每个输出图元素有4个（前面说了，每个卷积核会得到2*2的特征图，即l.out_h=l.out_w=2）。这个例子说到这，只说完了
**         前向传播部分，可以看出im2col_cpu()这个函数的重要性。而此处的col2im_cpu()是一个逆过程，主要用于反向传播中，由L层的敏感度图(sensitivity map，
**         可能每个地方叫的不一样，此处参考博客：https://www.zybuluo.com/hanbingtao/note/485480)反向求得第L-1层的敏感度图。顺承上面的例子，第L-1层的输出
**         是一个5*5*3（l.w=l.h=5,l.c=3）的矩阵，也就是敏感度图的维度为5*5*3（每个输出元素，对应一个敏感度值），第L层的输出是一个2*4的矩阵，敏感度图的维度为2*4，假设已经计算得到了
**         第L层的2*4的敏感度图，那么现在的问题是，如何由第L层的2*4敏感度图以及2个卷积核（2*27）反向获取第L-1层的敏感度图呢？上面给的博客链接给出了一种很好的求解方式，
**         但darknet并不是这样做的，为什么？因为前面有im2col_cpu()，im2col_cpu()函数中的重排方式，使得我们不再需要博客中提到的将sensitivity map还原为步长为1的sensitivity map，
**         只需再使用col2im_cpu()就可以了！过程是怎样的呢，看backward_convolutional_layer()函数中if(net.delta)中的语句就知道了，此处仅讨论col2im_cpu()的过程，
**         在backward_convolutional_layer()已经得到了data_col，这个矩阵含有了所有的第L-1层敏感度的信息，但遗憾的是，不能直接用，需要整理，因为此时的data_col还是一个
**         27*4的矩阵，而我们知道第L-1层的敏感度图是一个5*5*3的矩阵，如何将一个27*4的矩阵变换至一个5*5*3的矩阵是本函数要完成的工作，前面说到27*4元素个数多于5*5*3,
**         很显然要从27*4变换至5*5*3，肯定会将某些元素相加合并（下面col2im_add_pixel()函数就是干这个的），具体怎样，先不说，先来看看输入参数都代表什么意思吧：
** 输入：data_col    backward_convolutional_layer()中计算得到的包含上一层所有敏感度信息的矩阵，行数为l.n*l.size*l.size（l代表本层/当前层），列数为l.out_h*l.out_w（对于本例子，行数为27,列数为4,上一层为第L-1层，本层是第L层）
**       channels    当前层输入图像的通道数（对于本例子，为3）
**       height      当前层输入图像的行数（对于本例子，为5）
**       width       当前层输入图像的列数（对于本例子，为5）
**       ksize       当前层卷积核尺寸（对于本例子，为3）
**       stride      当前层卷积跨度（对于本例子，为2）
**       pad         当前层对输入图像做卷积时四周补0的长度
**       data_im     经col2im_cpu()重排恢复之后得到的输出矩阵，也即上一层的敏感度图，尺寸为l.c * l.h * l.w（刚好就是上一层的输出当前层输入的尺寸，对于本例子，5行5列3通道），
**                   注意data_im的尺寸，是在本函数之外就已经确定的，不是在本函数内部计算出来的，这与im2col_cpu()不同，im2col_cpu()计算得到的data_col的尺寸都是在函数内部计算得到的，
**                   并不是事先指定的。也就是说，col2im_cpu()函数完成的是指定尺寸的输入矩阵往指定尺寸的输出矩阵的转换。
** 原理：原理比较复杂，很难用文字叙述，博客：https://www.zybuluo.com/hanbingtao/note/485480中基本原理说得很详细了，但是此处的实现与博客中并不一样，所以具体实现的原理此处简要叙述一下，具体见个人博客。
**      第L-1层得到l.h*l.w*l.c输出，也是第L层的输入，经L层卷积及激活函数处理之后，得到l.out_h*l.out_w*l.out_c的输出，也就是由l.h*l.w*l.c-->l.out_h*l.out_w*l.out_c，
**      由于第L层有多个卷积核，所以第L-1层中的一个输出元素会流入到第L层多个输出中，除此之外，由于卷积核之间的重叠，也导致部分元素流入到第L层的多个输出中，这两种情况，都导致第L-1层中的某个敏感度会与第L层多个输出有关，
**      为清晰，还是用上面的例子来解释，第L-1层得到5*5*3(3*25)的输出，第L层得到2*2*2（2*4）的输出，在backward_convolutional_layer()已经计算得到的data_col实际是27*2矩阵与2*4矩阵相乘的结果，
**      为方便，我们记27*2的矩阵为a，记2*4矩阵为b，那么a中一行（2个元素）与b中一列（2个元素）相乘对应这什么呢？对应第一情况，因为有两个卷积核，使得L-1中一个输出至少与L层中两个输出有关系，经此矩阵相乘，得到27*4的矩阵，
**      已经考虑了第一种情况（27*4这个矩阵中的每一个元素都是两个卷积核影响结果的求和），那么接下来的就是要考虑第二种情况：卷积核重叠导致的一对多关系，具体做法就是将data_col中对应相同像素的值相加，这是由
**      im2col_cpu()函数决定的（可以配合im2col_cpu()来理解），因为im2col_cpu()将这些重叠元素也铺陈保存在data_col中，所以接下来，只要按照im2col_cpu()逆向将这些重叠元素的影响叠加就可以了，
**      大致就是这个思路，具体的实现细节可能得见个人博客了（这段写的有点罗嗦～）。
**
*/
void col2im_cpu(float* data_col,
         int channels,  int height,  int width,
         int ksize,  int stride, int pad, float* data_im)
{
    int c,h,w;
	// 当前层输出图的尺寸（对于上面的例子，height_col=2,width_col=2）
    int height_col = (height + 2*pad - ksize) / stride + 1;
    int width_col = (width + 2*pad - ksize) / stride + 1;
	// 当前层每个卷积核在所有输入图像通道上的总元素个数（对于上面的例子，channels_col=3*3*3=27）
    // 注意channels_col实际是data_col的行数
    int channels_col = channels * ksize * ksize;
	// 开始遍历：外循环遍历data_col的每一行（对于上面的例子，data_col共27行）
    for (c = 0; c < channels_col; ++c) {
		// 列偏移，卷积核是一个二维矩阵，并按行存储在一维数组中，利用求余运算获取对应在卷积核中的列数，比如对于
        // 3*3的卷积核，当c=0时，显然在第一列，当c=5时，显然在第2列，当c=9时，在第二通道上的卷积核的第一列
        int w_offset = c % ksize;
		// 行偏移，卷积核是一个二维的矩阵，且是按行（卷积核所有行并成一行）存储在一维数组中的，
        // 比如对于3*3的卷积核，处理3通道的图像，那么一个卷积核具有27个元素，每9个元素对应一个通道上的卷积核（互为一样），
        // 每当c为3的倍数，就意味着卷积核换了一行，h_offset取值为0,1,2
        int h_offset = (c / ksize) % ksize;
		// 通道偏移，channels_col是多通道的卷积核并在一起的，比如对于3通道，3*3卷积核，每过9个元素就要换一通道数，
        // 当c=0~8时，c_im=0;c=9~17时，c_im=1;c=18~26时，c_im=2
        // c_im是data_im的通道数（即上一层输出当前层输入的通道数），对于上面的例子，c_im取值为0,1,2
        int c_im = c / ksize / ksize;
		// 中循环与内循环和起来刚好遍历data_col的每一行（对于上面的例子，data_col的列数为4,height_col*width_col=4）
        for (h = 0; h < height_col; ++h) {
            for (w = 0; w < width_col; ++w) {
				// 获取在输出data_im中的行数im_row与列数im_col
                // 由上面可知，对于3*3的卷积核，h_offset取值为0,1,2,当h_offset=0时，会提取出所有与卷积核第一行元素进行运算的像素，
                // 依次类推；加上h*stride是对卷积核进行行移位操作，比如卷积核从图像(0,0)位置开始做卷积，那么最先开始涉及(0,0)~(3,3)
                // 之间的像素值，若stride=2，那么卷积核进行行移位一次时，下一行的卷积操作是从元素(2,0)（2为图像行号，0为列号）开始
                int im_row = h_offset + h * stride;
				// 对于3*3的卷积核，w_offset取值也为0,1,2，当w_offset取1时，会提取出所有与卷积核中第2列元素进行运算的像素，
                // 实际在做卷积操作时，卷积核对图像逐行扫描做卷积，加上w*stride就是为了做列移位，
                // 比如前一次卷积其实像素元素为(0,0)，若stride=2,那么下次卷积元素起始像素位置为(0,2)（0为行号，2为列号）
                int im_col = w_offset + w * stride;
				// 计算在输出data_im中的索引号
                // 对于上面的例子，im_row的取值范围为0~4,im_col从0~4，c从0~2（其中h_offset从0~2,w_offset从0~2, h从0~1,w从0~1）
                // 输出的data_im的尺寸为l.c * l.h * lw，对于上面的例子，为3*5*5,因此，im_row,im_col,c的取值范围刚好填满data_im

                // 获取data_col中索引为col_index的元素，对于上面的例子，data_col为27*4行，按行存储
                // col_index = c * height_col * width_col + h * width_col + w逐行读取data_col中的每一个元素。
                // 相同的im_row,im_col与c_im可能会对应多个不同的col_index，这就是卷积核重叠带来的影响，处理的方式是将这些val都加起来，
                // 存在data_im的第im_row - pad行，第im_col - pad列（c_im通道上）中。
                // 比如上面的例子，上面的例子，如果固定im_row = 0, im_col =2, c_im =0，由c_im = 0可以知道c在0~8之间，由im_row=0,可以确定h = 0, h_offset =0，
                // 可以得到两组：1)w_offset = 0, w = 1; 2) w_offset = 2, w =0，第一组，则可以完全定下：c=0,h=0,w=1，此时col_index=1，由第二组，可完全定下：c=2,h=0,w=0，
                // 此时col_index = 2*2*2=8
                int col_index = (c * height_col + h) * width_col + w;
                float val = data_col[col_index];
				// 从data_im找出c_im通道上第im_row - pad行im_col - pad列处的像素，使其加上val
                // height, width, channels都是上一层输出即当前层输入图像的尺寸，也是data_im的尺寸（对于本例子，三者的值分别为5,5,3）,
                // im_row - pad,im_col - pad,c_im都是某一具体元素在data_im中的行数与列数与通道数（因为im_row与im_col是根据卷积过程计算的，
                // 所以im_col和im_row中实际还包含了补零长度pad，需要减去之后，才是原本的没有补零矩阵data_im中的行列号）
                col2im_add_pixel(data_im, height, width, channels,
                        im_row, im_col, c_im, pad, val);
            }
        }
    }
}