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

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

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, поскольку a[[i]] в момент оценки не является списком.

Простое решение - просто игнорировать это сообщение или отключить его. Конечно, есть и другие способы его обойти. Например, можно обойти оценку a[[i]], заменив его внутри Compile-body перед компиляцией

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

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

Да, и это связано с причиной, по которой я не могу использовать здесь :=. Я использую 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 работает в этом случае, но оно оставляет все выражение без оценки. Позвольте мне объяснить, почему важно, чтобы выражение было расширено перед компиляцией. Должен признать, что мой игрушечный пример был не самым удачным. Так что давайте попробуем сделать лучше, используя 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 и мне нужна его производная внутри Compile. В приведенном выше коде я хочу, чтобы система Mathematica вычисляла производную для неприписанных корней a[[i]], чтобы я мог часто использовать формулу very в скомпилированном коде. Если посмотреть на скомпилированный код выше видно, что 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]
  ]
 ]

Edit: I put the Quiet to the place it belongs. Я поставил его перед блоком кода, чтобы всем было видно, что я использовал Quiet здесь для подавления предупреждения. Как уже отметил Рубенко, в реальном коде он всегда должен быть как можно ближе к месту, где ему место. При таком подходе вы, вероятно, не пропустите другие важные предупреждения/ошибки.

Обновление 2

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

Обновление 3

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

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

Тизер:

Dynamic Newton fractal

7
задан halirutan 6 January 2012 в 10:35
поделиться