wolfram-mathematica - a[[i]] 的未求值形式

Question

考虑以下简单的示例

cf = Block[{a, x, degree = 3},
  With[{expr = Product[x - a[[i]], {i, degree}]},
   Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]
   ]
  ]

这是在编译语句的主体中传输代码的一种可能方式。它会产生 Part::partd 错误，因为 a[[i]] 在评估时不是列表。

简单的解决方案是忽略此消息或将其关闭。当然还有其他方法。例如，可以通过在编译之前在 Compile-body 中替换它来规避对 a[[i]] 的评估

cf = ReleaseHold[Block[{a, x, degree = 3},
   With[{expr = Product[x - a[i], {i, degree}]},
    Hold[Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]] /. 
     a[i_] :> a[[i]]]
   ]
  ]

如果编译后的函数代码很多，最后的 Hold、Release 和替换有点违背我对漂亮代码的想法。有没有我还没有考虑过的简短而好的解决方案？

回复 Szabolcs 的帖子

Could you tell me though why you are using With here?

是的，这与我不能在这里使用 := 的原因有关。我使用 With，得到类似 #define 的东西在 C 中，这意味着在我需要的地方进行代码替换。在 With 中使用 := 会延迟评估，并且 Compile 的主体所看到的并不是它应该编译的最后一段代码。所以，

<< CompiledFunctionTools`
cf = Block[{a, x, degree = 3}, 
   With[{expr := Product[x - a[[i]], {i, degree}]}, 
    Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]]];
CompilePrint[cf]

向您展示，在编译函数中有对 Mathematica 内核的调用

I4 = MainEvaluate[ Function[{x, a}, degree][ R0, T(R1)0]]

这很糟糕，因为 Compile 应该只使用局部变量来计算结果。

更新

Szabolcs 解决方案在这种情况下有效，但它使整个表达式未求值。让我解释一下，为什么在编译之前扩展表达式很重要。我不得不承认，我的玩具示例并不是最好的。因此，让我们在 Szabolcs 的解决方案中尝试使用 With 和 SetDelayed 更好的方法

Block[{a, x}, With[
  {expr := D[Product[x - a[[i]], {i, 3}], x]}, 
  Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]
  ]
 ]

假设我有一个 3 次多项式，我需要它在 Compile 中的导数。在上面的代码中，我希望 Mathematica 计算未分配根 a[[i]] 的导数，因此我可以在编译代码中经常使用该公式。看上面的编译代码
可以看到，D[..] 不能像 Product 一样被编译，并且没有被评估

11  R1 = MainEvaluate[ Hold[D][ R5, R0]]

因此，我更新的问题是:是否可以在不评估其中的 Part[] 访问的情况下评估一段代码比使用更好/更好

Block[{a, x}, With[
  {expr = D[Quiet@Product[x - a[[i]], {i, 3}], x]}, 
  Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]
  ]
 ]

编辑:我把安静放在它所属的地方。我把它放在代码块前面，让每个人都可以看到我在这里使用 Quiet 来抑制警告。正如 Ruebenko 已经指出的那样，在实际代码中，它应该始终尽可能接近它所属的位置。使用这种方法，您可能不会错过其他重要的警告/错误。

更新 2

由于我们正在远离最初的问题，我们应该将这个讨论转移到一个新的线程。我不知道我应该给谁我的问题的最佳答案奖，因为我们讨论了 Mathematica 和范围界定，而不是如何抑制 a[[i]] 问题。

更新 3

给出最终的解决方案:我只是用 Quiet 压制(不幸的是，就像我一直做的那样)a[[i]] 警告。在下面的一个真实示例中，我必须在整个 Block 之外使用 Quiet 来抑制警告。

为了将所需的代码注入(inject)到 Compile 的主体中，我使用了一个纯函数并将代码作为参数内联。这与 Michael Trott 使用的方法相同，例如他的数字书。这有点像 where Haskell 中的子句，您可以在其中定义之后使用的内容。

newtonC = Function[{degree, f, df, colors},
   Compile[{{x0, _Complex, 0}, {a, _Complex, 1}},
    Block[{x = x0, xn = 0.0 + 0.0 I, i = 0, maxiter = 256, 
      eps = 10^(-6.), zeroId = 1, j = 1},
     For[i = 0, i < maxiter, ++i,
      xn = x - f/(df + eps);
      If[Abs[xn - x] < eps,
       Break[]
       ];
      x = xn;
      ];
     For[j = 1, j <= degree, ++j,
      If[Abs[xn - a[[j]]] < eps*10^2,
        zeroId = j + 1;
        Break[];
        ];
      ];
     colors[[zeroId]]*(1 - (i/maxiter)^0.3)*1.5
     ],
    CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
    RuntimeOptions -> "Speed", Parallelization -> True]]@@
    (Quiet@Block[{degree = 3, polynomial, a, x},
     polynomial = HornerForm[Product[x - a[[i]], {i, degree}]];
     {degree, polynomial, HornerForm[D[polynomial, x]], 
      List @@@ (ColorData[52, #] & /@ Range[degree + 1])}])

这个函数现在足够快，可以计算根位置不固定的多项式的牛顿分形。因此，我们可以动态调整根。
随意调整 n。在这里它流畅地运行到 n=756

(* ImageSize n*n, Complex plange from -b-I*b to b+I*b *)
With[{n = 256, b = 2.0},
 DynamicModule[{
   roots = RandomReal[{-b, b}, {3, 2}],
   raster = Table[x + I y, {y, -b, b, 2 b/n}, {x, -b, b, 2 b/n}]},
  LocatorPane[Dynamic[roots],
   Dynamic[
    Graphics[{Inset[
       Image[Reverse@newtonC[raster, Complex @@@ roots], "Real"],
       {-b, -b}, {1, 1}, 2 {b, b}]}, PlotRange -> {{-b, b}, {-
         b, b}}, ImageSize -> {n, n}]], {{-b, -b}, {b, b}}, 
   Appearance -> Style["\[Times]", Red, 20]
   ]
  ]
 ]

预告片:

Answer 1

好的，这是我用于各种目的的代码生成框架的非常简化的版本:

ClearAll[symbolToHideQ]
SetAttributes[symbolToHideQ, HoldFirst];
symbolToHideQ[s_Symbol, expandedSymbs_] :=! MemberQ[expandedSymbs, Unevaluated[s]];

ClearAll[globalProperties]
globalProperties[] := {DownValues, SubValues, UpValues (*,OwnValues*)};

ClearAll[getSymbolsToHide];
Options[getSymbolsToHide] = {
     Exceptions -> {List, Hold, HoldComplete, 
        HoldForm, HoldPattern, Blank, BlankSequence, BlankNullSequence, 
       Optional, Repeated, Verbatim, Pattern, RuleDelayed, Rule, True, 
       False, Integer, Real, Complex, Alternatives, String, 
       PatternTest,(*Note-  this one is dangerous since it opens a hole 
                    to evaluation leaks. But too good to be ingored *)
       Condition, PatternSequence, Except
      }
 };

getSymbolsToHide[code_Hold, headsToExpand : {___Symbol}, opts : OptionsPattern[]] :=
  Join @@ Complement[
       Cases[{
          Flatten[Outer[Compose, globalProperties[], headsToExpand]], code},
            s_Symbol /; symbolToHideQ[s, headsToExpand] :> Hold[s],
            Infinity,
            Heads -> True
       ],
       Hold /@ OptionValue[Exceptions]];

ClearAll[makeHidingSymbol]
SetAttributes[makeHidingSymbol, HoldAll];
makeHidingSymbol[s_Symbol] := 
    Unique[hidingSymbol(*Unevaluated[s]*) (*,Attributes[s]*)];

ClearAll[makeHidingRules]
makeHidingRules[symbs : Hold[__Symbol]] :=
     Thread[List @@ Map[HoldPattern, symbs] -> List @@ Map[makeHidingSymbol, symbs]];

ClearAll[reverseHidingRules];
reverseHidingRules[rules : {(_Rule | _RuleDelayed) ..}] :=
   rules /. (Rule | RuleDelayed)[Verbatim[HoldPattern][lhs_], rhs_] :> (rhs :> lhs);


FrozenCodeEval[code_Hold, headsToEvaluate : {___Symbol}] :=   
   Module[{symbolsToHide, hidingRules, revHidingRules,  result}, 
      symbolsToHide = getSymbolsToHide[code, headsToEvaluate];
      hidingRules = makeHidingRules[symbolsToHide];
      revHidingRules = reverseHidingRules[hidingRules];
      result = 
         Hold[Evaluate[ReleaseHold[code /. hidingRules]]] /. revHidingRules;
      Apply[Remove, revHidingRules[[All, 1]]];
      result];

该代码通过使用一些虚拟符号暂时隐藏大多数符号来工作，并允许评估某些符号。这是在这里的工作方式:

In[80]:= 
FrozenCodeEval[
  Hold[Compile[{{x,_Real,0},{a,_Real,1}},D[Product[x-a[[i]],{i,3}],x]]],
  {D,Product,Derivative,Plus}
]

Out[80]= 
Hold[Compile[{{x,_Real,0},{a,_Real,1}},
  (x-a[[1]]) (x-a[[2]])+(x-a[[1]]) (x-a[[3]])+(x-a[[2]]) (x-a[[3]])]]

因此，要使用它，您必须将代码包装在 Hold 中。并指出您要评估的磁头。剩下的就是申请 ReleseHold给它。请注意，上面的代码只是说明了这些想法，但仍然非常有限。我的方法的完整版本涉及其他步骤，使其功能更强大但也更复杂。

编辑

虽然上面的代码仍然太有限，无法容纳许多真正有趣的案例，但这里有一个额外的简洁示例，说明 rather hard to achieve使用传统的评估控制工具:

In[102]:= 
FrozenCodeEval[
  Hold[f[x_, y_, z_] := 
    With[Thread[{a, b, c} = Map[Sqrt, {x, y, z}]], 
       a + b + c]], 
  {Thread, Map}]

Out[102]= 
Hold[
  f[x_, y_, z_] := 
    With[{a = Sqrt[x], b = Sqrt[y], c = Sqrt[z]}, a + b + c]]