将随时间变化的加权平均向量化
我在尝试随着时间推移计算一个储罐的质量属性。质量是上一周期质量的加权平均值（按期末存货余额加权）与当前周期进入的材料及其质量的加权平均值。我已经推导出一个慢速的顺序算法，但在向量化方面遇到困难。我以为累积和可能有帮助，但在尝试中并未取得明显成效。
储罐质量属性的公式如下：
q_{n} = (q_{n-1} * v_{n-1} + qinc_{n} * vinc_{n}) / (v_{n-1} + vinc_{n})
q_{n}..... quality at time period n
v_{n-1}... carry over balance (volume) from previous period
vinc_{n}.. incoming new material quantity
qinc_{n}.. incoming new material quality
我的慢速顺序算法：
q = np.empty_like(qinc)
q.fill(np.nan)
if v[0] > 0:
q[0] = qinc[0]
for i in range(1, len(qinc)):
v_prev = v[i-1]
q_prev = np.nan
if v_prev > 0:
q_prev = q[i-1]
q_curr = qinc[i]
v_curr = vinc[i]
if v[i] > 0:
q[i] = (np.nan_to_num(q_prev * v_prev) + np.nan_to_num(q_curr * v_curr)) / (v_prev + v_curr)
一个示例
qinc = np.asarray([0.5, 2.0, np.nan, 4, 8])
vinc = np.asarray([100, 200, 0, 300, 50])
v = np.asarray([0, 10, 100, 400, 200])
结果为
q == [nan, 2., 2., 3.5, 4.]
^ ^ ^ ^ ^
| | | | |
| | | | +- ==((q[3] * v[3]) + (qinc[4] * vinc[4])) / (v[3] + vinc[4]) => previous quality weighted by closing stock and this period's incoming quality weighted by the amount of material coming in
| | | +- ((2 * 100) + (4 * 300)) / (100 + 300)
| | +- no new material vinc[2] == 0 (use carry over quality from q[1])
| +- incoming quality is qinc[1] == 2 (no carry over from q[0])
+- closing stock (v[0]) == 0 -> quality value is nan
任何线索将不胜感激。
解决方案
现实世界建模
我理解递推关系， 它当然可以有意义地被利用。 但我并不明白为什么我们一定要使用它。
当你写到“quality”时，我更倾向写成“solite”（溶质）。 你把数量表示为分数（或百分比）， 而我更愿意把所有量用升来表达。
一升可销售的“产品”由一定体积的“溶剂”和另一部分分数升的“溶质”组成。 把乙醇溶解在水中的情景很方便来理解， 或者在乳品场景里，溶解的脂肪在奶业背景下也类似。 我们假设一个封闭系统，零蒸发损失，在恒温条件下， 因此质量和体积的守恒都会成立。 重申以上：溶质+ 溶剂的体积将给出总的产品体积。 在这样的加法关系下，我们只需要跟踪其中三个变量中的两个； 用你的数据来看，我们可能是在跟踪乙醇和啤酒的升数。 （引入换算因子将使我们能跟踪乙醇和啤酒的克数，如果这确实是你关注的重点。）
采用这种表示法自然会引导我们去找出
进入储罐的所有物质的cumsum()。
请注意，这立刻会对
示例中的qinc向量里的NaN提出质疑。
是否有某个“未定义”的溶质量流入？
没有！
恰好流入了0.0 升的溶质。
有限的罐体容量
显然，部分“流入”量应为负值。 否则，系统运行一段时间后会溢出。
按比例提取
原始问题的设定显然假定“完美混合”。 取出一升产品（相当于加上 -1升产品）会使我们在提取时依据当前浓度提取出一定量的溶质。
因此，我们把“需要知道前一步的数字”这个恼人的细节
从“往罐里灌”阶段移到“排空罐子”阶段，推迟处理。
我们需要在提取时知道浓度，而这会干扰使用 cumsum()。
除非 我们选择把取出的量表示为当前产品量的分数。 这样我们就可以在下一步通过加法（添加的量）再乘以（排出的分数）来得到取出的量。 在最简单的情况下，我们每天进行批量处理，排空罐子100%；任何其他分数都同样可行。
我们也许可以通过假设一个容量无限的罐子来“精简”这一点（没有溢出），并将所有正量在一个总账中记录，负量（取出）记录在另一个总账中。
如果我们经常进行每日多次添加（也许来自每头奶牛一次），并且取出次数较少（每晚送货车到达），那么这可能是一个净收益。
备选方案
递归计算并不总是能仅靠Python向量化。但你可以在这里使用 for 循环配合 numba。for 循环在numba中的性能与普通的Python for循环不同。使用 numba，for循环的处理速度与C 语言相同。只需在函数前添加一个 @njit 装饰器。
from numba import njit
import numpy as np
@njit
def calculate_quality_fast(qinc, vinc, v):
q = np.empty_like(qinc)
q.fill(np.nan)
if v[0] > 0:
q[0] = qinc[0]
for i in range(1, len(qinc)):
v_prev = v[i-1]
q_prev = q[i-1] if v_prev > 0 else 0.0
if v[i] > 0:
# Replaced nan_to_num with simple inline checks for speed
q_curr = qinc[i] if not np.isnan(qinc[i]) else 0.0
v_curr = vinc[i]
q[i] = (q_prev * v_prev + q_curr * v_curr) / (v_prev + v_curr)
return q
@njit 表示 no_python，意思是你不应该在某些 numba 无法处理的地方使用 - 仅在你提到的那个特定函数中才允许使用 @njit。 在其他地方你可以使用。Python与 numba 之间的切换会自动完成。（例如，numba 不知道 dict。你不能在其中使用混合数据类型的list。） 请参阅numba官方文档中关于在 numba 中可以/不能使用的Python功能列表，点此查看
相关文章
