前言

近日心血来潮想做一个开源项目，目标是做一款可以适配多端、功能完备的模板工程，包含后台管理系统和前台系统，开发者基于此项目进行裁剪和扩展来完成自己的功能开发。本项目为前后端分离开发，后端基于
Java21
和
SpringBoot3
开发，后端使用
Spring Security
、
JWT
、
Spring Data JPA
等技术栈，前端提供了
vue
、
angular
、
react
、
uniapp
、
微信小程序
等多种脚手架工程。

项目地址：
https://gitee.com/breezefaith/fast-alden

在前后端分离的项目开发过程中，我们通常会对数据返回格式进行统一的处理，这样可以方便前端人员取数据。但如果定义好响应对象
R
后，Controller类中每一个方法的返回值类型都只能是这个响应对象类，会使代码显得很不优雅。

@RestController
@RequestMapping("/admin")
public class AdminController {
    @PostMapping(value = "/register")
    public R<UmsAdmin> register(@Validated @RequestBody UmsAdminParam umsAdminParam) {
        return R.success(new UmsAdmin());
    }

    @PostMapping(value = "/logout")
    public R logout() {
        return R.success(null);
    }

    @PostMapping(value = "/login")
    public R login() {
        return R.success(new UmsAdmin());
    }
}

为了能够实现统一的响应对象，又能优雅的定义Controller类的方法，使其每个方法的返回值是其应有的类型，可以参考本文，主要是借助
RestControllerAdvice
注解和
ResponseBodyAdvice
接口来实现。

实现步骤

定义统一响应对象类

/**
 * 响应结果类
 *
 * @param <T> 任意类型
 */
@Data
public class ResponseResult<T> {
    /**
     * 响应状态码，200是正常，非200表示异常
     */
    private int status;
    /**
     * 异常编号
     */
    private String errorCode;
    /**
     * 异常信息
     */
    private String message;
    /**
     * 响应数据
     */
    private T data;

    public static <T> ResponseResult<T> success() {
        return success(HttpServletResponse.SC_OK, null, null);
    }

    public static <T> ResponseResult<T> success(T data) {
        return success(HttpServletResponse.SC_OK, null, data);
    }

    public static <T> ResponseResult<T> fail(String message) {
        return fail(HttpServletResponse.SC_INTERNAL_SERVER_ERROR, null, message, null);
    }

    public static <T> ResponseResult<T> fail(String errorCode, String message) {
        return fail(HttpServletResponse.SC_INTERNAL_SERVER_ERROR, errorCode, message, null);
    }

    public static <T> ResponseResult<T> success(int status, String message, T data) {
        ResponseResult<T> r = new ResponseResult<>();
        r.setStatus(status);
        r.setMessage(message);
        r.setData(data);

        return r;
    }

    public static <T> ResponseResult<T> fail(int status, String errorCode, String message) {
        return fail(status, errorCode, message, null);
    }

    public static <T> ResponseResult<T> fail(int status, String errorCode, String message, T data) {
        ResponseResult<T> r = new ResponseResult<>();
        r.setStatus(status);
        r.setErrorCode(errorCode);
        r.setMessage(message);
        r.setData(data);
        return r;
    }

}

定义一个忽略响应封装的注解

有些场景下我们不希望Controller方法的返回值被包装为统一响应对象，可以先定义一个忽略响应封装的注解，配合后续代码实现。

/**
 * 忽略响应封装注解
 */
@Target({ElementType.METHOD})
@Retention(RetentionPolicy.RUNTIME)
public @interface IgnoreRestControllerResponseAdvice {
}

实现ResponseBodyAdvice接口

本步骤需要使用
@RestControllerAdvice
注解，它是一个组合注解，由
@ControllerAdvice
、
@ResponseBody
组成，而
@ControllerAdvice
继承了
@Component
，因此
@RestControllerAdvice
本质上是个
Component
，用于定义
@ExceptionHandler
，
@InitBinder
和
@ModelAttribute
方法，适用于所有使用
@RequestMapping
方法。

还要用到
ResponseBodyAdvice
，它是Spring框架提供的一个接口，用于对Controller方法返回的响应体进行全局处理。它可以在Controller方法执行完毕并且响应体已经生成之后，对响应体进行自定义的修改或者增强操作。它本质上就是使用Spring AOP定义的一个切面，作用于Controller方法执行完成后。

具体实现代码如下：

/**
 * 响应实体封装切面
 */
@RestControllerAdvice(basePackages = {"com.demo.controller"})
public class GlobalResponseAdvice implements ResponseBodyAdvice<Object> {
    @Override
    public boolean supports(MethodParameter returnType, Class<? extends HttpMessageConverter<?>> converterType) {
        // 方法没有IgnoreRestControllerResponseAdvice注解，且response不是ResponseResult类型时启用beforeBodyWrite
        return !returnType.hasMethodAnnotation(IgnoreRestControllerResponseAdvice.class)
        && !returnType.getParameterType().isAssignableFrom(ResponseResult.class);
    }

    @Override
    public Object beforeBodyWrite(Object body, MethodParameter returnType, MediaType selectedContentType, Class<? extends HttpMessageConverter<?>> selectedConverterType, ServerHttpRequest request, ServerHttpResponse response) {
        // 如果返回值是void类型，直接返回200状态信息
        if (returnType.getParameterType().isAssignableFrom(void.class)) {
            return ResponseResult.success();
        }
        if (!(body instanceof ResponseResult)) {
            // warning: RestController方法上返回值类型为String时，响应的Content-Type是text/plain，需要手动指定为application/json
            if (body instanceof String) {
                try {
                    return JsonUtils.toJSON(ResponseResult.success(body));
                } catch (JsonProcessingException e) {
                    throw new RuntimeException(e);
                }
            }
            return ResponseResult.success(body);
        }
        return body;
    }
}

上述代码会对
com.demo.controller
包下所有的含有
@RequestMapping
注解的方法进行拦截，如果方法上没有
IgnoreRestControllerResponseAdvice
注解且返回值类型不是
ResponseResult
时，执行
beforeBodyWrite
方法。在
beforeBodyWrite
中将方法返回值包装为
ResponseResult
对象。

定义Controller类

下面我们就可以定义一个Controller类来进行简单的开发和测试。

@RestController
@RequestMapping("/demo")
public class DemoController {
    @GetMapping("/method1")
    public ResponseResult<Integer> method1() {
        return ResponseResult.success(100);
    }

    @GetMapping("/method2")
    public void method2() {

    }

    @GetMapping(value = "/method3")
    @IgnoreRestControllerResponseAdvice
    public String method3() {
        return "不会被封装";
    }

    /**
     * RestController中返回值类型是String的方法默认响应类型是text/plain，需要手动指定为application/json方可对其进行包装
     */
    @GetMapping(value = "/method4", produces = MediaType.APPLICATION_JSON_VALUE)
    public String method4() {
        return "会被封装";
    }
}

总结

本文介绍了SpringBoot项目中优雅地实现统一响应对象，如有错误，还望批评指正。

在后续实践中我也是及时更新自己的学习心得和经验总结，希望与诸位看官一起进步。

相关性描述的是
⼀个⽂档和查询语句匹配的程度
。我们从搜索引擎召回时，肯定希望召回相关性高的数据，那么如何来量化相关度呢。

首先，我们定义，一个文档doc，由多个词语 term 组成。

最早，通过最简单的TF-IDF来衡量。

TF-IDF

朴素的思想，相关度应该是词语权重、文档权重的融合。

• 词频 TF（Term Frequency）朴素的理解，一个词语term在一个doc中出现的频率越高，那么doc的相关性也越高。

$$ \text{TF}(t, d) = \frac{\text{词t在文档d中的出现次数}}{\text{文档d的总词数}} $$

• 逆向文档频率 IDF（Inverse Document Frequency），通俗的解释
  每个检索词在索引中出现的频率，频率越高，相关性越低。
  比如"的“ 这样的词语，可能每个doc都有，相反一些专业术语出现的文档数很少，这些文档的权重就很高。

$$ \text{IDF}(t, D) = \log\left(\frac{\text{文档集合D的总文档数}}{\text{包含词t的文档数} + 1}\right) $$

• TF和IDF结合起来计算一个词的TF-IDF值，公式为：

$$ \text{TF-IDF}(t, d, D) = \text{TF}(t, d) \times \text{IDF}(t, D) $$

TF-IDF
的优点在于简单有效，但它也有一些局限性，例如不考虑词序和上下文信息。

BM25（Best Matching 25）

BM25是一种用于文本信息检索的算法，是BM（Best Matching）系列中的一员。
它在信息检索领域中广泛应用，特别是在搜索引擎中用于排序搜索结果。整体而言
BM25 就是对 TF-IDF 算法的改进**，以下是BM25算法的公式：

$$ \text{BM25}(D, Q) = \sum_{i=1}^{n} \frac{{(k+1) \cdot f_i}}{{f_i + k \cdot (1 - b + b \cdot \frac{{\lvert D \rvert}}{{\text{avgdl}}})}} \cdot \log\left(\frac{{N - n_i + 0.5}}{{n_i + 0.5}}\right) $$

• D：文档
• Q：查询
• n：查询中的term数
• N：文档集中的文档数
• fi​：文档中term i的频率
• ni​：包含术语 term i 的文档数
• DL：文档长度（term数）
• AVG_DL：文档集的平均长度(term数量)

对于 TF-IDF 算法，
TF(t) 部分的值越大，整个公式返回的值就会越大，如果一个doc文章很长，词语很多，tf频率就会很大。
BM25 针对这个问题做了优化，通过b参数，对文档长度进行打压，
随着TF(t) 的逐步加大，该算法的返回值会趋于一个数值。

BM25的优势在于它对于长文本和短文本的处理更为灵活，并且能够适应不同查询的特征。这些可调整的参数使得BM25能够通过调整来适应不同的信息检索场景。

• b 参数 ，b 默认 0.75（经验值），主要是对长文档做惩罚，如果不希望文档长度更大的相似度搜索更好，可以把 b 设置得更大，如果设置为 0，文档的长度将与分数无关。从下图可以看到，b=0时，L与分数无关，b=1时，L越大，分数打压越厉害。
  

• k 参数， 默认值 1.2，会影响词语在文档中出现的次数对于得分的重要性，如果希望词语出现次数越大，文档越相关，这个参数可以设置更大。

BM25+

BM只考虑了term和doc的维度，对query 里term的频率没有考虑，BM25+（Best Matching 25 Plus）正是基于这一点，来改进BM25算法，以下是BM25+的公式，以及每个参数的含义：

$$ \text{BM25+}(D, Q) = \sum_{i=1}^{n} \frac{{(k_1+1) \cdot f_i \cdot (k_3+1) \cdot qf}}{{(f_i + k_1 \cdot (1 - b + b \cdot \frac{{\lvert D \rvert}}{{\text{avgdl}}})) \cdot (k_3 + qf)}} \cdot \log\left(\frac{{N - n_i + 0.5}}{{n_i + 0.5}}\right) $$

相比BM25，增加k3 和qf，
其中：

• qf：查询中词项的频率（Query Term Frequency）
• k3：控制查询中词项频率对得分的影响程度的参数，下面是不同的k3，在不同的query weight情况下，对分数的影响

![[Pasted image 20240202151809.png]]

• 注意k1就是BM25里的k

这里的qf，我们可以当做term weight来使用，训练term weight模型，来对重点的term做激励。

谷歌的End-to-End Query Term Weighting

谷歌在2023年发布了一篇论文，介绍如何端到端学习term weighting。

背景是term based recall的方法相对于embedding recall来说很繁琐，好处就是时延低且对基建要求低，向量召回会遇到分不清楚query中哪个词是核心词的问题，导致召回出现了非核心词的结果。

谷歌论文

在左侧，展示了传统的信息检索（IR），所有的term都是默认的权重，在右侧，我们插入了一个BERT模型来评估term的权重，BM25+ 打分时将这个weighting作为query freq，就是上面说的qf。

可以看谷歌论文的打分：

上面的$f(T_i, T,w)$ 就是模型学习的权重。

附录

可视化不同参数变化，BM25分数的变化

import numpy as np
import matplotlib.pyplot as plt


def calculate_bm25(tf, term_weight, corpus_freq, total_docs,
                   k1=1.5, k3=8, doc_len=100, avg_doc_len=120, b=0.75):
    idf = np.log((total_docs - corpus_freq + 0.5) / (corpus_freq + 0.5) + 1)  # IDF计算
    K = k1 * ((1.0 - b) + b * doc_len / avg_doc_len + tf)
    tf_term = tf * (k3 + 1) * term_weight / (K * (k3 + term_weight))  # TF计算

    return idf * tf_term  # BM25计算


def visualize_bm25_scores_k3(term_weight_values, tf_values, corpus_freq, total_docs, k1=1.5):
    # Plotting
    plt.figure(figsize=(12, 8))

    for k3 in np.linspace(2, 10, 5):
        scores = [calculate_bm25(5, query_weight, corpus_freq, total_docs, k1=k1, k3=k3)
                  for query_weight in term_weight_values]
        plt.plot(term_weight_values, scores, label=f'k3={k3}')

    plt.title('BM25 Scores with Different Query_freq Values (b=0)')
    plt.xlabel('Query weight')
    plt.ylabel('BM25 Score')
    plt.legend()
    plt.show()

def visualize_bm25_scores_b(term_weight_values, tf_values, corpus_freq, total_docs, k1=1.5):
    # b 用于惩罚文档长度
    plt.figure(figsize=(12, 8))

    doc_lens = np.linspace(100, 1000, 10)
    avg_doc_len = 100
    L = [doc_len /avg_doc_len for doc_len in doc_lens]
    for b in [0, 0.5, 0.75, 1.0]:
        scores = [calculate_bm25(5, 1, corpus_freq, total_docs, k1=k1, b=b, doc_len=doc_len)
                  for doc_len in doc_lens]
        plt.plot(L, scores, label=f'b={b}')

    plt.title('BM25 Scores with Different b')
    plt.xlabel('L(doc_length/avg_doc_length')
    plt.ylabel('BM25 Score')
    plt.legend()
    plt.show()



term_weight_values = np.linspace(0.1, 10, 10)  # query weight
corpus_freq = 500  # 包含term的文档数量，文档频率
total_docs = 5000000  # 总文档数量
tf_values = list(range(1, 101))
visualize_bm25_scores_k3(term_weight_values, tf_values, corpus_freq, total_docs)

visualize_bm25_scores_b(term_weight_values, tf_values, corpus_freq, total_docs)

本篇是mygin这个系列的最后一篇。如果想自己动手实现一个类似Gin的Web框架，建议从
mgin第一篇
开始，
总代码行数有效行数只有600多行

github
源码
mygin

目的

• 实现错误处理机制

panic简介

在实现错误处理机制之前，我们知道在
Go
中，错误的处理方式依靠return返回，由调用者处理。如果是不可恢复的错误，可以手动抛出错误，当然在实际运行中，也会遇到不可处理的错误，比如说除数为0的时候
painc
也会触发。终止当前的程序。

手动触发错误

//main.go
func main() {
  fmt.Println("before ....")
  panic("some err message")
  fmt.Println("after ....")
}

• shell

~  go run cmd/main.go
before ....
panic: some err message

goroutine 1 [running]:
main.main()
        /var/www/gophp/cmd/main.go:12 +0x5f
exit status 2

程序发生错误

func main() {
  i := []int{0}
  j := 2
  fmt.Println(j / i[0])
}

• shell

go run cmd/main.go
panic: runtime error: integer divide by zero

goroutine 1 [running]:
main.main()
        /var/www/gophp/cmd/main.go:13 +0x194
exit status 2

painc
的介绍到此为止，也就是发生了
painc
错误，对程序来说就是比较严重的，就会中断，但是这个时候需要捕获到这个
painc
错误，
Go
中没有
try...Catch...
又想捕获错误，这个时候就该
defer
出场了。

defer简介

panic
会导致程序被中止，但是在程序终止前，会先把当前协程上已经
defer
的任务，执行完成后再终止。效果类似于其他语言的
try...catch
。

示例

//cmd/main.go

func main() {
  defer func() {
  fmt.Println("defer....")
  }()
  i := []int{0}
  j := 2
  fmt.Println(j / i[0])
}

• shell

~ go run cmd/main.go
defer....
panic: runtime error: integer divide by zero

goroutine 1 [running]:
main.main()
        /var/www/gophp/cmd/main.go:16 +0x1d3
exit status 2

可以看到，在程序退出之前先执行了
defer
中的函数。如果能在
defer
中捕获
painc
错误，那么就能实现其他语言的
try...catch

recover函数

介绍之前看官方文档怎么定义的

// The recover built-in function allows a program to manage behavior of a
// panicking goroutine. Executing a call to recover inside a deferred
// function (but not any function called by it) stops the panicking sequence
// by restoring normal execution and retrieves the error value passed to the
// call of panic. If recover is called outside the deferred function it will
// not stop a panicking sequence. In this case, or when the goroutine is not
// panicking, or if the argument supplied to panic was nil, recover returns
// nil. Thus the return value from recover reports whether the goroutine is
// panicking.
翻译后...
recover 内置函数允许程序管理恐慌的 goroutine 的行为。在延迟函数（但不是它调用的任何函数）
中执行恢复调用会通过恢复正常执行来停止 panic 序列，并检索传递给 panic 调用的错误值。如果
在延迟函数之外调用 recover，则不会停止紧急序列。在这种情况下，或者当 goroutine 没有 
panic 时，或者如果提供给 panic 的参数为 nil，则 recover 返回 nil。因此，recover 的
返回值报告 goroutine 是否处于恐慌状态。

应用

• cmd/main.go

func main() {
  defer func() {
    fmt.Println("defer....")
    if err := recover(); err != nil {
      fmt.Printf("recover 到错误信息：%s\n", err.(error).Error())
	  fmt.Println("recover success")
     }
  }()
  i := []int{0}
  j := 2
  fmt.Println(j / i[0])
  fmt.Println("after panic")

• shell

~ go run cmd/main.go
defer....
recover 到错误信息：runtime error: integer divide by zero
recover success

可以看到已经捕获到错误信息了，程序正常结束。 after panic 没有打印，这是正确的，当
panic
被触发时，程序的执行栈就到定义的
defer
函数。就像在try代码块中发生了异常，执行栈来到 catch，接下来执行 catch 代码块中的代码。而在 main() 中打印了 recover success，说明程序已经恢复正常，继续往下执行直到结束。

mygin的错误处理

对一个 Web 框架而言，错误处理机制是必要的。如果发生了painc错误，应当返回错误信息或告诉对方失败了，而不至于什么都不返回，对调用方十分不友好。
例如我有如下的逻辑

package main

import (
	"fmt"
	"github.com/scott-pb/mygin"
	"net/http"
	"strconv"
)

func main() {

	r := mygin.Default()
	group := r.Group("/api")
	group.GET("/recovery/:index", func(c *mygin.Context) {
		index := c.Params.ByName("index")
		i, _ := strconv.ParseInt(index, 10, 10)
		s := []int{1, 3, 5, 7, 9}

		c.String(http.StatusOK, fmt.Sprintf("index:%d result:%d success!\n", i, s[i]))
	})

	err := r.Run(":8088")
	if err != nil {
		fmt.Println(err)
	}
}

根据调用方传递的index下标函数数组中的值，当传递的index超过数组下标时，就会发生painc错误，如果执行/api/recovery/10 时，就会发生错误。调用方什么返回都没有。
这个时候，就需要添加一个错误的处理机制，当错误发生时，向调用方返回
Internal Server Error
，且打印必要的错误信息，方便进行错误定位。

实现中间件 Recovery

新增文件 mygin/recovery.go，在这个文件中实现中间件 Recovery

• mygin/recovery.go

package mygin

import (
	"fmt"
	"net/http"
)

// Recovery 发生错误时，恢复函数，且返回相应错误信息。
func Recovery() HandlerFunc {
	return func(c *Context) {
		// 使用defer延迟执行，以便在函数退出时进行recover
		defer func() {
			if err := recover(); err != nil {
				// 如果发生panic，打印错误信息并返回500 Internal Server Error响应
				fmt.Println(err.(error).Error())
				c.Writer.Write([]byte("Internal Server Error\n"))
				c.status = http.StatusInternalServerError
				c.Abort() // 终止后续中间件的执行
			}
		}()

		c.Next() // 调用下一个中间件或处理函数
	}
}

Recovery
方法很简单，使用
defer
挂载上错误恢复的函数，在这个函数中调用
recover方法
，捕获
panic
，打印错误信息，并且向调用方返回 Internal Server Error。

• mygin/engine.go

// Default 返回一个默认的引擎实例
func Default() *Engine {
	engine := New()

    //Logger Recovery 中间件
	engine.Use(Logger(), Recovery())

	// Group 保存 engine 的指针
	engine.RouterGroup.engine = engine

	return engine
}

接下来就是测试了

测试

• shell

~ curl http://127.0.0.1:8088/api/recovery/1
index:1 result:3 success!
~ curl http://127.0.0.1:8088/api/recovery/10
Internal Server Error

• 控制台输出
  
  可以看到第一次请求返回200成功，第二次请求，先打印了错误信息，然后对应的返回500错误，且调用方接收到了Internal Server Error。

技术背景

所谓的分子力场，就是用一些计算量较小的函数，来拟合并替代一部分传统第一性原理计算的结果。这个结果，包含了势能和作用力，再用朗之万动力学进行演化，这才使得我们可以在计算机上模拟一个分子动力学的过程。否则在第一性原理计算的框架下，要想获得动力学统计的信息，是非常困难的。

分子力场，常见的有成键相互作用、非成键相互作用以及多体相互作用。本文主要解释一下其中的成键相互作用的Bond Energy和Angle Energy这两项，并给出一些简单的计算演示。

Bond Energy

关于键能的定义，一般采用谐振势的形式，这里用
\(\textbf{r}_B\)
和
\(\textbf{r}_A\)
来分别表示原子
\(B\)
和
\(A\)
的空间坐标：

这里有个细节是，在很多非成键相互作用中，会采用无穷远处的势能为0。而成键相互作用，例如这里的键能，则是采用了平衡位置
\(|\textbf{r}_B-\textbf{r}_A|=b\)
处的势能为0。因为自然作用力下，物体会从高势能位置向低势能位置运动，因此这里
\(E_{bond}\)
采取的是正号。而计算作用力
\(\textbf{F}_{bond}=-\frac{d E_{bond}}{d l}\)
时，前面会带有一个负号。其表示的意义是：当键长大于平衡键长时，通过成键相互作用力把键长拉短。而势能表达式前的系数
\(\frac{1}{2}\)
主要是为了使得力的形式符合胡克定律，更方便理解，但实际上很多力场中把这个系数融入到了参数
\(k_b\)
中，且
\(k_b>0\)
。

由于我们在动力学模拟的过程中，保存的轨迹信息只有分子的空间坐标和速度的信息，因此我们可以把
\(\textbf{F}_{bond_i}\)
写成
\(\textbf{F}_{bond_i}(\textbf{r}_A)\)
和
\(\textbf{F}_{bond_i}(\textbf{r}_B)\)
的函数形式，其中
\(i\)
表示的是第
\(i\)
根键带来的作用力，而一个原子可能同时有多根键连接，因此每一个原子上的作用力最后需要把所有相关的键都加和起来：

因此，要求解力，首先要求解每一根键上的作用力。由于是同一根键相连的两个原子，所以对应的力应该是大小相等、方向相反：
\(F_{bond_i}(\textbf{r}_A)=-F_{bond_i}(\textbf{r}_B)\)
。那么理论上说，我们只要保存其中的一个力向量
\(F_{bond_i}\)
就够了。那么令
\(l=|\textbf{r}_B-\textbf{r}_A|\)
，有第
\(i\)
根键上的作用力：

这里有个问题是，上述公式计算出来的作用力的符号只表示用于拉伸键的方向，并不表示实际空间中对应单个原子的受力方向。但我们要计算
\(\textbf{F}_{bond_i}(\textbf{r}_A)\)
怎么办呢？其实拉伸的方向就是每一个原子相对于键的中点
\(\frac{\textbf{r}_B+\textbf{r}_A}{2}\)
的向量。以原子
\(A,B\)
和第
\(i\)
根键为例，原子
\(A\)
相对于
\(A\)
和
\(B\)
的中点
\(C\)
的向量为：
\(\textbf{r}_{AC}=\frac{\textbf{r}_B+\textbf{r}_A}{2}-\textbf{r}_A=\frac{\textbf{r}_B-\textbf{r}_A}{2}\)
，单位化之后就是
\(\textbf{e}_{AC}=\frac{\textbf{r}_B-\textbf{r}_A}{|\textbf{r}_B-\textbf{r}_A|}\)
。类似的，
\(\textbf{e}_{BC}=\frac{\textbf{r}_A-\textbf{r}_B}{|\textbf{r}_B-\textbf{r}_A|}=-\textbf{e}_{AC}\)
。

结合上述分析的关于作用力的符号的物理意义，当力为负数时，表示跟
\(\textbf{r}_{AC}\)
和
\(\textbf{r}_{BC}\)
同向，当力为正数时，表示跟
\(\textbf{r}_{AC}\)
和
\(\textbf{r}_{BC}\)
反向：

同理可以得到：

有了力的表达式之后，可以探讨一下计算Bond Force在CPU上的时间复杂度。假设有
\(N\)
个原子和
\(P\)
对键连原子，在一个3维空间下，首先要计算一个矢量
\(\textbf{r}=\textbf{r}_A-\textbf{r}_B\)
，这需要
\(3P\)
的加法计算量。然后是计算这些向量的模
\(|\textbf{r}|\)
，需要做
\(3P\)
次乘法运算和
\(2P\)
次加法运算以及
\(P\)
次指数运算，一共是
\(6P\)
次运算。再计算
\(k_b(\textbf{r}_B-\textbf{r}_A)\)
再计算
\(k_bb(\textbf{r}_B-\textbf{r}_A)\)
，一共是
\(6P\)
次乘法运算。最后要把得到的
Pair-Wise-Force
加和到
Atom-Wise-Force
上，需要做
\(6P\)
次加法运算。最后，我们统计出来计算Bond Force的复杂度为
\(\Omega(21P)\)
。

接下来要说的是一个关于代码实现上的问题。如果我们选择用CPU上（如果选择用GPU计算，那么评估标准会有所不同）的Python来计算一个bond force的话，那么就可能有这么几种思路：

1. 用现有的自动微分框架，只需要实现一个
   bond_energy
   的计算函数，就可以自动求得每个原子上的力；
2. 用for循环来实现一个
   bond_force
   函数，使得计算满足最低的复杂度要求；
3. 用Python的向量化运算来实现一个
   vector_bond_force
   函数，可以直接调用比较成熟的、经过多种优化的向量化运算函数。

自动微分实现

这里我们使用MindSpore框架的函数式编程来实现bond energy的自动微分：

import numpy as np
np.random.seed(1)
import mindspore as ms
from mindspore import context, Tensor, grad
from mindspore import numpy as msnp
# 动态图模式，在CPU上运行
context.set_context(mode=context.PYNATIVE_MODE, device_target="CPU")
# 键能函数，这里简化k为1，b为0
def bond_energy(crd, bonds, k=1., b=0.):
    bond_vector = crd[bonds]
    energy = 0.5 * k * (msnp.norm(bond_vector[:, 1] - bond_vector[:, 0], axis=-1) - b) ** 2
    return energy.sum()
# 初始化10个原子的坐标信息
N = 10
crd = Tensor(np.random.random((N, 3)), ms.float32)
bonds_np = np.random.randint(0, N, (N*2, 2))
bonds = Tensor(bonds_np, ms.int32)
bond_force = grad(bond_energy)
adf_force = -bond_force(crd, bonds)
print (adf_force)

运行输出结果如下所示：

[[-0.6740933  -0.59173465  0.8205794 ]
 [ 0.31077877  0.7709253   0.77850956]
 [ 1.2435249   0.717391   -0.49775103]
 [-1.4042054  -0.39493966 -0.07195999]
 [ 0.84755206 -1.7858155   2.0126157 ]
 [-0.9960958   0.24160932 -0.28924745]
 [ 2.0020847   1.6122136  -1.6734666666 ]
 [-3.06986    -0.14928058  0.0689815 ]
 [-2.0460608  -2.140861    1.438695  ]
 [ 3.7863746   1.7204924  -2.5869672 ]]

最后我们得到的是每一个原子的受力总和，用自动微分框架实现的最大好处就是：方便、简单。我们甚至都不需要去手动推导bond force的解析形式和各种正负符号，只需要把bond energy函数写对就可以了，这也是对科研人员的一大福音。而这其中的计算的效率，那就只能依赖于AI框架本身的实现和优化。

For循环实现

在前面的推导中我们已经得到了成键相互作用力的计算公式，那么我们只需要先计算出每一个键上的受力，然后用一个for循环加到对应的原子位置上即可：

import numpy as np
np.random.seed(1)
# for循环实现的力函数
def bond_force(crd, bonds):
    bond_crd = crd[bonds]
    bond_vector = bond_crd[:, 1] - bond_crd[:, 0]
    force = np.zeros_like(crd)
    for i, bond in enumerate(bonds):
        force[bond[0]] += bond_vector[i]
        force[bond[1]] -= bond_vector[i]
    return force
# 初始化10个原子的坐标
N = 10
crd = np.random.random((N, 3)).astype(np.float32)
bonds_np = np.random.randint(0, N, (N*2, 2))
loop_force = bond_force(crd, bonds_np)
print (loop_force)

得到的结果如下：

[[-0.6740933  -0.59173465  0.8205794 ]
 [ 0.31077877  0.7709253   0.77850956]
 [ 1.2435249   0.717391   -0.49775103]
 [-1.4042056  -0.39493972 -0.07195997]
 [ 0.84755206 -1.7858155   2.0126157 ]
 [-0.9960958   0.2416093  -0.28924745]
 [ 2.002085    1.6122137  -1.6734666666 ]
 [-3.06986    -0.14928058  0.0689815 ]
 [-2.046061   -2.140861    1.4386952 ]
 [ 3.786375    1.7204924  -2.5869677 ]]

我们可以看到，得到的结果是跟自动微分框架计算出来的保持一致的。

向量化运算

在这几个方案中，向量化运算在计算复杂度上是不占优势的，这是因为成键相互作用力的向量化运算要求
双向索引
和
静态Shape
。在这个问题中，比较难的一点是，每一个原子的键连数量都是不一样的，所以要做一个向量化的实现很有可能是吃力不讨好的，这里仅提供一个可以参考的实现方案：

import numpy as np
np.random.seed(1)
# for循环实现的力函数
def vector_bond_force(crd, bonds, atom_bond, atom_bond_mask, atom_bond_zeros_mask):
    bond_crd = crd[bonds]
    bond_vector = bond_crd[:, 1] - bond_crd[:, 0]
    force = (bond_vector[atom_bond] * atom_bond_zeros_mask[..., None] * atom_bond_mask[..., None]).sum(axis=-2)
    return force
# 初始化10个原子的坐标
N = 10
crd = np.random.random((N, 3)).astype(np.float32)
bonds_np = np.random.randint(0, N, (N*2, 2))
# 计算单原子最大成键数量
max_bonded_atom = np.bincount(bonds_np.reshape(-1)).max()
print ('The maximum bonds per atom is: {}'.format(max_bonded_atom))
# 初始化一个反向索引矩阵
atom_bond = -1 * np.ones((N, max_bonded_atom)).astype(np.int32)
# 初始化一个反向标记矩阵，用于区分力的方向
atom_bond_mask = np.zeros((N, max_bonded_atom))
for i, bond in enumerate(bonds_np):
    # 在每一轮的迭代中，要roll一个位置，把上一步更新的元素轮转到第一位，这样就可以一直轮转更新
    atom_bond[bond[0]] = np.roll(atom_bond[bond[0]], 1)
    atom_bond_mask[bond[0]] = np.roll(atom_bond_mask[bond[0]], 1)
    atom_bond[bond[0]][-1] = i
    # 正方向标记为正1
    atom_bond_mask[bond[0]][-1] = 1.
    atom_bond[bond[1]] = np.roll(atom_bond[bond[1]], 1)
    atom_bond_mask[bond[1]] = np.roll(atom_bond_mask[bond[1]], 1)
    atom_bond[bond[1]][-1] = i
    # 逆方向标记为负1
    atom_bond_mask[bond[1]][-1] = -1.
# 初始化一个冗余标记矩阵，用于标记反向索引矩阵中的冗余位置
atom_bond_zeros_mask = np.where(atom_bond >= 0, 1, 0)
vector_force = vector_bond_force(crd, bonds_np, atom_bond, atom_bond_mask, atom_bond_zeros_mask)
print (vector_force)

在这个向量化运算的代码实现中，有个比较显著的缺点是，我们需要给力的计算函数传入三个额外的参量：
反向索引矩阵、方向标记矩阵和冗余标记矩阵
。从形式上来说，不仅仅是算法变得更加复杂了，而且可读性也并不是很好，但是这里面包含了几个值得思考的点：

1. 对于分子动力学模拟中可能出现的大部分体系而言，成键关系总是稀疏的；
2. 额外的三个初始化操作是在力的计算函数外部完成的，并不参与力的迭代过程；
3. 向量化的运算，并且排去除了自动微分可能带来的一些性能上的损失，在个别大规模体系下，预期能够展现较好的性能增益效果。

上述代码的运行结果如下图所示：

The maximum bonds per atom is: 8
[[-0.67409331 -0.59173465  0.82057939]
 [ 0.31077877  0.77092531  0.77850953]
 [ 1.24352494  0.71739101 -0.49775103]
 [-1.40420556 -0.39493971 -0.07195997]
 [ 0.84755203 -1.78581548  2.01261576]
 [-0.99609581  0.24160931 -0.28924745]
 [ 2.00208497  1.6122137  -1.67346666664]
 [-3.0698598  -0.14928058  0.0689815 ]
 [-2.04606098 -2.14086115  1.43869516]
 [ 3.78637475  1.72049223 -2.58696735]]

在这个案例中，由于成键关系只是随机生成的，所以出现了一个非常尴尬的最大键连数量，这种键连数量在正常的分子动力学模拟中也是不太可能出现的。也正是因为各种模拟体系的特殊性，所以这里暂不对以上的几个算法进行进一步的性能测试，我们只是介绍一下这些计算方案的可行性，同时对于分子之间键连关系的计算进行一定的思考。

Angle Force

要计算键角相互作用力，我们首先回顾一下两个向量之间夹角的计算方法：

其中
\(\textbf{e}_{ji}=\frac{\textbf{v}_{ji}}{|\textbf{v}_{ji}|},\textbf{v}_{ji}=\textbf{r}_i-\textbf{r}_j\)
。那么，如果把
\(\theta_{ijk}\)
类比于前面介绍的键连关系中的
\(bond_{AB}\)
，我们也可以使用一个谐振势的力场形式来计算键角势能：

类似于前面键作用力的推导，这里我们可以得到角作用力的形式：

对于一个键相互作用力而言，我们比较容易可以理解，力的符号用于表示将键拉伸或者缩短。而这里键角作用力也是如此，如果计算出来的力为负数，就是往键角变小的方向演化，如果力为正数，则表示把键角拉大。而如果力为0，按照谐振势的定义，也就是处于平衡位置的时候，此时键角势能
\(E(\theta_{ijk})\)
的值也是0。接下来我们关注一个力的方向问题，对于Bond Force我们很容易可以理解，力的方向自然是沿着键的方向，那么
键角作用力是沿着什么方向
的呢？按照经验来说，只有三种可能性，一种是沿着
\(e_{ik}\)
的方向，一种是分别垂直于
\(e_{ji}\)
和
\(e_{jk}\)
的方向，还有一种是沿着对向成键的方向。接下来要通过计算来确认一下单个原子的受力方向，我们先计算一下作用在原子
\(i\)
上的作用力：

这里面还涉及到一个向量点乘的求导，例如
\(f=\textbf{v}_i\cdot\textbf{v}_j=x_ix_j+y_iy_j+z_iz_j\)
，那么对于其中一个向量的求导形式为：

所以向量求导也是符合我们日常所用到的求导法则的。另外还涉及到一个归一化的求导函数：

由于
\(\textbf{e}_i\)
的长度始终是保持不变的，按照导数的定义有

另外，由于键矢量的定义
\(\textbf{v}_{ji}=\textbf{r}_i-\textbf{r}_j\)
，因此可以得到关系：

那么，通过观察
\(i\)
原子的键角作用力，我们可以发现，
\(\left.\textbf{F}_{\theta_{ijk}}\right|_i\)
的方向仅仅与
\(\textbf{v}_{jk}\)
有关，也就是说，第
\(i\)
个原子受到的键角作用力是沿着
\(e_{jk}\)
方向的。类似的我们可以推导出来第
\(k\)
个原子的键角作用力：

可以看出，第
\(k\)
个原子受到的键角作用力是沿着
\(e_{ji}\)
方向的。而第
\(j\)
个原子的受力如何呢？已知
\(j\)
原子的受力跟
\(\textbf{v}_{ji}\)
和
\(\textbf{v}_{jk}\)
都相关，则有：

可见，每一个键角所对应的键角作用力都是保守的，并且
\(i,k\)
原子的受力方向沿着对向成键的方向，
\(j\)
原子的受力则是沿着平行四边形的对角线方向。得到了每一个键角所产生的作用力之后，因为一个原子可能同时对应多个键角，因此，我们计算单个原子在所有键角作用力下的作用时，需要对所有相关的键角进行加和：

类似的，我们也可以统计一下键角这部分的计算复杂度。假设有
\(S\)
组键角，每组键角由3个原子组成。需要注意的是，这里所有构成键角的键矢量和单位键矢量，都已经在计算Bond Energy的时候算过了，理论上说是可以不用再进行重复计算的。那么就可以直接计算单位键矢量的两两点乘，这里面使用到了
\(3S\)
次计算，然后计算
\(S\)
次反余弦函数，就可以得到
\(S\)
个键角的值。有了键角值，代入力场参数，计算
\(2S\)
次就可以得到每个键角上的作用力大小，再计算
\(S\)
次可以得到
\(1-\textbf{e}_{ji}\cdot\textbf{e}_{jk}\)
的值，再经过
\(S\)
次计算可以得到
\(\textbf{F}_{\theta_{ijk}}|_i\)
和
\(\textbf{F}_{\theta_{ijk}}|_k\)
的前置系数。将这些系数跟每一个键矢量做点乘得到
\(\textbf{F}_{\theta_{ijk}}|_i\)
和
\(\textbf{F}_{\theta_{ijk}}|_k\)
，需要做
\(6S\)
次运算，
\(\textbf{F}_{\theta_{ijk}}|_j\)
是两者加和的相反数，需要
\(3S\)
次运算。最后要把这些键角对应的三个原子的力矢量加回到每个原子的受力上，最少需要做
\(9S\)
次运算。整体统计下来，键角作用力的计算复杂度为
\(\Omega(26S)\)
。那么，对于键长和键角这两项力场作用项而言，不可减免的计算复杂度就有
\(\Omega(21P+26S)\)
。对应到具体实现中，根据不同的实现方案，计算量只会大于这个数量级。而且需要注意的是，这里只统计了计算复杂度，实际的分子动力学模拟过程中，还会涉及到大量的索引和通信复杂度，这里不做赘述。关于Angle Force部分的代码实现先暂时忽略，因为本质上跟计算Bond Force的原理是一样的，既可以选择简单方便的自动微分，也可以自定义一些反向索引和Mask用矢量化计算的形式来实现。

总结概要

本文介绍了在分子力场中经常有可能被使用到的键长和键角项的谐振势模型，并且分别从自动微分的Python代码实现以及解析形式的矢量化编程形式给出了初步的实现方案。虽然力场形式较为简单，但是在实际的计算中，我们统计出来，至少需要21P+26S的计算量，其中P指键的数量，S指键角的数量。这里提到的矢量化计算的实现方案，虽然从计算的角度来说有大量的冗余，但由于一般情况下，一个分子系统单个原子的成键数量都在4以内（比如C原子的sp3杂化），因此矢量化计算的实现方案也不失为一个很好的参考。

版权声明

本文首发链接为：
https://www.cnblogs.com/dechinphy/p/bond-angle.html

作者ID：DechinPhy

更多原著文章：
https://www.cnblogs.com/dechinphy/

请博主喝咖啡：
https://www.cnblogs.com/dechinphy/gallery/image/379634.html

1：推荐系统简易架构图

2：推荐引擎流程图

Mixer系统的用户推荐就是按照上面流程实现，接下来分别介绍各个功能

3：索引内容

1：索引目的是提供高效查询

索引内容就是将内容生产索引文件，目的是提高查询速度，不可能每个请求过来都把所有的数据查出来，然后过滤、计算特征、排序、最后将排在前面的数据推给用户，数据量实在是太多不能这么玩，所以会将一些关键数据生成索引文件，通过databus将文件推送到推荐引擎的机器上，推荐引擎加载索引文件，在内存中生成各种索引（正排索引、倒排索引等），提供快速查询。

2：怎么生成索引文件

我们把一个表的数据查出来，写入文件，那么这个文件就是索引文件，但是这个表的一些字段我们可能用不上，查出来只会增大索引文件，影响加载速度和占用多余的内存，所以我们只会将需要的字段加到索引文件中，mixer作为用户推荐系统，那么用户表是主表，还设计到用户资料表、用户兴趣表等多张表，将多张表信息和成一张表，一个用户一条记录，定义好对应的struct，填充每个用户struct中的字段，得到一个struct数组，将数组的内容以pb格式序列化写入文件，索引文件就生成了。

由于数据库中的数据会实时变化，那么我们就需要定时生成索引文件。

3：databus传输文件

索引文件我们会用一个独立的服务生成，然后通过databus服务传送到推荐引擎的机器上。这个功能可以很简单，如果你的推荐系统比较简单，生成索引的服务和推荐的服务在一台机器上，约定一个文件路径就行，定时检查文件是否变更，变更重新加载就行。

那么正常的databus怎么实现？

数据巴士就是一个文件下载服务，你可以提供ftp服务或是http服务，下面以http服务为例简要说明databus的实现，他真的太简单了。databus提供一个server端和一个client端。server定时判定索引文件有没有更新，如果有更新，更新文件路径和文件的md5，每个推荐引擎机器部署一个client，client定时请求server端，比对client和server的md5是否一致，如果不一致说明文件更新了，就下载文件。可能你觉得定时请求有些没必要，想减少没必要的http请求，那么server只需要发现文件变更就通知所有客户端，你可以用etcd作为通知工具，但我们也要考虑引入一个中间件带来的问题，有可能是中间的问题，或是你使用的问题导致client没有加载最新的索引文件。

client下载索引文件放到指定路径，推荐引擎检测到文件变更，加载最新文件，实现索引文件传输。

4：内存中建立索引

推荐引擎加载索引文件建立内存索引，这块是核心。

索引文件我们可以理解成一张表，然后我们在内存中给这张表建立索引，比如我想按照性别查询、按照学历查询、按照年龄范围查询、按照城市id查询，我们可以将这4个查询分成2类，一类是这个字段的值是固定的几个（性别：0无，1男，2女；学历：0~7），第二类字段的值是多个的（年龄：18~150；城市id可能有几百个）

第一类的倒排索引可以用BitMap实现。

以性别为例，性别有3个值，我们建立3个BitMap为[3]BitMap，第一个性别bitmap表示无，第二个BitMap表示男，第三个BitMap表示女。比如我们有100W个用户，那么一个BitMap就是一个长度为15625的[]uint64，一个bit位对应一个用户，那么这个bitmap占用15625*8=122K，3个bitmap占用366K，占用内存是非常少的。

type BitMap []uint64


func NewBitMap(length int) BitMap { bitLen := (length + 63) / 64 return make([]uint64, bitLen) }
bt := NewBitMap(1000000) // 15625 = (1000000 +63) / 64

BitMap怎么赋值？

如第100个用户是男性，男性对应第二个bitmap，将[]uint64的100位设置成1，是第100位不是数组下标100。遍历这100w个用户，按照性别和下标设置对应bitmap的二进制位的值，把bitmap看成二进制数组就行，它的值只有0和1。

BitMap怎么查询？

如果我要给一个用户推荐男性用户，那么第二个BitMap中二进制位为1的下标的集合就是100W个用户组数中男性的集合，当然我们不需要现在就取出所有男性的信息，因为查询条件可能是：男性+本科+25岁+北京。索引性别只需要取第二个BitMap就行，索引查询性别就是取数组下标为1的元素，还有比这还快的查询吗？

同样的学历定义为bt := [8]BitMap，本科只需要取bt[4]，快得飞起，然后把学历和性别的两个BitMap取交集，就是SQL中的and，得到一个BitMap，二进制位1的用户就是我们要找的用户。时间复杂度是O(1)

第一类的倒排索引可以用排序+BitMap实现。

像年龄有100多个，城市有几百个，我们不能建立几百个BitMap的数组，那样占用内存较多，而且这些值随时变多，用固定长度数组不好解决，我们只能提前排好序，通过二分查找+遍历的方式实现，二分查找时间复杂度O(log n)。

比如年龄我们可以定义一下结构

type IndexAge struct {
Index uint32 //数组下标值，数组为100W用户数组
Age uint32 //用户对应的年龄
}

100W用户就定义[1000000]IndexAge，然后Age排序，如何按照年龄范围查找，以查找年龄范围是[18,25]为例，通过二分查找找到第一个年龄=18的用户，然后往后遍历直到年龄>25为止，Index指定的下标用户几位所求。类似mysql的between and。年龄查询时间复杂度O(n log n)。

不用BitMap，用数组也行，但每个索引的内存占用会扩大到64倍

5：召回

召回就是通过内存索引进行查询，然后取交集，如查找：性别男 + 学历（本科/硕士）+年龄（18~25）+城市（北京/上海）。

性别男，直接取下标，得到一个BitMap

学历（本科/硕士），直接取2次下标，得到本科BitMap和硕士BitMap，两个BitMap取并集，得类似mysql中的or，到一个BitMap

年龄（18~25），二分查找+遍历，类似mysql的between and，得到一个BitMap

城市（北京/上海），两次二分查找+遍历，类似mysql的in，得到两BitMap后取并集得到一个BitMap

然后4个BitMap取交集得到一个BitMap，BitMap为1的就是我们要找的用户。

6：过滤

作为用户推荐系统

1：推荐过的用户一天总不能推2次吧

2：不能把自己推荐给自己

3：我拉黑的人不能推荐给我

4：把我拉黑的人不能推荐给我

5：我喜欢的人，就不要再推荐给我了

6：………..

这些数据都会存在数据库中，同时在redis中保存一份对应的uid，这些都是要过滤的用户，通过查询redis得到所有要过滤的uid，将这些uid存到map中，将召回的用户列表都到map中判断一下是否存在，存在就被过滤，golang的map查询时间复杂度是O(1)，遍历用户列表时间复杂度是O(n)，整个过滤的时间复杂度是O(n)

7：特征计算

过滤完后得到一个用户列表，根据每个用户的特征信息，每个特征一个分值，然后按照一个公式，计算得到一个最终的分值。就像学生每门课成绩有一个分值，这里的公式可能不是简单的相加，而是按照业务随时调整的公式

8：排序

特征计算已经得到一个分值，按照这个分值排序，取top30推荐给用户，一个简单的推荐系统就完成了。