Mathematica：使用簡化去做常見的子表達式消除和減少力度

Question

最近我一直在玩弄Mathematica的模式匹配和術語重寫如何在編譯器優化中得到很好的利用......試圖高度優化短期代碼塊是循環的內部部分。減少評估表達式所花費的工作量的兩種常用方法是識別多次出現的子表達式並存儲結果，然後在隨後的點上使用存儲的結果來保存工作。另一種方法是儘可能使用更便宜的操作。例如，我的理解是，採取平方根需要比加法和乘法更多的時鐘週期。清楚的是，我對評估表達式所需的浮點運算成本感興趣，而不是Mathematica評估它需要多長時間。

我的第一個想法是，我會解決使用Mathematica的simplify函數開發的問題。可以指定一個複雜度函數來比較兩個表達式的相對簡單性。我打算爲相關的算術運算使用權重創建一個，並添加到表達式的LeafCount中以說明所需的賦值操作。這解決了強度方面的減少，但它是消除了我絆倒了常見的子表達式。

我在考慮將常見的子表達式消除添加到簡化使用的可能轉換函數中。但是對於一個大型的表達，可能會有許多可能的子表達式被替代，直到看到表達式纔會知道它們是什麼。我已經寫了一個函數來給出可能的替換，但是您指定的轉換函數似乎需要返回一個可能的轉換，至少從文檔中的示例中返回。關於如何解決這個限制的任何想法？有沒有人對如何簡化使用可能暗示前進方向的轉換函數有更好的想法？

我想在Simplify的幕後做一些動態編程，嘗試對錶達式的不同部分進行不同的簡化並返回複雜度最低的那個。使用常用的代數簡化（例如因子和收集），我是否會更好地嘗試自己動態編程？

編輯：我添加產生可能的子表達式的代碼一旦公共子表達式從由CommonSubExpressions返回，做更換低於功能列表選擇刪除

(*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]; 
    ] 

。

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

在這個問題變得漫長的危險中，我會舉一個小例子代碼。我認爲一個體面的表達來嘗試優化將是經典的求解微分方程的方法。

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 

Answer 1

要找出重複的子表達式，你可以使用類似這樣

(*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 

Answer 2

也有一些在這裏通過例行筆者在這裏實現：http://stoney.sb.org/wordpress/2009/06/converting-symbolic-mathematica-expressions-to-c-code/

我打包成一個* .M文件，並修復了一個錯誤（如果表達式沒有重複的子表達式，它會死掉），我試圖找到作者的聯繫信息，看看我能否將他修改的代碼上傳到pastebin或任何地方。

編輯：我已經收到許可，筆者把它上傳，並粘貼在這裏吧：http://pastebin.com/fjYiR0B3