Determinant Equation \(W(s) = sI - H(s)\)¶
The function \(H(s)\) is an integral defined as:
It is possible the function \(H\) as well as its determinant using the DeterminantEq objects. This is the object used when solving for the approximate missed-events likelihood. The determinant equation is initialized in one of two ways, either from a matrix or QMatrix.
| python: | from dcprogs.likelihood import DeterminantEq, QMatrix
# Define parameters.
matrix = [ [-3050, 50, 3000, 0, 0],
[2./3., -1502./3., 0, 500, 0],
[ 15, 0, -2065, 50, 2000],
[ 0, 15000, 4000, -19000, 0],
[ 0, 0, 10, 0, -10] ]
qmatrix = QMatrix(matrix, 2)
tau = 1e-4
det0 = DeterminantEq(qmatrix, 1e-4);
det1 = DeterminantEq(matrix, 2, 1e-4);
print("{0!s}\n\n{1!s}".format(det0, det1))
|
|---|---|
| c++11: | #include <iostream>
#include <likelihood/determinant_equation.h>
#include <likelihood/root_finder.h>
int main() {
// Define parameters.
DCProgs::t_rmatrix matrix(5 ,5);
matrix << -3050, 50, 3000, 0, 0,
2./3., -1502./3., 0, 500, 0,
15, 0, -2065, 50, 2000,
0, 15000, 4000, -19000, 0,
0, 0, 10, 0, -10;
DCProgs::QMatrix qmatrix(matrix, /*nopen=*/2);
// Create determinant using a QMatrix and a matrix+nopen.
DCProgs::DeterminantEq det0(qmatrix, 1e-4);
DCProgs::DeterminantEq det1(matrix, 2, 1e-4);
std::cout << det0 << "\n\n" << det1 << "\n";
<< " * d[sI-H(s)]/ds for s=-1e2:\n" << DCProgs::numpy_io(det0.s_derivative(-1e2)) << "\n\n";
return 0;
}
|
With an object in hand, it is possible to compute \(\mathop{det}W(s)\) for any \(s\). In the following we demonstrate that the two initialization methods are equivalent and that the determinant is zero at the roots of \(W(s)\), per definition.
| python: | The python bindings accept both scalars and arrays as input. x = [0, -1, -1e2]
assert all(abs(det0(x) - det1(x)) < 1e-8)
roots = array([-3045.285776037674, -162.92946543451328])
assert all(abs(det0(roots)) < 1e-6 * abs(roots))
|
|---|---|
| c++11: | if( std::abs(det0(0) - det1(0)) > 1e-6
or std::abs(det0(-1) - det1(-1)) > 1e-6
or std::abs(det0(-1e2) - det1(-1e2)) > 1e-6)
throw DCProgs::errors::Runtime("instanciations differ.");
if( std::abs( det0(-3045.285776037674) ) > 1e-6 * 3e3
or std::abs( det0(-162.92946543451328) ) > 1e-6 * 2e2 )
throw DCProgs::errors::Runtime("Roots are not roots.");
|
There exists a convenience function to transform a determinant equation into its “transpose”, e.g. one where A states become F states and F states become A states:
| python: | transpose = det0.transpose()
roots = array( [-17090.192769236815, -2058.0812921673496, -0.24356535498785126],
dtype=internal_dtype )
assert all(abs(transpose(roots)) < 1e-6 * abs(roots))
Note Here we choose to create an input which has same internal type as the dcprogs package. This may result in faster code since no conversion are required. |
|---|---|
| c++11: | DCProgs::DeterminantEq transpose = det0.transpose();
if( std::abs( transpose(-17090.192769236815) ) > 1e-6 * 2e5
or std::abs( transpose(-2058.0812921673496) ) > 1e-6 * 2e3
or std::abs( transpose(-0.24356535498785126) ) > 1e-6 )
throw DCProgs::errors::Runtime("Roots are not roots.");
|
Finally, it is possible to compute \(H(s)\) directly, as well as \(\frac{\partial W(s)}{\partial s}\), as demonstrated below.
| python: | print(" * H({0}):\n{1}\n".format(0, det0.H(0)))
print(" * H({0}):\n{1}\n".format(-1e2, det0.H(-1e2)))
print(" * d[sI-H(s)]/ds for s={0}:\n{1}\n".format(0, det0.s_derivative(0)))
print(" * d[sI-H(s)]/ds for s={0}:\n{1}\n".format(-1e2, det0.s_derivative(-1e2)))
|
|---|---|
| c++11: | std::cout << " * H(0):\n" << DCProgs::numpy_io(det0.H(0)) << "\n\n"
<< " * H(-1e2):\n" << DCProgs::numpy_io(det0.H(-1e2)) << "\n\n";
std::cout << " * d[sI-H(s)]/ds for s=0:\n" << DCProgs::numpy_io(det0.s_derivative(0)) << "\n\n"
<< " * d[sI-H(s)]/ds for s=-1e2:\n" << DCProgs::numpy_io(det0.s_derivative(-1e2)) << "\n\n";
|