yolov3 -tiny 网络实现和源码分析

摘自:https://blog.csdn.net/alangaixiaoxiao/article/details/105533746

相关推荐

网上的资源很多,也有很多博主的原理讲解。在这里推荐几个资源:
darknet官网:https://pjreddie.com/darknet/yolo/ linux系统
github上基于VS的工程:https://github.com/AlexeyAB/darknet
github上带有注释的工程:https://github.com/hgpvision/darknet

yolov3-tiny 原理

Yolo算法采用一个单独的CNN模型实现end-to-end的目标检测,首先将输入图片resize到448×448,然后送入CNN网络,最后处理网络预测结果得到检测的目标。
YOLO 的核心思想就是利用整张图作为网络的输入,直接在输出层回归 bounding box(边界框) 的位置及其所属的类别。将一幅图像分成 SxS 个网格(grid cell),如果某个 object 的中心落在这个网格中,则这个网格就负责预测这个 object。

每个 bounding box 要预测 (x, y, w, h) 和 confidence 共5个值,每个网格还要预测一个类别信息,记为 C 类。则 SxS个 网格,每个网格要预测 B 个 bounding box, 每个box中都有 C 个 classes对应的概率值。输出就是 S x S x B x(5+C) 的一个 tensor。

注意:class 信息是针对每个网格的,confidence 信息是针对每个 bounding box 的。

yolov3-tiny中,共有两个输出层(yolo层),分别为13×13和26×26,每个网格可以预测3个bounding box,共有80个分类数。所以最后的yolo层的尺寸为:13x13x255和26x26x255。
yolov3-tiny网络层结构如下:

可以看出,yolov3-tiny共有23层网络,其中包含五种不同的网络层:卷积层convolutional(13个),池化层maxpool(6个),路由层route(2个),上采样层upsample(1个),输出层yolo(2个)。Yolov3-tiny中,除了Yolo层之前的那个卷积层,每个卷积层之后都有BN层,且每个卷积层之后都有激活函数LEAKY(yolo层之前是linear)。

yolov3-tiny 源码分析

配置网络结构

yolov3-tiny前向传播主要在detector.c中的test_detector函数中完成:

/** 本函数是检测模型的一个前向推理测试函数.
* @param datacfg       数据集信息文件路径(也即cfg/*.data文件),文件中包含有关数据集的信息,比如cfg/coco.data
* @param cfgfile       网络配置文件路径(也即cfg/*.cfg文件),包含一个网络所有的结构参数,比如cfg/yolo.cfg
* @param weightfile    已经训练好的网络权重文件路径,比如darknet网站上下载的yolo.weights文件
* @param filename      待进行检测的图片路径(单张图片)
* @param thresh        阈值,类别检测概率大于该阈值才认为其检测结果有效
* @param hier_thresh
* @param outfile
* @param fullscreen
* @details 该函数为一个前向推理测试函数,不包括训练过程,因此如果要使用该函数,必须提前训练好网络,并加载训练好的网络参数文件,
*          这些文件可以在作者网站上根据作者的提示下载到。本函数由darknet.c中的主函数调用,严格来说,本文件不应纳入darknet网络结构文件夹中,
*          其只是一个测试文件,或者说是一个example,应该放入到example文件夹中(新版的darknet已经这样做了,可以在github上查看)。
*          本函数的流程为:.
*/
void test_detector(char *datacfg, char *cfgfile, char *weightfile, char *filename, float thresh,
    float hier_thresh, int dont_show, int ext_output, int save_labels, char *outfile, int letter_box, int benchmark_layers)
{
	// 从指定数据文件datacfg(.data文件)中读入数据信息(测试、训练数据信息)到options中
	// options是list类型数据,其中的node包含的void指针具体是kvp数据类型,具有键值和值(类似C++中的Map)
    list *options = read_data_cfg(datacfg);
	// 获取数据集的名称(包括路径),第二个参数"names"表明要从options中获取所用数据集的名称信息(如names = data/coco.names)
    char *name_list = option_find_str(options, "names", "data/names.list");
    int names_size = 0;
	// 从data/**.names中读取物体名称/标签信息
    char **names = get_labels_custom(name_list, &names_size); //get_labels(name_list);
	
    // 加载data/labels/文件夹中所有的字符标签图片
    image **alphabet = load_alphabet();

    network net = parse_network_cfg_custom(cfgfile, 1, 1); // set batch=1  配置各网络层参数,重要

在parser.c中的parse_network_cfg_custom函数中,根据yolov3-tiny.cfg文件对网络结构进行配置,明确各层网络的类型、输入输出通道数、图像尺寸、卷积核大小等。

//配置各网络层参数
network parse_network_cfg_custom(char *filename, int batch, int time_steps)
{
	// 从神经网络结构参数文件中读入所有神经网络层的结构参数,存储到sections中,
	// sections的每个node包含一层神经网络的所有结构参数
    list *sections = read_cfg(filename);
	// 获取sections的第一个节点,可以查看一下cfg/***.cfg文件,其实第一块参数(以[net]开头)不是某层神经网络的参数,
	// 而是关于整个网络的一些通用参数,比如学习率,衰减率,输入图像宽高,batch大小等,
	// 具体的关于某个网络层的参数是从第二块开始的,如[convolutional],[maxpool]...,
	// 这些层并没有编号,只说明了层的属性,但层的参数都是按顺序在文件中排好的,读入时,
	// sections链表上的顺序就是文件中的排列顺序。
    node *n = sections->front;
    if(!n) error("Config file has no sections");
	// 创建网络结构并动态分配内存:输入网络层数为sections->size - 1,sections的第一段不是网络层,而是通用网络参数
    network net = make_network(sections->size - 1);
	// 所用显卡的卡号(gpu_index在cuda.c中用extern关键字声明)
	// 在调用parse_network_cfg()之前,使用了cuda_set_device()设置了gpu_index的值号为当前活跃GPU卡号
    net.gpu_index = gpu_index;
	// size_params结构体元素不含指针变量
    size_params params;

    if (batch > 0) params.train = 0;    // allocates memory for Detection only
    else params.train = 1;              // allocates memory for Detection & Training

    section *s = (section *)n->val;
    list *options = s->options;
    if(!is_network(s)) error("First section must be [net] or [network]");
    parse_net_options(options, &net);

#ifdef GPU
    printf("net.optimized_memory = %d \n", net.optimized_memory);
    if (net.optimized_memory >= 2 && params.train) {
        pre_allocate_pinned_memory((size_t)1024 * 1024 * 1024 * 8);   // pre-allocate 8 GB CPU-RAM for pinned memory
    }
#endif  // GPU

    params.h = net.h;
    params.w = net.w;
    params.c = net.c;
    params.inputs = net.inputs;
    if (batch > 0) net.batch = batch;
    if (time_steps > 0) net.time_steps = time_steps;
    if (net.batch < 1) net.batch = 1;
    if (net.time_steps < 1) net.time_steps = 1;
    if (net.batch < net.time_steps) net.batch = net.time_steps;
    params.batch = net.batch;
    params.time_steps = net.time_steps;
    params.net = net;
    printf("mini_batch = %d, batch = %d, time_steps = %d, train = %d \n", net.batch, net.batch * net.subdivisions, net.time_steps, params.train);

    int avg_outputs = 0;
    float bflops = 0;
    size_t workspace_size = 0;
    size_t max_inputs = 0;
    size_t max_outputs = 0;
    n = n->next;
    int count = 0;
    free_section(s);

	// 此处stderr不是错误提示,而是输出结果提示,提示网络结构
    fprintf(stderr, "   layer   filters  size/strd(dil)      input                output\n");
    while(n){
        params.index = count;
        fprintf(stderr, "%4d ", count);
        s = (section *)n->val;
        options = s->options;
		// 定义网络层
        layer l = { (LAYER_TYPE)0 };
		// 获取网络层的类别

        LAYER_TYPE lt = string_to_layer_type(s->type);
		
		//通过读取网络类型,从而配置各网络层的参数
        if(lt == CONVOLUTIONAL){//yolov3-tiny  卷积层  13层
            l = parse_convolutional(options, params);
        }else if(lt == LOCAL){
            l = parse_local(options, params);
        }else if(lt == ACTIVE){
            l = parse_activation(options, params);
        }else if(lt == RNN){
            l = parse_rnn(options, params);
        }else if(lt == GRU){
            l = parse_gru(options, params);
        }else if(lt == LSTM){
            l = parse_lstm(options, params);
        }else if (lt == CONV_LSTM) {
            l = parse_conv_lstm(options, params);
        }else if(lt == CRNN){
            l = parse_crnn(options, params);
        }else if(lt == CONNECTED){
            l = parse_connected(options, params);
        }else if(lt == CROP){
            l = parse_crop(options, params);
        }else if(lt == COST){
            l = parse_cost(options, params);
            l.keep_delta_gpu = 1;
        }else if(lt == REGION){
            l = parse_region(options, params);
            l.keep_delta_gpu = 1;
        }else if (lt == YOLO) {//yolov3-tiny YOLO层  两层
            l = parse_yolo(options, params);
            l.keep_delta_gpu = 1;
        }else if (lt == GAUSSIAN_YOLO) {
            l = parse_gaussian_yolo(options, params);
            l.keep_delta_gpu = 1;
        }else if(lt == DETECTION){
            l = parse_detection(options, params);
        }else if(lt == SOFTMAX){
            l = parse_softmax(options, params);
            net.hierarchy = l.softmax_tree;
            l.keep_delta_gpu = 1;
        }else if(lt == NORMALIZATION){
            l = parse_normalization(options, params);
        }else if(lt == BATCHNORM){
            l = parse_batchnorm(options, params);
        }else if(lt == MAXPOOL){//yolov3-tiny 池化层 maxpool  6层
            l = parse_maxpool(options, params);
        }else if (lt == LOCAL_AVGPOOL) {
            l = parse_local_avgpool(options, params);
        }else if(lt == REORG){
            l = parse_reorg(options, params);        }
        else if (lt == REORG_OLD) {
            l = parse_reorg_old(options, params);
        }else if(lt == AVGPOOL){
            l = parse_avgpool(options, params);
        }else if(lt == ROUTE){//yolov3-tiny 路由层 2层
            l = parse_route(options, params);
            int k;
            for (k = 0; k < l.n; ++k) {
                net.layers[l.input_layers[k]].use_bin_output = 0;
                net.layers[l.input_layers[k]].keep_delta_gpu = 1;
            }
        }else if (lt == UPSAMPLE) {//yolov3-tiny 上采样层 1层
            l = parse_upsample(options, params, net);
        }else if(lt == SHORTCUT){
            l = parse_shortcut(options, params, net);
            net.layers[count - 1].use_bin_output = 0;
            net.layers[l.index].use_bin_output = 0;
            net.layers[l.index].keep_delta_gpu = 1;
        }else if (lt == SCALE_CHANNELS) {
            l = parse_scale_channels(options, params, net);
            net.layers[count - 1].use_bin_output = 0;
            net.layers[l.index].use_bin_output = 0;
            net.layers[l.index].keep_delta_gpu = 1;
        }
        else if (lt == SAM) {
            l = parse_sam(options, params, net);
            net.layers[count - 1].use_bin_output = 0;
            net.layers[l.index].use_bin_output = 0;
            net.layers[l.index].keep_delta_gpu = 1;
        }else if(lt == DROPOUT){
            l = parse_dropout(options, params);
            l.output = net.layers[count-1].output;
            l.delta = net.layers[count-1].delta;
            .........

下载权重文件

在parser.c的load_weights_upto中,根据卷积层的网络配置,开始下载读取各层的权重文件。

//读取权重文件函数
void load_weights_upto(network *net, char *filename, int cutoff)//cutoff = net->n
{
#ifdef GPU
    if(net->gpu_index >= 0){
        cuda_set_device(net->gpu_index);
    }
#endif
    fprintf(stderr, "Loading weights from %s...\n", filename);
    fflush(stdout);
    FILE *fp = fopen(filename, "rb");
    if(!fp) file_error(filename);

    int major;
    int minor;
    int revision;
    fread(&major, sizeof(int), 1, fp);//读取一个4字节的数据
    fread(&minor, sizeof(int), 1, fp);//读取一个4字节的数据
    fread(&revision, sizeof(int), 1, fp);//读取一个4字节的数据
	printf("the size of int in x64 is %d bytes,attention!!!\n", sizeof(int));//x86 x64: 4
	printf("major ,minor,revision of weight is %d, %d ,%d\n", major, minor, revision);//0.2.0
    if ((major * 10 + minor) >= 2) {//运行这一部分
        printf("\n seen 64");
        uint64_t iseen = 0;
        fread(&iseen, sizeof(uint64_t), 1, fp);//读取一个8字节的数据
		printf("the size of uint64_t is %d\n", sizeof(uint64_t));
        *net->seen = iseen;
    }
    else {
        printf("\n seen 32");
        uint32_t iseen = 0;
        fread(&iseen, sizeof(uint32_t), 1, fp);
        *net->seen = iseen;
    }
    *net->cur_iteration = get_current_batch(*net);
    printf(", trained: %.0f K-images (%.0f Kilo-batches_64) \n", (float)(*net->seen / 1000), (float)(*net->seen / 64000));
    int transpose = (major > 1000) || (minor > 1000);

    int i;
    for(i = 0; i < net->n && i < cutoff; ++i){//cutoff = net->n
        layer l = net->layers[i];
        if (l.dontload) continue;//always 0		跳过之后的循环体,直接运行++i
        if(l.type == CONVOLUTIONAL && l.share_layer == NULL){ //只运行这一个分支的代码
            load_convolutional_weights(l, fp);
			//printf("network layer [%d] is CONVOLUTIONAL \n",i);
        }
        .......

在读取yolov3-tiny各层权重文件前,先读取4个和训练有关的参数:major,minor, revision和iseen。在前向传播的工程当中,并没有实际的应用。

parser.c中的load_convolutional_weights函数,具体执行对yolov3-tiny权重文件的下载,包括节点参数weight,偏置参数bias和批量归一化参数BN。

void load_convolutional_weights(layer l, FILE *fp)
{
	static int flipped_num;
    if(l.binary){
        //load_convolutional_weights_binary(l, fp);
        //return;
    }
    int num = l.nweights;
	//int num = l.n*l.c*l.size*l.size;//l.n 输出的层数 l.c输入的层数 
    int read_bytes;
    read_bytes = fread(l.biases, sizeof(float), l.n, fp);//读取偏置参数 l.n个float数据
    if (read_bytes > 0 && read_bytes < l.n) printf("\n Warning: Unexpected end of wights-file! l.biases - l.index = %d \n", l.index);
    //fread(l.weights, sizeof(float), num, fp); // as in connected layer
    if (l.batch_normalize && (!l.dontloadscales)){
        read_bytes = fread(l.scales, sizeof(float), l.n, fp);//读取batch normalize 参数  l.n个float数据
        if (read_bytes > 0 && read_bytes < l.n) printf("\n Warning: Unexpected end of wights-file! l.scales - l.index = %d \n", l.index);
        read_bytes = fread(l.rolling_mean, sizeof(float), l.n, fp);//读取batch normalize 参数  l.n个float数据
        if (read_bytes > 0 && read_bytes < l.n) printf("\n Warning: Unexpected end of wights-file! l.rolling_mean - l.index = %d \n", l.index);
        read_bytes = fread(l.rolling_variance, sizeof(float), l.n, fp);//读取batch normalize 参数  l.n个float数据
        if (read_bytes > 0 && read_bytes < l.n) printf("\n Warning: Unexpected end of wights-file! l.rolling_variance - l.index = %d \n", l.index);

将权重参数批量归一化

yolov3-tiny每个卷积层之后,激活函数之前,都要对结果进行Batch Normalization:在这里插入图片描述
由于BN层和卷积操作都是线性的,将权重文件进行批量归一化,可以代替卷积层之后的BN层:
在这里插入图片描述
在network.c的fuse_conv_batchnorm函数中实现权重文件和BN层的合并。

void fuse_conv_batchnorm(network net)
{
    int j;
    for (j = 0; j < net.n; ++j) {
        layer *l = &net.layers[j];
		    // printf("the %d layer batch_normalize is %d,   groups is %d \n", j, l->batch_normalize, l->groups);
        if (l->type == CONVOLUTIONAL) { //只运行这一分支   合并卷积层和batch_normal
             //printf(" Merges Convolutional-%d and batch_norm \n", j);

            if (l->share_layer != NULL) {//l->share_layer always is 0,不运行这个分支
                l->batch_normalize = 0;
            }

            if (l->batch_normalize) {//#15,22层卷积,卷积之后没有batch normalize,其他都要运行这一分支
                int f;
                for (f = 0; f < l->n; ++f)//该层神经网络 1->n 个输出层权重
                {
                    l->biases[f] = l->biases[f] - (double)l->scales[f] * l->rolling_mean[f] / (sqrt((double)l->rolling_variance[f] + .00001));

                    const size_t filter_size = l->size*l->size*l->c / l->groups;//kernel_size * kernel_size * c/分组  l->groups存在于卷积层always is 1
                    int i;
                    for (i = 0; i < filter_size; ++i) {
                        int w_index = f*filter_size + i;

                        l->weights[w_index] = (double)l->weights[w_index] * l->scales[f] / (sqrt((double)l->rolling_variance[f] + .00001));
                    }
                }

                free_convolutional_batchnorm(l);//no use
                l->batch_normalize = 0;
                ......

输入图像

yolov3-tiny输入神经网络的图像尺寸为416×416,对不符合该尺寸的图像,要进行裁剪。在image.c的resize_image函数中完成。这个可以说是整个yolo算法对输入图像唯一进行预处理的地方了。这也是yolo算法在工程应用中极好的地方,没有那么多类似于降噪、滤波之类的预处理,直接送到网络里就完事了。

//im:输入图片  w:416 h:416
//函数作用:将输入图片热size到416x416的尺寸,基本按照缩放/扩大的策略
image resize_image(image im, int w, int h)
{
    if (im.w == w && im.h == h) return copy_image(im);

    image resized = make_image(w, h, im.c);//416 x 416 x 3空的地址空间
    image part = make_image(w, im.h, im.c);//416 x im.h x im.c空的地址空间
    int r, c, k;
    float w_scale = (float)(im.w - 1) / (w - 1);//宽度缩放因子
    float h_scale = (float)(im.h - 1) / (h - 1);//高度缩放因子
    for(k = 0; k < im.c; ++k){
        for(r = 0; r < im.h; ++r){
            for(c = 0; c < w; ++c){//416
                float val = 0;
                if(c == w-1 || im.w == 1){//c =415 最后一列
                    val = get_pixel(im, im.w-1, r, k);//取原图片最后一列的像素
                } else {
                    float sx = c*w_scale;
                    int ix = (int) sx;
                    float dx = sx - ix;
                    val = (1 - dx) * get_pixel(im, ix, r, k) + dx * get_pixel(im, ix+1, r, k);
                }
                set_pixel(part, c, r, k, val);
            }
        }
    }
    for(k = 0; k < im.c; ++k){
        for(r = 0; r < h; ++r){
            float sy = r*h_scale;
            int iy = (int) sy;
            float dy = sy - iy;
            for(c = 0; c < w; ++c){
                float val = (1-dy) * get_pixel(part, c, iy, k);
                set_pixel(resized, c, r, k, val);
            }
            if(r == h-1 || im.h == 1) continue;
            for(c = 0; c < w; ++c){
                float val = dy * get_pixel(part, c, iy+1, k);
                add_pixel(resized, c, r, k, val);
            }
        }
    }

    free_image(part);
    return resized;
}

前向传播网络

network.c中的forward_network函数是整个神经网络的核心部分,各层的网络都在函数指针l.forward(l, state)中完成。

void forward_network(network net, network_state state)
{
    state.workspace = net.workspace;
    int i;
	   /// 遍历所有层,从第一层到最后一层,逐层进行前向传播(网络总共有net.n层)
    for(i = 0; i < net.n; ++i){		  
        state.index = i;/// 置网络当前活跃层为当前层,即第i层		  
        layer l = net.layers[i];/// 获取当前层		  
        if(l.delta && state.train){//不执行此分支的代码
			/// 如果当前层的l.delta已经动态分配了内存,则调用fill_cpu()函数,将其所有元素的值初始化为0			   
            scal_cpu(l.outputs * l.batch, 0, l.delta, 1);/// 第一个参数为l.delta的元素个数,第二个参数为初始化值,为0
			printf("forward_network scal_cpu of %d layer done!\n ", i);
        }
           //double time = get_time_point();
		l.forward(l, state);//进行卷积运算,激活函数,池化运算/
		   //if layer_type = convolutional ;   l.forward = forward_convolutional_layer;
		   //if layer_type = maxpool           l.forward = forward_maxpool_layer;
		   //if layer_type = yolo              l.forward = forward_yolo_layer;
		   //if layer_type = ROUTE             l.forward = forward_route_layer;其实就是数据的复制和搬移
		   //if layer_type = upsample          l.forward = forward_upsample_layer;;		  
           //printf("%d - Predicted in %lf milli-seconds.\n", i, ((double)get_time_point() - time) / 1000);
		   /// 完成某一层的推理时,置网络的输入为当前层的输出(这将成为下一层网络的输入),要注意的是,此处是直接更改指针变量net.input本身的值,
		   /// 也就是此处是通过改变指针net.input所指的地址来改变其中所存内容的值,并不是直接改变其所指的内容而指针所指的地址没变,
		   /// 所以在退出forward_network()函数后,其对net.input的改变都将失效,net.input将回到进入forward_network()之前时的值。	
		   ......

卷积层[convolution]

卷积层在convolutional_layer.c中的forward_convolutional_layer函数实现。

void forward_convolutional_layer(convolutional_layer l, network_state state)
{
    
	int out_h = convolutional_out_height(l);//获得本层卷积层输出特征图的高、宽
    int out_w = convolutional_out_width(l);
    int i, j;
	
	// l.outputs = l.out_h * l.out_w * l.out_c在make各网络层函数中赋值(比如make_convolutional_layer()),
	// 对应每张输入图片的所有输出特征图的总元素个数(每张输入图片会得到n也即l.out_c张特征图)
	// 初始化输出l.output全为0.0;输入l.outputs*l.batch为输出的总元素个数,其中l.outputs为batch
	// 中一个输入对应的输出的所有元素的个数,l.batch为一个batch输入包含的图片张数;0表示初始化所有输出为0;
    fill_cpu(l.outputs*l.batch, 0, l.output, 1);//将地址l.output,l.outputs*l.batch个float地址空间的数据初始化0
    .......

作者在进行卷积运算前,将输入特征图进行重新排序:


```c
void im2col_cpu(float* data_im,
     int channels,  int height,  int width,
     int ksize,  int stride, int pad, float* data_col)
{
    int c,h,w;
	// 计算该层神经网络的输出图像尺寸(其实没有必要再次计算的,因为在构建卷积层时,make_convolutional_layer()函数
	// 已经调用convolutional_out_width(),convolutional_out_height()函数求取了这两个参数,
	// 此处直接使用l.out_h,l.out_w即可,函数参数只要传入该层网络指针就可了,没必要弄这么多参数)
    int height_col = (height + 2*pad - ksize) / stride + 1;
    int width_col = (width + 2*pad - ksize) / stride + 1;
	
	/// 卷积核大小:ksize*ksize是一个卷积核的大小,之所以乘以通道数channels,是因为输入图像有多通道,每个卷积核在做卷积时,
	/// 是同时对同一位置处多通道的图像进行卷积运算,这里为了实现这一目的,将三通道上的卷积核并在一起以便进行计算,因此卷积核
	/// 实际上并不是二维的,而是三维的,比如对于3通道图像,卷积核尺寸为3*3,该卷积核将同时作用于三通道图像上,这样并起来就得
	/// 到含有27个元素的卷积核,且这27个元素都是独立的需要训练的参数。所以在计算训练参数个数时,一定要注意每一个卷积核的实际
	/// 训练参数需要乘以输入通道数。
    int channels_col = channels * ksize * ksize;//输入通道
	// 外循环次数为一个卷积核的尺寸数,循环次数即为最终得到的data_col的总行数
    for (c = 0; c < channels_col; ++c) {

		//行,列偏置都是对应着本次循环要操作的输出位置的像素而言的,通道偏置,是该位置像素所在的输出通道的绝对位置(通道数)

		// 列偏移,卷积核是一个二维矩阵,并按行存储在一维数组中,利用求余运算获取对应在卷积核中的列数,比如对于
		// 3*3的卷积核(3通道),当c=0时,显然在第一列,当c=5时,显然在第2列,当c=9时,在第二通道上的卷积核的第一列,
		// 当c=26时,在第三列(第三输入通道上)
        int w_offset = c % ksize;//0,1,2
		// 行偏移,卷积核是一个二维的矩阵,且是按行(卷积核所有行并成一行)存储在一维数组中的,
		// 比如对于3*3的卷积核,处理3通道的图像,那么一个卷积核具有27个元素,每9个元素对应一个通道上的卷积核(互为一样),
		// 每当c为3的倍数,就意味着卷积核换了一行,h_offset取值为0,1,2,对应3*3卷积核中的第1, 2, 3行
        int h_offset = (c / ksize) % ksize;//0,1,2
		// 通道偏移,channels_col是多通道的卷积核并在一起的,比如对于3通道,3*3卷积核,每过9个元素就要换一通道数,
		// 当c=0~8时,c_im=0;c=9~17时,c_im=1;c=18~26时,c_im=2,操作对象是排序后的像素位置
        int c_im = c / ksize / ksize;
		// 中循环次数等于该层输出图像行数height_col,说明data_col中的每一行存储了一张特征图,这张特征图又是按行存储在data_col中的某行中
        for (h = 0; h < height_col; ++h) {
			// 内循环等于该层输出图像列数width_col,说明最终得到的data_col总有channels_col行,height_col*width_col列
            for (w = 0; w < width_col; ++w) {
				// 由上面可知,对于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;//yolov3-tiny stride = 1
				// 对于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;
				// col_index为重排后图像中的像素索引,等于c * height_col * width_col + h * width_col +w(还是按行存储,所有通道再并成一行),
				// 对应第c通道,h行,w列的元素
                int col_index = (c * height_col + h) * width_col + w;//将重排后的图片像素,按照左上->右下的顺序,计算一维索引

				//im_col + width*im_row +  width*height*channel 重排前的特征图在内存中的位置索引
				// im2col_get_pixel函数获取输入图像data_im中第c_im通道,im_row,im_col的像素值并赋值给重排后的图像,
				// height和width为输入图像data_im的真实高、宽,pad为四周补0的长度(注意im_row,im_col是补0之后的行列号,
				// 不是真实输入图像中的行列号,因此需要减去pad获取真实的行列号)
                data_col[col_index] = im2col_get_pixel(data_im, height, width, channels,
                        im_row, im_col, c_im, pad);
				// return data_im[im_col + width*im_row +  width*height*channel)];
            }
        }
    }
}

通过gemm进行卷积乘加操作,通过add_bias添加偏置。

//进行卷积的乘加运算,没有bias偏置参数参与运算;
gemm(0, 0, m, n, k, 1, a, k, b, n, 1, c, n);

add_bias(l.output, l.biases, l.batch, l.n, out_h*out_w);//每个输出特征图的元素都加上对应通道的偏置参数

池化层[maxpool]

maxpool_layer.c中的forward_maxpool_layer函数完成池化操作。yolov3-tiny保留了池化层,并使用最大值池化,将尺寸为2×2的核中最大值保留下来。

void forward_maxpool_layer_avx(float *src, float *dst, int *indexes, int size, int w, int h, int out_w, int out_h, int c,
    int pad, int stride, int batch)
{

    const int w_offset = -pad / 2;
    const int h_offset = -pad / 2;
    int b, k;

    for (b = 0; b < batch; ++b) {
		// 对于每张输入图片,将得到通道数一样的输出图,以输出图为基准,按输出图通道,行,列依次遍历
		// (这对应图像在l.output的存储方式,每张图片按行铺排成一大行,然后图片与图片之间再并成一行)。
		// 以输出图为基准进行遍历,最终循环的总次数刚好覆盖池化核在输入图片不同位置进行池化操作。
        #pragma omp parallel for
        for (k = 0; k < c; ++k) {
            int i, j, m, n;
            for (i = 0; i < out_h; ++i) {
                //for (j = 0; j < out_w; ++j) {
                j = 0;
                for (; j < out_w; ++j) {
					// out_index为输出图中的索引
                    int out_index = j + out_w*(i + out_h*(k + c*b));//j + out_w * i + out_w * iout_h * k
                    float max = -FLT_MAX;// FLT_MAX为c语言中float.h定义的对大浮点数,此处初始化最大元素值为最小浮点数
                    int max_i = -1;// 最大元素值的索引初始化为-1
                    // 下面两个循环回到了输入图片,计算得到的cur_h以及cur_w都是在当前层所有输入元素的索引,内外循环的目的是找寻输入图像中,
                    // 以(h_offset + i*l.stride, w_offset + j*l.stride)为左上起点,尺寸为l.size池化区域中的最大元素值max及其在所有输入元素中的索引max_i
                    for (n = 0; n < size; ++n) {//2
                        for (m = 0; m < size; ++m) {//2
                            // cur_h,cur_w是在所有输入图像中第k通道中的cur_h行与cur_w列,index是在所有输入图像元素中的总索引。
                            // 为什么这里少一层对输入通道数的遍历循环呢?因为对于最大池化层来说输入与输出通道数是一样的,并在上面的通道数循环了!
                            int cur_h = h_offset + i*stride + n;
                            int cur_w = w_offset + j*stride + m;
                            int index = cur_w + w*(cur_h + h*(k + b*c));
							// 边界检查:正常情况下,是不会越界的,但是如果有补0操作,就会越界了,这里的处理方式是直接让这些元素值为-FLT_MAX
							// (注意虽然称之为补0操作,但实际不是补0),总之,这些补的元素永远不会充当最大元素值。
                            int valid = (cur_h >= 0 && cur_h < h &&
                                cur_w >= 0 && cur_w < w);
                            float val = (valid != 0) ? src[index] : -FLT_MAX;
							// 记录这个池化区域中的最大的元素值及其在所有输入元素中的总索引
                            max_i = (val > max) ? index : max_i;
                            max = (val > max) ? val : max;
                        }
                    }
					// 由此得到最大池化层每一个输出元素值及其在所有输入元素中的总索引。
					// 为什么需要记录每个输出元素值对应在输入元素中的总索引呢?因为在下面的反向过程中需要用到,在计算当前最大池化层上一层网络的敏感度时,
					// 需要该索引明确当前层的每个元素究竟是取上一层输出(也即上前层输入)的哪一个元素的值,具体见下面backward_maxpool_layer()函数的注释。
                    dst[out_index] = max;
                    if (indexes) indexes[out_index] = max_i;
                }
            }
        }
    }
}

路由层[route]

yolov3-tiny中共有两层路由层。第17层路由层(从0层开始),其实直接将第13层网络的输出结果输入。第20层路由层,将第19层和第8层网络结果合并在一起,19层在前,8层在后。在route_layer.c中的forward_route_layer函数中实现。

void forward_route_layer(const route_layer l, network_state state)
{
    int i, j;
    int offset = 0;
    for(i = 0; i < l.n; ++i){//l.n:  卷积层:输出特征图通道数 路由层:有几层网络层输入本层  17层:1(路由第13层)   20:2(路由第19、8层)
        int index = l.input_layers[i];//输入本网络层的网络层的索引:如13,19,8
        float *input = state.net.layers[index].output;//输入等于 之前网络层索引值得输出(.output)
        int input_size = l.input_sizes[i];//输入的网络层的数据量
        int part_input_size = input_size / l.groups;//未分组
        for(j = 0; j < l.batch; ++j){
            //copy_cpu(input_size, input + j*input_size, 1, l.output + offset + j*l.outputs, 1);
			//从首地址input处复制input_size 个数据到 l.output中
            copy_cpu(part_input_size, input + j*input_size + part_input_size*l.group_id, 1, l.output + offset + j*l.outputs, 1);//l.group_id = 0
			//其实就是copy_cpu(part_input_size, input, 1, l.output + offset, 1);
        }
        //offset += input_size;
        offset += part_input_size;
    }
}

上采样层[upsample]

yolov3-tiny中第19层是上采样层,将18层13x13x128的输入特征图转变为26x26x128的输出特征图。在upsample_layer.c中的forward_upsample_layer函数中完成。

void upsample_cpu(float *in, int w, int h, int c, int batch, int stride, int forward, float scale, float *out)
{
	
    int i, j, k, b;
    for (b = 0; b < batch; ++b) {
        for (k = 0; k < c; ++k) {
            for (j = 0; j < h*stride; ++j) {
                for (i = 0; i < w*stride; ++i) {
                    int in_index = b*w*h*c + k*w*h + (j / stride)*w + i / stride;
                    int out_index = b*w*h*c*stride*stride + k*w*h*stride*stride + j*w*stride + i;
                    if (forward) out[out_index] = scale*in[in_index];
                    else in[in_index] += scale*out[out_index];
                }
            }
        }
    }
}

上采样效果:
在这里插入图片描述

输出层[yolo]

yolo层完成了对13x13x255和26x26x255输入特诊图的logistic逻辑回归计算。每个box的预测宽度和高度不参与逻辑回归,在yolo_layer.c中的forward_yolo_layer函数中完成。

//两个yolo层 只对数据进行了logistic处理,并没有预测box的位置
//将0-1通道(x,y) 4-84(confidence+class)计算logistic,三个prior(预测框都是这样)
void forward_yolo_layer(const layer l, network_state state)
{
    int i, j, b, t, n;
	//从state.input复制数据到l.output
    memcpy(l.output, state.input, l.outputs*l.batch * sizeof(float));

#ifndef GPU
	printf("yolo v3 tiny l.n and l.batch of yolo layer is %d and %d  \n ",l.n,l.batch);
    for (b = 0; b < l.batch; ++b) {//l.batch = 1
        for (n = 0; n < l.n; ++n) {//l.n:3(yolo层)mask 0,1,2  表示每个网络单元预测三个box?
			
			//printf("l.coords is %d in yolov3 tiny yolo layer ,l.scale_x_y is %f \n", l.coords, l.scale_x_y);
            // l.coords 坐标:0  l.classes分类数量:80   l.scale_x_y:1
			//l.w:输入特征图宽度 l.h输出特征图高度  
			int index = entry_index(l, b, n*l.w*l.h, 0);//index = n*l.w*l.h*(4 + l.classes + 1)
			
		   //起始地址为:l.output + index 个数为:2 * l.w*l.h  计算逻辑回归值,并保存
            activate_array(l.output + index, 2 * l.w*l.h, LOGISTIC);  // x,y,

			//起始地址为:l.output + index 个数为:2 * l.w*l.h  计算方式为:x = x*l.scale_x_y + -0.5*(l.scale_x_y - 1) 简化后:x = x
			//yolov3-tiny l.scale_x_y = 1  实际上该函数没有参与任何的运算   scal_add_cpu
            scal_add_cpu(2 * l.w*l.h, l.scale_x_y, -0.5*(l.scale_x_y - 1), l.output + index, 1);    // scale x,y
            
			//
			index = entry_index(l, b, n*l.w*l.h, 4);//index = n*l.w*l.h*(4 + l.classes + 1)+ 4*l.w*l.h
            
			//起始地址为:l.output + index,个数为:(1+80)*l.w*l.h   计算器其逻辑回归值
			activate_array(l.output + index, (1 + l.classes)*l.w*l.h, LOGISTIC);
        }
    }

预测结果统计[detection ]


//w:输入图像宽度640,不一定是416 h:输入图像高度424,不一定是416  thresh:图像置信度阈值0.25   hier:0.5
//map:0   relative:1  num:0   letter:0
//函数作用:统计两个yolo层中 置信度大于阈值的box个数,并对这个box初始化一段地址空间 dets
//根据网络来填充该地址空间dets:
//根据yolo层 计算满足置信度阈值要求的box相对的预测坐标、宽度和高度,并将结果保存在dets[count].bbox结构体中
//每个box有80个类别,有一个置信度,该类别对应的可能性prob:class概率*置信度
///舍弃prob小于阈值0.25的box
//将满足阈值的box个数保存到num中
detection *get_network_boxes(network *net, int w, int h, float thresh, float hier, int *map, int relative, int *num, int letter)
{
	//printf("w、h、thresh、hier and letter is %d 、%d 、%f 、%f and %d\n", w, h, thresh, hier, letter);

	//函数作用:统计两个yolo层中 置信度大于阈值的box个数,并对这个box初始化一段地址空间 dets
	//将满足阈值的box个数保存到num中
    detection *dets = make_network_boxes(net, thresh, num);
	
	//根据网络来填充该地址空间dets:
	//根据yolo层 计算满足置信度阈值要求的box相对的预测坐标、宽度和高度,并将结果保存在dets[count].bbox结构体中
	//每个box有80个类别,有一个置信度,该类别对应的可能性prob:class概率*置信度
	///舍弃prob小于阈值0.25的box
    fill_network_boxes(net, w, h, thresh, hier, map, relative, dets, letter);
    return dets;
}

使用make_network_boxes来创建预测信息的指针变量:

// thresh:  置信度阈值
//num: 0
//函数作用:统计置信度大于阈值的box个数,并对这个box初始化一段地址空间
detection *make_network_boxes(network *net, float thresh, int *num)
{
    layer l = net->layers[net->n - 1];//应该是神经网络最后一层 net->n:24 最后一层yolo层
	//printf(" net->n  of network is %d\n " ,(net->n));
    int i;
	// -thresh 0.25
	//yolo层:yolov3-tiny中共有两层
	//三个prior预测框,对每个预测框中,置信度大于thresh 0.25,记为一次,将次数进行累加,并输出
	//nboxes:即为要保留的box的个数 两个yolo层中的置信度个数一起累加
    int nboxes = num_detections(net, thresh);//-thresh 0.25

	if (num) {
		printf("nbox = %d \n", num);
		*num = nboxes;//不执行该语句
	}
    //申请内存,个数为nboxes,每个内存大小为:sizeof(detection)
    detection* dets = (detection*)xcalloc(nboxes, sizeof(detection));

	//遍历每个box,每个dets.prob申请80个float类型的内存:
	//dets.uc,申请4个float类型的空间:位置信息
    for (i = 0; i < nboxes; ++i) {
        dets[i].prob = (float*)xcalloc(l.classes, sizeof(float));
        // tx,ty,tw,th uncertainty
        dets[i].uc = (float*)xcalloc(4, sizeof(float)); // Gaussian_YOLOv3
        
		if (l.coords > 4) {//不执行这个分支 l.coords:0
            dets[i].mask = (float*)xcalloc(l.coords - 4, sizeof(float));
        }
    }
    return dets;
}

使用get_yolo_detections来统计两层yolo层的预测信息:

//w,h:640,424    netw, neth:416,416 thresh:图像置信度阈值0.25   hier:0.5
//map:0   relative:1    letter:0
//根据yolo层 计算满足置信度阈值要求的box相对的预测坐标、宽度和高度,并将结果保存在dets[count].bbox结构体中
//每个box有80个类别,有一个置信度,该类别对应的可能性prob:class概率*置信度
///舍弃prob小于阈值0.25的box
int get_yolo_detections(layer l, int w, int h, int netw, int neth, float thresh, int *map, int relative, detection *dets, int letter)
{
    printf("\n l.batch = %d, l.w = %d, l.h = %d, l.n = %d ,netw = %d, neth = %d \n", l.batch, l.w, l.h, l.n, netw, neth);
    int i,j,n;
    float *predictions = l.output;//yolo层的输出
    // This snippet below is not necessary
    // Need to comment it in order to batch processing >= 2 images
    //if (l.batch == 2) avg_flipped_yolo(l);
    int count = 0;

	//printf("yolo layer l.mask[0] is %d, l.mask[1] is %d, l.mask[2] is %d\n", l.mask[0], l.mask[1], l.mask[2]);
	//printf("yolo layer l.biases[l.mask[0]*2] is %f, l.biases[l.mask[1]*2] is %f, l.biases[l.mask[2]*2] is %f\n", l.biases[l.mask[0] * 2], l.biases[l.mask[1] * 2], l.biases[l.mask[2] * 2]);
	//遍历yolo层
    for (i = 0; i < l.w*l.h; ++i){//该yolo层输出特征图的宽度、高度:13x13 26x26
        int row = i / l.w;
        int col = i % l.w;
        for(n = 0; n < l.n; ++n){//yolo层,l.n = 3
			
            //obj_index:置信度层索引
            int obj_index  = entry_index(l, 0, n*l.w*l.h + i, 4);//obj_index  = n*l.w*l.h*(4+l.classes+1) + 4*l.w*l.h + i;
            float objectness = predictions[obj_index];//获得对应的置信度
            //if(objectness <= thresh) continue;    // incorrect behavior for Nan values
            
			if (objectness > thresh) {//只有置信度大于阈值才开始执行该分支
                //printf("\n objectness = %f, thresh = %f, i = %d, n = %d \n", objectness, thresh, i, n);
                
				//box_index:yolo层每个像素点有三个box,表示每个box的索引值
				int box_index = entry_index(l, 0, n*l.w*l.h + i, 0);//box_index = n*l.w*l.h*(4+l.classes+1)+ i;

				//l.biases->偏置参数起始地址    l.mask[n]:分别为3,4,5,0,1,2,biases偏置参数偏移量
				//根据yolo层 计算满足置信度阈值要求的box相对的预测坐标、宽度和高度,并将结果保存在dets[count].bbox结构体中
                dets[count].bbox = get_yolo_box(predictions, l.biases, l.mask[n], box_index, col, row, l.w, l.h, netw, neth, l.w*l.h);

				//获取对应的置信度,该置信度经过了logistic
                dets[count].objectness = objectness;

				//获得分类数:80(int类型)
                dets[count].classes = l.classes;
                for (j = 0; j < l.classes; ++j) {
					//80个类别,每个类别对应的概率,class_index为其所在层的索引
                    int class_index = entry_index(l, 0, n*l.w*l.h + i, 4 + 1 + j);//class_index  = n*l.w*l.h*(4+l.classes+1) + (4+1+j)*l.w*l.h + i;
                    //每个box有80个类别,有一个置信度,该类别对应的可能性prob:class概率*置信度
					float prob = objectness*predictions[class_index];
					
					//舍弃prob小于阈值0.25的box
                    dets[count].prob[j] = (prob > thresh) ? prob : 0;
                }
                ++count;
            }
        }
    }
    correct_yolo_boxes(dets, count, w, h, netw, neth, relative, letter);
    return count;
}

非极大值抑制[NMS]

//dets:box结构体 nboxes:满足阈值的box个数   l.classe:80    thresh=0.45f
//两个box,同一类别进行非极大值抑制,遍历
void do_nms_sort(detection *dets, int total, int classes, float thresh)
{
    int i, j, k;
    k = total - 1;
    for (i = 0; i <= k; ++i) {//box个数
        if (dets[i].objectness == 0) {//置信度==0  不执行该分支,理论上没有objectness = 0
			printf("there is no objectness == 0 !!! \n");
            detection swap = dets[i];
            dets[i] = dets[k];
            dets[k] = swap;
            --k;
            --i;
        }
    }
    total = k + 1;
	//同一类别进行比较
    for (k = 0; k < classes; ++k) {//80个        
        //box预测的类别
		for (i = 0; i < total; ++i) {//box个数
            dets[i].sort_class = k;
        }
		//函数作用:将prob较大的box排列到前面
        qsort(dets, total, sizeof(detection), nms_comparator_v3);
        for (i = 0; i < total; ++i) {//两个box,同一类别进行非极大值抑制
            //printf("  k = %d, \t i = %d \n", k, i);
            if (dets[i].prob[k] == 0) continue;
            box a = dets[i].bbox;
            for (j = i + 1; j < total;++j){
				box b = dets[j].bbox;
				if( box_iou(a, b) > thresh) dets[j].prob[k] = 0;
            }
        }
    }
}

FPGA—YOLO

device = “xcvu3p-ffvc1517-2-e”

参考:

https://thedatabus.io/convolver

https://github.com/sumanth-kalluri/cnn_hardware_acclerator_for_fpga

https://thedatabus.io/archive

1、网络结构分析:

YOLOv3 的剪枝量化:

https://github.com/coldlarry/YOLOv3-complete-pruning

说明:最后一层75输出通道

输出为单个通道:上述过程中,每一个卷积核的通道数量,必须要求与输入通道数量一致,因为要对每一个通道的像素值要进行卷积运算,所以每一个卷积核的通道数量必须要与输入通道数量保持一致

preview

多通道

输入 416*416*3 = 519,168 个值

conv0 : input_height=418 原因:这里是做1的填充,所以 最终输入是 416+2 =418

 conv_0_param={.kernel_dim=3,
                          .pad=1,
                          .input_channel=3,
                          .input_width=418,
                          .input_height=418,
                          .output_channel=16,
                          .output_width=416,
                          .output_height=416};

16*3个卷积核:一共16*3*9个卷积核参数 +16个bias == 》 432+16 个参数

输出: 416 *416* 16 个值,即2768896个输出

conv1 :

conv_1_param={.kernel_dim=3,
                          .pad=1,
                          .input_channel=16,
                          .input_width=210,
                          .input_height=210,
                          .output_channel=32,
                          .output_width=208,
                          .output_height=208};

一共 32*16个卷积核 4608个卷积核权重参数,36个bias参数

该层输出: 208 *208* 32 = 1384448

conv2 :

conv_2_param={.kernel_dim=3,
                          .pad=1,
                          .input_channel=32,
                          .input_width=106,
                          .input_height=106,
                          .output_channel=64,
                          .output_width=104,
                          .output_height=104};

一共 64 *32个卷积核 18432个卷积核权重参数,64个bias参数

该层输出: 104 *104*64 = 692224

conv3:

conv_3_param={.kernel_dim=3,
                          .pad=1,
                          .input_channel=64,
                          .input_width=54,
                          .input_height=54,
                          .output_channel=128,
                          .output_width=52,
                          .output_height=52};

一共 128*64个卷积核 73728个卷积核权重参数,128个bias参数

该层输出: 52 *52*128 = 346112

conv4: 256*128*9= 294912 个卷积核权重参数 ,256个偏置, 输出:26*26*256=173056

conv_4_param={.kernel_dim=3,
                          .pad=1,
                          .input_channel=128,
                          .input_width=28,
                          .input_height=28,
                          .output_channel=256,
                          .output_width=26,
                          .output_height=26};

conv5:256*512*9 = 1179648个卷积核权重 512 个bias, 输出:13*13*512=86528

conv6:512*1024*9=4718592个卷积核权重,1024个bias,输出:13*13*1024=173056

conv7:1024*256*1 =262144个卷积层权重,256个bias,输出:13*13*256=43264

conv8:256*512*9 =1179648 个卷积层权重 ,512个bias 输出:13*13*512=86528

上面的卷积层:后面紧跟一个batchnormal层

CONV 9 and 12 have no batch norm layer

conv9 : 最后一层:512*255*1=130560个卷积权重,255个bias, 13*13*255= 43095

conv10: 256*128*1=32768个权重 128个bias 输出13*13*128=21632

conv11: 384*256*3 =305280 个参数 256个bias 输出26*26*256=173056

conv12: 256*255*1=65280个权重 255个bias 输出26*26*255 = 172380

wight= [448,4608,36,18432,64,73728,128,294912,256,1179648,512,4718592,1024,262144,256,1179648 ,512,130560,255,32768,128,305280,256,65280,255]
out: [ 519168,2768896,1384448, 692224,346112,173056,86528,173056,43264,86528,43095,21632,173056,172380]

汇总:8269730个wight((八百多万) 输出结果中最大 2768896,最小21632

数据位宽:

在定点实现中,首先确定整数和小数部分使用了多少位是很重要的。对于线性量化,整数的位宽度与极值以及是否会发生溢出有关。另一方面,分数部分的长度影响量化误差。此外,量化的步长将影响数据的分布。对于不同的网络类型,应在比特宽度和网络精度之间进行彻底的权衡。

在定点实现中,首先确定整数和小数部分使用了多少位是很重要的。对于线性量化,整数的位宽度与极值以及是否会发生溢出有关。另一方面,分数部分的长度影响量化误差。此外,量化的步长将影响数据的分布。对于不同的网络类型,应在比特宽度和网络精度之间进行彻底的权衡。

就总比特而言,2的幂次更可取。为了确定合适的比特宽度,在包含5000幅图像的COCO2014-val5k数据集上进行了实验。下表总结了所有卷积层的数据分布。很明显,所有输入和输出的绝对值都小于128。这意味着在不造成任何溢出的情况下,为有符号整数部分分配8位就足够了。

尝试使用8位表示整数部分,另外8位表示小数。建立了软件定点仿真系统:

因此选择16bit存储权重以及16bit的数据权重,两者都由16位定点数字表示。对于卷积窗口的结果,采用32位以减少可能的溢出,并提供更高的累积精度。通过实施定点优化,数据宽度从32位压缩到16位,从而实现更高效的DMA传输。此外,使用定点表示为延迟和资源带来了实实在在的好处。

激活函数实现:

官方 float leaky_activate(float x){return (x>0) ? x : .1*x;}

使用动态定点量化时,激活函数 Leaky ReLU 也要进行相应的量化操作。当位宽 bw 为 16 时,使用与 0xccc 相乘和右移 15 位的定点运算来拟合与 0.1 相乘的操作。量化后的 Leaky ReLU 如下所示:
𝑓′(𝑥) = (𝑥 < 0)? (𝑥 ∗ 0xccc) ≫ 15: 𝑥

合并 Batch Normalisation

OLOv3 tiny中的大多数卷积层在激活前都会进行批量归一化。使用一些数学技术消除批次标准化是安全的。可以很容易地得到以下方程式。

w’0ij和b’j是新的权重和偏差。因为这种转换可以在运行时之前完成,浮点计算过程中的一些可忽略的精度损失外。合并批次归一化也会影响权重的分布。

二维卷积:

2D 卷积是一种采用权重内核(或窗口或矩阵)并使用这些权重(此矩阵中的值)以某种方式修改输入图像(或特征图或激活图)的操作,pad 默认补0操作

激活函数

在神经网络中的各种数学线性运算之间引入了激活函数,目的是为整个网络引入非线性。没有它,整个神经网络将简化为单个线性函数,我们都知道这样的函数不足以模拟任何复杂的东西,忘记像图像识别这样的东西。

文献中有多种用于此目的的函数,但最常用(也是第一个)使用的函数是 ReLu(整流线性单元)和 Tanh(双曲正切)函数。

池化函数

池化本质上是一种“下采样”操作,旨在减少数据在网络中传播时的参数数量和复杂性。此过程涉及在输入上运行“池化窗口”并使用某种算法减小输入的大小。到目前为止,最常见的算法是 Max – Pooling 和 Average – Pooling。

最大池化——在这种方法中,我们只保留落入池化窗口的所有值中的最大值,并丢弃其他值。

平均池化 – 在此方法中,计算池化窗口内所有元素的平均值,并保留该值而不是所有值。

upsaming 上采样:

torch.nn.functional.upsample

nearest 采样

yolov3 中的上采样层使用的是2*2上采样:

import torch
from torch import nn
input = torch.arange(1, 5, dtype=torch.float32).view(1, 1, 2, 2)
input


tensor([[[[1., 2.],
          [3., 4.]]]])
 

m = nn.Upsample(scale_factor=2, mode='nearest')
m(input)
返回:

tensor([[[[1., 1., 2., 2.],
          [1., 1., 2., 2.],
          [3., 3., 4., 4.],
          [3., 3., 4., 4.]]]])

最邻近插值算法

首先假设原图是一个像素大小为W*H的图片,缩放后的图片是一个像素大小为w*h的图片,这时候我们是已知原图中每个像素点上的像素值(即灰度值等)的(⚠️像素点对应像素值的坐标都是整数)。这个时候已知缩放后有一个像素点为(x,y),想要得到该像素点的像素值,那么就要根据缩放比例去查看其对应的原图的像素点的像素值,然后将该像素值赋值给该缩放后图片的像素点(x,y)

缩放公式为:

  • 根据横轴,即宽可得:X/x = W/w
  • 根据纵轴,即高可得:Y/y = H/h
  • 那么能够得到 f(X,Y)= f( W/w * x, H/h *y)

因此这个时候缩放后的图片像素点(x,y)的像素值就对应着原图像素点( W/w * x, H/h *y)的像素值

但是这个时候会出现一个问题就是因为缩放比例的原因,会导致像素点( W/w * x, H/h *y)中的值不是整数,那么就不知道应该对应的是哪个像素点的像素值

这个时候最邻近插值算法使用的方法就是四舍五入法,表示为[.],所以像素值f(x,y) = f( [W/w * x],  [H/h *y])

举个例子,如果原图为5*5,缩放后的图为3*3,那么缩放后的图的像素点(1,1)对应的就是原图中([5/3 * 1], [5/3 * 1]) = ([0.6], [0.6]) = (1,1) 像素点对应的像素值

这种方法的好处就是简单,但是坏处就是太过粗暴,会缺失精度,造成缩放后的图像灰度上的不连续,在变化地方可能出现明显锯齿状,如下图所示:(左原图,右缩放后)

累加层accumulate:

在yolo中存在累加操作:conv层 的输出与另外一个conv层输出加和:

yolo输出层

yolo在经过多个卷积和上采样之后最终得到的是2个个卷积结果(13*13*255,26*26*255),每一个卷积结果的长和宽分别是(13×13,26×26),深度信息是 [4(box信息)+1(物体判别信息)+80(classNum置信度)] *3(每个点上面计算几个anchor)

对于具有Gh×Gw×Nin输入的YoLO层,它将原始图像划分为 Gh×Gw 网格。通道数n等于(4 +1+C)×B,其中B表示在一个网格中可以检测到多少个对象,C表示对象类别的数量。在YOLOv3 tiny中,B设置为3,COCO数据集中有80个类(c=80).

因此,YLO层的输入具有255的恒定值,这进一步分为3组( 表示在一个网格中可以检测 3种对象)。在每组中,4个通道提供关于边界框的信息,1个通道表示对象分数,其余80个通道表示单个类分数。Yolo层在所有通道上使用Sigmod激活函数,除了表示边界框宽度和高度的通道以外。关键在于指定哪些通道应该通过sigmoid函数,哪些通道应该不被触及。解决方案是提供一个表,其中每个位都表示是否应转换通道。

最终得到的边框坐标值是bx,by,bw,bh即边界框bbox相对于feature map的位置和大小,是我们需要的预测输出坐标。但网络实际上的学习目标是tx,ty,tw,th这4个偏移量(offsets,其中tx,ty是预测的坐标偏移值,tw,th是尺度缩放,有了这4个offsets,自然可以根据图2的公式去求得真正需要的bx,by,bw,bh4个坐标。bx,by是预测框中心坐标,bw,bh是预测框的宽高

其中 \(C_{x} C_{y}\)是 feature map上anchor box的宽和高, \(t_{y}, t_{w} t_{z} t_{x}\)是 4个通道提供关于边界框的信息

Cx,Cy是feature map中grid cell的左上角坐标,在yolov3中每个grid cell在feature map中的宽和高均为1。如图3的情形时,这个bbox边界框的中心属于第二行第二列的grid cell,它的左上角坐标为(1,1),故Cx=1,Cy=1。Pw、Ph是预设的anchor box映射到feature map中的宽和高,在yolov3.cfg文件中的anchor box原本设定是相对于416*416坐标系下的坐标,代码中是把cfg中读取的坐标除以stride如32映射到feature map坐标系中。

yolov3.cfg文件中的anchor box :

anchors = 10,14,  23,27,  37,58,  81,82,  135,169,  344,319

问题:anchor box作用的详细描述。
解答:YOLO3为每种FPN预测特征图(13*13,26*26,52*52)设定3种anchor box,总共聚类出9种尺寸的anchor box。在COCO数据集这9个anchor box是:(10×13),(16×30),(33×23),(30×61),(62×45),(59×119),(116×90),(156×198),(373×326)。分配上,在最小的13*13特征图上由于其感受野最大故应用最大的anchor box (116×90),(156×198),(373×326),(这几个坐标是针对416*416下的,当然要除以32把尺度缩放到13*13下),适合检测较大的目标。中等的26*26特征图上由于其具有中等感受野故应用中等的anchor box (30×61),(62×45),(59×119),适合检测中等大小的目标。较大的52*52特征图上由于其具有较小的感受野故应用最小的anchor box(10×13),(16×30),(33×23),适合检测较小的目标。同Faster-Rcnn一样,特征图的每个像素(即每个grid)都会有对应的三个anchor box,如13*13特征图的每个grid都有三个anchor box (116×90),(156×198),(373×326)(这几个坐标需除以32缩放尺寸)。

卷积:

yolo卷积默认补0

主要思想是构建一个高度流水线的流式架构,其中处理模块不必在任何时间点停止。即卷积器的任何部分都不会等待任何其他部分完成其工作并提供结果。每个阶段在每个时钟周期的输入的不同部分上连续执行整个工作的一小部分。这不仅仅是该设计的一个特点,它是一种称为“流水线”的一般原则,广泛用于将大型计算过程分解为更小的步骤,并提高整个电路可以运行的最高频率。

  • 设计使用 MAC(乘法和累加)单元,旨在将这些操作映射到 FPGA 的 DSP 模块。实现这一点将使乘法和加法运算更快,并且消耗更少的功率,因为​​ DSP 模块是在硬宏中实现的。也就是说,DSP 模块已经以最有效的方式合成、放置和路由到 FPGA 设备上的硅片中,这与一般 IP 模块不同,后者只向您提供经过试验和测试的 RTL 代码并且可以合成随心所欲。
  • 卷积器只不过是一组 MAC 单元和一些移位寄存器,当它们提供正确的输入时,在固定数量的时钟周期后输出卷积的结果
  • 输入特征图(图像)的大小是一个维度NxN,我们的内核(过滤器/窗口)是维度KxK。显然可以理解,即K < N.
  • s表示窗口在特征图上移动的步幅值。
  • 我们还有一些信号,例如valid_convend_conv信号,它们告诉外界这个卷积器模块的输出是否有效。但是为什么卷积器首先会产生无效的结果呢?窗口在输入的特定行上完成移动后,它会继续环绕到下一行,从而创建无效输出。可以避免在环绕期间进行计算,但这需要我们停止流水线,这违反了我们的流式设计原则。因此,我们只是在信号的帮助下丢弃我们认为无效的输出valid conv

具体流程:

下一次从a1开始输入,输出 输出 = w0*a01+ w1*a2 + w2*a3 + w3*a5 + w4*a6 + w5*a7 + w6*a9 + w7*a10 + w8*a11

卷积动画

输出 = w0*a0 + w1*a1 + w2*a2 + w3*a4 + w4*a5 + w5*a6 + w6*a8 + w7*a9 + w8*a10

这种架构的优点:

  • 输入特征图只需发送一次(即存储的激活必须从内存中访问一次),从而大大减少了内存访问时间和存储需求。移位寄存器为先前访问的值创建一种本地缓存。
  • 可以使用此架构计算任何大小输入的卷积,而无需由于计算能力低而中断输入或将其临时存储在其他地方
  • 对于一些需要更大步幅的奇异架构,步幅值也可以更改。

代码:

import numpy as np 
from scipy import signal

ksize = 3      															#ksize x ksize convolution kernel
im_size = 4   															#im_size x im_size input activation map

def strideConv(arr, arr2, s): 											  #the function that performs the 2D convolution
    return signal.convolve2d(arr, arr2[::-1, ::-1], mode='valid')[::s, ::s]

kernel =  np.arange(0,ksize*ksize,1).reshape((ksize,ksize))                   #the kernel is a matrix of increasing numbers
act_map = np.arange(0,im_size*im_size,1).reshape((im_size,im_size))           #the activation map is a matrix of increasin numbers

conv = strideConv(act_map,kernel,1)
print(kernel)
print(act_map)
print(conv) 

单个加法:

module mac_manual #(
	parameter N = 16,
    parameter Q = 12
    )(
    input clk,sclr,ce,
    input [N-1:0] a,
    input [N-1:0] b,
    input [N-1:0] c,
    output reg [N-1:0] p
    );

always@(posedge clk,posedge sclr)
 begin
    if(sclr)
    begin
        p<=0;
    end
    else if(ce)
    begin
        p <= (a*b+c);                   //performs the multiply accumulate operation
    end
 end
endmodule


卷积中的偏置:
preview

多卷积核的偏置:同一个通道的卷积核共用一个偏置。每个卷积核卷积后的结果+bias

https://cs231n.github.io/assets/conv-demo/index.html

定点表示:

定点数字是小数点位置保持固定的数字,与数字所代表的值无关。与浮点数相比,这使得定点数更易于理解以及在硬件中实现。与浮点运算相比,定点运算使用的资源也少得多。当然,所有这些都需要权衡。浮点运算可以为特定位宽提供比定点运算更高的精度。当我们使用二进制点位于固定位置的数字时,我们会在量化噪声方面受到打击,尽管该数字表示幅度。

Bellman-Ford 单源最短路径(回路)算法

Bellman-ford 算法比dijkstra算法更具普遍性,因为它对边没有要求,可以处理负权边与负权回路。缺点是时间复杂度过高,高达O(VE), V为顶点数,E为边数。

特点:支持有环,负权重,但不能存在负权重环

其主要思想:对所有的边进行n-1轮松弛操作,因为在一个含有n个顶点的图中,任意两点之间的最短路径最多包含n-1边。换句话说,第1轮在对所有的边进行松弛后,得到的是源点最多经过一条边到达其他顶点的最短距离;第2轮在对所有的边进行松弛后,得到的是源点最多经过两条边到达其他顶点的最短距离;第3轮在对所有的边进行松弛后,得到的是源点最多经过一条边到达其他顶点的最短距离……

Bellman-Ford 算法描述:

  1. 创建源顶点 v 到图中所有顶点的距离的集合 distSet,为图中的所有顶点指定一个距离值,初始均为 Infinite,源顶点距离为 0;
  2. 计算最短路径,执行 V – 1 次遍历;
    • 对于图中的每条边:如果起点 u 的距离 d 加上边的权值 w 小于终点 v 的距离 d,则更新终点 v 的距离值 d;
  3. 检测图中是否有负权边形成了环,遍历图中的所有边,计算 u 至 v 的距离,如果对于 v 存在更小的距离,则说明存在环;
for (var i = 0; i < n - 1; i++) {
    for (var j = 0; j < m; j++) {//对m条边进行循环
      var edge = edges[j];
      // 松弛操作
      if (distance[edge.to] > distance[edge.from] + edge.weight ){ 
        distance[edge.to] = distance[edge.from] + edge.weight;
      }
    }
}

其中, n为顶点的个数,m为边的个数,edges数组储存了所有边,distance数组是源点到所有点的距离估计值,循环结束后就是最小值。

算法过程:首先初始化,对所有的节点V来说,所有的边E进行松弛操作,再然后循环遍历每条边,如果d[v] > d[u] + w(u,v),表示有一个负权重的环路存在。

 最后如果没有负权重环路,那么d[v] 是最小路径值。

python 实现:


def bellman_ford(graph, source):
    dist = {}
    p = {}
    max = 10000
    for v in graph:
        dist[v] = max  #赋值为负无穷完成初始化
        p[v] = None
    dist[source] = 0
 
    for i in range(len( graph ) - 1):
        for u in graph:
            for v in graph[u]:
                if dist[v] > graph[u][v] + dist[u]:
                    dist[v] = graph[u][v] + dist[u]
                    p[v] = u    #完成松弛操作,p为前驱节点
 
    for u in graph:
        for v in graph[u]:
            if dist[v] > dist[u] + graph[u][v]:
                return None, None  #判断是否存在环路
 
    return dist, p
 
def test():
    graph = {
        'a': {'b': -1, 'c':  4},
        'b': {'c':  2, 'd':  3, 'e':  2},
        'c': {},
        'd': {'b':  3, 'c':  5},
        'e': {'d': -3}
    }
    dist, p = bellman_ford(graph, 'a')
    print(dist)
    print(p)
def testfail():
    graph = {
        'a': {'b': -1, 'c':  4},
        'b': {'c':  2, 'd':  3, 'e':  2},
        'c': {'d': -5},
        'd': {'b':  3},
        'e': {'d': -3}
    }
    dist, p = bellman_ford(graph, 'a')
    print(dist)
    print(p)
 
test()
testfail()

Verilog HDL语法

  1. 算数运算符:+ – * / %(取模)
  2. 赋值运算符 =(堵塞) <=(非堵塞)
  3. 关系运算符 > ,< , >= , <= , ==:逻辑相等, !=:逻辑不等,===:全等;!==:不全等 “===”和”!==”可以比较含有x和z的操作数,在模块的功能仿真中有着广泛的应用。
  4. 逻辑运算符: && || !
  5. 位运算符: 按位 & | !~ ,^按位异或 ,^~ 按位同或
  6. 条件运算符: ?: r=s?t:u
  7. 移位运算符: << >>
  8. 拼接运算符: {信号1的某几位,信号2的某几位,信号3的某几位….} 2{信号的某几位} 表示{}中的信号重复2次
  9. 缩减运算符: 单目运算符 &,|,~ : c= &b 指的是b的前一位&后一位,知道最终只剩下一个bit

Vector:

type [upper:lower] vector_name;

如果声明vector时候指定了vector方向 [8:0] 高到低,那么在使用该vector时,也要从高到低指定:[3:0] [4:1]

向量的字节顺序(或非正式地称为“方向”)是指最低有效位是具有较低的索引(较小的字节序,例如[3:0])还是具有较高的索引(较大的字节序,例如[[ 0:3])。在Verilog中,一旦以特定的字节序声明了向量,就必须始终以相同的方式使用它。例如,声明vec[0:3]时写是非法的。与字节序一致是一种好习惯,因为如果将不同字节序的向量一起分配或使用,则会发生奇怪的错误。

wire [7:0] w;         // 8-bit wire
reg  [4:1] x;         // 4-bit reg
output reg [0:0] y;   // 1-bit reg that is also an output port (this is still a vector)
input wire [3:-2] z;  // 6-bit wire input (negative ranges are allowed)
output [3:0] a;       // 4-bit output wire. Type is 'wire' unless specified otherwise.
wire [0:7] b;         // 8-bit wire where b[0] is the most-significant bit.
wire [2:0] a, c;   // Two vectors
assign a = 3'b101;  // a = 101
assign b = a;       // b =   1  implicitly-created wire
assign c = b;       // c = 001  <-- bug
my_module i1 (d,e); // d and e are implicitly one-bit wide if not declared.
                    // This could be a bug if the port was intended to be a vector.

Adding `default_nettype none would make the second line of code an error, which makes the bug more visible.
assign w = a;
takes the entire 4-bit vector a and assigns it to the entire 8-bit vector w (declarations are taken from above). If the lengths of the right and left sides don't match, it is zero-extended or truncated as appropriate.

wire 和 reg 的区别和使用:

wire表示直通,即输入有变化,输出马上无条件地反映(如与、非门的简单连接)。

reg表示一定要有触发,输出才会反映输入的状态。

reg相当于存储单元,wire相当于物理连线。reg表示一定要有触发,没有输入的时候可以保持原来的值,但不直接实际的硬件电路对应。

      两者的区别是:寄存器型数据保持最后一次的赋值,而线型数据需要持续的驱动。wire使用在连续赋值语句中,而reg使用在过程赋值语句(initial ,always)中。wire若无驱动器连接,其值为z,reg默认初始值为不定值 x 。

     在连续赋值语句中,表达式右侧的计算结果可以立即更新表达式的左侧。在理解上,相当于一个逻辑之后直接连了一条线,这个逻辑对应于表达式的右侧,而这条线就对应于wire。在过程赋值语句中,表达式右侧的计算结果在某种条件的触发下放到一个变量当中,而这个变量可以声明成reg类型的。根据触发条件的不同,过程赋值语句可以建模不同的硬件结构:如果这个条件是时钟的上升沿或下降沿,那么这个硬件模型就是一个触发器;如果这个条件是某一信号的高电平或低电平,那么这个硬件模型就是一个锁存器;如果这个条件是赋值语句右侧任意操作数的变化,那么这个硬件模型就是一个组合逻辑。

      对组合逻辑输出变量,可以直接用assign。即如果不指定为reg类型,那么就默认为1位wire类型,故无需指定1位wire类型的变量。当然专门指定出wire类型,可能是多位或为使程序易读。wire只能被assign连续赋值,reg只能在initial和always中赋值。

      输入端口可以由wire/reg驱动,但输入端口只能是wire;输出端口可以是wire/reg类型,输出端口只能驱动wire;若输出端口在过程块中赋值则为reg型,若在过程块外赋值则为net型(wire/tri)。用关键词inout声明一个双向端口, inout端口不能声明为reg类型,只能是wire类型。

      默认信号是wire类型,reg类型要申明。这里所说的默认是指输出信号申明成output时为wire。如果是模块内部信号,必须申明成wire或者reg.

      对于always语句而言,赋值要申明成reg,连续赋值assign的时候要用wire。

模块调用时 信号类型确定方法总结如下:

•信号可以分为端口信号和内部信号。出现在端口列表中的信号是端口信号,其它的信号为内部信号。

•对于端口信号,输入端口只能是net类型。输出端口可以是net类型,也可以是register类型。若输出端口在过程块中赋值则为register类型;若在过程块外赋值(包括实例化语句),则为net类型。

•内部信号类型与输出端口相同,可以是net或register类型。判断方法也与输出端口相同。若在过程块中赋值,则为register类型;若在过程块外赋值,则为net类型。

•若信号既需要在过程块中赋值,又需要在过程块外赋值。这种情况是有可能出现的,如决断信号。这时需要一个中间信号转换。

下面所列是常出的错误及相应的错误信息(error message)

•用过程语句给一个net类型的或忘记声明类型的信号赋值。

           信息:illegal …… assignment.

•将实例的输出连接到声明为register类型的信号上。

           信息:<name> has illegal output port specification.

•将模块的输入信号声明为register类型。

case 、casex和casez的区别

用法:

case(in)
 2'b00:out = 0;
 2'b01:out = 1;
 2'b10:out = 2;
 2'b11:out = 3;
endcase

casex: 该语句不考虑高阻值z以及不定值x,即在表达式进行 比较时,不会比较x、z所在位的状态:in = 001 c = xx1 那么 in在case语句中等于c,因为不考虑x位,只需要比较最后一位 1=1


//要实现 如果in的第一个1出现的位置 in == 000 则out = 0  in = 101则out =1  in =111  则out仍是1
case(in)
 2'bxx1:out = 3;
 2'bx1x:out = 2;
 2'b1xx:out = 1;
 2'b000:out = 0
endcase

同理 casez 该语句不考虑高阻值z , 即在表达式进行 比较时,不会比较z所在位的状态

Floyd算法:求多源解负权图最短路径问题

思想:动态规划

Floyd-Warshall算法(Floyd-Warshall algorithm)是解决任意两点间的最短路径的一种算法,可以正确处理有向图或负权的最短路径问题,同时也被用于计算有向图的传递闭包。

Floyd-Warshall算法的时间复杂度为O(N^3),空间复杂度为O(N^2)。

算法描述

1) 算法思想原理:

Floyd算法是一个经典的动态规划算法。用通俗的语言来描述的话,首先我们的目标是寻找从点i到点j的最短路径。

从动态规划的角度看问题,我们需要为这个目标重新做一个诠释(这个诠释正是动态规划最富创造力的精华所在)

从任意节点i到任意节点j的最短路径不外乎2种可能,1是直接从i到j,2是从i经过若干个节点k到j。所以,我们假设Dis(i,j)为节点u到节点v的最短路径的距离,对于每一个节点k,我们检查Dis(i,k) + Dis(k,j) < Dis(i,j)是否成立,如果成立,证明从i到k再到j的路径比i直接到j的路径短,我们便设置Dis(i,j) = Dis(i,k) + Dis(k,j),这样一来,当我们遍历完所有节点k,Dis(i,j)中记录的便是i到j的最短路径的距离。

2) 算法描述:

a. 从任意一条单边路径开始。所有两点之间的距离是边的权,如果两点之间没有边相连,则权为无穷大。   

b. 对于每一对顶点 u 和 v,看看是否存在一个顶点 w 使得从 u 到 w 再到 v 比己知的路径更短。如果是更新它。

Floyd算法过程矩阵的计算—-十字交叉法

方法:两条线,从左上角开始计算一直到右下角 如下所示

先来简单分析下,由于矩阵中对角线上的元素始终为0,因此以k为中间点时,从上一个矩阵到下一个矩阵变化时,矩阵的第k行,第k列和对角线上的元素是不发生改变的(对角线上都是0,因为一个顶点到自己的距离就是0,一直不变;而当k为中间点时,k到其他顶点(第k行)和其他顶点到k(第k列)的距离是不变的

因此每一步中我们只需要判断4*4-3*4+2=6个元素是否发生改变即可,也就是要判断既不在第k行第k列又不在对角线上的元素。具体计算步骤如下:以k为中间点(1)“三条线”:划去第k行,第k列,对角线(2)“十字交叉法”:对于任一个不在三条线上的元素x,均可与另外在k行k列上的3个元素构成一个2阶矩阵x是否发生改变与2阶矩阵中不包含x的那条对角线上2个元素的和有关,若二者之和小于x,则用它们的和替换x,对应的Path矩阵中的与x相对应的位置用k来替代。

给出矩阵,其中矩阵A是邻接矩阵,而矩阵Path记录u,v两点之间最短路径所必须经过的点(指的是u->….->该点->v)

分别以每一个点作为K点转接v

image
image
image
image

代码实现

java 版本

核心代码只有 4 行,实在令人赞叹。  [java]

private void floyd() {
    for (int k = 0; k < n; k++) {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                a[i][j] = Math.min(a[i][j], a[i][k] + a[k][j]);
            }
        }
    }

    // 打印
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            System.out.println(i + " " + j + ":" + a[i][j]);
        }
    }
}

正确性分析

核心代码只有四行,三行循环,一行更新,的确十分简洁优雅,可是这四行代码为什么就能求出任意两点的最短路径呢?

看代码的特点,很显然特别像是一种动态规划,确实,之前也说过该算法是基于动态规划的。

我们从动态规划的角度考虑,解动态规划题目的重点就是合理的定义状态,划分阶段,我们定义 f[k][i][j] 为经过前k的节点,从i到j所能得到的最短路径,f[k][i][j]可以从f[k-1][i][j]转移过来,即不经过第k个节点,也可以从f[k-1][i][k]+f[k-1][k][j]转移过来,即经过第k个节点。

观察一下这个状态的定义,满足不满足最优子结构和无后效性原则。

  • 最优子结构

图结构中一个显而易见的定理:

最短路径的子路径仍然是最短路径, 这个定理显而易见,比如一条从a到e的最短路a->b->c->d->e 那么 a->b->c 一定是a到c的最短路c->d->e一定是c到e的最短路,反过来,如果说一条最短路必须要经过点k,那么i->k的最短路加上k->j的最短路一定是i->j 经过k的最短路,因此,最优子结构可以保证。

  • 无后效性

状态f[k][i][j]由f[k-1][i][j],f[k-1][i,k] 以及f[k-1][k][j]转移过来,很容易可以看出,f[k] 的状态完全由f[k-1]转移过来,只要我们把k放倒最外层循环中,就能保证无后效性。

满足以上两个要求,我们即证明了Floyd算法是正确的。

我们最后得出一个状态转移方程:f[k][i][j] = min(f[k-1][i][k],f[k-1][i][k],f[k-1][k][j]) ,可以观察到,这个三维的状态转移方程可以使用滚动数组进行优化。

K 为什么要放在最外层

采用动态规划思想,f[k][i][j] 表示 i 和 j 之间可以通过编号为 1…k 的节点的最短路径。

f[0][i][j] 初值为原图的邻接矩阵。

f[k][i][j]则可以从f[k-1][i][j]转移来,表示 i 到 j 不经过 k 这个节点。

也可以 f[k-1][i][k]+f[k-1][k][j] 从转移过来,表示经过 k 这个点。

意思即 f[k][i][j]=min(f[k-1][i][j], f[k-1][i][k]+f[k-1][k][j])

然后你就会发现 f 最外层一维空间可以省略,因为 f[k] 只与 f[k-1] 有关。

虽然这个算法非常简单,但也需要找点时间理解这个算法,就不会再有这种问题啦。

Floyd算法的本质是DP,而k是DP的阶段,因此要写最外面。

想象一个图,讨论的是要从1点到达3点,是直接走还是经过中间点2,从而确定两点之间的最短路径

滚动数组优化

滚动数组是一种动态规划中常见的降维优化的方式,应用广泛(背包dp等),可以极大的减少空间复杂度。

在我们求出的状态转移方程中,我们在更新f[k]层状态的时候,用到f[k-1]层的值,f[k-2] f[k-3]层的值就直接废弃了。

所以我们大可让第一层的大小从n变成2,再进一步,我们在f[k]层更新f[k][i][j]的时候,f[k-1][i][k] 和 f[k-1][k][j] 我们如果能保证,在更新k层另外一组路径m->n的时候,不受前面更新过的f[k][i][j]的影响,就可以把第一维度去掉了。

我们现在去掉第一个维度,写成我们在代码中的那样,就是f[i][j] 依赖 f[i][k] + f[k][j] 我们在更新f[m][n]的时候,用到了f[m][k] + f[k][n] 假设f[i][j]的更新影响到了f[m][k] 或者 f[k][m] 即 {m=i,k=j} 或者 {k=i,n=j} 这时候有两种情况,j和k是同一个点,或者i和k是同一个点,那么 f[i][j] = f[i][j] + f[j][j],或者f[i][j] = f[i][i]+f[i][j] 这时候,我们所谓的“前面更新的点对”还是这两个点本来的路径,也就是说,唯一两种在某一层先更新的点,影响到后更新的点的情况,是完全合理的,所以说,我们即时把第一维去掉,也满足无后效性原则。

因此可以用滚动数组优化。

优化之后的状态转移方程即为:f[i][j] = min(f[i][j],f[i][k]+f[k][j])

求具体路径:

我们上面仅仅是知道了最短路径的长度,实际应用中我们常常是需要知道具体的路径,在Floyd算法中怎么求具体路径呢,很简单,我们只需要记录下来在更新f[i][j]的时候,用到的中间节点是哪个就可以了。

假设我们用path[i][j]记录从i到j松弛的节点k,那么从i到j,肯定是先从i到k,然后再从k到j, 那么我们在找出path[i][k] , path[k][j]即可。

即 i到k的最短路是 i -> path[i][k] -> k -> path[k][j] -> k q` 然后求path[i][k]和path[k][j] ,一直到某两个节点没有中间节点为止,代码如下:

在更新路径的时候:  [java]

if(a[i][k]>temp){
    a[i][j] = temp;
    path[i][j] = k;
}

求路径的时候:  [java]

public String getPath(int[][] path, int i, int j) {
    if (path[i][j] == -1) {
        return " "+i+" "+j;
    } else {
        int k = path[i][j];
        return getPath(path, i, k)+" "+getPath(path, k, j)+" ";
    }
}

总结

Floyd算法只能在不存在负权环的情况下使用,因为其并不能判断负权环,上面也说过,如果有负权环,那么最短路将无意义,因为我们可以不断走负权环,这样最短路径值便成为了负无穷。

但是其可以处理带负权边但是无负权环的情况。

以上就是求多源最短路的Floyd算法,基于动态规划,十分优雅。

但是其复杂度确实是不敢恭维,因为要维护一个路径矩阵,因此空间复杂度达到了O(n^2),时间复杂度达到了O(n^3),只有在数据规模很小的时候,才适合使用这个算法,但是在实际的应用中,求单源最短路是最常见的,比如在路由算法中,某个节点仅需要知道从这个节点出发到达其他节点的最短路径,而不需要知道其他点的情况,因此在路由算法中最适合的应该是单源最短路算法,Dijkstra算法。

Hierholzer 算法 —有向图求解欧拉路径(环路)

  • 欧拉路径
    如果在一张图中,可以从一点出发遍历所有的边,那么遍历过程中的这条路径就叫做欧拉路径。如果这条路径是闭合的,那就称为欧拉回路
    简单地说,如果一张图可以被“一笔画”,那么“一笔画”的那个轨迹就叫做欧拉路径。
  • 欧拉图判定定理
    含有欧拉回路的图称为欧拉图,含有欧拉路径的图称为半欧拉图。
    无向图中,如果所有顶点的度数都为偶数,则为欧拉图;如果有两个顶点的度数为奇数,其他的为偶数,则为半欧拉图。
    有向图中,如果所有顶点的入度等于出度,那么就是欧拉图。判定半欧拉图的方法也很简单,大家可以自行推理。

Hierholzer 算法

问题简述:现给出一个有向图,且为欧拉图。求欧拉回路。

Hierholzer 算法过程

  1. 选择任一顶点为起点,遍历所有相邻边。
  2. 深度搜索,访问相邻顶点。将经过的边都删除。
  3. 如果当前顶点没有相邻边,则将顶点入栈。
  4. 栈中的顶点倒序输出,就是从起点出发的欧拉回路。

从一个可能的起点出发,进行深度优先搜索,但是每次沿着辅助边从某个顶点移动到另外一个顶点的时候,都需要删除这个辅助边。如果没有可移动的路径,则将所在结点加入到栈中,并返回。

dfs(node, trace){
	while(!node.adj.isEmpty()){
		Node next = node.adj.removeLast();
		dfs(next, trace);
	}
	trace.addLast(node);
}

最后得到的栈中保存的就是整个欧拉闭迹中的顶点。(要恢复我们需要不断出栈,因此如果你用列表来存欧垃迹的话需要反转一次)。

重新安排行程:给你一份航线列表 tickets ,其中 tickets[i] = [fromi, toi] 表示飞机出发和降落的机场地点。请你对该行程进行重新规划排序。所有这些机票都属于一个从 JFK(肯尼迪国际机场)出发的先生,所以该行程必须从 JFK 开始。如果存在多种有效的行程,请你按字典排序返回最小的行程组合。例如,行程 [“JFK”, “LGA”] 与 [“JFK”, “LGB”] 相比就更小,排序更靠前。
假定所有机票至少存在一种合理的行程。且所有的机票 必须都用一次且只能用一次。

在这里插入图片描述
class Solution:
    ''' Hierholzer 算法
    '''
    def findItinerary(self, tickets):
        def dfs(cur, graph, res):
            while graph[cur]: # 遍历所有边
                dfs(graph[cur].pop(), graph, res) # 访问并删除
            if not graph[cur]: # 将顶点入栈
                res.append(cur)
        import collections
        # 建图
        graph = collections.defaultdict(list)
        for start, end in sorted(tickets)[::-1]:
            graph[start].append(end)
        res = []
        dfs("JFK", graph, res)
        return res[::-1] # 倒序输出

AlphaFold2论文

https://www.nature.com/articles/s41586-021-03819-2

https://github.com/deepmind/alphafold

沐神AlphaFold讲解

论文补充材料:https://www.biorxiv.org/content/10.1101/2021.10.04.463034v1

摘要

蛋白质对生命至关重要,了解它们的结构可以促进对其功能的系统理解。通过大量的实验,已经确定了大约 100,000 种独特的蛋白质结构,但这仅代表了数十亿已知蛋白质序列中的一小部分。

蛋白质结构是指蛋白质分子的空间结构。作为一类重要的生物大分子,蛋白质主要由化学元素组成。所有蛋白质都是由20种不同的L型α氨基酸连接形成的多聚体,在形成蛋白质后,这些氨基酸又被称为残基。
蛋白质一级结构:组成蛋白质多肽链的线性氨基酸序列。
蛋白质二级结构:依靠不同氨基酸之间的C=O和N-H基团间的氢键形成的稳定结构,主要为α螺旋和β折叠。
蛋白质三级结构:通过多个二级结构元素在三维空间的排列所形成的一个蛋白质分子的三维结构。
蛋白质四级结构:用于描述由不同多肽链(亚基)间相互作用形成具有功能的蛋白质复合物分子。

仅根据其氨基酸序列预测蛋白质的三级结构,这是“蛋白质折叠问题”, 50 多年来一直是一个重要的开放性研究问题。尽管最近取得了进展,但现有方法仍远未达到原子级准确度,尤其是当没有可用的同源结构时。

同源结构,那些不同物种因来自共同祖先而具有的相似性结构。 例如现代马经较长时间的修饰成为具有一个趾,鼹鼠及其他洞穴动物成为瘤状肢体,大象的肢体成为柱状,这样功能各异的前肢有一个共同来源,他们都来自原始陆生脊椎动物五趾型的肢体。

在这里,我们提供了第一种计算方法,即使在不知道相似结构的情况下,它也可以以原子精度定期预测蛋白质结构。我们在具有挑战性的第 14 次蛋白质结构预测关键评估 (CASP14) 中验证了我们基于神经网络模型的完全重新设计的AlphaFold,在大多数情况下表现出与实验相媲美的准确性,并且大大优于其他方法。支持最新版本的 AlphaFold 是一种新颖的机器学习方法,它将关于蛋白质结构的物理和生物学知识,利用多序列比对,融入深度学习算法的设计中。

序列比对指将两个或多个序列排列在一起,标明其相似之处。序列中可以插入间隔(通常用短横线“-”表示)。对应的相同或相似的符号(在核酸中是A, T(或U), C, G,在蛋白质中是氨基酸残基的单字母表示)排列在同一列上。
tcctctgcctctgccatcat—caaccccaaagt
|||| ||| ||||| ||||| ||||||||||||
tcctgtgcatctgcaatcatgggcaaccccaaagt
多序列比对是成对比对的延伸,是为了在一次比对里面处理多于两条的的序列。多序列比对方法试图比对一个指定序列集合里面的所有序列,这可以帮助确定这些序列的共同区段。进行多序列比对有几种方法,最常用的一种是Clustal程序集,它使用渐进多序列比对算法。Clustal在cladistics中被用来建立进化树,在PSI-BLAST和Hidden Markov model (HMM)中用来建立序列档案以在序列数据库中搜索更远的同源序列。

从蛋白质序列预测蛋白质3D结构的计算方法的发展沿着两条互补的路径前进,分别关注物理相互作用或进化历史。 物理相互作用方案将我们认知的分子驱动力(molecular driving forces)整合到物理热力学或动力学模拟或统计模型逼近中。 虽然理论上非常吸引人,但由于分子模拟的计算难度、蛋白质稳定性的上下文依赖性以及难以产生足够准确的蛋白质物理学模型,这种方法已被证明对即使是中等大小的蛋白质也极具挑战性。 近年来,进化方案提供了一种替代方案,其中蛋白质结构约束来自蛋白质进化历史的生物信息学分析、与已解决结构的同源性和成对进化相关性。

这种生物信息学方法极大地受益于蛋白质数据库 (PDB) 中存储的实验蛋白质结构的稳定增长、基因组测序的爆炸式增长以及用于解释这些相关性的深度学习技术的快速发展。 尽管取得了这些进展,基于物理和进化历史的方法产生的预测远低于实验准确性。

我们(DeepMind团队)开发了第一个能够在大多数情况下预测蛋白质结构接近实验准确性的计算方法。 我们开发的神经网络 AlphaFold 已进入 CASP14 评估(2020 年 5 月至 7 月)。 CASP 评估每两年进行一次,使用未在 PDB 中存放或公开披露的最近解决的结构,因此它是对参与方法的盲测,长期以来一直作为结构预测准确性的金标准评估。

蛋白质结构预测 (CASP) 实验的批判性评估旨在建立蛋白质结构预测的当前技术水平,确定已取得的进展,并突出未来可能最有成效的工作重点。

… 省略了部分内容

在图 2a 中证明,CASP14 中展示AlphaFold的高准确率扩展到最近大量PDB结构样本,其中所有结构在我们的训练数据截止后都存储在PDB中,并作为完整链进行分析。此外,当主链预测准确时,观察到高侧链准确性(图 2b),并且表明我们预测的局部距离差异测试(pLDDT)置信度可靠地预测了 Ca 局部距离差异测试(lDDT-Cα)准确度相应的预测(图 2c)。我们还发现可以准确估计全局叠加度量模板建模分数 (TM-score)(图 2d)。总体而言,这些验证了 AlphaFold 在 CASP14 蛋白上的高精度和可靠性也可以迁移到最近 PDB 提交的未经整理的数据集中。

Fig. 1 | AlphaFold produces highly accurate structures.

PDB蛋白质结构数据库(Protein Data Bank,简称PDB)是美国Brookhaven国家实验室于1971年创建的,由结构生物信息学研究合作组织(Research Collaboratory for Structural Bioinformatics,简称RCSB)维护。和核酸序列数据库一样,可以通过网络直接向PDB数据库提交数据。

AlphaFold网络

AlphaFold 通过结合基于蛋白质结构的进化、物理和几何约束的新型神经网络架构和训练程序,大大提高了结构预测的准确性。特别是,我们展示了一种联合嵌入多序列比对 (MSA) 和成对特征的新架构、一种新的输出表示和相关损失函数,可实现准确的端到端结构预测、新的等变注意架构、使用中间损失函数来实现预测的迭代改进,屏蔽 MSA 损失与结构联合训练,使用自蒸馏从未标记的蛋白质序列中学习,以及自我估计准确率。

蒸馏,就是知识蒸馏,将教师网络(teacher network)的知识迁移到学生网络(student network)上,使得学生网络的性能表现如教师网络一般;或者大型模型迁移到小型模型中,小型模型的参数规模小,运行速度快,但性能与大型模型参不多。

AlphaFold 网络使用一级氨基酸序列和同源物的比对序列作为输入,直接预测给定蛋白质的所有重原子的 3-D 坐标(图 1e,请参阅方法了解输入的详细信息,包括数据库、MSA 构建和使用的模板)。 最重要的方法和组件的描述如下在附件中提供了完整的网络架构和训练过程在附件的方法章节中。

该网络包括两个主要阶段。 首先,网络的主干通过我们称为 Evoformer 的新型神经网络块的重复层处理输入,以生成 Nseq × Nres 数组(Nseq:序列数,Nres:残基数),表示已处理的 MSA 和 Nres × Nres 数组,表示残基对。 MSA 表示是用原始 MSA 初始化的,但请参阅 附件-方法 1.2.7 了解处理非常深的 MSA 的详细信息。 Evoformer 块包含许多新颖的基于注意力和非基于注意力的组件。 我们在“可解释性”部分展示了证据,表明在 Evoformer 块中早期出现了具体的结构假设并不断完善。 Evoformer 模块的关键创新是在 MSA 内交换信息的新机制和允许直接推理空间和进化关系的配对表示。

网络的主干之后是结构模块,该模块以蛋白质的每个残基(全局刚体框架)的旋转和平移的形式引入了明确的 3-D 结构。这些表示在简单的状态下初始化,所有旋转设置为一致,所有位置设置为原点,但快速发展和完善,具有精确原子细节的高度准确的蛋白质结构。网络这一部分的关键创新包括打破链原子结构以允许同时对结构的所有部分进行局部细化,一种新颖的等变变换器允许网络隐式推理未表示的侧链原子,以及一个损失项代替残基的方向正确性的重要权重。在结构模块和整个网络中,我们通过反复将最终损失函数应用于输出,然后将输出递归地提供给相同的模块来强化迭代细化的概念。使用整个网络的迭代细化(我们称之为“循环”) 对准确性有显着贡献,而额外的训练时间很少(有关详细信息,请参阅附件-方法-1.8)。

Evoformer模块

名为 Evoformer(图 1e 和 3a)的网络构建块的关键原理是将蛋白质结构预测视为 3-D 空间中的图推理问题,其中图的边缘由邻近的残基定义。 配对表示的元素编码有关残基之间关系的信息(图 3b)。 MSA 表示的列编码输入序列的各个残基,而行表示这些残基出现的序列。 在这个框架内,我们定义了许多更新操作,这些更新操作应用于每个块中,其中不同的更新操作被串联应用。

Fig. 3 | Architectural details.

MSA 表示通过在 MSA 序列维度上求和的逐元素外积更新配对表示。 与之前的工作不同,此操作应用于每个块中,而不是在网络中应用一次,这使得从不断发展的 MSA 表示到配对表示的连续通信成为可能。

在配对表示中,有两种不同的更新模式。两者都受到配对表示一致性必要性的启发——为了将氨基酸的配对描述表示为单个 3-D 结构,必须满足许多约束,包括距离上的三角不等式。基于这种直觉,我们根据涉及三个不同节点的边三角形来安排对表示的更新操作(图 3c)。特别是,我们向轴向注意力添加了一个额外的 logit 偏差,以包括三角形的“缺失边”,并且我们定义了一个非注意力更新操作“三角形乘法更新”,它使用两条边来更新缺失的第三条边(参见 附件-方法-1.6.5 了解详情)。三角形乘法更新最初是作为注意力更对称且更便宜的替代品而开发的,仅使用注意力或乘法更新的网络都能够产生高精度结构。然而,两个更新的组合更准确。

我们还在 MSA 表示中使用了一种轴向注意力的变体。 在 MSA 中的 per-sequence attention 期间,我们从 pair stack 中投射额外的 logits 以偏置 MSA attention。 这通过提供从配对表示返回到 MSA 表示的信息流来关闭循环,确保整个 Evoformer 模块能够完全混合对和 MSA 表示之间的信息,并为结构模块中的结构生成做好准备。

端到端结构预测

结构模块(图 3d)使用对表示和来自主干的 MSA 表示的原始序列行(“single representation”,“单一表示”)在具体的 3-D 主干结构上运行。 3-D 主干结构表示为 Nres 独立的旋转和平移,每个旋转和平移相对于全局框架(residue gas,“残余气”?,图 3e)。这些旋转和平移,代表 N-Cα-C 原子的几何形状,优先考虑蛋白质骨架的方向,以便每个残基的侧链位置在该框架内受到高度限制。相反,肽键几何形状完全不受约束,并且在应用结构模块期间观察到网络经常违反链约束,因为打破此约束允许对链的所有部分进行局部细化,而无需解决复杂的闭环问题。在微调期间通过违规损失项鼓励满足肽键几何形状。只有在 Amber力场中的梯度下降结构的预测后松弛中才能实现肽键几何形状的精确执行。根据经验,这种最终松弛不会提高模型的准确性,如通过全局距离测试 (GDT) 或 IDDT-Cα34 测量的,但确实消除了分散注意力的立体化学违规而不会损失准确性。

AMBER力场是在生物大分子的模拟计算领域有着广泛应用的一个分子力场。AMBER力场的优势在于对生物大分子的计算,其对小分子体系的计算结果常常不能令人满意。

residue gas表示分两个阶段迭代更新(图 3d)。首先,我们称为不变点注意力 ( Point Attention) 的新型几何感知注意操作用于更新 Nres 神经激活集(single representation,“单一表示”)而不改变 3-D 位置,然后对residue gas使用更新的激活。不变点注意力通过在每个残基的局部框架中产生的 3-D 点来增强每个通常的注意力查询、键和值,这样最终值对全局旋转和平移是不变的(参见方法“不变点注意力(IPA)”了解详情)。 3-D 查询和键也对注意力施加了强烈的空间/局部性偏差,这非常适合蛋白质结构的迭代细化。在每个注意力操作和逐元素转换块之后,该模块计算每个主干帧的旋转和平移的更新。这些更新在每个残差的局部框架内的应用使得整体注意力和更新块成为对residue gas的等变操作。

侧链 chi 角的预测以及结构的最终每个残基精度 (pLDDT) 是在网络末端的最终激活上使用小的每个残基网络计算的。 TM 分数 (pTM) 的估计值是从成对错误预测中获得的,该预测被计算为最终对表示的线性投影。最后的损失(我们称之为帧对齐点误差(FAPE)(图 3f))将预测的原子位置与许多不同对齐下的真实位置进行比较。对于每个对齐,通过将预测帧 (Rk,tk) 对齐到相应的真实帧来定义,我们计算所有预测原子位置 xi 与真实原子位置的距离。由此产生的 Nframes × Natoms 距离受到限制的 L1 损失的惩罚。这对原子相对于每个残基的局部框架是正确的产生了强烈的偏见,因此在其侧链相互作用方面是正确的,并为 AlphaFold 提供了手性的主要来源(增刊-方法 1.9.3 和增刊-图 9)

使用标记和未标记数据进行训练

AlphaFold 架构能够仅使用对 PDB 数据的监督学习来训练到高精度,但我们能够使用类似于noisy-student自我蒸馏的方法来提高准确性(见图 4a)。 在这个过程中,我们使用一个训练有素的网络来预测来自 Uniclust30的约 350,000 个不同序列的结构,并将预测结构的新数据集过滤为高置信度子集。 然后,我们使用 PDB 和这个新的预测结构数据集的混合作为训练数据从头开始训练相同的架构,其中各种训练数据增强(例如裁剪和 MSA 子采样)使网络难以重现先前预测的结构。 这种自蒸馏过程有效地利用了未标记的序列数据,并显着提高了所得网络的准确性。

Self-training是最简单的半监督方法之一,其主要思想是找到一种方法,用未标记的数据集来扩充已标记的数据集。算法流程如下:
(1)首先,利用已标记的数据来训练一个好的模型,然后使用这个模型对未标记的数据进行标记。
(2)然后,进行伪标签的生成,因为我们知道,已训练好的模型对未标记数据的所有预测都不可能都是好的,因此对于经典的Self-training,通常是使用分数阈值过滤部分预测,以选择出未标记数据的预测标签的一个子集。
(3)其次,将生成的伪标签与原始的标记数据相结合,并在合并后数据上进行联合训练。
(4)整个过程可以重复n次,直到达到收敛。

此外,我们随机屏蔽或突变 MSA 中的单个残基,并具有来自 Transformers (BERT) 式目标的双向编码器表示来预测 MSA 序列的屏蔽元素。 这个目标鼓励网络学习解释系统发育和协变关系,而无需将特定的相关统计量硬编码到特征中。 与最近的独立工作相比,BERT 目标是在相同的训练示例上与正常的 PDB 结构损失联合训练的,并且没有进行预训练。

解释神经网络

为了了解 AlphaFold 如何预测蛋白质结构,我们为网络中的 48 个 Evoformer 块中的每一个训练了一个单独的结构模块,同时保持主网络的所有参数保持不变(补充方法 1.14)。包括我们的回收阶段,这提供了 192 个中间结构的轨迹,每个完整的 Evoformer 块一个,其中每个中间体代表网络对该块最可能结构的信念。在前几个块之后产生的轨迹出奇地平滑,表明 AlphaFold 对结构进行了不断的增量改进,直到它不能再改进为止(见图 4b 的准确度轨迹)。这些轨迹也说明了网络深度的作用。对于非常具有挑战性的蛋白质,如 SARS-CoV-2 Orf8 (T1064),网络搜索并重新排列多层的二级结构元素,然后再确定一个好的结构。对于 LmrP (T1024) 等其他蛋白质,网络会在前几层内找到最终结构。请参阅增刊-视频 1-4 的 CASP14 目标 T1024、T1044、T1064 和 T1091 的结构轨迹显示了一系列蛋白质大小和难度的清晰迭代构建过程。在补充-方法 1.16 和增刊-图12-13,我们解释了 AlphaFold 层产生的注意力图。

Fig. 4 | Interpreting the neural network

图 4a 包含 AlphaFold 组件的详细消融,表明各种不同的机制有助于 AlphaFold 的准确性。 请参阅增刊-方法 1.13 详细描述了每个消融模型、它们的训练细节、消融结果的扩展讨论以及 MSA 深度对每次消融的影响(补充-图 10)。

MSA 深度和跨链接触

虽然 AlphaFold 在绝大多数沉积的 PDB 结构中具有很高的准确性,但我们注意到仍然存在影响准确性或限制模型适用性的因素。当平均比对深度小于约 30 个序列时,该模型使用多个序列比对,准确度大幅下降(详见图 5a)。我们观察到阈值效应,其中 MSA 深度超过约 100 个序列的改进导致小增益。我们假设需要 MSA 信息在网络的早期阶段粗略地找到正确的结构,但是将该预测细化为高精度模型并不关键取决于 MSA 信息。我们观察到的另一个实质性限制是,与异型接触的数量相比,AlphaFold 对于链内或同型接触很少的蛋白质要弱得多。这通常发生在较大复合物中的桥接结构域,其中蛋白质的形状几乎完全由与复合物中其他链的相互作用产生。相反,AlphaFold 通常能够为同聚体提供高精度预测,即使链基本上交织在一起(例如图 5b)。我们希望 AlphaFold 的想法很容易适用于预测未来系统中的完整异质复合物,并且这将消除具有大量异质接触的蛋白质链的困难。

Fig. 5 | Effect of MSA depth and cross-chain contacts.

相关工作

蛋白质结构预测有一个漫长而多样的发展,在许多优秀的评论中得到了广泛的介绍。尽管将神经网络应用于结构预测的历史悠久,但它们最近才开始改进结构预测。这些方法通过将蛋白质结构预测问题处理为将进化耦合的“图像” 转换为蛋白质距离矩阵的“图像”,然后将距离预测集成到启发式系统中,从而有效地利用了计算机视觉系统的快速改进。产生最终的 3-D 坐标预测。最近开发了一些工作来直接预测 3-D 坐标 ,但这些方法的准确性与传统的手工结构预测管道不匹配 。同时,基于注意力的语言处理网络的成功 和最近的计算机视觉激发了对解释蛋白质序列的基于注意力的方法的探索。

讨论

我们在设计 AlphaFold 时采用的方法是生物信息学和物理方法的结合:我们使用物理和几何归纳偏差来构建组件,这些组件可以从 PDB 数据中学习,并最大限度地减少手工制作的特征(例如,AlphaFold 在没有氢的情况下有效地构建氢键债券评分函数)。这使得网络从 PDB 中的有限数据中更有效地学习,但能够应对结构数据的复杂性和多样性。特别是,AlphaFold 能够处理缺失的物理环境,并在具有挑战性的情况下生成准确的模型,例如交织的同源异构体或仅在未知血红素组存在时才会折叠的蛋白质。处理未指定结构条件的能力对于从 PDB 结构中学习至关重要,因为 PDB 代表了结构已被求解的所有条件。一般来说,AlphaFold 被训练产生最有可能作为 PDB 结构的一部分出现的蛋白质结构。在特定化学计量或配体/离子可单独从序列预测的情况下,AlphaFold 可能会产生一种隐式遵守这些约束的结构。

AlphaFold 已经向实验界展示了它的实用性,包括分子置换和解释低温电子显微镜 (cryo-EM) 图 。 此外,由于 AlphaFold 直接输出蛋白质坐标,因此 AlphaFold 会根据蛋白质序列的长度以图形处理单元 (GPU) 分钟到 GPU 小时生成预测(例如,对于 384 个残基,每个模型大约 1 个 GPU 分钟,请参阅方法了解详细信息 )。 这开辟了在蛋白质组规模及其他范围内预测结构的令人兴奋的可能性。

可用基因组测序技术和数据的爆炸式增长彻底改变了生物信息学,但实验结构测定的内在挑战阻止了我们结构知识的类似扩展。 通过开发准确的蛋白质结构预测算法,再加上由实验社区组装的现有大型且精心准备的结构和序列数据库,我们希望加速结构生物信息学的进步,以跟上基因组学革命的步伐。 我们希望 AlphaFold 以及将其技术应用于其他生物物理问题的计算方法将成为现代生物学的重要工具。

在线内容

任何方法、附加参考资料、自然研究报告摘要、源数据、扩展数据、补充信息、致谢、同行评审信息; 作者贡献和竞争利益的详细信息; 数据和代码可用性声明可在 Highly accurate protein structure prediction with AlphaFold – Nature 上获得。

 Highly accurate protein structure prediction with AlphaFold – Nature

dijkstra算法求最短路径

基本思想

dijkstra的算法思想是从以上最短距离数组中每次选择一个最近的点,将其作为下一个点,然后重新计算从起始点经过该点到其他所有点的距离,更新最短距离数据。已经选取过的点就是确定了最短路径的点,不再参与下一次计算。

  1. 通过Dijkstra计算图G中的最短路径时,需要指定起点s(即从顶点s开始计算)。
  2. 此外,引进两个集合S和U。S的作用是记录已求出最短路径的顶点(以及相应的最短路径长度),而U则是记录还未求出最短路径的顶点(以及该顶点到起点s的距离)。
  3. 初始时,S中只有起点s;U中是除s之外的顶点,并且U中顶点的路径是”起点s到该顶点的路径”。然后,从U中找出路径最短的顶点,并将其加入到S中;接着,更新U中的顶点和顶点对应的路径。 然后,再从U中找出路径最短的顶点,并将其加入到S中;接着,更新U中的顶点和顶点对应的路径。 … 重复该操作,直到遍历完所有顶点。

操作步骤

  1. 初始时,S只包含起点s;U包含除s外的其他顶点,且U中顶点的距离为”起点s到该顶点的距离”[例如,U中顶点v的距离为(s,v)的长度,然后s和v不相邻,则v的距离为∞]。
  2. 从U中选出”距离最短的顶点k”,并将顶点k加入到S中;同时,从U中移除顶点k。
  3. 更新U中各个顶点到起点s的距离。之所以更新U中顶点的距离,是由于上一步中确定了k是求出最短路径的顶点,从而可以利用k来更新其它顶点的距离;例如,(s,v)的距离可能大于(s,k)+(k,v)的距离。
  4. 重复步骤(2)和(3),直到遍历完所有顶点。

第1步:初始化距离,其实指与D直接连接的点的距离。dis[c]代表D到C点的最短距离,因而初始dis[C]=3,dis[E]=4,dis[D]=0,其余为无穷大。设置集合S用来表示已经找到的最短路径。此时,S={D}。现在得到D到各点距离{D(0),C(3),E(4),F(*),G(*),B(*),A(*)},其中*代表未知数也可以说是无穷大,括号里面的数值代表D点到该点的最短距离。

第2步:不考虑集合S中的值,因为dis[C]=3,是当中距离最短的,所以此时更新S,S={D,C}。接着我们看与C连接的点,分别有B,E,F,已经在集合S中的不看,dis[C-B]=10,因而dis[B]=dis[C]+10=13,dis[F]=dis[C]+dis[C-F]=9,dis[E]=dis[C]+dis[C-E]=3+5=8>4(初始化时的dis[E]=4)不更新。此时{D(0),C(3),E(4),F(9),G(*),B(13),A(*)}。

第3步:在第2步中,E点的值4最小,更新S={D,C,E},此时看与E点直接连接的点,分别有F,G。dis[F]=dis[E]+dis[E-F]=4+2=6(比原来的值小,得到更新),dis[G]=dis[E]+dis[E-G]=4+8=12(更新)。此时{D(0),C(3),E(4),F(6),G(12),B(13),A(*)}

第4步:在第3步中,F点的值6最小,更新S={D,C,E,F},此时看与F点直接连接的点,分别有B,A,G。dis[B]=dis[F]+dis[F-B]=6+7=13,dis[A]=dis[F]+dis[F-A]=6+16=22,dis[G]=dis[F]+dis[F-G]=6+9=15>12(不更新)。此时{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第5步:在第4步中,G点的值12最小,更新S={D,C,E,F,G},此时看与G点直接连接的点,只有A。dis[A]=dis[G]+dis[G-A]=12+14=26>22(不更新)。{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第6步:在第5步中,B点的值13最小,更新S={D,C,E,F,G,B},此时看与B点直接连接的点,只有A。dis[A]=dis[B]+dis[B-A]=13+12=25>22(不更新)。{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第6步:最后只剩下A值,直接进入集合S={D,C,E,F,G,B,A},此时所有的点都已经遍历结束,得到最终结果{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

python实现:

将以上的过程使用python来实现。
首先总结一个Dijkstra算法的核心思想,分成两步走:

  1. 构造一个最短路径数组,每次找到数组中未访问的节点里最小的点
  2. 以上一步的节点为最新节点,更新起始点到所有点的距离

使用python就是实现这两步即可



MAX= float('inf')

matrix = [
    [0,10,MAX,4,MAX,MAX],
    [10,0,8,2,6,MAX],
    [MAX,8,10,15,1,5],
    [4,2,15,0,6,MAX],
    [MAX,6,1,6,0,12],
    [MAX,MAX,5,MAX,12,0]
    ]


def dijkstra(matrix, start_node):
    
    #矩阵一维数组的长度,即节点的个数
    matrix_length = len(matrix)

    #访问过的节点数组
    used_node = [False] * matrix_length

    #最短路径距离数组
    distance = [MAX] * matrix_length

    #初始化,将起始节点的最短路径修改成0
    distance[start_node] = 0
    
    #将访问节点中未访问的个数作为循环值,其实也可以用个点长度代替。
    while used_node.count(False):
        min_value = float('inf')
        min_value_index = 999
        
        #在最短路径节点中找到最小值,已经访问过的不在参与循环。
        #得到最小值下标,每循环一次肯定有一个最小值
        for index in range(matrix_length):
            if not used_node[index] and distance[index] < min_value:
                min_value = distance[index]
                min_value_index = index
        
        #将访问节点数组对应的值修改成True,标志其已经访问过了
        used_node[min_value_index] = True

        #更新distance数组。
        #以B点为例:distance[x] 起始点达到B点的距离,
        #distance[min_value_index] + matrix[min_value_index][index] 是起始点经过某点达到B点的距离,比较两个值,取较小的那个。
        for index in range(matrix_length):
            distance[index] = min(distance[index], distance[min_value_index] + matrix[min_value_index][index])

    return distance




start_node = int(input('请输入起始节点:'))
result = dijkstra(matrix,start_node)
print('起始节点到其他点距离:%s' % result)

数据结构 — 单调队列

顾名思义,单调队列的重点分为 “单调” 和 “队列”

“单调” 指的是元素的的 “规律”——递增(或递减)

“队列” 指的是元素只能从队头和队尾进行操作

  • 单调递增队列:保证队列头元素一定是当前队列的最小值,用于维护区间的最小值。
  • 单调递减队列:保证队列头元素一定是当前队列的最大值,用于维护区间的最大值。

1、对于单调递增队列,对于一个元素如果它大于等于队尾元素, 那么直接把元素进队, 如果这个元素小于队尾元素,则将队尾元素出队列,直到满足这个元素大于等于队尾元素。

2、对于单调递减队列,对于一个元素如果它小于等于队尾元素, 那么直接把元素进队, 如果这个元素大于队尾元素,则将队尾元素出队列,直到满足这个元素小于等于队尾元素。

3、判断队头元素是否在待求解的区间之内,如果不在,就将其删除。(这个很好理解呀,因为单调队列的队头元素就是待求解区间的极值)

4、经过上面两个操作,取出 队列的头元素 ,就是 当前区间的极值

可以发现队列中的元素都是单调递减的(不然也就不叫单调递减队列啦),同时既有入队列的操作、也有出队列的操作。

实现单调队列,主要分为三个部分:

  • 去尾操作队尾元素出队列。当队列有新元素待入队,需要从队尾开始,删除影响队列单调性的元素,维护队列的单调性。(删除一个队尾元素后,就重新判断新的队尾元素)

去尾操作结束后,将该新元素入队列。

  • 删头操作队头元素出队列。判断队头元素是否在待求解的区间之内,如果不在,就将其删除。(这个很好理解呀,因为单调队列的队头元素就是待求解区间的极值)
  • 取解操作 :经过上面两个操作,取出 队列的头元素 ,就是 当前区间的极值

参考代码

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#define maxn 1000100
using namespace std;
int q[maxn], a[maxn];
int n, k;

void getmin() {  // 得到这个队列里的最小值,直接找到最后的就行了
  int head = 0, tail = 0;
  for (int i = 1; i < k; i++) {
    while (head <= tail && a[q[tail]] >= a[i]) tail--;
    q[++tail] = i;
  }
  for (int i = k; i <= n; i++) {
    while (head <= tail && a[q[tail]] >= a[i]) tail--;
    q[++tail] = i;
    while (q[head] <= i - k) head++;
    printf("%d ", a[q[head]]);
  }
}

void getmax() {  // 和上面同理
  int head = 0, tail = 0;
  for (int i = 1; i < k; i++) {
    while (head <= tail && a[q[tail]] <= a[i]) tail--;
    q[++tail] = i;
  }
  for (int i = k; i <= n; i++) {
    while (head <= tail && a[q[tail]] <= a[i]) tail--;
    q[++tail] = i;
    while (q[head] <= i - k) head++;
    printf("%d ", a[q[head]]);
  }
}

int main() {
  scanf("%d%d", &n, &k);
  for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
  getmin();
  printf("\n");
  getmax();
  printf("\n");
  return 0;
}