Skip to content

Commit 4d86b03

Browse files
bhallerbhaller
committed
add matrixPow() function to Eidos
1 parent 20a06af commit 4d86b03

5 files changed

Lines changed: 62 additions & 3 deletions

File tree

‎EidosScribe/EidosHelpFunctions.rtf‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4736,7 +4736,26 @@ Both
47364736
\f3\i0 matrix). Vectors will not be promoted to matrices by this function, even if such promotion would lead to a conformable matrix.\
47374737
\pard\pardeftab720\li720\fi-446\ri720\sb180\sa60\partightenfactor0
47384738

4739-
\f1\fs18 \cf0 \kerning1\expnd0\expndtw0 (integer$)nrow(*\'a0x)
4739+
\f1\fs18 \cf2 \kerning1\expnd0\expndtw0 (numeric)matrixPow(numeric\'a0x, integer$\'a0power)\
4740+
\pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0
4741+
4742+
\f0\b\fs20 \cf2 Returns the result of raising matrix
4743+
\f11\fs18 x
4744+
\f0\fs20 to an
4745+
\f11\fs18 integer
4746+
\f0\fs20
4747+
\f11\fs18 power
4748+
\f0\fs20 .
4749+
\f3\b0 The parameter x must be a square matrix (or an error will be raised). This operation is performed by repeated matrix multiplication with
4750+
\f1\fs18 matrixMult()
4751+
\f3\fs20 , and uses
4752+
\f1\fs18 inverse()
4753+
\f3\fs20 to compute the inverse of the matrix if
4754+
\f1\fs18 power
4755+
\f3\fs20 is negative.\
4756+
\pard\pardeftab720\li720\fi-446\ri720\sb180\sa60\partightenfactor0
4757+
4758+
\f1\fs18 \cf0 (integer$)nrow(*\'a0x)
47404759
\f2 \
47414760
\pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0
47424761

‎QtSLiM/help/EidosHelpFunctions.html‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -378,6 +378,8 @@
378378
<p class="p2">(numeric)matrixMult(numeric x, numeric y)</p>
379379
<p class="p5"><span class="s5"><b>Returns the result of matrix multiplication</b> of </span><span class="s6">x</span><span class="s5"> with </span><span class="s6">y</span><span class="s5">.<span class="Apple-converted-space">  </span>In Eidos (as in R), with two matrices </span><span class="s6">A</span><span class="s5"> and </span><span class="s6">B</span><span class="s5"> the simple product </span><span class="s6">A * B</span><span class="s5"> multiplies the corresponding elements of the matrices; in other words, if </span><span class="s6">X</span><span class="s5"> is the result of </span><span class="s6">A * B</span><span class="s5">, then </span><span class="s6">X</span><span class="s21"><i><sub>ij</sub></i></span><span class="s5"> = </span><span class="s6">A</span><span class="s21"><i><sub>ij</sub></i></span><span class="s5"> * </span><span class="s6">B</span><span class="s21"><i><sub>ij</sub></i></span><span class="s5">.<span class="Apple-converted-space">  </span>This is parallel to the definition of other operators; A + B adds the corresponding elements of the matrices (</span><span class="s6">X</span><span class="s21"><i><sub>ij</sub></i></span><span class="s5"> = </span><span class="s6">A</span><span class="s21"><i><sub>ij</sub></i></span><span class="s5"> + </span><span class="s6">B</span><span class="s21"><i><sub>ij</sub></i></span><span class="s5">), etc.<span class="Apple-converted-space">  </span>In R, true matrix multiplication is achieved with a special operator, </span><span class="s6">%*%</span><span class="s5">; in Eidos, the </span><span class="s6">matrixMult()</span><span class="s5"> function is used instead.</span></p>
380380
<p class="p5"><span class="s5">Both </span><span class="s6">x</span><span class="s5"> and </span><span class="s6">y</span><span class="s5"> must be matrices, and must be conformable according to the standard definition of matrix multiplication (i.e., if </span><span class="s6">x</span><span class="s5"> is an <i>n</i> × <i>m</i> matrix then </span><span class="s6">y</span><span class="s5"> must be a <i>m</i> × <i>p</i> matrix, and the result will be a <i>n</i> × <i>p</i> matrix).<span class="Apple-converted-space">  </span>Vectors will not be promoted to matrices by this function, even if such promotion would lead to a conformable matrix.</span></p>
381+
<p class="p4">(numeric)matrixPow(numeric x, integer$ power)</p>
382+
<p class="p5"><b>Returns the result of raising matrix </b><span class="s2"><b>x</b></span><b> to an </b><span class="s2"><b>integer</b></span><b> </b><span class="s2"><b>power</b></span><b>.</b><span class="Apple-converted-space">  </span>The parameter x must be a square matrix (or an error will be raised).<span class="Apple-converted-space">  </span>This operation is performed by repeated matrix multiplication with <span class="s2">matrixMult()</span>, and uses <span class="s2">inverse()</span> to compute the inverse of the matrix if <span class="s2">power</span> is negative.</p>
381383
<p class="p2">(integer$)nrow(* x)</p>
382384
<p class="p5"><span class="s5"><b>Returns the number of rows</b> in matrix or array </span><span class="s6">x</span><span class="s5">.<span class="Apple-converted-space">  </span>For vector </span><span class="s6">x</span><span class="s5">, </span><span class="s6">nrow()</span><span class="s5"> returns </span><span class="s6">NULL</span><span class="s5">; </span><span class="s6">size()</span><span class="s5"> should be used.<span class="Apple-converted-space">  </span>An equivalent of R’s </span><span class="s6">NROW()</span><span class="s5"> function, which treats vectors as </span><span class="s6">1</span><span class="s5">-column matrices, is not provided but would be trivial to implement as a user-defined function.</span></p>
383385
<p class="p2">(integer$)ncol(* x)</p>

‎VERSIONS‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ development head (in the master branch):
3131
partial fix for #530, add function (numeric)outerProduct(numeric x, numeric y) to obtain the outer product of vectors x and y
3232
partial fix for #530, add function (numeric$)det(numeric x) to obtain the determinant of a square matrix
3333
partial fix for #530, add function (float)inverse(numeric x) to obtain the inverse of a square non-singular matrix
34+
partial fix for #530, add function (numeric)matrixPow(numeric x, integer$ power) to raise a matrix to an integer power
3435

3536

3637
version 5.0 (Eidos version 4.0):

‎eidos/eidos_functions.cpp‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,6 +62,29 @@ R"V0G0N({
6262
return matrixMult(matrix(x), t(matrix(y)));
6363
})V0G0N";
6464

65+
// - (numeric)matrixPow(numeric x, integer$ power)
66+
const char *gEidosSourceCode_matrixPow =
67+
R"V0G0N({
68+
d = dim(x);
69+
70+
if (size(d) != 2)
71+
stop("ERROR (matrixPow): matrixPow() requires x to be a matrix");
72+
if (d[0] != d[1])
73+
stop("ERROR (matrixPow): matrixPow() requires x to be a square matrix");
74+
75+
if (power == 0)
76+
return diag(d[0]);
77+
else if (power == 1)
78+
return x;
79+
else if (power < 0) {
80+
// check for singularity; we could just let inverse() raise an error, but it would be more confusing for the user
81+
if (det(x) == 0)
82+
stop("ERROR (matrixPow): in matrixPow() the vector x is singular and thus non-invertible, so it cannot be raised to a negative power.");
83+
return matrixPow(inverse(x), -power);
84+
} else // (power >= 2)
85+
return matrixMult(x, matrixPow(x, power - 1));
86+
})V0G0N";
87+
6588

6689
//
6790
// Construct our built-in function map
@@ -355,6 +378,8 @@ const std::vector<EidosFunctionSignature_CSP> &EidosInterpreter::BuiltInFunction
355378
//
356379
signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature(gEidosStr_source, gEidosSourceCode_source, kEidosValueMaskVOID))->AddString_S(gEidosStr_filePath)->AddLogical_OS("chdir", gStaticEidosValue_LogicalF));
357380

381+
signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("matrixPow", gEidosSourceCode_matrixPow, kEidosValueMaskNumeric))->AddNumeric("x")->AddInt_S("power"));
382+
358383
signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("outerProduct", gEidosSourceCode_outerProduct, kEidosValueMaskNumeric))->AddNumeric("x")->AddNumeric("y"));
359384

360385

‎eidos/eidos_test_functions_other.cpp‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -61,9 +61,9 @@ void _RunFunctionMatrixArrayTests(void)
6161
EidosAssertScriptSuccess_I("det(matrix(c(1, 2, 3, 6, 5, 4, 8, 9, 7), nrow=3));", 21);
6262
EidosAssertScriptSuccess_L("abs(det(matrix(c(1.0, 2, 3, 6, 5, 4, 8, 9, 7), nrow=3)) - 21.0) < 1e-10;", true);
6363
EidosAssertScriptSuccess_L("x = matrix(rdunif(4, -100, 100), ncol=2); identical(det(x), drop(x[0,0]*x[1,1] - x[0,1]*x[1,0]));", true);
64-
EidosAssertScriptSuccess_L("x = matrix(runif(4, -100, 100), ncol=2); abs(det(x) - (x[0,0]*x[1,1] - x[0,1]*x[1,0])) < 1e-10;", true);
64+
EidosAssertScriptSuccess_L("x = matrix(runif(4, -100, 100), ncol=2); abs(det(x) - drop(x[0,0]*x[1,1] - x[0,1]*x[1,0])) < 1e-10;", true);
6565
EidosAssertScriptSuccess_L("x = matrix(rdunif(9, -100, 100), ncol=3); identical(det(x), drop(x[0,0]*x[1,1]*x[2,2] + x[0,1]*x[1,2]*x[2,0] + x[0,2]*x[1,0]*x[2,1] - x[0,2]*x[1,1]*x[2,0] - x[0,1]*x[1,0]*x[2,2] - x[0,0]*x[1,2]*x[2,1]));", true);
66-
EidosAssertScriptSuccess_L("x = matrix(runif(9, -100, 100), ncol=3); abs(det(x) - (x[0,0]*x[1,1]*x[2,2] + x[0,1]*x[1,2]*x[2,0] + x[0,2]*x[1,0]*x[2,1] - x[0,2]*x[1,1]*x[2,0] - x[0,1]*x[1,0]*x[2,2] - x[0,0]*x[1,2]*x[2,1])) < 1e-10;", true);
66+
EidosAssertScriptSuccess_L("x = matrix(runif(9, -100, 100), ncol=3); abs(det(x) - drop(x[0,0]*x[1,1]*x[2,2] + x[0,1]*x[1,2]*x[2,0] + x[0,2]*x[1,0]*x[2,1] - x[0,2]*x[1,1]*x[2,0] - x[0,1]*x[1,0]*x[2,2] - x[0,0]*x[1,2]*x[2,1])) < 1e-10;", true);
6767
EidosAssertScriptSuccess_L("x = matrix(rdunif(16, -100, 100), ncol=4); det(x); T;", true);
6868
EidosAssertScriptSuccess_L("x = matrix(runif(16, -100, 100), ncol=4); det(x); T;", true);
6969

@@ -236,6 +236,18 @@ void _RunFunctionMatrixArrayTests(void)
236236
EidosAssertScriptSuccess_L("A = matrix(1.0:6.0, nrow=2); B = matrix(1.0:6.0, ncol=2); identical(matrixMult(A, B), matrix(c(22.0, 28.0, 49.0, 64.0), nrow=2));", true);
237237
EidosAssertScriptSuccess_L("A = matrix(1.0:6.0, ncol=2); B = matrix(1.0:6.0, nrow=2); identical(matrixMult(A, B), matrix(c(9.0, 12.0, 15.0, 19.0, 26.0, 33.0, 29.0, 40.0, 51.0), nrow=3));", true);
238238

239+
// matrixPow()
240+
EidosAssertScriptRaise("matrixPow(5.0, 3);", 0, "requires x to be a matrix");
241+
EidosAssertScriptRaise("matrixPow(matrix(1:6, nrow=2), 3);", 0, "requires x to be a square matrix");
242+
EidosAssertScriptSuccess_L("A = matrix(1:9, nrow=3); B = matrixPow(A, 0); identical(B, diag(3));", true);
243+
EidosAssertScriptSuccess_L("A = matrix(1:9, nrow=3); B = matrixPow(A, 1); identical(B, A);", true);
244+
EidosAssertScriptSuccess_L("A = matrix(1:9, nrow=3); B = matrixPow(A, 2); identical(B, matrix(c(30, 36, 42, 66, 81, 96, 102, 126, 150), nrow=3));", true);
245+
EidosAssertScriptSuccess_L("A = matrix(1:9, nrow=3); B = matrixPow(A, 3); identical(B, matrix(c(468, 576, 684, 1062, 1305, 1548, 1656, 2034, 2412), nrow=3));", true);
246+
EidosAssertScriptRaise("A = matrix(1:9, nrow=3); B = matrixPow(A, -1);", 29, "singular and thus non-invertible");
247+
EidosAssertScriptSuccess_L("A = matrix(c(1, 3, 2, 6, 4, 5, 10, 2, 7), nrow=3); B = matrixPow(A, -1); identical(B, inverse(A));", true);
248+
EidosAssertScriptSuccess_L("A = matrix(c(1, 3, 2, 6, 4, 5, 10, 2, 7), nrow=3); B = matrixPow(A, -2); identical(B, matrixPow(inverse(A), 2));", true);
249+
EidosAssertScriptSuccess_L("A = matrix(c(1, 3, 2, 6, 4, 5, 10, 2, 7), nrow=3); B = matrixPow(A, -3); identical(B, matrixPow(inverse(A), 3));", true);
250+
239251
// ncol()
240252
EidosAssertScriptSuccess_NULL("ncol(NULL);");
241253
EidosAssertScriptSuccess_NULL("ncol(T);");

0 commit comments

Comments
 (0)