Home

今天看Gilbert Strang的公开课，他提到了矩阵的应用中说了化学，
这种方法在我之前写kinetic模型的时候就用了，核心代码很少，就能把一连串基元反应配平并合并了哈哈哈。
来来来，上code(主要代码)：

1
2
3
4
5
6
7
8
9

def null(A, eps=1e-10):
    "get null space of transposition of site_matrix"
     u, s, vh = np.linalg.svd(A, full_matrices=1, compute_uv=1)
     null_space = np.compress(s <= eps, vh, axis=0)
     return null_space.T
 x = null(site_matrix.T)  # basis of null space
 if not x.any():  # x is not empty
     raise ValueError('Failed to get basis of nullspace.')
 x = map(abs, x.T.tolist()[0])

之前的动力学代码的数据文件我是在catmap的基础上进行改进，写了个tablemaker来生成一个csv文件，然后让要计算动力学的人把每个需要的物种的绝对能量输入到表格中，模型会在处理数据的时候自动计算其所有基元反应的能垒以及反应能量变化。之所有这样进行能量输入，是方便在后面能够基于每个物种的能量进行计算，比如在做出是猜测的时候进行Boltzmann分布，在进行敏感度分析时候能够针对每个中间态和过渡态进行分析。
但是之前有一次帮师姐计算它的反应，他只有相对能量。如下，

决定空闲的时候把第五版《C-Primer-Plus》的编程练习敲一边来复习下C，在GitHub(https://github.com/PytLab/C-Primer-Plus)上同步更新。

Table of Content

• 第三章 数据和C
• 第四章 字符串和格式化输入/输出
• 第五章 运算符、表达式和语句
• 第六章 C控制语句: 循环
• 第七章 C控制语句: 分支和跳转
• 第八章 字符输入/输出和输入确认
• 第九章 函数

在用sympy符号运算将动力学模型重新实现一遍以后，出了能进行latex输出等其他好处外，最初使用符号运算的目的是符号求解敏感度 $X_{TRC}$的值，不过由于完全用solve暴力解析求解稳态覆盖度，稍微复杂点的机理就会计算很久。于是我打算用平衡态近似来进行求解，并且写完平衡态solver后不仅仅是求$X_{TRC}$，单独用QuasiEquilibriumSolver求解动力学模型也是可以的。

为了避免直接用sympy.solve()函数直接求解，分析了平衡态近似的求解方法后，自己用python先对求解过程进行了一下”预处理”，这样从一定程度上减小了对solve()依赖。由于这部分还是有点麻烦的，所以写在这里记录下来，方便自己以及其他人学习和改进。

由于平衡态近似求解的方法还是比较单一的，大概过程就是在决定Rate Determining Step后，通过联立其他基元反应的平衡条件，也就是$\prod_{j} \theta_{i,j}^{-} \prod_{j} p_{i,j}^{-} =K_{i}\prod_{j} \theta_{i,j}^{+}\prod_{j} p_{i,j}^{+}$用$\theta^{*}$将其他的$\theta_{i}$表示出来，在回代到归一化条件中$\sum_{i}{\theta_{i}} = 1$,将$\theta^{*}$的表达式求出来，然后分别回代求出不同中间态的覆盖度。
计算方法很简单，但是这里涉及到对sympy.Symbol对象的引用和substitution以及后面基元反应的特殊情况，所以写起来略复杂。下面以甲酸裂解的一条路径的基元反应为例说明一下。

这两天想学习下优化代码，既然手底下有自己写完的动力学代码，那正好可以拿着练手。按我现在的认识程度，能对现在写好的动力学代码优化提升效率的方法除了改进算法外，单纯在编程方面一个是使用多线程并行，一个是用c api把循环部分重写。
把python多线程的一点点皮毛看了看就现学现卖了下。
把迭代的部分看了看，每步迭代里面都不是各自独立不相关的运算，真心不能给每个部分的运算分配子线程。在整个动力学模型里面最适合用多线程的地方就是作图了，因为作图中的阴影以及循环添加注释是可以分别分配子线程来进行的。于是我就把plotter里面画图的部分改了下:
分别添加了两个threading.Thread的子类ShadowThread，NoteThread用来实例化阴影和添加注释的子进程。

在这里把C中有关内存管理的部分稍微总结下，方便以后回头浏览复习。