# error when defining a compiled function with ReplacePart

I always have difficulty in writing a non-trivial compiled function. I'm using Mathematica 9. Please see the following code

Clear[longtermCompiled];
longtermCompiled =
Compile[{{M, _Real, 2}},
Module[{len = Length[M]},
ReplacePart[Transpose[M] - IdentityMatrix[len], {len, i_} -> Table[1, {len}]]]];


where M is expected to be a matrix (two-dimensional). Upon executing the above definition, it complained:

ReplacePart::argrx: ReplacePart called with 2 arguments; 3 arguments are expected.


I don't understand this error message, because the following code works without any problems:

Clear[longtermX];
longtermX[M_] :=
Module[{len = Length[M]},
ReplacePart[Transpose[M] - IdentityMatrix[len], {len, i_} -> Table[1, {len}]]];


What's wrong?

-
That's really strange. ReplacePart does not look like it accepts 3 arguments, so I'm wondering if this is a bug. – rcollyer Jan 6 '13 at 7:20
Related thread at MathGroup. – István Zachar Oct 9 '13 at 15:38

It's not you, it's Mathematica. You are not expected to know this, but basically, in compiled code, ReplacePart merely acts as syntactic sugar for setting a part, i.e.:

l = Range[3]; ReplacePart[l, 2 -> 0]
(* -> {1, 0, 3} *)


would be compiled (but see below) into exactly the same bytecode as

l = Range[3]; Block[{l = l}, l[[2]] = 0; l]
(* -> {1, 0, 3} *)


(which you can easily verify if you so desire).

It should be obvious from this equivalence that the pattern syntax allowed by ReplacePart in top-level code cannot be used in compiled code; indeed, no patterns or even Rules are accepted by the compiler. (As such, the example given above actually will not compile as written.) To get around the latter limitation it is possible to use an undocumented (since v6, but see version 5.2) variant syntax of ReplacePart which works both at the top level and in compiled code: replace the Rule with a Sequence, and reverse the order of the arguments. That is, the first example now looks like this:

l = Range[3]; ReplacePart[l, 0, 2]
(* -> {1, 0, 3} *)


So, this should suffice to explain (a) that ReplacePart can only be used in a strictly limited sense within compiled code and (b) the meaning of the strange message that you receive as a result of your attempt to use it in an unsupported way.

It happens in this case that you can re-write your ReplacePart call in terms of Part and Set:

Clear[longtermX];
longtermX[M_?MatrixQ] := Module[{len = Length[M], tmp},
tmp = Transpose[M] - IdentityMatrix[len];
tmp[[len, All]] = Sequence@ConstantArray[1, len];
tmp
];


However, this does not really help in practice because a ragged array (which the VM does not support) is produced as output, and this usage of Sequence together with Set/Part is not properly understood by the compiler anyway. What one might try instead, and which will work in compiled code, is something along the lines of tmp[[len, All]] = 1 in place of tmp[[len, All]] = Sequence@ConstantArray[1, len], then tmp2 = ConstantArray[Last[tmp], len]; tmp = Most[tmp] to have the two (full-rank) pieces of the array separately.

-
I think that the three argument form of ReplacePart was documented behavior for older versions (certainly before version 6). So I guess it might still be there for compatibility reasons... – Albert Retey Jan 6 '13 at 11:27
I don't understand the purpose of Sequence in your version of longtermX, which gives some wrong results on Mathematica 9. I found that replacing Sequence@ConstantArray[1, len] with ConstantArray[1, len] will give me correct results. Also, the compiled version of longtermX (by direct translating your version of longtermX) seems to work. At least I don't see any complaints from Mathematica. – wdg Jan 15 '13 at 13:28
@wdg My version gives the same results as yours (top-level code, not compiled). I inserted the Sequence to replicate what you are doing with ReplacePart, but the result of this seemed odd to me so possibly I replicated a mistake. Compile does not generate a warning but silently gives incorrect output. If the output it gives is actually "correct" then your ReplacePart call contains an error and Sequence is not needed after all. – Oleksandr R. Jan 15 '13 at 14:26