代码之家  ›  专栏  ›  技术社区  ›  Samsdram

Mathematica:使用Simplify进行常用子表达式消除和强度降低

  •  8
  • Samsdram  · 技术社区  · 15 年前

    所以最近我一直在玩弄Mathematica的模式匹配和术语重写如何在编译器优化中得到很好的应用…尝试高度优化循环内部的短代码块。减少计算表达式所需工作量的两种常见方法是标识出现多次的子表达式并存储结果,然后在随后的点使用存储的结果来保存工作。另一种方法是尽可能使用更便宜的操作。例如,我的理解是取平方根比加和乘需要更多的时钟周期。很明显,我对浮点运算的成本感兴趣,因为计算表达式需要花费多少时间,而不是mathematica计算表达式需要多长时间。

    我的第一个想法是用Mathematica解决这个问题 simplify 功能。可以指定一个比较两个表达式相对简单性的复杂度函数。我准备为相关的算术运算创建一个权重,并在此基础上添加表达式的leafcount,以说明所需的赋值运算。这解决了强度方面的减少,但正是消除了常见的子表达式,我才绊倒了。

    我正在考虑将公共子表达式消除添加到可能简化使用的转换函数中。但是对于一个大型表达式,可能会有许多可能被替换的子表达式,在看到表达式之前,不可能知道它们是什么。我已经编写了一个函数,它给出了可能的替换,但您指定的转换函数似乎只需要返回一个可能的转换,至少从文档中的示例中是这样的。有什么想法可以绕过这个限制吗?是否有人对简化如何使用可能暗示未来方向的转换函数有更好的了解?

    我想,在简化的幕后,正在做一些动态编程,尝试对表达式的不同部分进行不同的简化,并返回复杂度得分最低的部分。我是否最好自己使用常见的代数简化(如factor和collect)来进行这种动态编程?

    编辑:我添加了生成可能子表达式的代码来删除

    (*traverses entire expression tree storing each node*)
    AllSubExpressions[x_, accum_] := Module[{result, i, len},
      len = Length[x];
      result = Append[accum, x];
      If[LeafCount[x] > 1,
       For[i = 1, i <= len, i++,
         result = ToSubExpressions2[x[[i]], result];
         ];
       ];
      Return[Sort[result, LeafCount[#1] > LeafCount[#2] &]]
      ]
    
    CommonSubExpressions[statements_] := Module[{common, subexpressions},
    subexpressions = AllSubExpressions[statements, {}];
      (*get the unique set of sub expressions*)
      common = DeleteDuplicates[subexpressions];
      (*remove constants from the list*)
      common = Select[common, LeafCount[#] > 1 &];
      (*only keep subexpressions that occur more than once*)
      common = Select[common, Count[subexpressions, #] > 1 &];
      (*output the list of possible subexpressions to replace with the \
    number of occurrences*)
      Return[common];
      ]
    

    一旦从commonSubExpressions返回的列表中选择了一个公共子表达式,则执行替换的函数将在下面。

    eliminateCSE[statements_, expr_] := Module[{temp},
    temp = Unique["r"];
    Prepend[ReplaceAll[statements, expr -> temp], temp[expr]]
    ]
    

    冒着这个问题变长的风险,我将编写一个小示例代码。我认为一个合适的表达来优化将是经典的 Runge-Kutta 解微分方程的方法。

    Input:
    nextY=statements[y + 1/6 h (f[t, n] + 2 f[0.5 h + t, y + 0.5 h f[t, n]] + 
        2 f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]] + 
        f[h + t, 
         y + h f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]]])];
    possibleTransformations=CommonSubExpressions[nextY]
    transformed=eliminateCSE[nextY, First[possibleTransformations]]
    
    Output:
    {f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]], 
    y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]], 
    0.5 h f[0.5 h + t, y + 0.5 h f[t, n]], 
    f[0.5 h + t, y + 0.5 h f[t, n]], y + 0.5 h f[t, n], 0.5 h f[t, n], 
    0.5 h + t, f[t, n], 0.5 h}
    
    statements[r1[f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]]], 
    y + 1/6 h (2 r1 + f[t, n] + 2 f[0.5 h + t, y + 0.5 h f[t, n]] + 
     f[h + t, h r1 + y])]
    

    最后,给出了判断不同表达式相对成本的代码。在这一点上,权重是概念上的,因为这仍然是我正在研究的领域。

    Input:
    cost[e_] := 
     Total[MapThread[
       Count[e, #1, Infinity, Heads -> True]*#2 &, {{Plus, Times, Sqrt, 
         f}, {1, 2, 5, 10}}]]
    cost[transformed]
    
    Output:
    100
    
    2 回复  |  直到 15 年前
        1
  •  4
  •   Yaroslav Bulatov    15 年前

    要标识重复的子表达式,可以使用类似这样的内容

    (*helper functions to add Dictionary-like functionality*)
    
    index[downvalue_, 
       dict_] := (downvalue[[1]] /. HoldPattern[dict[x_]] -> x) // 
       ReleaseHold;
    value[downvalue_] := downvalue[[-1]];
    indices[dict_] := 
      Map[#[[1]] /. {HoldPattern[dict[x_]] -> x} &, DownValues[dict]] // 
       ReleaseHold;
    values[dict_] := Map[#[[-1]] &, DownValues[dict]];
    items[dict_] := Map[{index[#, dict], value[#]} &, DownValues[dict]];
    indexQ[dict_, index_] := 
      If[MatchQ[dict[index], HoldPattern[dict[index]]], False, True];
    
    
    (*count number of times each sub-expressions occurs *)
    expr = Cos[x + Cos[Cos[x] + Sin[x]]] + Cos[Cos[x] + Sin[x]];
    Map[(counts[#] = If[indexQ[counts, #], counts[#] + 1, 1]; #) &, expr, 
      Infinity];
    items[counts] // Column
    
        2
  •  5
  •   eacousineau    13 年前

    作者在这里还实现了一些例程: http://stoney.sb.org/wordpress/2009/06/converting-symbolic-mathematica-expressions-to-c-code/

    我把它打包成一个*.m文件,并修复了一个bug(如果表达式没有重复的子表达式,它就会死掉),我试图找到作者的联系信息,看看是否可以将修改后的代码上传到Pastebin或任何地方。

    编辑:我已收到作者的上载权限,并已将其粘贴到此处: http://pastebin.com/fjYiR0B3