Неоцененная форма [[i]] - PullRequest
       36

Неоцененная форма [[i]]

7 голосов
/ 05 января 2012

Рассмотрим следующий простой, иллюстрирующий пример

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

Это один из возможных способов передачи кода в теле оператора Compile. Он выдает ошибку Part :: partd, так как [[i]] на момент оценки не является списком.

Самое простое решение - просто проигнорировать это сообщение или выключить его. Конечно, есть и другие способы. Например, можно обойти оценку [[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 и замена в конце идут вразрез с моим представлением о прекрасном коде. Есть ли короткое и приятное решение, которое я еще не рассмотрел?

Ответ на пост сабольца

Не могли бы вы рассказать мне, почему вы используете здесь?

Да, и это связано с причиной, по которой я не могу использовать: = здесь. Я использую With, чтобы иметь что-то вроде #define в C, что означает замену кода в том месте, где оно мне нужно. Использование: = in 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 работает, но оно оставляет все выражение без оценки. Позвольте мне объяснить, почему важно, чтобы выражение было расширено до его компиляции. Должен признать, мой игрушечный пример был не лучшим. Итак, давайте попробуем лучше, используя With и SetDelayed, как в решении Szabolcs

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

Скажем, у меня есть многочлен степени 3, и мне нужна его производная внутри компиляции. В приведенном выше коде я хочу, чтобы Mathematica вычисляла производную для неназначенных корней a [[i]], чтобы я могла использовать формулу очень часто в скомпилированном коде. Глядя на скомпилированный код выше каждый видит, что D [..] не может быть скомпилирован так же хорошо, как Продукт, и остается неоцененным

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 здесь, чтобы подавить предупреждение. Как уже указывал Рубенко, в реальном коде он всегда должен быть как можно ближе к месту, к которому он принадлежит. При таком подходе вы, вероятно, не пропустите другие важные предупреждения / ошибки.

Обновление 2

Поскольку мы отходим от первоначального вопроса, нам следует перенести это обсуждение, возможно, в новую ветку. Я не знаю, кому мне следует дать лучший ответ на мой вопрос, поскольку мы обсудили Mathematica и Scoping больше, чем способы подавления проблемы [[i]].

Обновление 3

Чтобы дать окончательное решение: я просто подавляю (к сожалению, как делал все время) предупреждение [[i]] с помощью Quiet. В настоящем примере ниже я должен использовать Quiet за пределами полного блока, чтобы подавить предупреждение.

Чтобы вставить требуемый код в тело компиляции, я использую чистую функцию и даю код в качестве аргумента. Это тот же подход, который использует Майкл Тротт, например его книга чисел. Это немного похоже на предложение 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]
   ]
  ]
 ]

Teaser:

Dynamic Newton fractal

Ответы [ 4 ]

10 голосов
/ 05 января 2012

Хорошо, вот очень упрощенная версия инфраструктуры генерации кода, которую я использую для различных целей:

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 к нему. Обратите внимание, что приведенный выше код просто иллюстрирует идеи, но все еще довольно ограничен. Полная версия моего метода включает в себя другие шаги, которые делают его намного более мощным, но и более сложным.

Редактировать

Хотя приведенный выше код все еще слишком ограничен, чтобы вместить много действительно интересных случаев, вот еще один изящный пример того, что было бы довольно трудно достичь с помощью традиционных инструментов контроля оценки:

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]]
4 голосов
/ 05 января 2012

РЕДАКТИРОВАТЬ - Большое предупреждение !! Внедрение кода, использующего With или Function в Compile, который использует некоторые из локальных переменных Compile, не является надежным ! Учтите следующее:

In[63]:= With[{y=x},Compile[x,y]]
Out[63]= CompiledFunction[{x$},x,-CompiledCode-]

In[64]:= With[{y=x},Compile[{{x,_Real}},y]]
Out[64]= CompiledFunction[{x},x,-CompiledCode-]

Обратите внимание на переименование x в x$ в первом случае. Я рекомендую вам прочитать о локализации здесь и здесь . (Да, это сбивает с толку!) Мы можем догадаться, почему это происходит только в первом случае, а не во втором, но я хочу сказать, что такое поведение может быть непреднамеренным (назовите это ошибкой, темным углом или неопределенным поведением), поэтому полагаться на это хрупко ...

Решения на основе

Replace, такие как моя withRules функция , работают, хотя (это не было моим предназначением для этой функции, но она хорошо подходит здесь ...)

In[65]:= withRules[{y->x},Compile[x,y]]
Out[65]= CompiledFunction[{x},x,-CompiledCode-]

Оригинальные ответы

Вы можете использовать := в With, вот так:

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

Это позволит избежать оценки expr и ошибки от Part.

Обычно = и := работают, как и ожидалось, во всех With, Module и Block.


Не могли бы вы сказать мне, почему вы используете With здесь? (Уверен, у вас есть веская причина, я просто не вижу этого в этом упрощенном примере.)


Дополнительный ответ

Устранение проблемы @ halirutan по поводу того, что degree не указывается во время компиляции

Я вижу, что это точно такая же ситуация, как если бы у нас была определена глобальная переменная, которую мы используем в Compile. Взять, к примеру:

In[18]:= global=1
Out[18]= 1

In[19]:= cf2=Compile[{},1+global]
Out[19]= CompiledFunction[{},1+global,-CompiledCode-]

In[20]:= CompilePrint[cf2]
Out[20]= 
        No argument
        3 Integer registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = 1
        Result = I2

1   I1 = MainEvaluate[ Function[{}, global][ ]]
2   I2 = I0 + I1
3   Return

Это общая проблема. Решение состоит в том, чтобы указать Compile встроенным глобальным переменным , например:

cf = Block[{a, x, degree = 3}, 
   With[{expr := Product[x - a[[i]], {i, degree}]}, 
    Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr, 
     CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];

CompilePrint[cf]

Вы можете проверить, что теперь нет обратного вызова для главного оценщика.


В качестве альтернативы вы можете ввести значение degree, используя дополнительный слой With вместо Block. Это заставит вас очень сильно захотеть что-то вроде этого .


Расширение макроса в Mathematica

Это несколько не связано, но вы упоминаете в своем посте, что вы используете With для расширения макроса. Вот мой первый (возможно, ошибочный) подход к реализации расширения макросов в Mathematica. Это не очень хорошо проверено, не стесняйтесь попробовать сломать его и оставить комментарий.

Clear[defineMacro, macros, expandMacros]

macros = Hold[];

SetAttributes[defineMacro, HoldAllComplete]
defineMacro[name_Symbol, value_] := (AppendTo[macros, name]; name := value)

SetAttributes[expandMacros, HoldAllComplete]
expandMacros[expr_] := Unevaluated[expr] //. Join @@ (OwnValues /@ macros)

Описание:

macros - это (удерживаемый) список всех символов, которые должны быть расширены. defineMacro создаст новый макрос. expandMacros развернет определения макросов в выражении.

Осторожно: я не реализовал переопределение макросов, это не будет работать, если расширение используется $Pre. Также остерегайтесь рекурсивных определений макросов и бесконечных циклов.

Использование:

Расширение макроса для всех входных данных путем определения $Pre:

$Pre = expandMacros;

Определите a, чтобы иметь значение 1:

defineMacro[a, 1]

Установить отложенное определение для b:

b := a + 1

Обратите внимание, что определение b не полностью оценено, но a расширено.

?b

Global`b

b:=1+1

Отключить расширение макроса ($Pre может быть опасно, если в моем коде есть ошибка):

$Pre =.
2 голосов
/ 05 января 2012

В одну сторону:

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

Будьте осторожны, если вы действительно этого хотите.

0 голосов
/ 08 февраля 2017

Оригинальный код:

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, 
...
RuntimeOptions -> "Speed", Parallelization -> True]]@@
(Quiet@Block[{degree = 3, polynomial, a, x},
 polynomial = HornerForm[Product[x - a[[i]], {i, degree}]];
...

Модифицированный код:

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, 
...
RuntimeOptions -> "Speed", Parallelization -> True],HoldAllComplete]@@
( (( (HoldComplete@@#)/.a[i_]:>a[[i]] )&)@Block[{degree = 3, polynomial, a, x},
 polynomial = HornerForm[Product[x - a[i], {i, degree}]];
...

Добавить атрибут HoldAllComplete к функции.

Напишите a[i] вместо a[[i]].

Заменить Quiet на (( (HoldComplete@@#)/.a[i_]:>a[[i]] )&)

Создает идентичный код, без Quiet, и все вещи Hold находятся в одном месте.

...